Changed to use the newest version of MoL. Used its support for evolving

grid arrays to introduce tracking of the generators of the horizon along
with the horizon itself. Still needs some work on setting up the initial
generators (now it only does one) but first it needs to be tested.


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

9 files changed, 276 insertions, 24 deletions

diff --git a/interface.ccl b/interface.ccl
index b061a46..ebb03d0 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -5,7 +5,11 @@ implements: ehfinder
 inherits: grid admbase coordgauge staticconformal spacemask boundary
 
 USES INCLUDE: Boundary.h
-USES INCLUDE: MoL.h
+#USES INCLUDE: MoL.h
+
+CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, \
+                                     CCTK_INT IN RHSIndex)
+USES FUNCTION MoLRegisterEvolved
 
 public:
 
@@ -163,3 +167,18 @@ CCTK_REAL eh_circumference2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_n
 
 CCTK_REAL eh_group_test[maximum_surface_number] TYPE=SCALAR
 CCTK_REAL eh_group_array_test[maximum_surface_number] TYPE=ARRAY DIM=1 SIZE=100 GHOSTSIZE=0 DISTRIB=DEFAULT
+
+CCTK_REAL generators TYPE=ARRAY DIM=1 TIMELEVELS=2 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
+{
+  xg, yg, zg
+} "The position of the generators of the event horizon"
+
+CCTK_REAL generators_rhs TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
+{
+  dxg, dyg, dzg
+}
+
+CCTK_REAL generators_arrays TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
+{
+  alpg, betaxg, betayg, betazg, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg, dfxg, dfyg, dfzg, psig
+}
diff --git a/param.ccl b/param.ccl
index 0eefed4..f54d5e8 100644
--- a/param.ccl
+++ b/param.ccl
@@ -198,19 +198,28 @@ BOOLEAN use_user_center "Should the user prescribed center be used"
 
 CCTK_REAL center_x "The x-coordinate of the center"
 {
-*:*              :: "Anything"
+  *:*            :: "Anything"
 } 0.0
 
 CCTK_REAL center_y "The y-coordinate of the center"
 {
-*:*              :: "Anything"
+  *:*            :: "Anything"
 } 0.0
 
 CCTK_REAL center_z "The z-coordinate of the center"
 {
-*:*              :: "Anything"
+  *:*            :: "Anything"
 } 0.0
 
+BOOLEAN evolve_generators "Should the generators be evolved"
+{
+} "no"
+
+CCTK_INT number_of_generators "How many generators should be evolved"
+{
+  1:*            :: "Postive please"
+} 1
+
 shares: grid
 
 USES KEYWORD domain
@@ -239,3 +248,28 @@ EXTENDS KEYWORD initial_shift
 shares: spacemask
 
 USES KEYWORD use_mask
+
+shares: MethodOfLines
+
+USES CCTK_INT MoL_Num_Evolved_Vars
+USES CCTK_INT MoL_Max_Evolved_Array_Size
+USES CCTK_INT MoL_Num_ArrayEvolved_Vars
+
+restricted:
+
+CCTK_INT EHFinder_MaxNumEvolvedVars "The maximum number of evolved variables used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Num_Evolved_Vars
+{
+  1:1		:: "Only evolve the level set function"
+} 1
+
+CCTK_INT EHFinder_Max_Evolved_Array_Size "The maximum size of evolved grid arrays used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Max_Evolved_Array_Size
+{
+  1:*		:: "The size of the generator grid arrays"
+} 1
+
+CCTK_INT EHFinder_Num_ArrayEvolved_Vars "The maximum number of evolved grid arrays used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Num_ArrayEvolved_Vars
+{
+  0:3		:: "Should be exactly zero or three"
+} 3
+
+
diff --git a/schedule.ccl b/schedule.ccl
index 3e5f227..c234ae5 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -21,6 +21,12 @@ else
   STORAGE: eh_area2, eh_centroid2, eh_circumference2
   STORAGE: find_surface_status
 #  STORAGE: surface_index
+  if ( evolve_generators )
+  {
+    STORAGE: generators[2]
+    STORAGE: generators_rhs
+    STORAGE: generators_arrays
+  }
 }
 
 # Check for metric_state
@@ -297,6 +303,14 @@ if (CCTK_Equals(mode,"normal"))
     LANG: Fortran
   } "Calculate the source terms"
 
+  if ( evolve_generators)
+  {
+    schedule EHFinder_Generator_Sources in MoL_CalcRHS
+    {
+      LANG: Fortran
+    } "Calculate the source terms for the generator evolution"
+  }
+
   schedule GROUP EHFinder_PostStep in MoL_PostStep
   {
   } "Schedule group for symmetry boundaries and syncing"
diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
new file mode 100644
index 0000000..924588f
--- /dev/null
+++ b/src/EHFinder_Generator_Sources.F90
@@ -0,0 +1,157 @@
+! Calculation of the sources for the level set function.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i
+  CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+  CCTK_INT, dimension(1) :: 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, dfux, dfuy, dfuz, factor
+  CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
+
+  out_types = CCTK_VARIABLE_REAL
+
+  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    end if
+!    print*,lsh 
+
+    call CCTK_InterpHandle ( interp_handle, &
+                                     "Lagrange polynomial interpolation" )
+    if ( interp_handle .lt. 0 ) then
+      call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+    end if
+    ! For now order=2 is hard wired.
+    call Util_TableCreateFromString ( table_handle, "order=3" )
+    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
+      call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
+    endif
+
+    interp_coords(1) = CCTK_PointerTo(xg)
+    interp_coords(2) = CCTK_PointerTo(yg)
+    interp_coords(3) = CCTK_PointerTo(zg)
+
+    out_arrays(1) = CCTK_PointerTo(alpg)
+    out_arrays(2) = CCTK_PointerTo(betaxg)
+    out_arrays(3) = CCTK_PointerTo(betayg)
+    out_arrays(4) = CCTK_PointerTo(betazg)
+    out_arrays(5) = CCTK_PointerTo(gxxg)
+    out_arrays(6) = CCTK_PointerTo(gxyg)
+    out_arrays(7) = CCTK_PointerTo(gxzg)
+    out_arrays(8) = CCTK_PointerTo(gyyg)
+    out_arrays(9) = CCTK_PointerTo(gyzg)
+    out_arrays(10) = CCTK_PointerTo(gzzg)
+    out_arrays(11) = CCTK_PointerTo(dfxg)
+    out_arrays(12) = CCTK_PointerTo(dfyg)
+    out_arrays(13) = CCTK_PointerTo(dfzg)
+
+    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 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
+
+    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 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
+    do i = 1, lsh(1)
+      alp2 = alpg(i)**2
+
+      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
+      
+      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) + factor * dfux
+      dyg(i) = - betayg(i) + factor * dfuy
+      dzg(i) = - betazg(i) + factor * dfuz
+!      print*, alpg
+!      print*
+!      print*, betaxg, betayg, betazg
+!      print*
+!      print*, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg
+!      print*
+!      print*, dfxg, dfyg, dfzg
+!      print*
+!      print*, alp2
+!      print*
+!      print*, dfux, dfuy, dfuz
+!      print*
+!      print*, factor
+!      print*
+!      print*, dxg(i), dyg(i), dzg(i)
+!      stop
+    end do
+  end if
+  return
+end subroutine EHFinder_Generator_Sources
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index cf700f3..ebb7cb7 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -41,6 +41,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
     f = sqrt( ( x - translate_x )**2 + &
               ( y - translate_y )**2 + &
               ( z - translate_z )**2 ) - initial_rad
+    if ( evolve_generators .gt. 0 ) then
+      xg(1) = initial_rad + translate_x
+      yg(1) = translate_y
+      zg(1) = translate_z
+    end if
   end if
 
   ! If an ellipsoid is requested...
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index 81745ea..48bd609 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -35,7 +35,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
                                                       CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL /)
   
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -212,7 +212,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   character(len=30) :: info_message
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -274,7 +274,7 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
   CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -382,7 +382,7 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS)
   CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha
   CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90
index 679c8fb..bccc870 100644
--- a/src/MoL_Init.F90
+++ b/src/MoL_Init.F90
@@ -10,6 +10,7 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
   implicit none
 
   CCTK_INT :: ierr, ierr_cum, varindex, rhsindex
+  CCTK_INT :: MoLRegisterEvolved
 
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
@@ -18,8 +19,28 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
   ierr_cum = 0
   call CCTK_VarIndex(varindex, 'ehfinder::f')
   call CCTK_VarIndex(rhsindex, 'ehfinder::sf')
-  call MoL_RegisterVar(ierr, varindex, rhsindex)
-  ierr_cum = ierr_cum + ierr
+!  call MoL_RegisterVar(ierr, varindex, rhsindex)
+!  ierr_cum = ierr_cum + ierr
+  ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+  if ( evolve_generators .gt. 0 ) then
+
+!    call CCTK_GroupIndex (varindex, 'ehfinder::generators')
+!    call CCTK_GroupIndex(rhsindex, 'ehfinder::generators_rhs')
+
+    call CCTK_VarIndex(varindex, 'ehfinder::xg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dxg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+    call CCTK_VarIndex(varindex, 'ehfinder::yg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dyg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+    call CCTK_VarIndex(varindex, 'ehfinder::zg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dzg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+  end if
+ 
 !  call CCTK_VarIndex(varindex, 'admbase::gxx')
 !  call MoL_RegisterPrimitive(ierr, varindex)
 !  ierr_cum = ierr_cum + ierr
diff --git a/src/make.code.defn b/src/make.code.defn
index c54b696..c8e8ce4 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -7,6 +7,7 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_Init.F90 \
        MoL_Init.F90 \
        EHFinder_Sources.F90 \
+       EHFinder_Generator_Sources.F90 \
        EHFinder_ReParametrize.F90 \
        EHFinder_Check.F90 \
        EHFinder_Sources2.F90 \
diff --git a/src/make.code.deps b/src/make.code.deps
index 1c09018..d8e8918 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -1,17 +1,18 @@
 # Dependencies file for thorn EHFinder.
 # $Header$
 
-EHFinder_mod.F90.o:		EHFinder_Constants.F90.o
-EHFinder_Init.F90.o:		EHFinder_mod.F90.o
-EHFinder_Sources.F90.o:		EHFinder_mod.F90.o
-EHFinder_ReParametrize.F90.o:	EHFinder_mod.F90.o
-EHFinder_Check.F90.o:		EHFinder_mod.F90.o
-EHFinder_Sources2.F90.o:	EHFinder_mod.F90.o
-EHFinder_Sources3.F90.o:	EHFinder_mod.F90.o
-EHFinder_SetMask.F90.o:		EHFinder_mod.F90.o
-EHFinder_SetSym.F90.o:		EHFinder_mod.F90.o
-EHFinder_ReadData.F90.o:	EHFinder_mod.F90.o
-EHFinder_Check.F90.o:		EHFinder_mod.F90.o
-EHFinder_Integrate.F90.o:	EHFinder_mod.F90.o
-EHFinder_Integrate2.F90.o:	EHFinder_mod.F90.o
-EHFinder_FindSurface.F90.o:	EHFinder_mod.F90.o
+EHFinder_mod.F90.o:			EHFinder_Constants.F90.o
+EHFinder_Init.F90.o:			EHFinder_mod.F90.o
+EHFinder_Sources.F90.o:			EHFinder_mod.F90.o
+EHFinder_Generator_Sources.F90.o:	EHFinder_mod.F90.o
+EHFinder_ReParametrize.F90.o:		EHFinder_mod.F90.o
+EHFinder_Check.F90.o:			EHFinder_mod.F90.o
+EHFinder_Sources2.F90.o:		EHFinder_mod.F90.o
+EHFinder_Sources3.F90.o:		EHFinder_mod.F90.o
+EHFinder_SetMask.F90.o:			EHFinder_mod.F90.o
+EHFinder_SetSym.F90.o:			EHFinder_mod.F90.o
+EHFinder_ReadData.F90.o:		EHFinder_mod.F90.o
+EHFinder_Check.F90.o:			EHFinder_mod.F90.o
+EHFinder_Integrate.F90.o:		EHFinder_mod.F90.o
+EHFinder_Integrate2.F90.o:		EHFinder_mod.F90.o
+EHFinder_FindSurface.F90.o:		EHFinder_mod.F90.o