First attempt at using 2d grid arrays to hold the generator information.

Something fishy is going on. Don't use yet.


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

6 files changed, 535 insertions, 71 deletions

diff --git a/interface.ccl b/interface.ccl
index f52afdc..9d887be 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -186,3 +186,35 @@ CCTK_REAL generator_gf TYPE=GF TIMELEVELS=1
 {
   xgf, ygf, zgf
 } "Temporary grid function used in calculating the right hand side of the generator evolution equation"
+
+# The following is for a 2-d distribution of generators on the surface.
+
+# The X-position of the generators of the event horizons.
+CCTK_REAL xg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+# The Y-position of the generators of the event horizons.
+CCTK_REAL yg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+# The Z-position of the generators of the event horizons.
+CCTK_REAL zg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+# The right hand side of the X-position of the generators.
+CCTK_REAL dxg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+# The right hand side of the Y-position of the generators.
+CCTK_REAL dyg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+# The right hand side of the Z-position of the generators.
+CCTK_REAL dzg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+
+
+CCTK_REAL generator_arrays2 TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT
+{
+  alpg2, betaxg2, betayg2, betazg2, gxxg2, gxyg2, gxzg2, gyyg2, gyzg2, gzzg2, dfxg2, dfyg2, dfzg2, psig2
+} "Arrays to hold the interpolated metric, gauge and level set data"
+
+CCTK_REAL generator_gf2 TYPE=GF TIMELEVELS=1
+{
+  xgf2, ygf2, zgf2
+} "Temporary grid function used in calculating the right hand side of the generator evolution equation"
+
diff --git a/param.ccl b/param.ccl
index b08888f..0007865 100644
--- a/param.ccl
+++ b/param.ccl
@@ -215,9 +215,20 @@ CCTK_INT number_of_generators "How many generators should be evolved"
   1:*            :: "Postive please"
 } 1
 
+CCTK_INT number_of_generators_theta "How many generators in the theta direction"
+{
+  1:*           :: "Positive please"
+} 1
+
+CCTK_INT number_of_generators_phi "How many generators in the phi direction"
+{
+  1:*           :: "Positive please"
+} 1
+
 KEYWORD generator_distribution "What initial distribution should be used"
 {
   "line"	:: "Put the generators on a line in the xz-plane"
+  "2D array"    :: "Put the generators on a surface with spherical topology"
 } "line"
 
 STRING surface_interpolator "What interpolator should be used to locate the surface"
diff --git a/schedule.ccl b/schedule.ccl
index 3e9f6ab..f522f5c 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -19,8 +19,16 @@ else
   STORAGE: center_arrays
   if ( evolve_generators )
   {
-    STORAGE: xg[2], yg[2], zg[2]
-    STORAGE: dxg, dyg, dzg
+    if (CCTK_Equals(generator_distribution,"line"))
+    {
+      STORAGE: xg[2], yg[2], zg[2]
+      STORAGE: dxg, dyg, dzg
+    }
+    if (CCTK_Equals(generator_distribution,"2D array"))
+    {
+      STORAGE: xg2[2], yg2[2], zg2[2]
+      STORAGE: dxg2, dyg2, dzg2
+    }
   }
 }
 
@@ -341,11 +349,22 @@ if (CCTK_Equals(mode,"normal"))
   {
     if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
     {
-      schedule EHFinder_Generator_Sources in MoL_CalcRHS
+      if (CCTK_Equals(generator_distribution,"line"))
       {
-        LANG: Fortran
-        STORAGE: generator_arrays
-      } "Calculate the source terms for the generator evolution"
+        schedule EHFinder_Generator_Sources in MoL_CalcRHS
+        {
+          LANG: Fortran
+          STORAGE: generator_arrays
+        } "Calculate the source terms for the generator evolution"
+      }
+      if (CCTK_Equals(generator_distribution,"2D array"))
+      {
+        schedule EHFinder_Generator_Sources_2D in MoL_CalcRHS
+        {
+          LANG: Fortran
+          STORAGE: generator_arrays2
+        } "Calculate the source terms for the 2D generator evolution"
+      }
     }
     if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
     {
@@ -353,7 +372,7 @@ if (CCTK_Equals(mode,"normal"))
       {
         LANG: Fortran
         STORAGE: generator_arrays
-      } "Calculate the source terms for the generator evolution"
+      } "Calculate the source terms2 for the generator evolution"
     }
   }
 
@@ -533,23 +552,23 @@ if (CCTK_Equals(mode,"normal"))
     LANG: Fortran
   } "Apply symmetry boundaries"
 
-  if ( evolve_generators)
-  {
-    if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
-    {
-      schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS
-      {
-        LANG: Fortran
-        STORAGE: generator_arrays
-      } "Calculate the source terms for the generator evolution"
-    }
-    if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
-    {
-      schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS
-      {
-        LANG: Fortran
-        STORAGE: generator_arrays
-      } "Calculate the source terms for the generator evolution"
-    }
-  }
+#  if ( evolve_generators)
+#  {
+#    if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
+#    {
+#      schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS
+#      {
+#        LANG: Fortran
+#        STORAGE: generator_arrays
+#      } "Calculate the source terms for the generator evolution"
+#    }
+#    if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
+#    {
+#      schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS
+#      {
+#        LANG: Fortran
+#        STORAGE: generator_arrays
+#      } "Calculate the source terms for the generator evolution"
+#    }
+#  }
 }
diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index 984d80f..3dee327 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -53,9 +53,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   call CCTK_InterpHandle ( interp_handle, gen_interp )
 
   if ( interp_handle .lt. 0 ) then
-    warn_message = 'Cannot get handle for interpolation. '
-    warn_message = trim(warn_message)//'Forgot to activate an implementation '
-    warn_message = trim(warn_message)//'providing interpolation operators?'
+    warn_message = 'Cannot get handle for interpolation.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing interpolation operators?'
     call CCTK_WARN( 0, trim(warn_message) )
   end if
 
@@ -73,9 +73,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
 ! Get the 3D coordinate system handle.
   call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
   if ( coord_system_handle .lt. 0) then
-    warn_message = 'Cannot get handle for cart3d coordinate system. '
-    warn_message = trim(warn_message)//'Forgot to activate an implementation '
-    warn_message = trim(warn_message)//'providing coordinates?'
+    warn_message = 'Cannot get handle for cart3d coordinate system.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing coordinates?'
     call CCTK_WARN( 0, trim(warn_message) )
   endif
 
@@ -194,6 +194,18 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
         dxg(i,l) = - betaxg(i) + ssign * factor * dfux
         dyg(i,l) = - betayg(i) + ssign * factor * dfuy
         dzg(i,l) = - betazg(i) + ssign * factor * dfuz
+        if ( i == 5 ) then
+          print*,'Debug : ', xg(i,l), yg(i,l), zg(i,l), &
+                             dxg(i,l), dyg(i,l), dzg(i,l)
+            print*,'lapse : ', alp2
+            print*,'shift : ', betaxg(i), betayg(i), betazg(i)
+            print*,'metric : ', gxxg(i), gxyg(i), gxzg(i), &
+                                gyyg(i), gyzg(i), gzzg(i)
+            print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz
+            print*,'df : ', dfxg(i), dfyg(i), dfzg(i)
+            print*,'dfu : ', dfux, dfuy, dfuz
+            print*,'factor : ', factor
+        end if
       end do
     else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
       call Util_TableSetIntArray ( status, table_handle, 14, &
@@ -260,3 +272,324 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
 
   return
 end subroutine EHFinder_Generator_Sources
+
+
+subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, l
+  CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+
+  character(len=200) :: gen_interp
+  character(len=128) :: warn_message
+  CCTK_INT :: gen_interp_len
+  character(len=7) :: gen_order
+  character(len=15) :: vname
+
+  CCTK_INT, dimension(2) :: lsh
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(14) :: out_arrays
+  CCTK_INT, dimension(12) :: in_arrays
+  CCTK_INT, dimension(14), parameter :: op_indices = (/ 0, 1, 2, 3, 4, &
+                                                        5, 6, 7, 8, 9, &
+                                                        10, 10, 10, 11 /), &
+                                        op_codes = (/ 0, 0, 0, 0, 0, &
+                                                      0, 0, 0, 0, 0, &
+                                                      1, 2, 3, 0 /)
+  CCTK_INT, dimension(14) :: out_types
+  CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
+  CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
+  CCTK_REAL, dimension(:,:), allocatable :: tmp, postmp
+  CCTK_INT, dimension(1) :: shape
+
+! Set the sign correctly depending on the surface direction.
+  if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
+  if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
+
+  out_types = CCTK_VARIABLE_REAL
+
+! Convert the generator_interpolator string parameter to a Fortran string.
+  call CCTK_FortranString ( gen_interp_len, generator_interpolator, &
+                                            gen_interp )
+
+! Get the corresponding interpolator handle.
+  call CCTK_InterpHandle ( interp_handle, gen_interp )
+
+  if ( interp_handle .lt. 0 ) then
+    warn_message = 'Cannot get handle for interpolation.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing interpolation operators?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  end if
+
+! Convert the interpolation order parameter to a Fortran string to be placed
+! in the interpolator table. Note that the order is assumed to contain only
+! 1 digit.
+  write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order
+
+! Create the table directly from the string.
+  call Util_TableCreateFromString ( table_handle, gen_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
+  end if
+
+! Get the 3D coordinate system handle.
+  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
+  if ( coord_system_handle .lt. 0) then
+    warn_message = 'Cannot get handle for cart3d coordinate system.'
+    warn_message = trim(warn_message)//' Forgot to activate an implementation'
+    warn_message = trim(warn_message)//' providing coordinates?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  endif
+
+! Find out how many interpolation points are located on this processor.
+  call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, 'ehfinder::xg2' )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
+  end if
+
+  shape(1) = lsh(1)*lsh(2)
+  allocate ( tmp(lsh(1)*lsh(2),14), postmp(lsh(1)*lsh(2),3) )
+
+  do l = 1, 14
+    out_arrays(l) = CCTK_PointerTo(tmp(:,l))
+  end do
+  out_arrays(1) = CCTK_PointerTo(alpg2)
+  out_arrays(2) = CCTK_PointerTo(betaxg2)
+  out_arrays(3) = CCTK_PointerTo(betayg2)
+  out_arrays(4) = CCTK_PointerTo(betazg2)
+  out_arrays(5) = CCTK_PointerTo(gxxg2)
+  out_arrays(6) = CCTK_PointerTo(gxyg2)
+  out_arrays(7) = CCTK_PointerTo(gxzg2)
+  out_arrays(8) = CCTK_PointerTo(gyyg2)
+  out_arrays(9) = CCTK_PointerTo(gyzg2)
+  out_arrays(10) = CCTK_PointerTo(gzzg2)
+  out_arrays(11) = CCTK_PointerTo(dfxg2)
+  out_arrays(12) = CCTK_PointerTo(dfyg2)
+  out_arrays(13) = CCTK_PointerTo(dfzg2)
+  out_arrays(14) = CCTK_PointerTo(psig2)
+
+!! Set the pointers to the output arrays.
+
+! Loop over the level sets
+  do l = 1, eh_number_level_sets
+
+!   Set the pointers to the points to be interpolated to.
+!    postmp(:,1) = reshape(xg2(:,:,l),shape)
+!    postmp(:,2) = reshape(yg2(:,:,l),shape)
+!    postmp(:,3) = reshape(zg2(:,:,l),shape)
+!    interp_coords(1) = CCTK_PointerTo(postmp(:,1))
+!    interp_coords(2) = CCTK_PointerTo(postmp(:,2))
+!    interp_coords(3) = CCTK_PointerTo(postmp(:,3))
+
+    interp_coords(1) = CCTK_PointerTo(xg2(1,1,l))
+    interp_coords(2) = CCTK_PointerTo(yg2(1,1,l))
+    interp_coords(3) = CCTK_PointerTo(zg2(1,1,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
+        warn_message = 'Cannot set operand indices array in parameter table' 
+        call CCTK_WARN ( 0, trim(warn_message) )
+      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
+        warn_message = 'Cannot set operation codes array in parameter table'
+        call CCTK_WARN ( 0, trim(warn_message) )
+      endif
+
+!     Call the interpolator.
+      call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                   table_handle, coord_system_handle, &
+                                   lsh(1)*lsh(2), 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
+
+!      alpg2 = reshape(tmp(:,1),lsh)
+!      betaxg2 = reshape(tmp(:,2),lsh)
+!      betayg2 = reshape(tmp(:,3),lsh)
+!      betazg2 = reshape(tmp(:,4),lsh)
+!      gxxg2 = reshape(tmp(:,5),lsh)
+!      gxyg2 = reshape(tmp(:,6),lsh)
+!      gxzg2 = reshape(tmp(:,7),lsh)
+!      gyyg2 = reshape(tmp(:,8),lsh)
+!      gyzg2 = reshape(tmp(:,9),lsh)
+!      gzzg2 = reshape(tmp(:,10),lsh)
+!      dfxg2 = reshape(tmp(:,11),lsh)
+!      dfyg2 = reshape(tmp(:,12),lsh)
+!      dfzg2 = reshape(tmp(:,13),lsh)
+
+!     For each point on this processor calculate the right hand side of the
+!     characteristic evolution equation.
+      do j = 1, lsh(2)
+        do i = 1, lsh(1)
+
+!         calculate the square of the lapse.
+          alp2 = alpg2(i,j)**2
+
+!         Calculate the inverse of the 3-metric.
+          guxx = gyyg2(i,j) * gzzg2(i,j) - gyzg2(i,j)**2
+          guxy = gxzg2(i,j) * gyzg2(i,j) - gxyg2(i,j) * gzzg2(i,j)
+          guxz = gxyg2(i,j) * gyzg2(i,j) - gxzg2(i,j) * gyyg2(i,j)
+
+          idetg = one / ( gxxg2(i,j) * guxx + gxyg2(i,j) * guxy &
+                                            + gxzg2(i,j) * guxz )
+
+          guxx = idetg * guxx
+          guxy = idetg * guxy
+          guxz = idetg * guxz
+
+          guyy = ( gxxg2(i,j) * gzzg2(i,j) - gxzg2(i,j)**2 ) * idetg
+          guyz = ( gxyg2(i,j) * gxzg2(i,j) - gxxg2(i,j) * gyzg2(i,j) ) * idetg
+          guzz = ( gxxg2(i,j) * gyyg2(i,j) - gxyg2(i,j)**2 ) * idetg
+      
+!         Raise the index of the partial derivatives of f.
+          dfux = guxx * dfxg2(i,j) + guxy * dfyg2(i,j) + guxz * dfzg2(i,j)
+          dfuy = guxy * dfxg2(i,j) + guyy * dfyg2(i,j) + guyz * dfzg2(i,j)
+          dfuz = guxz * dfxg2(i,j) + guyz * dfyg2(i,j) + guzz * dfzg2(i,j)
+
+!         Calculate the overall multiplication factor.
+          factor = alp2 / sqrt ( alp2 * ( dfux * dfxg2(i,j) + &
+                                          dfuy * dfyg2(i,j) + &
+                                          dfuz * dfzg2(i,j) ) )
+
+!         Finally obtain dx^i/dt.
+          dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux
+          dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy
+          dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz
+          if ( j == 1 .and. i == 5 ) then
+              
+            print*,'debug2 : ', xg2(i,j,l), yg2(i,j,l), zg2(i,j,l), &
+                                dxg2(i,j,l), dyg2(i,j,l), dzg2(i,j,l)
+            print*,'lapse : ', alp2
+            print*,'shift : ', betaxg2(i,j), betayg2(i,j), betazg2(i,j)
+            print*,'metric : ', gxxg2(i,j), gxyg2(i,j), gxzg2(i,j), &
+                                gyyg2(i,j), gyzg2(i,j), gzzg2(i,j)
+            print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz
+            print*,'df : ', dfxg2(i,j), dfyg2(i,j), dfzg2(i,j)
+            print*,'dfu : ', dfux, dfuy, dfuz
+            print*,'factor : ', factor
+          end if
+        end do
+      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
+        warn_message = 'Cannot set operand indices array in parameter table'
+        call CCTK_WARN ( 0, trim(warn_message) )
+      endif
+
+      call Util_TableSetIntArray ( status, table_handle, 14, &
+                                 op_codes, 'operation_codes' )
+      if ( status .lt. 0 ) then
+        warn_message = 'Cannot set operation codes array in parameter table'
+        call CCTK_WARN ( 0, trim(warn_message) )
+      endif
+
+      call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                   table_handle, coord_system_handle, &
+                                   lsh(1)*lsh(2), 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
+
+!      alpg2 = reshape(tmp(:,1),lsh)
+!      betaxg2 = reshape(tmp(:,2),lsh)
+!      betayg2 = reshape(tmp(:,3),lsh)
+!      betazg2 = reshape(tmp(:,4),lsh)
+!      gxxg2 = reshape(tmp(:,5),lsh)
+!      gxyg2 = reshape(tmp(:,6),lsh)
+!      gxzg2 = reshape(tmp(:,7),lsh)
+!      gyyg2 = reshape(tmp(:,8),lsh)
+!      gyzg2 = reshape(tmp(:,9),lsh)
+!      gzzg2 = reshape(tmp(:,10),lsh)
+!      dfxg2 = reshape(tmp(:,11),lsh)
+!      dfyg2 = reshape(tmp(:,12),lsh)
+!      dfzg2 = reshape(tmp(:,13),lsh)
+!      psig2 = reshape(tmp(:,14),lsh)
+
+      do j = 1, lsh(2)
+        do i = 1, lsh(1)
+          alp2 = alpg2(i,j)**2
+
+!         The inverse of psi^4
+          psi4 = one / psig2(i,j)**4
+
+          guxx = gyyg2(i,j) * gzzg2(i,j) - gyzg2(i,j)**2
+          guxy = gxzg2(i,j) * gyzg2(i,j) - gxyg2(i,j) * gzzg2(i,j)
+          guxz = gxyg2(i,j) * gyzg2(i,j) - gxzg2(i,j) * gyyg2(i,j)
+
+!         The inverse of the determinant divided by psi^4.
+          idetg = psi4 / ( gxxg2(i,j) * guxx + &
+                           gxyg2(i,j) * guxy + &
+                           gxzg2(i,j) * 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 = ( gxxg2(i,j) * gzzg2(i,j) - gxzg2(i,j)**2 ) * idetg
+          guyz = ( gxyg2(i,j) * gxzg2(i,j) - gxxg2(i,j) * gyzg2(i,j) ) * idetg
+          guzz = ( gxxg2(i,j) * gyyg2(i,j) - gxyg2(i,j)**2 ) * idetg
+
+          dfux = guxx * dfxg2(i,j) + guxy * dfyg2(i,j) + guxz * dfzg2(i,j)
+          dfuy = guxy * dfxg2(i,j) + guyy * dfyg2(i,j) + guyz * dfzg2(i,j)
+          dfuz = guxz * dfxg2(i,j) + guyz * dfyg2(i,j) + guzz * dfzg2(i,j)
+          factor = alp2 / sqrt ( alp2 * ( dfux * dfxg2(i,j) + &
+                                          dfuy * dfyg2(i,j) + &
+                                          dfuz * dfzg2(i,j) ) )
+          dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux
+          dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy
+          dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz
+
+        end do
+      end do
+    end if
+  end do
+
+  deallocate ( tmp, postmp )
+  return
+end subroutine EHFinder_Generator_Sources_2D
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index b7f2a25..c9dd60c 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -22,7 +22,9 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
   CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc
   CCTK_REAL :: last_time
   CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el
+  CCTK_REAL :: phi, dphi, phimin, phimax
   CCTK_INT, dimension(1) :: lsh, lbnd
+  CCTK_INT, dimension(2) :: lsh2, lbnd2
 
 ! Get the size of the local grid.
   nx = cctk_lsh(1)
@@ -52,41 +54,76 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
   if ( evolve_generators .gt. 0 ) then
 
-    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'ehfinder::xg' )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' )
-    end if
-    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' )
-    if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, 'cannot get local size for generator arrays' )
-    end if
-
-    if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
-
-      if ( CCTK_EQUALS( domain, 'full' ) ) then
+    if ( CCTK_EQUALS( domain, 'full' ) ) then
+      thetamin = zero; thetamax = pi
+      phimin = zero; phimax = 2*pi
+    else if ( CCTK_EQUALS( domain, 'bitant') ) then
+      if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then
+        thetamin = zero; thetamax = half * pi
+        phimin = zero; phimax = 2*pi
+      else  if ( CCTK_EQUALS( bitant_plane, 'xz' ) ) then
         thetamin = zero; thetamax = pi
-      else if ( CCTK_EQUALS( domain, 'bitant') ) then
-        if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then
-          thetamin = zero; thetamax = half * pi
-        else 
-          thetamin = zero; thetamax = pi
-        end if
-      else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then
-        if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. &
-             CCTK_EQUALS( quadrant_direction, 'y' ) ) then
-          thetamin = zero; thetamax = half * pi
-        else
-          thetamin = zero; thetamax = pi
-        end if
-      else if ( CCTK_EQUALS( domain, 'octant' ) ) then
+        phimin = zero; phimax = pi
+      else
+        thetamin = zero; thetamax = pi
+        phimin = -half * pi; phimax = half * pi
+      end if
+    else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then
+      if ( CCTK_EQUALS( quadrant_direction, 'x' ) ) then
         thetamin = zero; thetamax = half * pi
+        phimin = zero; phimax = pi
+      else if ( CCTK_EQUALS( quadrant_direction, 'y' ) ) then
+        thetamin = zero; thetamax = half * pi
+        phimin = -half * pi; phimax = half * pi
+      else
+        thetamin = zero; thetamax = pi
+        phimin = zero; phimax = half * pi
       end if
+    else if ( CCTK_EQUALS( domain, 'octant' ) ) then
+      thetamin = zero; thetamax = half * pi
+      phimin = zero; phimax = half * pi
+    end if
                    
+    if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
+
+      call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'ehfinder::xg' )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' )
+      end if
+      call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, 'cannot get local size for generator arrays' )
+      end if
+
       if ( number_of_generators .eq. 1 ) then
         theta = half * ( thetamax - thetamin ) + thetamin
       else
         dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 )
       end if
+
+    else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then
+
+      call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd2, 'ehfinder::xg2' )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' )
+      end if
+      call CCTK_GrouplshGN ( status, cctkGH, 2, lsh2, 'ehfinder::xg2' )
+      if ( status .lt. 0 ) then
+        call CCTK_WARN ( 0, 'cannot get local size for generator arrays' )
+      end if
+
+      if ( number_of_generators_theta .eq. 1 ) then
+        theta = half * ( thetamax - thetamin ) + thetamin
+      else
+        dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta - 1 )
+      end if
+
+      if ( number_of_generators_phi .eq. 1 ) then
+        phi = half * ( phimax - phimin ) + phimin
+      else
+        dphi = ( phimax - phimin ) / ( number_of_generators_phi - 1 )
+      end if
+
     end if
   end if
 
@@ -103,10 +140,21 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       if ( evolve_generators .gt. 0 ) then
         if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
           do i = 1, lsh(1)
-             theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-             xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l)
-             yg(i,l) = translate_y(l)
-             zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l)
+            theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
+            xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l)
+            yg(i,l) = translate_y(l)
+            zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l)
+          end do
+        else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then
+          do j = 1, lsh2(2)
+            phi = phimin + dphi * ( j + lbnd2(2) - 1 )
+            do i = 1, lsh2(1)
+              theta = thetamin + dtheta * ( i + lbnd2(1) - 1 )
+              xg2(i,j,l) = initial_rad(l) * sin(theta) * cos(phi) + translate_x(l)
+              yg2(i,j,l) = initial_rad(l) * sin(theta) * sin(phi) + translate_x(l)
+              zg2(i,j,l) = initial_rad(l) * cos(theta) + translate_z(l)
+              print*,xg2(i,j,l),yg2(i,j,l),zg2(i,j,l)
+            end do
           end do
         end if
       end if
diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90
index 61a84e6..d01eec1 100644
--- a/src/MoL_Init.F90
+++ b/src/MoL_Init.F90
@@ -26,20 +26,41 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
 
   if ( evolve_generators .gt. 0 ) then
 
-    call CCTK_GroupIndex (varindex, 'ehfinder::xg')
-    call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg')
+    if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
 
-    ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+      call CCTK_GroupIndex (varindex, 'ehfinder::xg')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg')
 
-    call CCTK_GroupIndex (varindex, 'ehfinder::yg')
-    call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg')
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
 
-    ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+      call CCTK_GroupIndex (varindex, 'ehfinder::yg')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg')
 
-    call CCTK_GroupIndex (varindex, 'ehfinder::zg')
-    call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg')
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
 
-    ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+      call CCTK_GroupIndex (varindex, 'ehfinder::zg')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg')
+
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+
+    else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then
+
+      call CCTK_GroupIndex (varindex, 'ehfinder::xg2')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg2')
+
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+
+      call CCTK_GroupIndex (varindex, 'ehfinder::yg2')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg2')
+
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+
+      call CCTK_GroupIndex (varindex, 'ehfinder::zg2')
+      call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg2')
+
+      ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex)
+
+    end if
 
   end if