Major modification of almost all source files to use vector groups for various

variables. Not tested in all detail, but the standard features seem to work.
The version before all these changes was tagged with PRE_MULTI.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@126 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 153 insertions, 146 deletions

diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index a89a548..bb8e350 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -15,12 +15,13 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i
+  CCTK_INT :: i, l
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
 
   character(len=200) :: gen_interp
   CCTK_INT :: gen_interp_len
   character(len=7) :: gen_order
+  character(len=15) :: vname
 
   CCTK_INT, dimension(1) :: lsh
   CCTK_POINTER, dimension(3) :: interp_coords
@@ -70,16 +71,11 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   endif
 
   ! Find out how many interpolation points are located on this processor.
-  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" )
   if ( status .lt. 0 ) then
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
 
-  ! Set the pointers to the points to be interpolated to.
-  interp_coords(1) = CCTK_PointerTo(xg)
-  interp_coords(2) = CCTK_PointerTo(yg)
-  interp_coords(3) = CCTK_PointerTo(zg)
-
   ! Set the pointers to the output arrays.
   out_arrays(1) = CCTK_PointerTo(alpg)
   out_arrays(2) = CCTK_PointerTo(betaxg)
@@ -96,147 +92,158 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   out_arrays(13) = CCTK_PointerTo(dfzg)
   out_arrays(14) = CCTK_PointerTo(psig)
 
-  ! Set the indices to the input grid functions.
-  call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
-  call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
-  call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
-  call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
-  call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
-  call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
-  call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
-  call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
-  call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
-  call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
-  call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" )
-  call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" )
-
-  ! Check the metric type. At present physical and static_conformal are
-  ! supported.
-  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
-
-    ! Set the operand indices table entry, corresponding
-    ! to interpolation of admbase::alp (1), admbase::shift (3),
-    ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for
-    ! a total of 13 output arrays.
-    call Util_TableSetIntArray ( status, table_handle, 13, &
-                               op_indices(1:13), "operand_indices" )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
-    endif
-
-    ! Set the corresponding table entry for the operation codes.
-    call Util_TableSetIntArray ( status, table_handle, 13, &
-                               op_codes(1:13), "operation_codes" )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
-    endif
-
-    ! Call the interpolator.
-    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
-                                 table_handle, coord_system_handle, &
-                                 lsh(1), CCTK_VARIABLE_REAL, &
-                                 interp_coords, 11, in_arrays(1:11), &
-                                 13, out_types(1:13), out_arrays(1:13) )
-
-    if ( status .lt. 0 ) then
-      call CCTK_INFO ( 'Interpolation failed.' )
-    end if
-
-    ! For each point on this processor calculate the right hand side of the
-    ! characteristic evolution equation.
-    do i = 1, lsh(1)
-
-      ! calculate the square of the lapse.
-      alp2 = alpg(i)**2
-
-      ! Calculate the inverse of the 3-metric.
-      guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
-      guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
-      guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
-
-      idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz )
-
-      guxx = idetg * guxx
-      guxy = idetg * guxy
-      guxz = idetg * guxz
-
-      guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
-      guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
-      guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
+  ! Loop over the level sets
+  do l = 1, eh_number_level_sets
+
+    ! Set the pointers to the points to be interpolated to.
+    interp_coords(1) = CCTK_PointerTo(xg(:,l))
+    interp_coords(2) = CCTK_PointerTo(yg(:,l))
+    interp_coords(3) = CCTK_PointerTo(zg(:,l))
+
+    ! Set the indices to the input grid functions.
+    call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
+    call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
+    call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
+    call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
+    call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
+    call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
+    call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
+    call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
+    call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
+    call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
+    write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1,']'
+    call CCTK_VarIndex ( in_arrays(11), vname )
+    call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" )
+
+    ! Check the metric type. At present physical and static_conformal are
+    ! supported.
+    if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
+
+      ! Set the operand indices table entry, corresponding
+      ! to interpolation of admbase::alp (1), admbase::shift (3),
+      ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for
+      ! a total of 13 output arrays.
+      call Util_TableSetIntArray ( status, table_handle, 13, &
+                                 op_indices(1:13), "operand_indices" )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+      endif
+
+      ! Set the corresponding table entry for the operation codes.
+      call Util_TableSetIntArray ( status, table_handle, 13, &
+                                 op_codes(1:13), "operation_codes" )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+      endif
+
+      ! Call the interpolator.
+      call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                   table_handle, coord_system_handle, &
+                                   lsh(1), CCTK_VARIABLE_REAL, &
+                                   interp_coords, 11, in_arrays(1:11), &
+                                   13, out_types(1:13), out_arrays(1:13) )
+
+      if ( status .lt. 0 ) then
+        call CCTK_INFO ( 'Interpolation failed.' )
+      end if
+
+      ! For each point on this processor calculate the right hand side of the
+      ! characteristic evolution equation.
+      do i = 1, lsh(1)
+
+        ! calculate the square of the lapse.
+        alp2 = alpg(i)**2
+
+        ! Calculate the inverse of the 3-metric.
+        guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
+        guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
+        guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
+
+        idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz )
+
+        guxx = idetg * guxx
+        guxy = idetg * guxy
+        guxz = idetg * guxz
+
+        guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
+        guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
+        guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
       
-      ! Raise the index of the partial derivatives of f.
-      dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
-      dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
-      dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
-
-      ! Calculate the overall multiplication factor.
-      factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
-                                      dfuy * dfyg(i) + &
-                                      dfuz * dfzg(i) ) )
-
-      ! Finally obtain dx^i/dt.
-      dxg(i) = - betaxg(i) + ssign * factor * dfux
-      dyg(i) = - betayg(i) + ssign * factor * dfuy
-      dzg(i) = - betazg(i) + ssign * factor * dfuz
-    end do
-  else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
-    call Util_TableSetIntArray ( status, table_handle, 14, &
-                               op_indices, "operand_indices" )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
-    endif
-
-    call Util_TableSetIntArray ( status, table_handle, 14, &
-                               op_codes, "operation_codes" )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
-    endif
-
-    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
-                                 table_handle, coord_system_handle, &
-                                 lsh(1), CCTK_VARIABLE_REAL, &
-                                 interp_coords, 12, in_arrays, &
-                                 14, out_types, out_arrays )
-
-    if ( status .lt. 0 ) then
-      call CCTK_INFO ( 'Interpolation failed.' )
+        ! Raise the index of the partial derivatives of f.
+        dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
+        dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
+        dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
+
+        ! Calculate the overall multiplication factor.
+        factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
+                                        dfuy * dfyg(i) + &
+                                        dfuz * dfzg(i) ) )
+
+        ! Finally obtain dx^i/dt.
+        dxg(i,l) = - betaxg(i) + ssign * factor * dfux
+        dyg(i,l) = - betayg(i) + ssign * factor * dfuy
+        dzg(i,l) = - betazg(i) + ssign * factor * dfuz
+      end do
+    else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
+      call Util_TableSetIntArray ( status, table_handle, 14, &
+                                 op_indices, "operand_indices" )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+      endif
+
+      call Util_TableSetIntArray ( status, table_handle, 14, &
+                                 op_codes, "operation_codes" )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+      endif
+
+      call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                   table_handle, coord_system_handle, &
+                                   lsh(1), CCTK_VARIABLE_REAL, &
+                                   interp_coords, 12, in_arrays, &
+                                   14, out_types, out_arrays )
+
+      if ( status .lt. 0 ) then
+        call CCTK_INFO ( 'Interpolation failed.' )
+      end if
+
+      do i = 1, lsh(1)
+        alp2 = alpg(i)**2
+
+!       The inverse of psi^4
+        psi4 = one / psig(i)**4
+
+        guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
+        guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
+        guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
+
+!       The determinant divided by psi^4.
+        idetg = psi4 / ( gxxg(i) * guxx + &
+                         gxyg(i) * guxy + &
+                         gxzg(i) * guxz )
+
+!       The inverse metric. Since the determinant is already divided
+!       by psi^4, this gives the inverse of the physical metric.
+        guxx = idetg * guxx
+        guxy = idetg * guxy
+        guxz = idetg * guxz
+
+        guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
+        guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
+        guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
+      
+        dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
+        dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
+        dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
+        factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
+                                        dfuy * dfyg(i) + &
+                                        dfuz * dfzg(i) ) )
+        dxg(i,l) = - betaxg(i) + ssign * factor * dfux
+        dyg(i,l) = - betayg(i) + ssign * factor * dfuy
+        dzg(i,l) = - betazg(i) + ssign * factor * dfuz
+      end do
     end if
+  end do
 
-    do i = 1, lsh(1)
-      alp2 = alpg(i)**2
-
-!     The inverse of psi^4
-      psi4 = one / psig(i)**4
-
-      guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
-      guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
-      guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
-
-!     The determinant divided by psi^4.
-      idetg = psi4 / ( gxxg(i) * guxx + &
-                       gxyg(i) * guxy + &
-                       gxzg(i) * guxz )
-
-!     The inverse metric. Since the determinant is already divided
-!     by psi^4, this gives the inverse of the physical metric.
-      guxx = idetg * guxx
-      guxy = idetg * guxy
-      guxz = idetg * guxz
-
-      guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
-      guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
-      guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
-      
-      dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
-      dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
-      dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
-      factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
-                                      dfuy * dfyg(i) + &
-                                      dfuz * dfzg(i) ) )
-      dxg(i) = - betaxg(i) + ssign * factor * dfux
-      dyg(i) = - betayg(i) + ssign * factor * dfuy
-      dzg(i) = - betazg(i) + ssign * factor * dfuz
-    end do
-  end if
   return
 end subroutine EHFinder_Generator_Sources