Dont use initial guess origins from parameters when tracking from gridscalar

Problem: My horizon will appear sometimes during the simulation at some position. I don't know this position in advance, so I set origin_* and AHFinderDirect::initial_guesscoord_sphere*_center to zero initially in my par-file. However, I have a grid scalar which tracks the coordinate location where the horizon will eventually appear.
When the horizon finder starts to search for a horizon, it will first setup the
coordinate ellipsoid. This routine is executed *before* the new origin is set
from the grid scalar.
The tracking occurs in the routine Newton(...). The ellipsoid is set before
Newton(...) gets executed. Hence, the ellipsoid uses the value that got set via
parameters (which would be zero in my case).

Patch by Christian Reisswig



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@1573 f88db872-0e4f-0410-b76b-b9085cfa78c5

4 files changed, 188 insertions, 73 deletions

diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 26460d2..44ed611 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -83,6 +83,49 @@ void Newton_step(patch_system& ps,
 		 const struct verbose_info& verbose_info);
 	  }
 
+void track_origin(const cGH* const cctkGH, patch_system& ps, 
+                  struct AH_data* const AH_data_ptr, 
+                  const int hn, const bool print_algorithm_highlights);
+
+
+//******************************************************************************
+//
+// This function reads a coordinate origin from a grid scalar, 
+// and sets the patch system's origin to that new value
+//
+//
+void track_origin(const cGH* const cctkGH, patch_system& ps, 
+                  struct AH_data* const AH_data_ptr, 
+                  const int hn, const bool print_algorithm_highlights)
+{
+DECLARE_CCTK_ARGUMENTS;
+DECLARE_CCTK_PARAMETERS;
+
+if (AH_data_ptr->depends_on == 0) {
+        // move the origin as specified in the grid scalars
+        fp const * const ox =
+          static_cast<CCTK_REAL const *>
+          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_x[hn]));
+        assert (ox);
+        fp const * const oy =
+          static_cast<CCTK_REAL const *>
+          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_y[hn]));
+        assert (oy);
+        fp const * const oz =
+          static_cast<CCTK_REAL const *>
+          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_z[hn]));
+        assert (oz);
+        if (print_algorithm_highlights) {
+          std::cout << "AHF tracked position ox " << *ox << " oy " << *oy << " oz " << *oz << std::endl;
+        }
+        ps.origin_x (*ox);
+        ps.origin_y (*oy);
+        ps.origin_z (*oz);
+}
+
+}
+
+
 //******************************************************************************
 
 //
@@ -331,27 +374,8 @@ if (hs.has_genuine_horizons())
           ps.synchronize();
         }
       }
-      if (track_origin_from_grid_scalar[hn] && AH_data_ptr->depends_on == 0) {
-        // move the origin as specified in the grid scalars
-        fp const * const ox =
-          static_cast<CCTK_REAL const *>
-          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_x[hn]));
-        assert (ox);
-        fp const * const oy =
-          static_cast<CCTK_REAL const *>
-          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_y[hn]));
-        assert (oy);
-        fp const * const oz =
-          static_cast<CCTK_REAL const *>
-          (CCTK_VarDataPtr (cctkGH, 0, track_origin_source_z[hn]));
-        assert (oz);
-        if (verbose_info.print_algorithm_highlights) {
-          std::cout << "AHF tracked position ox " << *ox << " oy " << *oy << " oz " << *oz << std::endl;
-        }
-        patch_system& ps = *ps_ptr;
-        ps.origin_x (*ox);
-        ps.origin_y (*oy);
-        ps.origin_z (*oz);
+      if (track_origin_from_grid_scalar[hn]) {
+         track_origin(cctkGH, *ps_ptr, AH_data_ptr, hn, verbose_info.print_algorithm_highlights);
       }
       
       // modify the initial guess
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 6569ba8..c998ce9 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -431,6 +431,10 @@ void setup_initial_guess(patch_system& ps,
 enum initial_guess_method
   decode_initial_guess_method(const char initial_guess_method_string[]);
 
+
+void set_initial_guess_parameters(struct AH_data& AH_data, const int hn,
+                                  const fp origin_x, const fp origin_y, const fp origin_z);
+
 // Newton.cc
 // returns true for success, false for failure to converge
 void Newton(const cGH* GH,
@@ -447,6 +451,11 @@ void Newton(const cGH* GH,
 	    const struct verbose_info& verbose_info,
 	    struct iteration_status_buffers& isb);
 
+// Tracks coordinate origin
+void track_origin(const cGH* const cctkGH, patch_system& ps, 
+                  struct AH_data* const AH_data_ptr, 
+                  const int hn, const bool print_algorithm_highlights);
+
 // io.cc
 void input_gridfn(patch_system& ps, int unknown_gfn,
 		  const struct IO_info& IO_info, const char base_file_name[],
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 45cc780..114f825 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -255,8 +255,13 @@ IO_info.output_mean_curvature
                         if (verbose_info.print_algorithm_highlights) {
                           printf ("AHF find_horizons[%d] setup_initial_guess\n", hn);
                         }
+                        if (track_origin_from_grid_scalar[hn] && state.method == method__find_horizons) {
+                           track_origin(cctkGH, ps, &AH_data, hn, verbose_info.print_algorithm_highlights);
+                           set_initial_guess_parameters(AH_data, hn, 
+                                                        ps.origin_x(), ps.origin_y(), ps.origin_z());
+                        }
         		setup_initial_guess(ps,
-        				    AH_data.initial_guess_info,
+        		          	    AH_data.initial_guess_info,
         				    IO_info,
         				    hn, N_horizons, verbose_info);
         		if (active_flag && IO_info.output_initial_guess)
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index c8ce8cc..79f34b3 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -85,6 +85,133 @@ enum patch_system::patch_system_type
 			   fp origin_x, fp origin_y, fp origin_z);
 	  }
 
+void set_initial_guess_parameters(struct AH_data& AH_data, const int hn,
+                                  const fp ini_origin_x, const fp ini_origin_y, const fp ini_origin_z);
+
+
+
+
+
+//******************************************************************************
+
+//
+// This function is used to setup initial guesses either from 
+// the corresponding parameters, or from a custom ini_origin_*.
+// The latter is necessary when tracking with a grid scalar.
+//
+//
+
+void set_initial_guess_parameters(struct AH_data& AH_data, const int hn,
+                                  const fp ini_origin_x, const fp ini_origin_y, const fp ini_origin_z)
+{
+DECLARE_CCTK_PARAMETERS;
+
+                AH_data.initial_guess_info.method
+			= decode_initial_guess_method(initial_guess_method[hn]);
+                AH_data.initial_guess_info.reset_horizon_after_not_finding
+                        = reset_horizon_after_not_finding[hn];
+		// ... read from named file
+		AH_data.initial_guess_info.read_from_named_file_info.file_name
+			= initial_guess__read_from_named_file__file_name[hn];
+
+if (!track_origin_from_grid_scalar[hn]) {
+
+		// ... Kerr/Kerr
+		AH_data.initial_guess_info.Kerr_Kerr_info.x_posn
+			= initial_guess__Kerr_Kerr__x_posn[hn];
+		AH_data.initial_guess_info.Kerr_Kerr_info.y_posn
+			= initial_guess__Kerr_Kerr__y_posn[hn];
+		AH_data.initial_guess_info.Kerr_Kerr_info.z_posn
+			= initial_guess__Kerr_Kerr__z_posn[hn];
+		AH_data.initial_guess_info.Kerr_Kerr_info.mass
+			= initial_guess__Kerr_Kerr__mass[hn];
+		AH_data.initial_guess_info.Kerr_Kerr_info.spin
+			= initial_guess__Kerr_Kerr__spin[hn];
+		// ... Kerr/Kerr-Schild
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.x_posn
+			= initial_guess__Kerr_KerrSchild__x_posn[hn];
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.y_posn
+			= initial_guess__Kerr_KerrSchild__y_posn[hn];
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.z_posn
+			= initial_guess__Kerr_KerrSchild__z_posn[hn];
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.mass
+			= initial_guess__Kerr_KerrSchild__mass[hn];
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.spin
+			= initial_guess__Kerr_KerrSchild__spin[hn];
+		// ... coordinate sphere
+		AH_data.initial_guess_info.coord_sphere_info.x_center
+			= initial_guess__coord_sphere__x_center[hn];
+		AH_data.initial_guess_info.coord_sphere_info.y_center
+			= initial_guess__coord_sphere__y_center[hn];
+		AH_data.initial_guess_info.coord_sphere_info.z_center
+			= initial_guess__coord_sphere__z_center[hn];
+		AH_data.initial_guess_info.coord_sphere_info.radius
+			= initial_guess__coord_sphere__radius[hn];
+		// ... coordinate ellipsoid
+		AH_data.initial_guess_info.coord_ellipsoid_info.x_center
+			= initial_guess__coord_ellipsoid__x_center[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.y_center
+			= initial_guess__coord_ellipsoid__y_center[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.z_center
+			= initial_guess__coord_ellipsoid__z_center[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.x_radius
+			= initial_guess__coord_ellipsoid__x_radius[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.y_radius
+			= initial_guess__coord_ellipsoid__y_radius[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.z_radius
+			= initial_guess__coord_ellipsoid__z_radius[hn];
+
+} else {
+
+		// ... Kerr/Kerr
+		AH_data.initial_guess_info.Kerr_Kerr_info.x_posn
+			= ini_origin_x;
+		AH_data.initial_guess_info.Kerr_Kerr_info.y_posn
+			= ini_origin_y;
+		AH_data.initial_guess_info.Kerr_Kerr_info.z_posn
+			= ini_origin_z;
+		AH_data.initial_guess_info.Kerr_Kerr_info.mass
+			= initial_guess__Kerr_Kerr__mass[hn];
+		AH_data.initial_guess_info.Kerr_Kerr_info.spin
+			= initial_guess__Kerr_Kerr__spin[hn];
+		// ... Kerr/Kerr-Schild
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.x_posn
+			= ini_origin_x;
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.y_posn
+			= ini_origin_y;
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.z_posn
+			= ini_origin_z;
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.mass
+			= initial_guess__Kerr_KerrSchild__mass[hn];
+		AH_data.initial_guess_info.Kerr_KerrSchild_info.spin
+			= initial_guess__Kerr_KerrSchild__spin[hn];
+		// ... coordinate sphere
+		AH_data.initial_guess_info.coord_sphere_info.x_center
+			= ini_origin_x;
+		AH_data.initial_guess_info.coord_sphere_info.y_center
+			= ini_origin_y;
+		AH_data.initial_guess_info.coord_sphere_info.z_center
+			= ini_origin_z;
+		AH_data.initial_guess_info.coord_sphere_info.radius
+			= initial_guess__coord_sphere__radius[hn];
+		// ... coordinate ellipsoid
+		AH_data.initial_guess_info.coord_ellipsoid_info.x_center
+			= ini_origin_x;
+		AH_data.initial_guess_info.coord_ellipsoid_info.y_center
+			= ini_origin_y;
+		AH_data.initial_guess_info.coord_ellipsoid_info.z_center
+			= ini_origin_z;
+		AH_data.initial_guess_info.coord_ellipsoid_info.x_radius
+			= initial_guess__coord_ellipsoid__x_radius[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.y_radius
+			= initial_guess__coord_ellipsoid__y_radius[hn];
+		AH_data.initial_guess_info.coord_ellipsoid_info.z_radius
+			= initial_guess__coord_ellipsoid__z_radius[hn];
+
+}
+
+}
+
 //******************************************************************************
 
 //
@@ -591,57 +718,7 @@ if (strlen(surface_interpolator_name) > 0)
 		if (verbose_info.print_algorithm_details)
 		   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "      setting initial guess parameters etc");
-		AH_data.initial_guess_info.method
-			= decode_initial_guess_method(initial_guess_method[hn]);
-                AH_data.initial_guess_info.reset_horizon_after_not_finding
-                        = reset_horizon_after_not_finding[hn];
-		// ... read from named file
-		AH_data.initial_guess_info.read_from_named_file_info.file_name
-			= initial_guess__read_from_named_file__file_name[hn];
-		// ... Kerr/Kerr
-		AH_data.initial_guess_info.Kerr_Kerr_info.x_posn
-			= initial_guess__Kerr_Kerr__x_posn[hn];
-		AH_data.initial_guess_info.Kerr_Kerr_info.y_posn
-			= initial_guess__Kerr_Kerr__y_posn[hn];
-		AH_data.initial_guess_info.Kerr_Kerr_info.z_posn
-			= initial_guess__Kerr_Kerr__z_posn[hn];
-		AH_data.initial_guess_info.Kerr_Kerr_info.mass
-			= initial_guess__Kerr_Kerr__mass[hn];
-		AH_data.initial_guess_info.Kerr_Kerr_info.spin
-			= initial_guess__Kerr_Kerr__spin[hn];
-		// ... Kerr/Kerr-Schild
-		AH_data.initial_guess_info.Kerr_KerrSchild_info.x_posn
-			= initial_guess__Kerr_KerrSchild__x_posn[hn];
-		AH_data.initial_guess_info.Kerr_KerrSchild_info.y_posn
-			= initial_guess__Kerr_KerrSchild__y_posn[hn];
-		AH_data.initial_guess_info.Kerr_KerrSchild_info.z_posn
-			= initial_guess__Kerr_KerrSchild__z_posn[hn];
-		AH_data.initial_guess_info.Kerr_KerrSchild_info.mass
-			= initial_guess__Kerr_KerrSchild__mass[hn];
-		AH_data.initial_guess_info.Kerr_KerrSchild_info.spin
-			= initial_guess__Kerr_KerrSchild__spin[hn];
-		// ... coordinate sphere
-		AH_data.initial_guess_info.coord_sphere_info.x_center
-			= initial_guess__coord_sphere__x_center[hn];
-		AH_data.initial_guess_info.coord_sphere_info.y_center
-			= initial_guess__coord_sphere__y_center[hn];
-		AH_data.initial_guess_info.coord_sphere_info.z_center
-			= initial_guess__coord_sphere__z_center[hn];
-		AH_data.initial_guess_info.coord_sphere_info.radius
-			= initial_guess__coord_sphere__radius[hn];
-		// ... coordinate ellipsoid
-		AH_data.initial_guess_info.coord_ellipsoid_info.x_center
-			= initial_guess__coord_ellipsoid__x_center[hn];
-		AH_data.initial_guess_info.coord_ellipsoid_info.y_center
-			= initial_guess__coord_ellipsoid__y_center[hn];
-		AH_data.initial_guess_info.coord_ellipsoid_info.z_center
-			= initial_guess__coord_ellipsoid__z_center[hn];
-		AH_data.initial_guess_info.coord_ellipsoid_info.x_radius
-			= initial_guess__coord_ellipsoid__x_radius[hn];
-		AH_data.initial_guess_info.coord_ellipsoid_info.y_radius
-			= initial_guess__coord_ellipsoid__y_radius[hn];
-		AH_data.initial_guess_info.coord_ellipsoid_info.z_radius
-			= initial_guess__coord_ellipsoid__z_radius[hn];
+		set_initial_guess_parameters(AH_data, hn, /* irrelevant here; leave at zero */0, 0, 0);
 		}
 
 	AH_data.search_flag = false;