<p><b>dwj07@fsu.edu</b> 2012-10-19 09:08:29 -0600 (Fri, 19 Oct 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Creating missing bit-reproducible version of del4 on momentum.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-10-19 14:51:58 UTC (rev 2223)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-10-19 15:08:29 UTC (rev 2224)
@@ -112,22 +112,22 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, cell1, cell2, vertex1, vertex2, k
+ integer :: iEdge, cell1, cell2, vertex1, vertex2, k, i
integer :: iCell, iVertex
- integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve
+ integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve, vertexDegree
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexBot, &
- maxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+ maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnVertex, edgesOnCell, edgeSignOnVertex, edgeSignOnCell
real (kind=RKIND) :: u_diffusion, invAreaCell1, invAreaCell2, invAreaTri1, &
- invAreaTri2, invDcEdge, invDvEdge, r_tmp, delsq_u
+ invAreaTri2, invDcEdge, invDvEdge, r_tmp
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
meshScalingDel4, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence, &
- delsq_circulation, delsq_vorticity
+ delsq_circulation, delsq_vorticity, delsq_u
err = 0
@@ -138,6 +138,8 @@
nEdgesSolve = grid % nEdgessolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
maxLevelCell => grid % maxLevelCell % array
@@ -149,46 +151,95 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ allocate(delsq_u(nVertLEvels, nEdges+1))
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_vorticity(nVertLevels, nVertices+1))
+ delsq_u(:,:) = 0.0
delsq_vorticity(:,:) = 0.0
delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
+ ! DWJ 10/19/12 Bit Reproducible Version
+ !Compute delsq_u
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
invDcEdge = 1.0 / dcEdge(iEdge)
invDvEdge = 1.0 / dvEdge(iEdge)
do k=1,maxLevelEdgeTop(iEdge)
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- delsq_u = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
+ delsq_u(k, iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
-viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
+ end do
+ end do
- ! vorticity using </font>
<font color="red">abla^2 u
- r_tmp = dcEdge(iEdge) * delsq_u
- delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
- delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
+ ! Compute delsq_vorticity
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexBot(iVertex)
+ delsq_vorticity(k, iVertex) = delsq_vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * dcEdge(iEdge) * delsq_u(k, iEdge) * invAreaTri1
+ end do
+ end do
+ end do
- ! Divergence using </font>
<font color="red">abla^2 u
- r_tmp = dvEdge(iEdge) * delsq_u
- delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
- delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
+ ! Compute delsq_divergence
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - edgeSignOnCell(i, iCell) * dvEdge(iEdge) * delsq_u(k, iEdge) * invAreaCell1
+ end do
end do
end do
+! ! DWJ 10/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+
+! vertex1 = verticesOnEdge(1,iEdge)
+! vertex2 = verticesOnEdge(2,iEdge)
+
+! invAreaTri1 = 1.0 / areaTriangle(vertex1)
+! invAreaTri2 = 1.0 / areaTriangle(vertex2)
+
+! invAreaCell1 = 1.0 / areaCell(cell1)
+! invAreaCell2 = 1.0 / areaCell(cell2)
+
+! invDcEdge = 1.0 / dcEdge(iEdge)
+! invDvEdge = 1.0 / dvEdge(iEdge)
+
+! do k=1,maxLevelEdgeTop(iEdge)
+! ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+! delsq_u(k, iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
+! -viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
+
+! ! vorticity using </font>
<font color="blue">abla^2 u
+! r_tmp = dcEdge(iEdge) * delsq_u(k, iEdge)
+! delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
+! delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
+
+! ! Divergence using </font>
<font color="blue">abla^2 u
+! r_tmp = dvEdge(iEdge) * delsq_u(k, iEdge)
+! delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
+! delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
+! end do
+! end do
+
! Compute - \kappa </font>
<font color="black">abla^4 u
! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
do iEdge=1,nEdgesSolve
@@ -209,6 +260,7 @@
end do
end do
+ deallocate(delsq_u)
deallocate(delsq_divergence)
deallocate(delsq_vorticity)
</font>
</pre>