<p><b>mpetersen@lanl.gov</b> 2012-11-01 17:04:28 -0600 (Thu, 01 Nov 2012)</p><p>branch commit, strain_rate: Added test cases and solutions for power functions u=x^n, y^n, and any arbitrary rotation. Strain rate and divergence of the stress tensor all agree between computed and analytic solution, with error due to grid resolution.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F
===================================================================
--- branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F 2012-11-01 18:36:57 UTC (rev 2291)
+++ branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F 2012-11-01 23:04:28 UTC (rev 2292)
@@ -379,7 +379,7 @@
type (block_type), pointer :: block
- integer :: nCells, nEdges, nVertices, nVertLevels, iCell, iEdge, iVertex, k, i
+ integer :: nCells, nEdges, nVertices, nVertLevels, iCell, iEdge, iVertex, k, i, p
integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnEdge
@@ -455,10 +455,14 @@
cn = 1.555
cs = 2.555
rot = 3.0 !pii/4.0 !0.0*pii/2
+ p = 3 ! power for polynomial function
! 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 = 'sin_arb_rot'
+! test_function_type = 'power_x'
+! test_function_type = 'power_y'
+ test_function_type = 'power_arb_rot'
ld = maxval(xEdge) ! note: this will not work in parallel
pi2l = 2*pii/ld
@@ -529,6 +533,101 @@
divStressTensorSolution = 0.0
+ elseif (test_function_type.eq.'power_x') then
+
+ ! power function in x: x^p
+ do iEdge = 1,nEdges
+ ux = cn*xEdge(iEdge)**p
+ uy = cs*xEdge(iEdge)**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = cn *p*(p-1)*xEdge(iEdge)**(p-2)
+ divSigma_y = 0.5*cs*p*(p-1)*xEdge(iEdge)**(p-2)
+
+ divStressTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ strainRateNormalSolution(1,:,iCell) = cn *p*xCell(iCell)**(p-1)
+ strainRateNormalSolution(2,:,iCell) = 0.0
+ strainRateNormalSolution(3,:,iCell) = 0.5*cs*p*xCell(iCell)**(p-1)
+ end do
+
+ do iVertex = 1,nVertices
+ strainRateShearSolution(1,:,iVertex) = cn *p*xVertex(iVertex)**(p-1)
+ strainRateShearSolution(2,:,iVertex) = 0.0
+ strainRateShearSolution(3,:,iVertex) = 0.5*cs*p*xVertex(iVertex)**(p-1)
+ end do
+
+ elseif (test_function_type.eq.'power_y') then
+ ! power function in y: y^n
+ do iEdge = 1,nEdges
+ ux = -cs*yEdge(iEdge)**p
+ uy = cn*yEdge(iEdge)**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = -0.5*cs*p*(p-1)*yEdge(iEdge)**(p-2)
+ divSigma_y = cn *p*(p-1)*yEdge(iEdge)**(p-2)
+
+ divStressTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ strainRateNormalSolution(1,:,iCell) = 0.0
+ strainRateNormalSolution(2,:,iCell) = cn *p*yCell(iCell)**(p-1)
+ strainRateNormalSolution(3,:,iCell) = -0.5*cs*p*yCell(iCell)**(p-1)
+ end do
+
+ do iVertex = 1,nVertices
+ strainRateShearSolution(1,:,iVertex) = 0.0
+ strainRateShearSolution(2,:,iVertex) = cn *p*yVertex(iVertex)**(p-1)
+ strainRateShearSolution(3,:,iVertex) = -0.5*cs*p*yVertex(iVertex)**(p-1)
+ 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)
+
+ ux = g1*f**p
+ uy = g2*f**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = p*(p-1)*f**(p-2) *(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+ divSigma_y = p*(p-1)*f**(p-2) *(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+
+ 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))
+ f = r*cos(theta-rot)
+
+ strainRateNormalSolution(1,:,iCell) = p *f**(p-1) *cos(rot)*g1
+ strainRateNormalSolution(2,:,iCell) = p *f**(p-1) *sin(rot)*g2
+ strainRateNormalSolution(3,:,iCell) = p *f**(p-1) *(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))
+ f = r*cos(theta-rot)
+
+ strainRateShearSolution(1,:,iVertex) = p *f**(p-1) *cos(rot)*g1
+ strainRateShearSolution(2,:,iVertex) = p *f**(p-1) *sin(rot)*g2
+ strainRateShearSolution(3,:,iVertex) = p *f**(p-1) *(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+
elseif (test_function_type.eq.'sin_arb_rot') then
! sin function, arbitrary rotation
@@ -536,19 +635,17 @@
! arbitrary rotation
r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ f = sin(pi2l*r*cos(theta-rot))
- f = sin(pi2l*r*cos(theta-rot))
ux = f*g1
uy = f*g2
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+! v(:,iEdge) = ux*sin(angleEdge(iEdge)) - uy*cos(angleEdge(iEdge))
! 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
</font>
</pre>