<p><b>mpetersen@lanl.gov</b> 2012-11-01 07:38:08 -0600 (Thu, 01 Nov 2012)</p><p>branch commit, strain_rate: Added test cases and solutions. Strain rate computation is tested and correct for all test functions. Divergence of stress tensor is not yet fully tested.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/strain_rate/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/strain_rate/src/core_ocean/Registry 2012-10-31 15:24:12 UTC (rev 2289)
+++ branches/ocean_projects/strain_rate/src/core_ocean/Registry 2012-11-01 13:38:08 UTC (rev 2290)
@@ -192,34 +192,34 @@
var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
-var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
+var persistent real referenceBottomDepth ( nVertLevels ) 0 ir referenceBottomDepth mesh - -
var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zstarWeight ( nVertLevels ) 0 - zstarWeight mesh - -
% Boundary conditions: read from input, saved in restart and written to output
-var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
-var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
-var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
-var persistent integer edgeMask ( nVertLevels nEdges ) 0 o edgeMask mesh - -
-var persistent integer vertexMask ( nVertLevels nVertices ) 0 o vertexMask mesh - -
-var persistent integer cellMask ( nVertLevels nCells ) 0 o cellMask mesh - -
-var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
-var persistent real temperatureRestore ( nCells ) 0 iro temperatureRestore mesh - -
-var persistent real salinityRestore ( nCells ) 0 iro salinityRestore mesh - -
+var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 ir boundaryEdge mesh - -
+var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 ir boundaryVertex mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 ir boundaryCell mesh - -
+var persistent integer edgeMask ( nVertLevels nEdges ) 0 - edgeMask mesh - -
+var persistent integer vertexMask ( nVertLevels nVertices ) 0 - vertexMask mesh - -
+var persistent integer cellMask ( nVertLevels nCells ) 0 - cellMask mesh - -
+var persistent real u_src ( nVertLevels nEdges ) 0 ir u_src mesh - -
+var persistent real temperatureRestore ( nCells ) 0 ir temperatureRestore mesh - -
+var persistent real salinityRestore ( nCells ) 0 ir salinityRestore mesh - -
% mrp trying to figure out why these do not appear
-var persistent real windStressMonthly ( nMonths nEdges ) 0 iro windStressMonthly mesh - -
-var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 iro temperatureRestoreMonthly mesh - -
-var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 iro salinityRestoreMonthly mesh - -
+var persistent real windStressMonthly ( nMonths nEdges ) 0 ir windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 ir temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 ir salinityRestoreMonthly mesh - -
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
-var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
-var persistent real rho ( nVertLevels nCells Time ) 2 iro rho state - -
-var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
-var persistent real salinity ( nVertLevels nCells Time ) 2 iro salinity state tracers dynamics
-var persistent real tracer1 ( nVertLevels nCells Time ) 2 iro tracer1 state tracers testing
+var persistent real h ( nVertLevels nCells Time ) 2 ir h state - -
+var persistent real rho ( nVertLevels nCells Time ) 2 ir rho state - -
+var persistent real temperature ( nVertLevels nCells Time ) 2 ir temperature state tracers dynamics
+var persistent real salinity ( nVertLevels nCells Time ) 2 ir salinity state tracers dynamics
+var persistent real tracer1 ( nVertLevels nCells Time ) 2 ir tracer1 state tracers testing
% Tendency variables: neither read nor written to any files
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
@@ -231,7 +231,7 @@
% state variables for Split Explicit timesplitting
var persistent real uBtr ( nEdges Time ) 2 - uBtr state - -
-var persistent real ssh ( nCells Time ) 2 o ssh state - -
+var persistent real ssh ( nCells Time ) 2 - ssh state - -
var persistent real uBtrSubcycle ( nEdges Time ) 2 - uBtrSubcycle state - -
var persistent real sshSubcycle ( nCells Time ) 2 - sshSubcycle state - -
var persistent real FBtr ( nEdges Time ) 2 - FBtr state - -
@@ -246,32 +246,32 @@
var persistent real uBolusGMX ( nVertLevels nEdges Time ) 2 - uBolusGMX state - -
var persistent real uBolusGMY ( nVertLevels nEdges Time ) 2 - uBolusGMY state - -
var persistent real uBolusGMZ ( nVertLevels nEdges Time ) 2 - uBolusGMZ state - -
-var persistent real uBolusGMZonal ( nVertLevels nEdges Time ) 2 o uBolusGMZonal state - -
-var persistent real uBolusGMMeridional ( nVertLevels nEdges Time ) 2 o uBolusGMMeridional state - -
+var persistent real uBolusGMZonal ( nVertLevels nEdges Time ) 2 - uBolusGMZonal state - -
+var persistent real uBolusGMMeridional ( nVertLevels nEdges Time ) 2 - uBolusGMMeridional state - -
var persistent real hEddyFlux ( nVertLevels nEdges Time ) 2 - hEddyFlux state - -
var persistent real h_kappa ( nVertLevels nEdges Time ) 2 - h_kappa state - -
var persistent real h_kappa_q ( nVertLevels nEdges Time ) 2 - h_kappa_q state - -
-var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
+var persistent real divergence ( nVertLevels nCells Time ) 2 - divergence state - -
+var persistent real vorticity ( nVertLevels nVertices Time ) 2 - vorticity state - -
var persistent real Vor_edge ( nVertLevels nEdges Time ) 2 - Vor_edge state - -
var persistent real h_edge ( nVertLevels nEdges Time ) 2 - h_edge state - -
var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
-var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real kev ( nVertLevels nVertices Time ) 2 o kev state - -
-var persistent real kevc ( nVertLevels nCells Time ) 2 o kevc state - -
+var persistent real ke ( nVertLevels nCells Time ) 2 - ke state - -
+var persistent real kev ( nVertLevels nVertices Time ) 2 - kev state - -
+var persistent real kevc ( nVertLevels nCells Time ) 2 - kevc state - -
var persistent real ke_edge ( nVertLevels nEdges Time ) 2 - ke_edge state - -
var persistent real Vor_vertex ( nVertLevels nVertices Time ) 2 - Vor_vertex state - -
-var persistent real Vor_cell ( nVertLevels nCells Time ) 2 o Vor_cell state - -
+var persistent real Vor_cell ( nVertLevels nCells Time ) 2 - Vor_cell state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 - uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 - uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 - uReconstructZ state - -
-var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
-var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
+var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 - uReconstructZonal state - -
+var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 - uReconstructMeridional state - -
var persistent real uSrcReconstructX ( nVertLevels nCells Time ) 2 - uSrcReconstructX state - -
var persistent real uSrcReconstructY ( nVertLevels nCells Time ) 2 - uSrcReconstructY state - -
var persistent real uSrcReconstructZ ( nVertLevels nCells Time ) 2 - uSrcReconstructZ state - -
-var persistent real uSrcReconstructZonal ( nVertLevels nCells Time ) 2 o uSrcReconstructZonal state - -
-var persistent real uSrcReconstructMeridional ( nVertLevels nCells Time ) 2 o uSrcReconstructMeridional state - -
+var persistent real uSrcReconstructZonal ( nVertLevels nCells Time ) 2 - uSrcReconstructZonal state - -
+var persistent real uSrcReconstructMeridional ( nVertLevels nCells Time ) 2 - uSrcReconstructMeridional state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
@@ -279,6 +279,12 @@
var persistent real strainRateNormal ( R3 nVertLevels nCells Time ) 2 o strainRateNormal state - -
var persistent real strainRateShear ( R3 nVertLevels nVertices Time ) 2 o strainRateShear state - -
var persistent real divStressTensor ( nVertLevels nEdges Time ) 2 o divStressTensor state - -
+var persistent real strainRateNormalSolution ( R3 nVertLevels nCells Time ) 2 o strainRateNormalSolution state - -
+var persistent real strainRateShearSolution ( R3 nVertLevels nVertices Time ) 2 o strainRateShearSolution state - -
+var persistent real divStressTensorSolution ( nVertLevels nEdges Time ) 2 o divStressTensorSolution state - -
+var persistent real strainRateNormalDiff ( R3 nVertLevels nCells Time ) 2 o strainRateNormalDiff state - -
+var persistent real strainRateShearDiff ( R3 nVertLevels nVertices Time ) 2 o strainRateShearDiff state - -
+var persistent real divStressTensorDiff ( nVertLevels nEdges Time ) 2 o divStressTensorDiff state - -
% Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -302,11 +308,11 @@
var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 - vertDiffTopOfCell diagnostics - -
var persistent real nAccumulate ( Time ) 2 o nAccumulate state - -
-var persistent real acc_ssh ( nCells Time ) 2 o acc_ssh state - -
-var persistent real acc_sshVar ( nCells Time ) 2 o acc_sshVar state - -
-var persistent real acc_uReconstructZonal ( nVertLevels nCells Time ) 2 o acc_uReconstructZonal state - -
-var persistent real acc_uReconstructMeridional ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridional state - -
-var persistent real acc_uReconstructZonalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructZonalVar state - -
-var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
-var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
-var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+var persistent real acc_ssh ( nCells Time ) 2 - acc_ssh state - -
+var persistent real acc_sshVar ( nCells Time ) 2 - acc_sshVar state - -
+var persistent real acc_uReconstructZonal ( nVertLevels nCells Time ) 2 - acc_uReconstructZonal state - -
+var persistent real acc_uReconstructMeridional ( nVertLevels nCells Time ) 2 - acc_uReconstructMeridional state - -
+var persistent real acc_uReconstructZonalVar ( nVertLevels nCells Time ) 2 - acc_uReconstructZonalVar state - -
+var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 - acc_uReconstructMeridionalVar state - -
+var persistent real acc_u ( nVertLevels nEdges Time ) 2 - acc_u state - -
+var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 - acc_uVar state - -
Modified: branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-31 15:24:12 UTC (rev 2289)
+++ branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F 2012-11-01 13:38:08 UTC (rev 2290)
@@ -116,7 +116,7 @@
print *, 'mpas_test_stress_tensor'
call mpas_test_stress_tensor(domain)
print *, 'mpas_test_stress_tensor complete'
- stop
+!!!!!!!! stop
! mrp added just to test stress tensor ops end
if (config_enforce_zstar_at_restart) then
Modified: branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F
===================================================================
--- branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F 2012-10-31 15:24:12 UTC (rev 2289)
+++ branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F 2012-11-01 13:38:08 UTC (rev 2290)
@@ -192,57 +192,53 @@
end do
end do
+! print '(a,100f8.3)', 'xCell ',grid % xCell % array(6:7)
+! print '(a,100f8.3)', 'yCell ',grid % yCell % array(6:7)
+! print '(a,100f8.3)', 'xEdge ',grid % xEdge % array(1:nEdges)
+! print '(a,100f8.3)', 'yEdge ',grid % yEdge % array(1:nEdges)
+! print '(a,100f8.3)', 'angleEdge ',angleEdge(1:nEdges)
+! print '(a,100f8.3)', 'u ',u(:,1:nEdges)
+! print '(a,100f8.3)', 'v ',v(:,1:nEdges)
+! print '(a,100f8.3)', 'eps Nor 1 ',strainRateNormal(1,:,6:7)
+! print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,6:7)
+! print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(3,:,6:7)
+! print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(10:11)
+! print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(10:11)
+! print '(a,100f8.3)', 'eps Shr 1 ',strainRateShear(1,:,10:11)
+! print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,10:11)
+! print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,10:11)
-
-
-
- print '(a,100f8.3)', 'xCell ',grid % xCell % array(6:7)
- print '(a,100f8.3)', 'yCell ',grid % yCell % array(6:7)
+if (1==2) then
+ print '(a,100f8.3)', 'xCell ',grid % xCell % array(1:nCells)
+ print '(a,100f8.3)', 'yCell ',grid % yCell % array(1:nCells)
print '(a,100f8.3)', 'xEdge ',grid % xEdge % array(1:nEdges)
print '(a,100f8.3)', 'yEdge ',grid % yEdge % array(1:nEdges)
print '(a,100f8.3)', 'angleEdge ',angleEdge(1:nEdges)
print '(a,100f8.3)', 'u ',u(:,1:nEdges)
print '(a,100f8.3)', 'v ',v(:,1:nEdges)
- print '(a,100f8.3)', 'eps Nor 1 ',strainRateNormal(1,:,6:7)
- print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,6:7)
- print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(3,:,6:7)
- print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(10:11)
- print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(10:11)
- print '(a,100f8.3)', 'eps Shr 1 ',strainRateShear(1,:,10:11)
- print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,10:11)
- print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,10:11)
+ print '(a,100f8.3)', 'eps Nor 1 ',strainRateNormal(1,:,1:nCells)
+ print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,1:nCells)
+ print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(3,:,1:nCells)
+ print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(1:nVertices)
+ print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 1 ',strainRateShear(1,:,1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,1:nVertices)
+endif
- ! print '(a,100f8.3)', 'xCell ',grid % xCell % array(1:nCells)
- ! print '(a,100f8.3)', 'yCell ',grid % yCell % array(1:nCells)
- ! print '(a,100f8.3)', 'xEdge ',grid % xEdge % array(1:nEdges)
- ! print '(a,100f8.3)', 'yEdge ',grid % yEdge % array(1:nEdges)
- ! print '(a,100f8.3)', 'angleEdge ',angleEdge(1:nEdges)
- ! print '(a,100f8.3)', 'u ',u(:,1:nEdges)
- ! print '(a,100f8.3)', 'v ',v(:,1:nEdges)
- ! print '(a,100f8.3)', 'eps Nor 1 ',strainRateNormal(1,:,1:nCells)
- ! print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,1:nCells)
- ! print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(3,:,1:nCells)
- ! print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(1:nVertices)
- ! print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(1:nVertices)
- ! print '(a,100f8.3)', 'eps Shr 1 ',strainRateShear(1,:,1:nVertices)
- ! print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,1:nVertices)
- ! print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,1:nVertices)
-
-
end subroutine mpas_strain_rate!}}}
!***********************************************************************
!
! routine mpas_div_stress_tensor
!
-!> \brief Computes normal and shear strain rate
+!> \brief Computes divergence of the stress tensor
!> \author Mark Petersen
!> \date 25 September 2012
!> \version SVN:$Id$
!> \details
-!> This routine computes the normal strain rate at the cell centers
-!> and the shear strain rate at vertices using the weak derivative.
+!> This routine computes the divergence of the stress tensor
!
!-----------------------------------------------------------------------
@@ -390,13 +386,17 @@
real (kind=RKIND) :: invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, &
cos1, sin1, cos2, sin2, cossin, uDvEdge, uDcEdge, eps11, eps22, eps12, &
- invdcEdge, invdvEdge, ux, uy, a, eoe, cn, cs, r, theta, rot
+ invdcEdge, invdvEdge, ux, uy, eoe, 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, xVertex, yVertex
real (kind=RKIND), dimension(:,:), pointer :: u, v, divStressTensor, weightsOnEdge
real (kind=RKIND), dimension(:,:,:), pointer :: strainRateNormal, strainRateShear
+! real (kind=RKIND), dimension(:,:,:), allocatable :: strainRateNormalSolution, strainRateShearSolution, divStressTensorSolution
+ real (kind=RKIND), dimension(:,:), pointer :: divStressTensorSolution
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRateNormalSolution, strainRateShearSolution
+ character(len=16) :: test_function_type
-
!-----------------------------------------------------------------
!
! call relevant routines for computing tendencies
@@ -446,31 +446,137 @@
strainRateNormal => block % state % time_levs(1) % state % strainRateNormal % array
strainRateShear => block % state % time_levs(1) % state % strainRateShear % array
divStressTensor => block % state % time_levs(1) % state % divStressTensor % array
+ strainRateNormalSolution => block % state % time_levs(1) % state % strainRateNormalSolution % array
+ strainRateShearSolution => block % state % time_levs(1) % state % strainRateShearSolution % array
+ divStressTensorSolution => block % state % time_levs(1) % state % divStressTensorSolution % array
+
+! allocate(strainRateNormalSolution(3,nVertLevels,nCells),strainRateShearSolution(3,nVertLevels,nVertices),divStressTensorSolution(3,nVertLevels,nEdges))
- cn = 1.5
- cs = 0.7
+ cn = 1.555
+ cs = 2.555
+ rot = 3.0 !pii/4.0 !0.0*pii/2
+! test_function_type = 'linear_x'
+! test_function_type = 'linear_y'
+! test_function_type = 'linear_arb_rot'
+ test_function_type = 'sin_arb_rot'
- do iEdge = 1,nEdges
- a = angleEdge(iEdge)
+ ld = maxval(xEdge) ! note: this will not work in parallel
+ pi2l = 2*pii/ld
+ g1 = cn*cos(rot) - cs*sin(rot)
+ g2 = cn*sin(rot) + cs*cos(rot)
- !ux = cn*xEdge(iEdge)
- !uy = cs*xEdge(iEdge)
- !ux = -cs*yEdge(iEdge)
- !uy = cn*yEdge(iEdge)
- r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
- theta = atan(yEdge(iEdge)/xEdge(iEdge))
- rot = 4*pii/5
- ux = r*cos(theta-rot)*(cn*cos(rot) - cs*sin(rot))
- uy = r*cos(theta-rot)*(cn*sin(rot) + cs*cos(rot))
- do k=1,nVertLevels
- !do k = 1,maxLevelEdgeTop(iEdge)
- u(k,iEdge) = ux*cos(a) + uy*sin(a)
- end do
- end do
+ if (test_function_type.eq.'linear_x') then
+ ! linear function in x
+ do iEdge = 1,nEdges
+ ux = cn*xEdge(iEdge)
+ uy = cs*xEdge(iEdge)
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+ enddo
+ strainRateNormalSolution(1,:,:) = cn
+ strainRateNormalSolution(2,:,:) = 0.0
+ strainRateNormalSolution(3,:,:) = 0.5*cs
+ strainRateShearSolution(1,:,:) = cn
+ strainRateShearSolution(2,:,:) = 0.0
+ strainRateShearSolution(3,:,:) = 0.5*cs
+
+ divStressTensorSolution = 0.0
+ elseif (test_function_type.eq.'linear_y') then
+ ! linear function in y
+ do iEdge = 1,nEdges
+ ux = -cs*yEdge(iEdge)
+ uy = cn*yEdge(iEdge)
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+ enddo
+
+ strainRateNormalSolution(1,:,:) = 0.0
+ strainRateNormalSolution(2,:,:) = cn
+ strainRateNormalSolution(3,:,:) = -0.5*cs
+
+ strainRateShearSolution(1,:,:) = 0.0
+ strainRateShearSolution(2,:,:) = cn
+ strainRateShearSolution(3,:,:) = -0.5*cs
+
+ divStressTensorSolution = 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)
+ ux = f*g1
+ uy = f*g2
+
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ enddo
+
+ strainRateNormalSolution(1,:,:) = cos(rot)*g1
+ strainRateNormalSolution(2,:,:) = sin(rot)*g2
+ strainRateNormalSolution(3,:,:) = 0.5*(cos(rot)*g2 + sin(rot)*g1)
+
+ strainRateShearSolution(1,:,:) = cos(rot)*g1
+ strainRateShearSolution(2,:,:) = sin(rot)*g2
+ strainRateShearSolution(3,:,:) = 0.5*(cos(rot)*g2 + sin(rot)*g1)
+
+ divStressTensorSolution = 0.0
+
+ 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))
+ ux = f*g1
+ uy = f*g2
+
+ ! solution for div(sigma)
+ ! normal to edge:
+ divSigma_x = -pi2l**2*f*(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+ ! tangent to edge:
+ divSigma_y = -pi2l**2*f*(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+! v(:,iEdge) = ux*sin(angleEdge(iEdge)) - uy*cos(angleEdge(iEdge))
+ divStressTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ 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))
+
+ strainRateNormalSolution(1,:,iCell) = pi2l*fcos*cos(rot)*g1
+ strainRateNormalSolution(2,:,iCell) = pi2l*fcos*sin(rot)*g2
+ strainRateNormalSolution(3,:,iCell) = pi2l*fcos*(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+ do iVertex = 1,nVertices
+ r = sqrt(xVertex(iVertex)**2 + yVertex(iVertex)**2)
+ theta = atan(yVertex(iVertex)/xVertex(iVertex))
+ fcos = cos(pi2l*r*cos(theta-rot))
+
+ strainRateShearSolution(1,:,iVertex) = pi2l*fcos*cos(rot)*g1
+ strainRateShearSolution(2,:,iVertex) = pi2l*fcos*sin(rot)*g2
+ strainRateShearSolution(3,:,iVertex) = 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
+
do iEdge = 1,nEdges
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
@@ -486,10 +592,12 @@
end do
+
call mpas_strain_rate(block % mesh, u,v, strainRateNormal, strainRateShear)
call mpas_div_stress_tensor(block % mesh, strainRateNormal, strainRateShear, divStressTensor)
+if (1==2) then
print *, 'nedges, size(u)',nEdges, size(u)
print '(a,3f20.10)', 'min max mean u ',minval(u(:,1:nEdges)),maxval(u(:,1:nEdges)),sum(u(:,1:nEdges))/nEdges/nVertLevels
print '(a,3f20.10)', 'min max mean v ',minval(v(:,1:nEdges)),maxval(v(:,1:nEdges)),sum(v(:,1:nEdges))/nEdges/nVertLevels
@@ -500,7 +608,23 @@
print '(a,3f20.10)', 'min max mean strainRateShear2 ',minval(strainRateShear(2,:,1:nVertices)),maxval(strainRateShear(2,:,1:nVertices)),sum(strainRateShear(2,:,1:nVertices))/nVertices/nVertLevels
print '(a,3f20.10)', 'min max mean strainRateShear3 ',minval(strainRateShear(3,:,1:nVertices)),maxval(strainRateShear(3,:,1:nVertices)),sum(strainRateShear(3,:,1:nVertices))/nVertices/nVertLevels
print '(a,3f20.10)', 'min max mean divStressTensor ',minval(divStressTensor(:,1:nEdges)),maxval(divStressTensor(:,1:nEdges)),sum(divStressTensor(:,1:nEdges))/nEdges/nVertLevels
+endif
+ block % state % time_levs(1) % state % strainRateNormalDiff % array = abs(strainRateNormal- strainRateNormalSolution)
+ block % state % time_levs(1) % state % strainRateShearDiff % array = abs(strainRateShear- strainRateShearSolution)
+ block % state % time_levs(1) % state % divStressTensorDiff % array = abs(divStressTensor- divStressTensorSolution)
+
+ print '(a,3f20.10)', 'max abs diff strainRateNormal1',maxval(abs(strainRateNormal(1,:,1:nCells) - strainRateNormalSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff strainRateNormal2',maxval(abs(strainRateNormal(2,:,1:nCells) - strainRateNormalSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff strainRateNormal3',maxval(abs(strainRateNormal(3,:,1:nCells) - strainRateNormalSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff strainRateShear1 ',maxval(abs(strainRateShear(1,:,1:nVertices) - strainRateShearSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff strainRateShear2 ',maxval(abs(strainRateShear(2,:,1:nVertices) - strainRateShearSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff strainRateShear3 ',maxval(abs(strainRateShear(3,:,1:nVertices) - strainRateShearSolution(1,:,1:nCells)))
+ print '(a,3f20.10)', 'max abs diff divStressTensor ',maxval(abs(divStressTensor(:,1:nEdges) - divStressTensorSolution(:,1:nEdges)))
+
+
+! deallocate(strainRateNormalSolution,strainRateShearSolution,divStressTensorSolution)
+
block => block % next
end do
@@ -530,3 +654,4 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
! vim: foldmethod=marker
+
</font>
</pre>