bugfix in pack and unpack when using RRM

fixes #684 indexing bug when a corner node was coupled to more than 1 element (for RRM grids)
E3SM-Project · Feb 10, 2016 · a5bccbf · a5bccbf
1 parent 47bd009
commit a5bccbf
Show file tree

Hide file tree

Showing 5 changed files with 63 additions and 56 deletions.
diff --git a/components/homme/src/share/cube_mod.F90 b/components/homme/src/share/cube_mod.F90
@@ -2301,13 +2301,6 @@ subroutine CubeSetupEdgeIndex(Edge)
     integer,allocatable                :: backwardV(:), backwardP(:)
     integer                            :: i,ii
 
-    ii=Edge%tail_face
-
-    !map to correct location - for now all on same nbr side have same wgt, so take the first one
-    ii = Edge%tail%nbrs_ptr(ii)
-
-    np0 = Edge%tail%nbrs_wgt(ii)
-
     sFace = Edge%tail_face
     dFace = Edge%head_face
     ! Do not reverse the indices

diff --git a/components/homme/src/share/edge_mod.F90 b/components/homme/src/share/edge_mod.F90
@@ -604,7 +604,6 @@ subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
     integer :: i,k,ir,ll,iptr
 
     integer :: is,ie,in,iw
-    integer :: nce
 
     !call t_adj_detailf(+2)
     !call t_startf('edgeVpack')
@@ -688,12 +687,15 @@ subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
     endif
 
     !set the max_corner_elem
-    nce = max_corner_elem
 ! SWEST
     do ll=swest,swest+max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                edge%buf(nce*(kptr+k-1)+edge%putmap(ll,ielem)+1)=v(1  ,1 ,k)
+                iptr = (kptr+k-1)+edge%putmap(ll,ielem)+1
+                if (iptr > size(edge%buf)) then
+                    print *,'ERROR 1: ',size(edge%buf),iptr,edge%putmap(ll,ielem)
+                endif
+                edge%buf(iptr)=v(1  ,1 ,k)
             end do
         end if
     end do
@@ -702,7 +704,11 @@ subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                edge%buf(nce*(kptr+k-1)+edge%putmap(ll,ielem)+1)=v(np ,1 ,k)
+                iptr = (kptr+k-1)+edge%putmap(ll,ielem)+1
+                if (iptr > size(edge%buf)) then
+                    print *,'ERROR 2: ',size(edge%buf),iptr
+                endif
+                edge%buf(iptr)=v(np ,1 ,k)
 !                edge%buf(kptr+k,edge%putmap(ll,ielem)+1)=v(np ,1 ,k)
             end do
         end if
@@ -712,7 +718,11 @@ subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                edge%buf(nce*(kptr+k-1)+edge%putmap(ll,ielem)+1)=v(np ,np,k)
+                iptr = (kptr+k-1)+edge%putmap(ll,ielem)+1
+                if (iptr > size(edge%buf)) then
+                    print *,'ERROR 3: ',size(edge%buf),iptr
+                endif
+                edge%buf(iptr)=v(np ,np,k)
 !                edge%buf(kptr+k,edge%putmap(ll,ielem)+1)=v(np ,np,k)
             end do
         end if
@@ -722,7 +732,11 @@ subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                edge%buf(nce*(kptr+k-1)+edge%putmap(ll,ielem)+1)=v(1  ,np,k)
+                iptr = (kptr+k-1)+edge%putmap(ll,ielem)+1
+                if (iptr > size(edge%buf)) then
+                    print *,'ERROR 4: ',size(edge%buf),iptr
+                endif
+                edge%buf(iptr)=v(1  ,np,k)
 !                edge%buf(kptr+k,edge%putmap(ll,ielem)+1)=v(1  ,np,k)
             end do
         end if
@@ -748,7 +762,6 @@ subroutine edgeSpack(edge,v,vlyr,kptr,ielem)
     integer :: i,k,ir,ll,iptr
 
     integer :: is,ie,in,iw
-    integer :: nce
     real (kind=real_kind) :: tmp
 
 !pw call t_adj_detailf(+2)
@@ -769,13 +782,11 @@ subroutine edgeSpack(edge,v,vlyr,kptr,ielem)
        edge%buf(iptr+in+1)   = v(k) ! North
        edge%buf(iptr+iw+1)   = v(k) ! West
     enddo
-    !set the max_corner_elem
-    nce = max_corner_elem
 ! SWEST
     do ll=swest,swest+max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 edge%buf(iptr+edge%putmap(ll,ielem)+1)=v(k)
             end do
         end if
@@ -785,7 +796,7 @@ subroutine edgeSpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 edge%buf(iptr+edge%putmap(ll,ielem)+1)=v(k)
             end do
         end if
@@ -795,7 +806,7 @@ subroutine edgeSpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 edge%buf(iptr+edge%putmap(ll,ielem)+1)=v(k)
             end do
         end if
@@ -805,7 +816,7 @@ subroutine edgeSpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if (edge%putmap(ll,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 edge%buf(iptr+edge%putmap(ll,ielem)+1)=v(k)
             end do
         end if
@@ -972,7 +983,7 @@ subroutine edgeVunpack(edge,v,vlyr,kptr,ielem)
     !type (EdgeDescriptor_t)            :: desc
 
     ! Local
-    integer :: i,k,ll,iptr,nce
+    integer :: i,k,ll,iptr
     integer :: is,ie,in,iw
     integer :: ks,ke,kblock
     logical :: done
@@ -1014,11 +1025,10 @@ subroutine edgeVunpack(edge,v,vlyr,kptr,ielem)
     enddo
 
 ! SWEST
-    nce = max_corner_elem
     do ll=swest,swest+max_corner_elem-1
         if(edge%getmap(ll,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,1 ,k)=v(1 ,1 ,k)+edge%receive(nce*(kptr+k-1)+edge%getmap(ll,ielem)+1)
+                v(1  ,1 ,k)=v(1 ,1 ,k)+edge%receive((kptr+k-1)+edge%getmap(ll,ielem)+1)
             enddo
         endif
     end do
@@ -1027,7 +1037,7 @@ subroutine edgeVunpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
         if(edge%getmap(ll,ielem) /= -1) then 
             do k=1,vlyr
-                v(np ,1 ,k)=v(np,1 ,k)+edge%receive(nce*(kptr+k-1)+edge%getmap(ll,ielem)+1)
+                v(np ,1 ,k)=v(np,1 ,k)+edge%receive((kptr+k-1)+edge%getmap(ll,ielem)+1)
             enddo
         endif
     end do
@@ -1036,7 +1046,7 @@ subroutine edgeVunpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if(edge%getmap(ll,ielem) /= -1) then 
             do k=1,vlyr
-                v(np ,np,k)=v(np,np,k)+edge%receive(nce*(kptr+k-1)+edge%getmap(ll,ielem)+1)
+                v(np ,np,k)=v(np,np,k)+edge%receive((kptr+k-1)+edge%getmap(ll,ielem)+1)
             enddo
         endif
     end do
@@ -1045,7 +1055,7 @@ subroutine edgeVunpack(edge,v,vlyr,kptr,ielem)
     do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if(edge%getmap(ll,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,np,k)=v(1 ,np,k)+edge%receive(nce*(kptr+k-1)+edge%getmap(ll,ielem)+1)
+                v(1  ,np,k)=v(1 ,np,k)+edge%receive((kptr+k-1)+edge%getmap(ll,ielem)+1)
             enddo
         endif
     end do
@@ -1338,6 +1348,7 @@ subroutine edgeDGVunpack(edge,v,vlyr,kptr,ielem)
     end do
 
     nce = max_corner_elem
+!   this is probably broken.  nce should be 1?  MT 2016/2/9
     i = swest
     if(edge%getmap(i,ielem) /= -1) then
       do k=1,vlyr
@@ -1387,7 +1398,7 @@ subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
 
     ! Local
 
-    integer :: i,k,l,iptr,nce
+    integer :: i,k,l,iptr
     integer :: is,ie,in,iw
 
     threadsafe=.false.
@@ -1406,12 +1417,11 @@ subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
        end do
     end do
 
-    nce = max_corner_elem
 ! SWEST
     do l=swest,swest+max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,1 ,k)=MAX(v(1 ,1 ,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(1  ,1 ,k)=MAX(v(1 ,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1420,7 +1430,7 @@ subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+max_corner_elem,swest+2*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(np ,1 ,k)=MAX(v(np,1 ,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(np ,1 ,k)=MAX(v(np,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1429,7 +1439,7 @@ subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then
             do k=1,vlyr
-                v(np ,np,k)=MAX(v(np,np,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(np ,np,k)=MAX(v(np,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1438,7 +1448,7 @@ subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,np,k)=MAX(v(1 ,np,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(1  ,np,k)=MAX(v(1 ,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1458,7 +1468,7 @@ subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
 
     ! Local
 
-    integer :: i,k,l,iptr,nce
+    integer :: i,k,l,iptr
     integer :: is,ie,in,iw
 
 !pw call t_startf('edgeSunpack')
@@ -1473,12 +1483,11 @@ subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
        v(k) = MAX(v(k),edge%receive(iptr+is+1),edge%receive(iptr+ie+1),edge%receive(iptr+in+1),edge%receive(iptr+iw+1))
     end do
 
-    nce = max_corner_elem
 ! SWEST
     do l=swest,swest+max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MAX(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1488,7 +1497,7 @@ subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+max_corner_elem,swest+2*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MAX(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1498,7 +1507,7 @@ subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MAX(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1508,7 +1517,7 @@ subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
     do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MAX(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1530,7 +1539,7 @@ subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
 
     ! Local
 
-    integer :: i,k,l,iptr,nce
+    integer :: i,k,l,iptr
     integer :: is,ie,in,iw
 
 !pw call t_startf('edgeSunpack')
@@ -1545,12 +1554,11 @@ subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
        v(k) = MIN(v(k),edge%receive(iptr+is+1),edge%receive(iptr+ie+1),edge%receive(iptr+in+1),edge%receive(iptr+iw+1))
     end do
 
-    nce = max_corner_elem
 ! SWEST
     do l=swest,swest+max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MiN(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1560,7 +1568,7 @@ subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+max_corner_elem,swest+2*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MIN(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1570,7 +1578,7 @@ subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MIN(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1580,7 +1588,7 @@ subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                iptr = nce*(kptr+k-1)
+                iptr = (kptr+k-1)
                 v(k)=MIN(v(k),edge%receive(iptr+edge%getmap(l,ielem)+1))
             enddo
         endif
@@ -1602,7 +1610,7 @@ subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
 
     ! Local
 
-    integer :: i,k,l,iptr,nce
+    integer :: i,k,l,iptr
     integer :: is,ie,in,iw
 
     threadsafe=.false.
@@ -1622,11 +1630,10 @@ subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
     end do
 
 ! SWEST
-    nce = max_corner_elem
     do l=swest,swest+max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,1 ,k)=MIN(v(1 ,1 ,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(1  ,1 ,k)=MIN(v(1 ,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1635,7 +1642,7 @@ subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+max_corner_elem,swest+2*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(np ,1 ,k)=MIN(v(np,1 ,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(np ,1 ,k)=MIN(v(np,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1644,7 +1651,7 @@ subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(np ,np,k)=MIN(v(np,np,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(np ,np,k)=MIN(v(np,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do
@@ -1653,7 +1660,7 @@ subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
     do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
         if(edge%getmap(l,ielem) /= -1) then 
             do k=1,vlyr
-                v(1  ,np,k)=MIN(v(1 ,np,k),edge%receive(nce*(kptr+k-1)+edge%getmap(l,ielem)+1))
+                v(1  ,np,k)=MIN(v(1 ,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
             enddo
         endif
     end do

diff --git a/components/homme/src/share/unit_tests_mod.F90 b/components/homme/src/share/unit_tests_mod.F90
@@ -431,7 +431,8 @@ subroutine test_subcell_div_fluxes(elem,deriv,nets,nete)
 
       do i=1,intervals
       do j=1,intervals
-        t = ABS(test(i,j)-values(i,j))/MAX(ABS(test(i,j)),ABS(values(i,j)))
+!        t = ABS(test(i,j)-values(i,j))/MAX(ABS(test(i,j)),ABS(values(i,j)))
+        t = ABS(test(i,j)-values(i,j))/MAX(ABS(MAXVAL(test(:,:))),ABS(MAXVAL(values(:,:))))
         if (.0000001<t) then
           print *,__FILE__,__LINE__,ie,i,j,test(i,j),values(i,j),t
           success = .false.