Changed to using vector groups of the level set function and related

grid functions. More changes are needed to make everything work as before
but all of the standard features seem to be working.


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

3 files changed, 316 insertions, 285 deletions

diff --git a/interface.ccl b/interface.ccl
index 09ad541..fd1f3b7 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -17,41 +17,47 @@ USES FUNCTION MoLRegisterEvolvedGroup
 
 public:
 
-CCTK_REAL level_set TYPE=GF TIMELEVELS=2
-{
-  f
-} "The scalar level set function that defines the null surface"
+# The scalar level set functions that defines the null surfaces.
+CCTK_REAL f[eh_number_level_sets] TYPE=GF TIMELEVELS=2
 
-CCTK_REAL slevel_set TYPE=GF TIMELEVELS=1
-{
-  sf
-} "Source for the level set function"
+# Sources for the level set functions.
+CCTK_REAL sf[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
-CCTK_REAL dlevel_set TYPE=GF TIMELEVELS=1
-{
-  dfx, dfy, dfz, dfsq
-} "Derivatives of the level set function"
+# X-derivatives of the level set functions.
+CCTK_REAL dfx[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
-CCTK_REAL ftmp_set TYPE=GF TIMELEVELS=1
-{
-  ftmp, sftmp
-} "temporary variables used in pde re-parametrization"
+# Y-derivatives of the level set functions.
+CCTK_REAL dfy[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
-CCTK_REAL level_bak TYPE=GF TIMELEVELS=1
-{
-  fbak
-} "Temporary placeholder for the level set during re-parametrization"
+# Z-derivatives of the level set functions.
+CCTK_REAL dfz[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
-CCTK_INT eh_mask_all TYPE=GF TIMELEVELS=1
-{
-  eh_mask, tm_mask
-} "Masks to define active cells"
+# Squares of the derivatives of the level set functions.
+CCTK_REAL dfsq[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
-CCTK_INT mask_bak TYPE=GF TIMELEVELS=1
-{
-  eh_mask_bak
-} "Temporary placeholder for the mask during re-parametrization"
+# Temporary variables used in pde re-parametrization.
+CCTK_REAL ftmp[eh_number_level_sets] TYPE=GF TIMELEVELS=1
 
+# Temporary variables for the right hand side used in pde re-parametrization.
+CCTK_REAL sftmp[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+# Temporary placeholder for the level set during re-parametrization.
+CCTK_REAL fbak[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+# Masks to define active cells.
+CCTK_INT eh_mask[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+# Temporary masks to define active cells.
+CCTK_INT tm_mask[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+# Temporary placeholders for the masks during re-parametrization.
+CCTK_INT eh_mask_bak[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+#CCTK_INT rep_mask[eh_number_level_sets] TYPE=GF TIMELEVELS=1
+
+CCTK_INT re_param_control_pde TYPE=SCALAR
+CCTK_INT re_param_control_approx TYPE=SCALAR
+ 
 CCTK_INT surface_index TYPE=GF TIMELEVELS=1
 {
   sc
@@ -59,6 +65,12 @@ CCTK_INT surface_index TYPE=GF TIMELEVELS=1
 
 CCTK_INT find_surface_status TYPE=SCALAR
 
+CCTK_INT levelset_integers TYPE=SCALAR
+{
+  levelset_counter
+  more_levelsets
+} "Integer variables used to loop over the level sets"
+
 CCTK_INT surface_integers TYPE=SCALAR
 {
   surface_counter
@@ -106,89 +118,73 @@ CCTK_REAL integrate_tmp_array TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHO
   int_tmp
 }
 
-CCTK_INT rep_mask TYPE=GF TIMELEVELS=1
+#CCTK_INT rep_mask TYPE=GF TIMELEVELS=1
 
-CCTK_INT re_param_control_pde TYPE=SCALAR
-CCTK_INT re_param_control_approx TYPE=SCALAR
+# Storage for the event horizon areas. This is used to trigger output.
+CCTK_REAL eh_area[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
+# Temporary storage for the event horizon areas.
+CCTK_REAL eh_area2[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT 
 
-#CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+# Storage for the event horizon x-centroids. This is used to trigger output.
+CCTK_REAL eh_centroid_x[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_area[maximum_surface_number] TYPE=SCALAR
+# Storage for the event horizon y-centroids. This is used to trigger output.
+CCTK_REAL eh_centroid_y[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT #"The area of the event horizon(s)"
+# Storage for the event horizon z-centroids. This is used to trigger output.
+CCTK_REAL eh_centroid_z[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT #"The area of the event horizon(s)"
+# Temporary storage for the event horizon x-centroids.
+CCTK_REAL eh_centroid2_x[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_area2[maximum_surface_number] TYPE=SCALAR
-#CCTK_REAL eh_area3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+# Temporary storage for the event horizon y-centroids.
+CCTK_REAL eh_centroid2_y[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_centroid TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
-#{
-#  eh_centroid_x, eh_centroid_y, eh_centroid_z
-#} "The centroid of the event horizon(s)"
+# Temporary storage for the event horizon z-centroids.
+CCTK_REAL eh_centroid2_z[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_centroid TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
-{
-  eh_centroid_x, eh_centroid_y, eh_centroid_z
-} "The centroid of the event horizon(s)"
+# Storage for the event horizon equatorial circumferences.
+# This is used to trigger output.
+CCTK_REAL eh_circ_eq[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_centroid2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
-#{
-#  eh_centroid2_x, eh_centroid2_y, eh_centroid2_z
-#} "The centroid of the event horizon(s)"
+# Storage for the event horizon polar circumferences.
+#This is used to trigger output.
+CCTK_REAL eh_circ_pol[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_centroid2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
-{
-  eh_centroid2_x, eh_centroid2_y, eh_centroid2_z
-} "The centroid of the event horizon(s)"
+# Temporary storage for the event horizon equatorial circumferences.
+# This is used to trigger output.
+CCTK_REAL eh_circ_eq2[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_circumference TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
-#{
-#  eh_circ_eq, eh_circ_pol
-#} "The circumferences of the event horizon(s)"
+# Temporary storage for the event horizon polar circumferences.
+# This is used to trigger output.
+CCTK_REAL eh_circ_pol2[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_circumference TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
-{
-  eh_circ_eq, eh_circ_pol
-} "The circumferences of the event horizon(s)"
+# The X-position of the generators of the event horizons.
+CCTK_REAL xg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=2 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_circumference2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
-#{
-#  eh_circ_eq2, eh_circ_pol2
-#} "The circumferences of the event horizon(s)"
+# The Y-position of the generators of the event horizons.
+CCTK_REAL yg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=2 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
 
-CCTK_REAL eh_circumference2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
-{
-  eh_circ_eq2, eh_circ_pol2
-} "The circumferences of the event horizon(s)"
+# The Z-position of the generators of the event horizons.
+CCTK_REAL zg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=2 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
 
-#CCTK_REAL eh_centroid3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
-#{
-#  eh_centroid3_x, eh_centroid3_y, eh_centroid3_z
-#} "The centroid of the event horizon(s)"
+# The right hand side of the X-position of the generators.
+CCTK_REAL dxg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
 
-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
+# The right hand side of the Y-position of the generators.
+CCTK_REAL dyg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators 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"
+# The right hand side of the Z-position of the generators.
+CCTK_REAL dzg[eh_number_level_sets] TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT
 
-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
+CCTK_REAL generator_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
-}
+} "Arrays to hold the interpolated metric, gauge and level set data"
 
 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"
-
diff --git a/param.ccl b/param.ccl
index 8f62563..530cb71 100644
--- a/param.ccl
+++ b/param.ccl
@@ -3,6 +3,11 @@
 
 private:
 
+INT eh_number_level_sets "How many level set functions should we evolve"
+{
+  1:          :: "Larger than or equal to one"
+} 1
+ 
 KEYWORD eh_metric_type "Do we use numerical or analytic metric information"
 {
   "numerical" :: "Read in metric from numerical data"
@@ -25,14 +30,14 @@ BOOLEAN read_conformal_factor_once "Should the conformal factor only be read onc
 {
 } "yes"
 
-KEYWORD initial_f "Initial surface choice"
+KEYWORD initial_f[10] "Initial surface choice"
 {
   "sphere"     :: "spherical surface"
   "ellipsoid"  :: "ellipsoidal surface"
   "cassini"    :: "ovaloid of cassini"
 } "sphere"
 
-REAL initial_rad "Initial radius of surface"
+REAL initial_rad[10] "Initial radius of surface"
 {
 0.0:            :: "Positive please"
 } 1.0
@@ -126,7 +131,7 @@ KEYWORD upwind_type "Type of upwinding used in evolving the ehfinder equations"
   "intrinsic" :: "Use the values of f itself to determine upwind direction"
   "shift"     :: "Use the shift to determine upwind direction"
   "characteristic" :: "Use characteristic information"
-} "intrinsic"
+} "characteristic"
 
 KEYWORD surface_direction "Should we track outward or inward moving surfaces"
 {
@@ -311,7 +316,7 @@ 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:10		:: "Only evolve the level set functions"
 } 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
@@ -321,7 +326,5 @@ CCTK_INT EHFinder_Max_Evolved_Array_Size "The maximum size of evolved grid array
 
 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"
+  0:30		:: "Should be exactly zero or a multiple of three"
 } 3
-
-
diff --git a/schedule.ccl b/schedule.ccl
index e6b3fd4..b3bff78 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -3,41 +3,36 @@
 
 if (CCTK_Equals(mode,"analysis"))
 {
-  STORAGE: level_set[1]
+  STORAGE: f[1]
 }
 else
 {
-  STORAGE: level_set[2]
-  STORAGE: slevel_set
-  STORAGE: dlevel_set
-  STORAGE: eh_mask_all
-  STORAGE: rep_mask
-  STORAGE: ftmp_set
-  STORAGE: level_bak
-  STORAGE: mask_bak
+  STORAGE: f[2]
+  STORAGE: sf
+  STORAGE: ftmp
+  STORAGE: eh_mask
   STORAGE: surface_arrays
-  STORAGE: center_arrays
-#  STORAGE: eh_area, eh_centroid
-  STORAGE: eh_area2, eh_centroid2, eh_circumference2
+  STORAGE: eh_area2, eh_centroid2_x, eh_centroid2_y, eh_centroid2_x
+  STORAGE: eh_circ_eq2, eh_circ_pol2
   STORAGE: find_surface_status
-#  STORAGE: surface_index
+  STORAGE: levelset_integers
+  STORAGE: center_arrays
   if ( evolve_generators )
   {
-    STORAGE: generators[2]
-    STORAGE: generators_rhs
-    STORAGE: generators_arrays
+    STORAGE: xg[2], yg[2], zg[2]
+    STORAGE: dxg, dyg, dzg
   }
 }
 
-if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
-{
-  STORAGE: generator_gf
-}
+#if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
+#{
+#  STORAGE: generator_gf
+#}
 
 # Check for metric_state
 
-if (!CCTK_EQUALS(mode,"analysis"))
-{
+#if (!CCTK_EQUALS(mode,"analysis"))
+#{
 
 schedule EHFinder_ParamCheck at CCTK_PARAMCHECK
 {
@@ -90,141 +85,163 @@ if (use_mask)
 schedule EHFinder_Init at CCTK_INITIAL
 {
   LANG: Fortran
-  SYNC: level_set
+  SYNC: f
 } "Setup local variables"
 
 schedule EHFinder_Init_F at CCTK_INITIAL
 {
   LANG: Fortran
-  SYNC: level_set
+  SYNC: f
 } "Setup the initial surface"
 
 schedule EHFinder_InitWeights at CCTK_INITIAL
 {
-  LANG: FORTRAN
+  LANG: Fortran
 } "Setup weights for Simpson integration"
 
 schedule EHFinder_Init at CCTK_POST_RECOVER_VARIABLES
 {
   LANG: Fortran
-  SYNC: level_set
+  SYNC: f
 } "Setup local variables"
 
 schedule EHFinder_InitWeights at CCTK_POST_RECOVER_VARIABLES
 {
-  LANG: FORTRAN
+  LANG: Fortran
 } "Setup weights for Simpson integration"
 
-schedule GROUP EHFinder_Surfaces at CCTK_ANALYSIS
+schedule GROUP EHFinder_Level_Sets at CCTK_ANALYSIS
+{
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
+} "Loop over the level set functions"
+
+schedule EHFinder_LevelSetLoopInit in EHFinder_Level_Sets
+{
+  LANG: Fortran
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
+} "Initialise the loop counter over the level set functions"
+
+schedule GROUP EHFinder_Surfaces in EHFinder_Level_Sets while ehfinder::more_levelsets
 {
   STORAGE: surface_integers
   STORAGE: surface_reals
   STORAGE: surface_index
-#  STORAGE: eh_area2, eh_centroid2
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Count the number of surfaces and integrate over them"
 
 
 schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces
 {
   LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Initialise while loop control"
 
 schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces
 {
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Counting and mark surfaces"
 
 schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces
 {
   LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
   SYNC: surface_index
 } "Check if there are more surfaces"
 
 schedule EHFinder_MarkSurfaces in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points
 {
   LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
   SYNC: surface_index
 } "Mark points inside the current surface"
 
 schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces
 {
   LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Output info about found surfaces"
 
 schedule group EHFinder_Integration in EHFinder_Surfaces after EHFinder_InfoSurfaces while ehfinder::integrate_counter
 {
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Find and integrate over surfaces"
 
 schedule EHFinder_FindSurface in EHFinder_Integration
 {
   LANG: Fortran
   STORAGE: surface_tmp_arrays, surface_int_array
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
+  SYNC: surface_arrays
 } "Find Surface"
 
 schedule EHFinder_FindSurfaceElement in EHFinder_Integration after EHFinder_FindSurface
 {
   LANG: Fortran
   STORAGE: surface_tmp_arrays, interp_metric_arrays 
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Find Surface Area Element"
 
 schedule EHFinder_IntegrateArea in EHFinder_Integration after EHFinder_FindSurfaceElement
 {
   LANG: Fortran
   STORAGE: integrate_tmp_array 
-  TRIGGERS: eh_area, eh_centroid, eh_circumference
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Calculate area integrals"
 
 schedule EHFinder_IntegrateCentroid in EHFinder_Integration after EHFinder_IntegrateArea
 {
   LANG: Fortran
   STORAGE: surface_tmp_arrays, integrate_tmp_array
-  TRIGGERS: eh_centroid
+  TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
 } "Calculate centroid integrals"
 
 schedule EHFinder_IntegrateEquatorial in EHFinder_Integration after EHFinder_IntegrateArea
 {
   LANG: Fortran
   STORAGE: surface_tmp_arrays, integrate_tmp_array
-  TRIGGERS: eh_circumference
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Calculate circumferences"
 
-schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Surfaces
+schedule EHFinder_UpdateCounter in EHFinder_Surfaces
+{
+  LANG: Fortran
+  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_circ_eq, eh_circ_pol
+} "Update the loop variables"
+
+schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Level_Sets
 {
   LANG: Fortran
   STORAGE: eh_area
   TRIGGERS: eh_area
 } "Copy areas to output variable"
 
-schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Surfaces
+schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Level_Sets
 {
   LANG: Fortran
-  STORAGE: eh_centroid
-  TRIGGERS: eh_centroid
+  STORAGE: eh_centroid_x, eh_centroid_y, eh_centroid_z
+  TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
 } "Copy centroids to output variable"
 
-schedule EHFinder_CopyCircumference at CCTK_ANALYSIS after EHFinder_Surfaces
+schedule EHFinder_CopyCircumference at CCTK_ANALYSIS after EHFinder_Level_Sets
 {
   LANG: Fortran
-  STORAGE: eh_circumference
-  TRIGGERS: eh_circumference
+  STORAGE: eh_circ_eq, eh_circ_pol
+  TRIGGERS: eh_circ_eq, eh_circ_pol
 } "Copy circumferences to output variable"
 
-#schedule EHFinder_Integrate at CCTK_ANALYSIS
-#schedule EHFinder_Integrate at CCTK_ANALYSIS
-#{
-#  LANG: Fortran
-#  STORAGE: eh_area
-#  TRIGGERS: eh_area
-#} "Calculate integrals"
-
 # Read in the data used in reconstructing the 4-metric if necessary
 
 if (CCTK_Equals(eh_metric_type,"numerical"))
@@ -288,7 +305,7 @@ schedule EHFinder_SetSym at CCTK_BASEGRID
 } "Register the symmetries for the level set function"
 
 
-# Set the initial eh_mask to at the physical outer boundaries.
+# Set the initial eh_mask at the physical outer boundaries.
 
 schedule EHFinder_MaskInit at CCTK_POSTINITIAL after EHFinder_Init
 {
@@ -306,6 +323,9 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_Sources in MoL_CalcRHS
   {
     LANG: Fortran
+    STORAGE: dfx
+    STORAGE: dfy
+    STORAGE: dfz
   } "Calculate the source terms"
 
   if ( evolve_generators)
@@ -315,6 +335,7 @@ if (CCTK_Equals(mode,"normal"))
       schedule EHFinder_Generator_Sources in MoL_CalcRHS
       {
         LANG: Fortran
+        STORAGE: generator_arrays
       } "Calculate the source terms for the generator evolution"
     }
     if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
@@ -322,6 +343,7 @@ if (CCTK_Equals(mode,"normal"))
       schedule EHFinder_Generator_Sources2 in MoL_CalcRHS after EHFinder_Sources
       {
         LANG: Fortran
+        STORAGE: generator_arrays
       } "Calculate the source terms for the generator evolution"
     }
   }
@@ -334,150 +356,160 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_ApplySymF in EHFinder_PostStep
   {
     LANG: Fortran
-    SYNC: level_set
+    SYNC: f
   } "Apply symmetry boundaries"
 
 
-  # Set up the schedule group for re-parametrization
-
-  schedule GROUP EHFinder_ReParametrize at CCTK_POSTSTEP
-  {
-  } "Re-parametrize the level set function"
+# Set up the schedule group for re-parametrization
 
-  if (CCTK_Equals(re_param_method,"pde") || CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed"))
-  {
-    # Schedule the control routine, that initialises re_param_control_approx
-    # and re_param_control_pde if it is time for re-parametrization
+schedule GROUP EHFinder_ReParametrize at CCTK_POSTSTEP
+{
+  STORAGE: sftmp
+  STORAGE: fbak
+  STORAGE: eh_mask_bak
+  STORAGE: re_param_control_pde
+} "Re-parametrize the level set function"
 
-    schedule EHFinder_ReParametrizeControl in EHFinder_ReParametrize
-    {
-      LANG: Fortran
-    } "Initializes the re-parametrization control"
-  }
+if (CCTK_Equals(re_param_method,"pde") || CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed"))
+{
+  # Schedule the control routine, that initialises re_param_control_approx
+  # and re_param_control_pde if it is time for re-parametrization
 
-  if (CCTK_Equals(re_param_method,"pde") ||  CCTK_Equals(re_param_method,"mixed"))
+  schedule EHFinder_ReParametrizeControl in EHFinder_ReParametrize
   {
-    if (CCTK_Equals(re_param_int_method,"euler"))
-    {
-
-      # Set up the schedule group for euler re-parametrization using 
-      # re_param_control_pde to control the exit from the while loop.
-
-      schedule GROUP Euler_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_pde
-      {
-      } "Schedule group for Euler re-parametrization"
-
-      # Schedule the routine that does one Euler step 
-
-      schedule EHFinder_ReParametrizeEuler in Euler_ReParametrize 
-      {
-        LANG: Fortran
-#        SYNC: level_set
-      } "Euler scheme"
+    LANG: Fortran
+  } "Initializes the re-parametrization control"
+}
 
-      # Then synchronize the level set and apply symmetry boundary conditions.
-     
-      schedule  EHFinder_ApplySymF in Euler_ReParametrize after EHFinder_ReParametrizeEuler
-      {
-        LANG: Fortran
-        SYNC: level_set
-      } "Apply symmetry boundaries and sync"
+if (CCTK_Equals(re_param_method,"pde") ||  CCTK_Equals(re_param_method,"mixed"))
+{
+  if (CCTK_Equals(re_param_int_method,"euler"))
+  {
 
-      schedule EHFinder_ReParametrize_Check in EHFinder_ReParametrize AFTER Euler_ReParametrize
-      {
-        LANG: Fortran
-      } "Check to see if re-parametrization has to be undone"
-    }
+    # Set up the schedule group for euler re-parametrization using 
+    # re_param_control_pde to control the exit from the while loop.
 
-    if (CCTK_Equals(re_param_int_method,"rk2"))
+    schedule GROUP Euler_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_pde
     {
+    } "Schedule group for Euler re-parametrization"
 
-      # Set up the schedule group for rk2 re-parametrization using 
-      # re_param_control_pde to control the exit from the while loop.
-
-      schedule GROUP RK2_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control_pde
-      {
-      } "Schedule group for RK2 re-parametrization"
-
-      # Schedule the routine that does the first RK2 step
-
-      schedule EHFinder_ReParametrizeRK2_1 in RK2_ReParametrize 
-      {
-        LANG: Fortran
-      } "Euler scheme"
-
-      # Then synchronize the level set and apply symmetry boundary conditions.
-     
-      schedule  EHFinder_ApplySymF AS Sym_RK2_1 in RK2_ReParametrize after EHFinder_ReParametrizeRK2_1
-      {
-        LANG: Fortran
-        SYNC: level_set
-      } "Apply symmetry boundaries and sync"
-
-      # Schedule the routine that does the second RK2 step
+    # Schedule the routine that does one Euler step 
 
-      schedule EHFinder_ReParametrizeRK2_2 in RK2_ReParametrize after Sym_RK2_1
-      {
-        LANG: Fortran
-      } "Euler scheme"
+    schedule EHFinder_ReParametrizeEuler in Euler_ReParametrize 
+    {
+      LANG: Fortran
+      STORAGE: dfx
+      STORAGE: dfy
+      STORAGE: dfz
+      STORAGE: dfsq
+      SYNC: f
+    } "Euler scheme"
+
+     # Then synchronize the level set and apply symmetry boundary conditions.
+   
+    schedule  EHFinder_ApplySymF in Euler_ReParametrize after EHFinder_ReParametrizeEuler
+    {
+      LANG: Fortran
+      SYNC: f
+    } "Apply symmetry boundaries and sync"
 
-      # Then synchronize the level set and apply symmetry boundary conditions.
-     
-      schedule  EHFinder_ApplySymF AS Sym_RK2_2 in RK2_ReParametrize after EHFinder_ReParametrizeRK2_2
-      {
-        LANG: Fortran
-        SYNC: level_set
-      } "Apply symmetry boundaries and sync"
-    }
+    schedule EHFinder_ReParametrize_Check in EHFinder_ReParametrize AFTER Euler_ReParametrize
+    {
+      LANG: Fortran
+    } "Check to see if re-parametrization has to be undone"
   }
+#
+#    if (CCTK_Equals(re_param_int_method,"rk2"))
+#    {
+#
+#      # Set up the schedule group for rk2 re-parametrization using 
+#      # re_param_control_pde to control the exit from the while loop.
+#
+#      schedule GROUP RK2_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control_pde
+#      {
+#      } "Schedule group for RK2 re-parametrization"
+#
+#      # Schedule the routine that does the first RK2 step
+#
+#      schedule EHFinder_ReParametrizeRK2_1 in RK2_ReParametrize 
+#      {
+#        LANG: Fortran
+#      } "RK2 scheme step 1"
+#
+#      # Then synchronize the level set and apply symmetry boundary conditions.
+#     
+#      schedule  EHFinder_ApplySymF AS Sym_RK2_1 in RK2_ReParametrize after EHFinder_ReParametrizeRK2_1
+#      {
+#        LANG: Fortran
+#        SYNC: f
+#      } "Apply symmetry boundaries and sync"
+#
+#      # Schedule the routine that does the second RK2 step
+#
+#      schedule EHFinder_ReParametrizeRK2_2 in RK2_ReParametrize after Sym_RK2_1
+#      {
+#        LANG: Fortran
+#      } "RK2 scheme step 2"
+#
+#      # Then synchronize the level set and apply symmetry boundary conditions.
+#     
+#      schedule  EHFinder_ApplySymF AS Sym_RK2_2 in RK2_ReParametrize after EHFinder_ReParametrizeRK2_2
+#      {
+#        LANG: Fortran
+#        SYNC: f
+#      } "Apply symmetry boundaries and sync"
+#    }
+#  }
 
-  if (CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed"))
-  {
+#  if (CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed"))
+#  {
 
-    # Schedule a routine to update the points closest to the zero-level surface.
+#    # Schedule a routine to update the points closest to the zero-level surface.
 
-    schedule EHFinder_ReParametrize5 in EHFinder_ReParametrize after EHFinder_ReParametrizeControl
-    {
-      LANG: Fortran
-    } "Update the points closest to the surface"
+#    schedule EHFinder_ReParametrize5 in EHFinder_ReParametrize after EHFinder_ReParametrizeControl
+#    {
+#      LANG: Fortran
+#    } "Update the points closest to the surface"
 
-   schedule EHFinder_ApplySymFRep in EHFinder_ReParametrize after EHFinder_ReParametrize5
-    {
-      LANG: Fortran
-      SYNC: level_set
-      SYNC: rep_mask
-    } "Apply symmetry boundaries and sync"
+#   schedule EHFinder_ApplySymFRep in EHFinder_ReParametrize after EHFinder_ReParametrize5
+#    {
+#      LANG: Fortran
+#      SYNC: level_set
+#      SYNC: rep_mask
+#    } "Apply symmetry boundaries and sync"
 
-    # Set up the schedule group for approx re-parametrization using 
-    # re_param_control_approx to control the exit from the while loop.
+#    # Set up the schedule group for approx re-parametrization using 
+#    # re_param_control_approx to control the exit from the while loop.
 
-    schedule GROUP Approx_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_approx
-    {
-    } "Schedule group for Approximate re-parametrization"
+#    schedule GROUP Approx_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_approx
+#    {
+#    } "Schedule group for Approximate re-parametrization"
 
-    schedule EHFinder_ReParametrize6 in Approx_ReParametrize                   
-    {
-      LANG: Fortran
-    } "Update cells with neighbours that have already been updated"
+#    schedule EHFinder_ReParametrize6 in Approx_ReParametrize                   
+#    {
+#      LANG: Fortran
+#    } "Update cells with neighbours that have already been updated"
 
-    schedule EHFinder_ReParametrize7 in Approx_ReParametrize after ReParametrize6
-    {
-      LANG: Fortran
-    } "Update rep_mask"
+#    schedule EHFinder_ReParametrize7 in Approx_ReParametrize after ReParametrize6
+#    {
+#      LANG: Fortran
+#    } "Update rep_mask"
 
-    schedule EHFinder_ApplySymFRep in Approx_ReParametrize after EHFinder_ReParametrize7
-    {
-      LANG: Fortran
-      SYNC: level_set
-      SYNC: rep_mask
-    } "Apply symmetry boundaries and sync"
-  }
+#    schedule EHFinder_ApplySymFRep in Approx_ReParametrize after EHFinder_ReParametrize7
+#    {
+#      LANG: Fortran
+#      SYNC: level_set
+#      SYNC: rep_mask
+#    } "Apply symmetry boundaries and sync"
+
+}
 
 
   # Schedule modifications of the mask in group EHFinder_SetMask
 
   schedule GROUP EHFinder_SetMask at CCTK_POSTSTEP after EHFinder_ReParametrize
   {
+    STORAGE: tm_mask
   } "Set the mask"
 
   # First schedule the routine that adds or removes cells from the mask
@@ -485,8 +517,8 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_SetMask1 in EHFinder_SetMask
   {
     LANG: Fortran
-    SYNC: level_set
-    SYNC: eh_mask_all
+    SYNC: f
+    SYNC: eh_mask, tm_mask
   } "Start modifying the mask"
 
   # Then apply the symmetry boundary conditions.
@@ -502,7 +534,7 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_SetMask2 in EHFinder_SetMask after EHFinder_ApplySymAll
   {
     LANG: Fortran
-    SYNC: eh_mask_all
+    SYNC: eh_mask, tm_mask
   } "Find excision boundaries"
 
   schedule  EHFinder_ApplySymMask as ASM1 in EHFinder_SetMask after EHFinder_SetMask2
@@ -513,7 +545,7 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_SetMask3 in EHFinder_SetMask after EHFinder_ASM1
   {
     LANG: Fortran
-    SYNC: eh_mask_all
+    SYNC: eh_mask, tm_mask
   } "Check to see if the mask needs to be modified"
 
   schedule  EHFinder_ApplySymMask as ASM2 in EHFinder_SetMask after EHFinder_SetMask3
@@ -524,7 +556,7 @@ if (CCTK_Equals(mode,"normal"))
   schedule EHFinder_SetMask2 as SM2 in EHFinder_SetMask after ASM2
   {
     LANG: Fortran
-    SYNC: eh_mask_all
+    SYNC: eh_mask, tm_mask
   } "Find excision boundaries"
 
   schedule  EHFinder_ApplySymMask as ASM3 in EHFinder_SetMask after SM2
@@ -534,27 +566,27 @@ if (CCTK_Equals(mode,"normal"))
  
 }
 
-  if ( evolve_generators)
+if ( evolve_generators)
+{
+  if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
   {
-    if (CCTK_Equals(generator_tracking_method,"interpolate_before"))
+    schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS
     {
-      schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS
-      {
-        LANG: Fortran
-      } "Calculate the source terms for the generator evolution"
-    }
-    if (CCTK_Equals(generator_tracking_method,"interpolate_after"))
+      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
     {
-      schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS after EHFinder_Sources
-      {
-        LANG: Fortran
-      } "Calculate the source terms for the generator evolution"
-    }
+      LANG: Fortran
+      STORAGE: generator_arrays
+    } "Calculate the source terms for the generator evolution"
   }
-
-
 }
 
+
 #schedule  EHFinder_ApplySym at CCTK_POSTINITIAL after EHFinder_MaskInit
 #{
 #  LANG: Fortran