<p><b>mpetersen@lanl.gov</b> 2013-04-19 11:49:14 -0600 (Fri, 19 Apr 2013)</p><p>branch commit, tensor_operations: tested strain rate and divergence of tensor subroutines.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/tensor_operations/namelist.input.ocean
===================================================================
--- branches/ocean_projects/tensor_operations/namelist.input.ocean 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/namelist.input.ocean 2013-04-19 17:49:14 UTC (rev 2779)
@@ -2,7 +2,7 @@
config_do_restart = .false.
config_start_time = '0000-01-01_00:00:00'
config_stop_time = 'none'
- config_run_duration = '0_06:00:00'
+ config_run_duration = '0_00:00:00'
config_calendar_type = '360day'
/
&io
@@ -10,11 +10,11 @@
config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
config_restart_interval = '0_06:00:00'
- config_output_interval = '0_06:00:00'
+ config_output_interval = '0_00:00:01'
config_stats_interval = '0_01:00:00'
config_write_stats_on_startup = .true.
config_write_output_on_startup = .true.
- config_frames_per_outfile = 1000
+ config_frames_per_outfile = 1
config_pio_num_iotasks = 0
config_pio_stride = 1
/
Modified: branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml
===================================================================
--- branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml 2013-04-19 17:49:14 UTC (rev 2779)
@@ -795,6 +795,27 @@
<var name="divTensor3DCell" type="real" dimensions="R3 nVertLevels nEdges Time" streams="o" units="s^{-2}"
description="divergence of the tensor at cell center, as a 3-vector in x,y,z space"
/>
+ <var name="outerProductEdge" type="real" dimensions="SIX nVertLevels nEdges Time" streams="o" units="m^2 s^{-1}"
+ description="Outer produce, u_e \cross n_e, at each edge. {\color{red} change to scratch}"
+ />
+ <var name="strainRate3DCellSolution" type="real" dimensions="SIX nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="strain rate solution tensor at cell center, 3D, in symmetric 6-index form"
+ />
+ <var name="strainRate3DEdgeSolution" type="real" dimensions="SIX nVertLevels nEdges Time" streams="o" units="s^{-1}"
+ description="strain rate solution tensor at edge, 3D, in symmetric 6-index form"
+ />
+ <var name="divTensor3DCellSolution" type="real" dimensions="R3 nVertLevels nEdges Time" streams="o" units="s^{-2}"
+ description="divergence of the tensor solution at cell center, as a 3-vector in x,y,z space"
+ />
+ <var name="strainRate3DCellDiff" type="real" dimensions="SIX nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="strain rate difference between solution and computed tensor at cell center, 3D, in symmetric 6-index form"
+ />
+ <var name="strainRate3DEdgeDiff" type="real" dimensions="SIX nVertLevels nEdges Time" streams="o" units="s^{-1}"
+ description="strain rate difference between solution and computed tensor at edge, 3D, in symmetric 6-index form"
+ />
+ <var name="divTensor3DCellDiff" type="real" dimensions="R3 nVertLevels nEdges Time" streams="o" units="s^{-2}"
+ description="divergence of the tensor difference between solution and computed at cell center, as a 3-vector in x,y,z space"
+ />
</var_struct>
<var_struct name="mesh" time_levs="0">
<var name="latCell" type="real" dimensions="nCells" streams="iro" units="radians"
@@ -908,9 +929,6 @@
<var name="cellTangentPlane" type="real" dimensions="R3 TWO nCells" streams="o" units="unitless"
description="The two vectors that define a tangent plane at a cell center."
/>
- <var name="outerProductEdgeSym6" type="real" dimensions="SIX nVertLevels nEdges" streams="o" units="m^2 s^{-1}"
- description="Outer produce, u_e \cross n_e, at each edge. {\color{red} change to scratch}"
- />
<var name="cellsOnCell" type="integer" dimensions="maxEdges nCells" streams="iro" units="unitless"
description="List of cells that neighbor each cell."
/>
Modified: branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F 2013-04-19 17:49:14 UTC (rev 2779)
@@ -5,6 +5,8 @@
use mpas_timekeeping
use mpas_dmpar
use mpas_timer
+ ! mrp added to test stress tensor:
+ use mpas_tensor_operations
use ocn_global_diagnostics
use ocn_time_integration
@@ -163,6 +165,13 @@
block => block % next
end do
+ ! mrp added just to test stress tensor ops:
+ print *, 'mpas_test_stress_tensor'
+ call mpas_test_tensor(domain)
+ print *, 'mpas_test_stress_tensor complete'
+! stop
+ ! mrp added just to test stress tensor ops end
+
if (config_write_stats_on_startup) then
call mpas_timer_start("global diagnostics", .false., globalDiagTimer)
call ocn_compute_global_diagnostics(domain, 1 , 0, dt)
Modified: branches/ocean_projects/tensor_operations/src/operators/mpas_matrix_operations.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/mpas_matrix_operations.F 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/src/operators/mpas_matrix_operations.F 2013-04-19 17:49:14 UTC (rev 2779)
@@ -364,6 +364,9 @@
matrixEdge(:,k,iEdge) = 0.5*(matrixCell(:,k,cell1) + matrixCell(:,k,cell2))
enddo
+!print '(/,a,30es10.2)','cel1',matrixCell(:,1,cell1)
+!print '(a,30es10.2)','cel2',matrixCell(:,1,cell2)
+!print '(a,30es10.2)','edge',matrixEdge(:,1,iEdge)
enddo
end subroutine mpas_matrix_cell_to_edge!}}}
Modified: branches/ocean_projects/tensor_operations/src/operators/mpas_tensor_operations.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/mpas_tensor_operations.F 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/src/operators/mpas_tensor_operations.F 2013-04-19 17:49:14 UTC (rev 2779)
@@ -66,7 +66,8 @@
!
!-----------------------------------------------------------------------
- subroutine mpas_strain_rate(grid, normalVelocity, tangentialVelocity, strainRate3DCell)!{{{
+ subroutine mpas_strain_rate(grid, normalVelocity, tangentialVelocity, &
+ strainRate3DCell, outerProductEdge)!{{{
!-----------------------------------------------------------------
!
@@ -90,6 +91,9 @@
real (kind=RKIND), dimension(:,:,:), intent(out) :: &
strainRate3DCell !< Output: strain rate tensor at cell center, 3D, in symmetric 6-index form
+ real (kind=RKIND), dimension(:,:,:), intent(out) :: outerProductEdge
+
+
!-----------------------------------------------------------------
!
! local variables
@@ -103,25 +107,38 @@
real (kind=RKIND) :: invAreaCell
real (kind=RKIND), dimension(3,3) :: outerProductEdge3x3
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, angleEdge
+ real (kind=RKIND), dimension(:), pointer :: xCell,yCell,xEdge,yEdge ! temp
real (kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors, edgeTangentVectors
- real (kind=RKIND), dimension(:,:,:), pointer :: outerProductEdgeSym6
nEdges = grid % nEdges
nCells = grid % nCells
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelCell => grid % maxLevelCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnCell => grid % edgesOnCell % array
edgeSignOnCell => grid % edgeSignOnCell % array
dvEdge => grid % dvEdge % array
+ angleEdge => grid % angleEdge % array
areaCell => grid % areaCell % array
edgeNormalVectors => grid % edgeNormalVectors % array
edgeTangentVectors => grid % edgeTangentVectors % array
- ! mrp change outerProductEdgeSym6 to a scratch variable later
- outerProductEdgeSym6 => grid % outerProductEdgeSym6 % array
+! temp:
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ xEdge => grid % xEdge % array
+ yEdge => grid % yEdge % array
+if (1==2) then
+ print '(a,3f20.10)', 'min max mean edgeNormalVectors1 ',minval(edgeNormalVectors(1,1:nEdges)),maxval(edgeNormalVectors(1,1:nEdges)),sum(edgeNormalVectors(1,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean edgeNormalVectors2 ',minval(edgeNormalVectors(2,1:nEdges)),maxval(edgeNormalVectors(2,1:nEdges)),sum(edgeNormalVectors(2,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean edgeNormalVectors3 ',minval(edgeNormalVectors(3,1:nEdges)),maxval(edgeNormalVectors(3,1:nEdges)),sum(edgeNormalVectors(3,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean edgeTangentVectors1 ',minval(edgeTangentVectors(1,1:nEdges)),maxval(edgeTangentVectors(1,1:nEdges)),sum(edgeTangentVectors(1,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean edgeTangentVectors2 ',minval(edgeTangentVectors(2,1:nEdges)),maxval(edgeTangentVectors(2,1:nEdges)),sum(edgeTangentVectors(2,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean edgeTangentVectors3 ',minval(edgeTangentVectors(3,1:nEdges)),maxval(edgeTangentVectors(3,1:nEdges)),sum(edgeTangentVectors(3,1:nEdges))/nEdges/grid%nVertLevels
+endif
do iEdge=1,nEdges
! mrp question: for this to be general across cores, we need to use nVertLevels
do k=1,maxLevelEdgeBot(iEdge)
@@ -132,15 +149,31 @@
outerProductEdge3x3(i,j) = edgeNormalVectors(i,iEdge) &
*( normalVelocity(k,iEdge) *edgeNormalVectors(j,iEdge) &
+ tangentialVelocity(k,iEdge)*edgeTangentVectors(j,iEdge) &
- ) * dvEdge(iEdge)
+ )
enddo
enddo
- call mpas_matrix_3x3_to_sym6index(outerProductEdge3x3,outerProductEdgeSym6(:,iEdge,k))
+ call mpas_matrix_3x3_to_sym6index(outerProductEdge3x3,outerProductEdge(:,k,iEdge))
+if (1==2) then
+ print '(a,i5,30f8.2)', '*********** iEdge, angleEdge',iEdge, angleEdge(iEdge)
+ print '(a,30es10.2)', 'nor',edgeNormalVectors(:,iEdge)
+ print '(a,30es10.2)', 'tan',edgeTangentVectors(:,iEdge)
+ print '(a,30es10.2)', 'vel',normalVelocity(k,iEdge),tangentialVelocity(k,iEdge)
+ print '(a,30es10.2)', '3x3',outerProductEdge3x3
+ print '(a,30es10.2)', 'sm6', outerProductEdge(:,k,iEdge)
+endif
enddo
enddo
+if (1==2) then
+ print '(a,3f20.10)', 'min max mean outerProductEdge 1 ',minval(outerProductEdge(1,1,1:nEdges)),maxval(outerProductEdge(1,1,1:nEdges)),sum(outerProductEdge(1,1,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean outerProductEdge 2 ',minval(outerProductEdge(2,1,1:nEdges)),maxval(outerProductEdge(2,1,1:nEdges)),sum(outerProductEdge(2,1,1:nEdges))/nEdges/grid%nVertLevels
+ print '(a,3f20.10)', 'min max mean outerProductEdge 4 ',minval(outerProductEdge(4,1,1:nEdges)),maxval(outerProductEdge(4,1,1:nEdges)),sum(outerProductEdge(4,1,1:nEdges))/nEdges/grid%nVertLevels
+print *, 'iCell, i, k, edgeSignOnCell(i,iCell),invAreaCell,outerProductEdge(:,k,iEdge)'
+endif
strainRate3DCell = 0.0
do iCell = 1, nCells
+! print '(a,i5,30f8.2)', '***************************************************** iCell,xCell,yCell',iCell,xCell(iCell),yCell(iCell)
+! print '(a,30i5,30f8.2)', '** nCells,cellsOnCell',nCells,grid % cellsOnCell % array(:,iCell)
invAreaCell = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
@@ -148,9 +181,26 @@
! mrp edgeSignOnCell is to get outward unit normal on edgeNormalVectors
! minus sign in front is to match form on divergence, below
strainRate3DCell(:,k,iCell) = strainRate3DCell(:,k,iCell) &
- - edgeSignOnCell(i,iCell)*outerProductEdgeSym6(:,iEdge,k)*invAreaCell
+ - edgeSignOnCell(i,iCell)*outerProductEdge(:,k,iEdge)*invAreaCell*dvEdge(iEdge)
+if (1==2) then
+!if (iCell>=15.and.iCell<=17) then
+ print '(a,2i5,30f8.2)', '*********** iEdge, edgeSignOnCell, angleEdge,xEdge,yEdge',iEdge, edgeSignOnCell(i,iCell), angleEdge(iEdge),xEdge(iEdge),yEdge(iEdge)
+ print '(a,2i5,30f8.2)', 'cell1, cell2',grid%cellsOnEdge%array(:,iEdge)
+ print '(a,30es10.2)', 'nor',edgeNormalVectors(:,iEdge)
+ print '(a,30es10.2)', 'tan',edgeTangentVectors(:,iEdge)
+ print '(a,30es10.2)', 'vel',normalVelocity(k,iEdge),tangentialVelocity(k,iEdge)
+ print '(a,30es10.2)', 'sm6', outerProductEdge(:,k,iEdge)
+ print '(a,30es10.2)', 'src', strainRate3DCell(:,k,iCell)
+endif
end do
end do
+if (1==2) then
+ print '(a,30es10.2)', 'nor',edgeNormalVectors(:,iEdge)
+ print '(a,30es10.2)', 'tan',edgeTangentVectors(:,iEdge)
+ print '(a,30es10.2)', 'vel',normalVelocity(k,iEdge),tangentialVelocity(k,iEdge)
+ print '(a,30es10.2)', '3x3',outerProductEdge3x3
+ print '(a,30es10.2)', 'sm6', outerProductEdge(:,k,iEdge)
+endif
end do
end subroutine mpas_strain_rate!}}}
@@ -214,6 +264,7 @@
maxLevelCell => grid % maxLevelCell % array
edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
edgeNormalVectors => grid % edgeNormalVectors % array
@@ -229,7 +280,7 @@
edgeNormalDotTensor(:) = 0.0
do q=1,3
do p=1,3
- edgeNormalDotTensor(q) = edgeNormalVectors(p,iEdge)*strainRate3DEdge3x3(p,q)
+ edgeNormalDotTensor(q) = edgeNormalDotTensor(q) + edgeNormalVectors(p,iEdge)*strainRate3DEdge3x3(p,q)
enddo
enddo
divTensor3DCell(:,k,iCell) = divTensor3DCell(:,k,iCell) &
@@ -477,9 +528,38 @@
type (block_type), pointer :: block
+ integer :: nCells, nEdges, nVertices, nVertLevels, iCell, iEdge, k, i, p, minCell, maxCell
+
+ real (kind=RKIND) :: xVelocity, yVelocity, cn, cs, r, theta, rot, f, g1, g2, fcos, pi2l, ld, &
+ divSigma_x, divSigma_y
+ real (kind=RKIND), dimension(:), pointer :: angleEdge, xCell, yCell, xEdge, yEdge
real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, tangentialVelocity
+ real (kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors, edgeTangentVectors
real (kind=RKIND), dimension(:,:,:), pointer :: strainRate3DCell, strainRate3DEdge, divTensor3DCell
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRate3DCellSolution, divTensor3DCellSolution
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRate3DCellDiff, divTensor3DCellDiff, outerProductEdge
+ character(len=16) :: test_function_type
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+
+ xCell => block % mesh % xCell % array
+ yCell => block % mesh % yCell % array
+ xEdge => block % mesh % xEdge % array
+ yEdge => block % mesh % yEdge % array
+ angleEdge => block % mesh % angleEdge % array
+
+ edgeNormalVectors => block % mesh % edgeNormalVectors % array
+ edgeTangentVectors => block % mesh % edgeTangentVectors % array
+
+
normalVelocity => block % state % time_levs(1) % state % normalVelocity % array
tangentialVelocity => block % state % time_levs(1) % state % tangentialVelocity % array
@@ -487,23 +567,247 @@
strainRate3DCell => block % state % time_levs(1) % state % strainRate3DCell % array
strainRate3DEdge => block % state % time_levs(1) % state % strainRate3DEdge % array
divTensor3DCell => block % state % time_levs(1) % state % divTensor3DCell % array
+ outerProductEdge => block % state % time_levs(1) % state % outerProductEdge % array
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
+ strainRate3DCellSolution => block % state % time_levs(1) % state % strainRate3DCellSolution % array
+ divTensor3DCellSolution => block % state % time_levs(1) % state % divTensor3DCellSolution % array
+ strainRate3DCellDiff => block % state % time_levs(1) % state % strainRate3DCellDiff % array
+ divTensor3DCellDiff => block % state % time_levs(1) % state % divTensor3DCellDiff % array
-! see old branch for testing code.
+ ! create test functions for normalVelocity and tangentialVelocity
+ normalVelocity = 0.0
+ tangentialVelocity = 0.0
+ strainRate3DCellSolution = 0.0
+ divTensor3DCellSolution = 0.0
- ! create test functions for normalVelocity and tangentialVelocity
- normalVelocity = 0.0
- tangentialVelocity = 0.0
+ cn = 15.0
+ cs = 20.0
+ rot = 3.0 !pii/4.0 !0.0*pii/2
+ p = 2 ! power for polynomial function
+! test_function_type = 'constant'
+ test_function_type = 'linear_x'
+! test_function_type = 'linear_y'
+! test_function_type = 'linear_arb_rot'
+! test_function_type = 'sin_arb_rot'
+! test_function_type = 'power_x'
+! test_function_type = 'power_y'
+! test_function_type = 'power_arb_rot'
- call mpas_strain_rate(block % mesh, normalVelocity, tangentialVelocity, strainRate3DCell)
+ ld = 100.0e3 ! wavelength in meters
+ pi2l = 2*pii/ld
+ g1 = cn*cos(rot) - cs*sin(rot)
+ g2 = cn*sin(rot) + cs*cos(rot)
- call mpas_matrix_cell_to_edge(block % mesh, strainRate3DCell, strainRate3DEdge)
+ if (test_function_type.eq.'constant') then
+ ! linear function in x
+ do iEdge = 1,nEdges
+ xVelocity = cn
+ yVelocity = cs
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
- call mpas_divergence_of_tensor(block % mesh, strainRate3DEdge, divTensor3DCell)
+ strainRate3DCellSolution = 0.0
+ divTensor3DCellSolution = 0.0
+ elseif (test_function_type.eq.'linear_x') then
+ ! linear function in x
+ do iEdge = 1,nEdges
+ xVelocity = cn*xEdge(iEdge)
+ yVelocity = cs*xEdge(iEdge)
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
+
+ strainRate3DCellSolution(1,:,:) = cn
+ strainRate3DCellSolution(2,:,:) = 0.0
+ strainRate3DCellSolution(4,:,:) = 0.5*cs
+
+ divTensor3DCellSolution = 0.0
+
+ elseif (test_function_type.eq.'linear_y') then
+
+ ! linear function in y
+ do iEdge = 1,nEdges
+ xVelocity = -cs*yEdge(iEdge)
+ yVelocity = cn*yEdge(iEdge)
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
+
+ strainRate3DCellSolution(1,:,:) = 0.0
+ strainRate3DCellSolution(2,:,:) = cn
+ strainRate3DCellSolution(4,:,:) = -0.5*cs
+
+ divTensor3DCellSolution = 0.0
+
+ elseif (test_function_type.eq.'linear_arb_rot') then
+
+ ! linear function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+
+ f = r*cos(theta-rot)
+ xVelocity = f*g1
+ yVelocity = f*g2
+
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+
+ enddo
+
+ strainRate3DCellSolution(1,:,:) = cos(rot)*g1
+ strainRate3DCellSolution(2,:,:) = sin(rot)*g2
+ strainRate3DCellSolution(4,:,:) = 0.5*(cos(rot)*g2 + sin(rot)*g1)
+
+ divTensor3DCellSolution = 0.0
+
+ elseif (test_function_type.eq.'power_x') then
+
+ ! power function in x: x^p
+ do iEdge = 1,nEdges
+ xVelocity = cn*xEdge(iEdge)**p
+ yVelocity = cs*xEdge(iEdge)**p
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
+
+ do iCell = 1,nCells
+ strainRate3DCellSolution(1,:,iCell) = cn *p*xCell(iCell)**(p-1)
+ strainRate3DCellSolution(2,:,iCell) = 0.0
+ strainRate3DCellSolution(4,:,iCell) = 0.5*cs*p*xCell(iCell)**(p-1)
+
+ divTensor3DCellSolution(1,:,iCell) = cn *p*(p-1)*xCell(iCell)**(p-2)
+ divTensor3DCellSolution(2,:,iCell) = 0.5*cs*p*(p-1)*xCell(iCell)**(p-2)
+ end do
+
+ elseif (test_function_type.eq.'power_y') then
+
+ ! power function in y: y^n
+ do iEdge = 1,nEdges
+ xVelocity = -cs*yEdge(iEdge)**p
+ yVelocity = cn*yEdge(iEdge)**p
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
+
+ do iCell = 1,nCells
+ strainRate3DCellSolution(1,:,iCell) = 0.0
+ strainRate3DCellSolution(2,:,iCell) = cn *p*yCell(iCell)**(p-1)
+ strainRate3DCellSolution(4,:,iCell) = -0.5*cs*p*yCell(iCell)**(p-1)
+
+ divTensor3DCellSolution(1,:,iCell) = -0.5*cs*p*(p-1)*yCell(iCell)**(p-2)
+ divTensor3DCellSolution(2,:,iCell) = cn *p*(p-1)*yCell(iCell)**(p-2)
+ end do
+
+ elseif (test_function_type.eq.'power_arb_rot') then
+
+ ! power function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ f = r*cos(theta-rot)
+
+ xVelocity = g1*f**p
+ yVelocity = g2*f**p
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+ enddo
+
+ do iCell = 1,nCells
+ r = sqrt(xCell(iCell)**2 + yCell(iCell)**2)
+ theta = atan(yCell(iCell)/xCell(iCell))
+ f = r*cos(theta-rot)
+
+ strainRate3DCellSolution(1,:,iCell) = p *f**(p-1) *cos(rot)*g1
+ strainRate3DCellSolution(2,:,iCell) = p *f**(p-1) *sin(rot)*g2
+ strainRate3DCellSolution(4,:,iCell) = p *f**(p-1) *(cos(rot)*g2+sin(rot)*g1)/2.0
+
+ divTensor3DCellSolution(1,:,iCell) = p*(p-1)*f**(p-2) *(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+ divTensor3DCellSolution(2,:,iCell) = p*(p-1)*f**(p-2) *(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+ end do
+
+ elseif (test_function_type.eq.'sin_arb_rot') then
+
+ ! sin function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ f = sin(pi2l*r*cos(theta-rot))
+
+ xVelocity = f*g1
+ yVelocity = f*g2
+ normalVelocity(:,iEdge) = xVelocity*edgeNormalVectors(1,iEdge) + yVelocity*edgeNormalVectors(2,iEdge)
+ tangentialVelocity(:,iEdge) = xVelocity*edgeTangentVectors(1,iEdge) + yVelocity*edgeTangentVectors(2,iEdge)
+
+ ! solution for div(sigma)
+! divTensor3DCellSolution(1,:,iEdge) = -pi2l**2*f*(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+! divTensor3DCellSolution(2,:,iEdge) = -pi2l**2*f*(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+
+ enddo
+
+ do iCell = 1,nCells
+ r = sqrt(xCell(iCell)**2 + yCell(iCell)**2)
+ theta = atan(yCell(iCell)/xCell(iCell))
+ fcos = cos(pi2l*r*cos(theta-rot))
+
+ strainRate3DCellSolution(1,:,iCell) = pi2l*fcos*cos(rot)*g1
+ strainRate3DCellSolution(2,:,iCell) = pi2l*fcos*sin(rot)*g2
+ strainRate3DCellSolution(3,:,iCell) = pi2l*fcos*(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+ else
+ print *, 'bad choice of test_function_type:',test_function_type
+ stop
+ endif
+
+! call mpas_tangential_velocity(block % mesh, normalVelocity, tangentialVelocity)
+
+ call mpas_strain_rate(block % mesh, normalVelocity, tangentialVelocity, strainRate3DCell, outerProductEdge)
+
+ call mpas_matrix_cell_to_edge(block % mesh, strainRate3DCell, strainRate3DEdge)
+
+ call mpas_divergence_of_tensor(block % mesh, strainRate3DEdge, divTensor3DCell)
+
+
+minCell=17
+maxCell=17
+if (1==1) then
+! print '(a,3es20.10)', 'min max mean normalVelocity ',minval(normalVelocity(:,1:nEdges)),maxval(normalVelocity(:,1:nEdges)),sum(normalVelocity(:,1:nEdges))/nEdges/nVertLevels
+! print '(a,3es20.10)', 'min max mean tangentialVelocity ',minval(tangentialVelocity(:,1:nEdges)),maxval(tangentialVelocity(:,1:nEdges)),sum(tangentialVelocity(:,1:nEdges))/nEdges/nVertLevels
+ print '(a,3es20.10)', 'min max mean strainRate3DCell1 ',minval(strainRate3DCell(1,:,minCell:maxCell)),maxval(strainRate3DCell(1,:,minCell:maxCell)),sum(strainRate3DCell(1,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean strainRate3DCellSol1',minval(strainRate3DCellSolution(1,:,minCell:maxCell)),maxval(strainRate3DCellSolution(1,:,minCell:maxCell)),sum(strainRate3DCellSolution(1,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(/,a,3es20.10)', 'min max mean strainRate3DCell2 ',minval(strainRate3DCell(2,:,minCell:maxCell)),maxval(strainRate3DCell(2,:,minCell:maxCell)),sum(strainRate3DCell(2,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean strainRate3DCell2Sol',minval(strainRate3DCellSolution(2,:,minCell:maxCell)),maxval(strainRate3DCellSolution(2,:,minCell:maxCell)),sum(strainRate3DCellSolution(2,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(/,a,3es20.10)', 'min max mean strainRate3DCell4 ',minval(strainRate3DCell(4,:,minCell:maxCell)),maxval(strainRate3DCell(4,:,minCell:maxCell)),sum(strainRate3DCell(4,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean strainRate3DCell4Sol',minval(strainRate3DCellSolution(4,:,minCell:maxCell)),maxval(strainRate3DCellSolution(4,:,minCell:maxCell)),sum(strainRate3DCellSolution(4,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(/,a,3es20.10)', 'min max mean divTensor3DCell1 ',minval(divTensor3DCell(1,:,minCell:maxCell)),maxval(divTensor3DCell(1,:,minCell:maxCell)),sum(divTensor3DCell(1,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean divTensor3DCellSol1 ',minval(divTensor3DCellSolution(1,:,minCell:maxCell)),maxval(divTensor3DCellSolution(1,:,minCell:maxCell)),sum(divTensor3DCellSolution(1,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+! print '(/,a,3es20.10)', 'min max mean strainRate3DEdge1 ',minval(strainRate3DEdge(1,:,minCell:maxCell)),maxval(strainRate3DEdge(1,:,minCell:maxCell)),sum(strainRate3DEdge(1,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+! print '(a,3es20.10)', 'min max mean strainRate3DEdge2 ',minval(strainRate3DEdge(2,:,minCell:maxCell)),maxval(strainRate3DEdge(2,:,minCell:maxCell)),sum(strainRate3DEdge(2,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+! print '(a,3es20.10)', 'min max mean strainRate3DEdge3 ',minval(strainRate3DEdge(3,:,minCell:maxCell)),maxval(strainRate3DEdge(3,:,minCell:maxCell)),sum(strainRate3DEdge(3,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+
+ print '(/,a,3es20.10)', 'min max mean divTensor3DCell2 ',minval(divTensor3DCell(2,:,minCell:maxCell)),maxval(divTensor3DCell(2,:,minCell:maxCell)),sum(divTensor3DCell(2,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean divTensor3DCellSol2 ',minval(divTensor3DCellSolution(2,:,minCell:maxCell)),maxval(divTensor3DCellSolution(2,:,minCell:maxCell)),sum(divTensor3DCellSolution(2,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(/,a,3es20.10)', 'min max mean divTensor3DCell3 ',minval(divTensor3DCell(3,:,minCell:maxCell)),maxval(divTensor3DCell(3,:,minCell:maxCell)),sum(divTensor3DCell(3,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+ print '(a,3es20.10)', 'min max mean divTensor3DCellSol3 ',minval(divTensor3DCellSolution(3,:,minCell:maxCell)),maxval(divTensor3DCellSolution(3,:,minCell:maxCell)),sum(divTensor3DCellSolution(3,:,minCell:maxCell))/(maxCell-minCell+1)/nVertLevels
+endif
+
+ block % state % time_levs(1) % state % strainRate3DCellDiff % array = abs(strainRate3DCell- strainRate3DCellSolution)
+ block % state % time_levs(1) % state % divTensor3DCellDiff % array = abs(divTensor3DCell- divTensor3DCellSolution)
+
+ print '(/,a,3es20.10)', 'max abs diff strainRate3DCell1',maxval(abs(strainRate3DCell(1,:,minCell:maxCell) - strainRate3DCellSolution(1,:,minCell:maxCell)))
+ print '(a,3es20.10)', 'max abs diff strainRate3DCell2',maxval(abs(strainRate3DCell(2,:,minCell:maxCell) - strainRate3DCellSolution(2,:,minCell:maxCell)))
+ print '(a,3es20.10)', 'max abs diff strainRate3DCell4',maxval(abs(strainRate3DCell(4,:,minCell:maxCell) - strainRate3DCellSolution(4,:,minCell:maxCell)))
+ print '(a,3es20.10)', 'max abs diff divTensor3DCell ',maxval(abs(divTensor3DCell(:,:,minCell:maxCell) - divTensor3DCellSolution(:,:,minCell:maxCell)))
+
block => block % next
end do
Modified: branches/ocean_projects/tensor_operations/src/operators/mpas_vector_operations.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/mpas_vector_operations.F 2013-04-19 17:48:27 UTC (rev 2778)
+++ branches/ocean_projects/tensor_operations/src/operators/mpas_vector_operations.F 2013-04-19 17:49:14 UTC (rev 2779)
@@ -35,6 +35,7 @@
public :: mpas_initialize_vectors, &
mpas_unit_vec_in_r3, &
mpas_cross_product_in_r3, &
+ mpas_tangential_velocity, &
mpas_vector_3DCell_to_2DEdge
!--------------------------------------------------------------------
@@ -166,6 +167,73 @@
!***********************************************************************
!
+! routine mpas_tangential_velocity
+!
+!> \brief compute tangential velocity at an edge from the normal velocities
+!> \author Mark Petersen
+!> \date 17 April 2013
+!> \version SVN:$Id$
+!> \details
+!> Compute tangential velocity at an edge from the normal velocities on the
+!> edges of the two neighboring cells.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_tangential_velocity(grid, normalVelocity, tangentialVelocity)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ normalVelocity !< Input: Horizontal velocity normal to edge
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ tangentialVelocity !< Output: Horizontal velocity tangent to edge
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges,i,k,nVertLevels, eoe
+ integer, dimension(:), pointer :: nEdgesOnEdge
+ integer, dimension(:,:), pointer :: edgesOnEdge
+
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
+
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ weightsOnEdge => grid % weightsOnEdge % array
+
+ do iEdge = 1,nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ do k=1,nVertLevels
+ tangentialVelocity(k,iEdge) = tangentialVelocity(k,iEdge) + weightsOnEdge(i,iEdge) * normalVelocity(k, eoe)
+ end do
+ end do
+ end do
+
+ end subroutine mpas_tangential_velocity!}}}
+
+!***********************************************************************
+!
! routine mpas_initialize_vectors
!
!> \brief MPAS RBF interpolation initialization routine
@@ -260,6 +328,7 @@
else ! this is not a boundary cell
! the normal points from the cell 1 to cell2
+ ! mrp problem: on periodic domains, vectors on edges of domain point the wrong way.
edgeNormalVectors(1,iEdge) = xCell(cell2) - xCell(cell1)
edgeNormalVectors(2,iEdge) = yCell(cell2) - yCell(cell1)
edgeNormalVectors(3,iEdge) = zCell(cell2) - zCell(cell1)
</font>
</pre>