fix a bug which made us use the -1 ("past") time level every 2nd time

step of an evolution:

Previously the code found the data pointers for the Cactus g_ij and K_ij
grid arrays once in its setup phase, then used those all through the
evolution.  The problem is... the "current time level" data pointers
change at each time step, because Cactus rotates the time-level arrays
rather than copying all the (big) grid array data.

This commit changes the coee to just lookup the Cactus variable indices
in setup, then get a fresh set of data pointers each time we look for
an apparent horizon (i.e. at each time step).


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

1 files changed, 98 insertions, 11 deletions

diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 29eabfc..4781372 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -4,9 +4,14 @@
 // <<<access to persistent data>>>
 // <<<prototypes for functions local to this file>>>
 // AHFinderDirect_find_horizons - top-level driver to find apparent horizons
-// find_horizon - find a horizon
-// BH_diagnostics - compute BH diagnostics for a horizon
-// surface_integral - compute surface integral of a gridfn over the horizon
+///
+/// setup_Cactus_gridfn_data_ptrs - get all data pointers given variable indices
+/// Cactus_gridfn_data_ptr - get a single data pointer from a variable index
+///
+/// find_horizon - find a horizon
+///
+/// BH_diagnostics - compute BH diagnostics for a horizon
+/// surface_integral - compute surface integral of a gridfn over the horizon
 //
 
 #include <stdio.h>
@@ -57,6 +62,11 @@ extern struct state state;
 // ***** prototypes for functions local to this file
 //
 namespace {
+void setup_Cactus_gridfn_data_ptrs(const cGH *GH, struct cactus_grid_info& cgi);
+const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
+					const char gridfn_name[],
+					bool check_for_NULL = true);
+
 bool find_horizon(enum method method,
 		  const struct verbose_info& verbose_info, int timer_handle,
 		  struct IO_info& IO_info,
@@ -65,6 +75,7 @@ bool find_horizon(enum method method,
 		  struct cactus_grid_info& cgi, struct geometry_info& gi,
 		  patch_system& ps, Jacobian* Jac_ptr,
 		  int hn, int N_horizons);
+
 void BH_diagnostics(enum patch::integration_method surface_integral_method,
 		    const struct verbose_info& verbose_info,
 		    struct AH_info& AH_info);
@@ -91,6 +102,8 @@ const struct solver_info&  solver_info  = state.solver_info;
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
 
+IO_info.time_iteration = cctk_iteration;
+
 // what are the semantics of the Cactus gxx variables?
 if      (CCTK_Equals(metric_type, "physical"))
    then state.cgi.Cactus_conformal_metric = false;
@@ -100,15 +113,12 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
 		   metric_type);				/*NOTREACHED*/
 
-// if Cactus is using a conformal metric,
-//    we need to re-fetch the data pointer for the conformal factor,
-//    since the data array may have moved in memory since the last
-//    time we were called due to storage being switched off/on
-if (state.cgi.Cactus_conformal_metric)
-   then state.cgi.psi_data = Cactus_gridfn_data_ptr(cctkGH,
-						    "StaticConformal::psi");
+//
+// we need to re-fetch the Cactus data pointers at least each time step,
+// because they change each time Cactus rotates the time levels
+//
+setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 
-IO_info.time_iteration = cctk_iteration;
 	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
 	{
 	struct AH_info& AH_info = * state.AH_info_ptrs[hn];
@@ -176,6 +186,81 @@ if (state.timer_handle >= 0)
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function sets up the geometry data pointers in a
+//  struct cactus_grid_info .
+//
+
+namespace {
+void setup_Cactus_gridfn_data_ptrs(const cGH *GH, struct cactus_grid_info& cgi)
+{
+cgi.g_dd_11_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_11_varindex, "g_11");
+cgi.g_dd_12_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_12_varindex, "g_12");
+cgi.g_dd_13_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_13_varindex, "g_13");
+cgi.g_dd_22_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_22_varindex, "g_22");
+cgi.g_dd_23_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_23_varindex, "g_23");
+cgi.g_dd_33_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_33_varindex, "g_33");
+cgi.K_dd_11_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_11_varindex, "K_11");
+cgi.K_dd_12_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_12_varindex, "K_12");
+cgi.K_dd_13_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_13_varindex, "K_13");
+cgi.K_dd_22_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_22_varindex, "K_22");
+cgi.K_dd_23_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_23_varindex, "K_23");
+cgi.K_dd_33_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_33_varindex, "K_33");
+cgi.psi_data     = Cactus_gridfn_data_ptr(GH, cgi.psi_varindex,     "psi",
+					  cgi.Cactus_conformal_metric);
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function gets the Cactus data pointer for a single gridfn, and
+// optionally checks to make sure this is non-NULL.
+//
+// Arguments:
+// gridfn_name[] = The character-string name of the grid function;
+//		   this is used only for formatting error messages.
+// check_for_NULL = true ==> check to make sure the data pointer is non-NULL
+//		    false ==> skip this check (presumably because a NULL
+//			      pointer is ok)
+//
+namespace {
+const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
+					const char gridfn_name[],
+					bool check_for_NULL /* = true */)
+{
+const int time_level = 0;
+
+//
+// CCTK_VarDataPtrI() returns a  void * , but we need a  const CCTK_REAL*;
+// since  static_cast<>  won't change const-ness, we need a 2-stage cast
+// for this:
+//
+const CCTK_REAL* data_ptr = static_cast<const fp*>(
+			       const_cast<const void *>(
+				  CCTK_VarDataPtrI(GH, time_level, varindex)
+						       )
+						  );
+
+if (check_for_NULL && (data_ptr == NULL))
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   Cactus_gridfn_data_ptr(): got unexpected NULL data pointer\n"
+"                             for Cactus geometry gridfn!\n"
+"                             name=\"%s\" varindex=%d"
+		   ,
+		   gridfn_name, varindex);			/*NOTREACHED*/
+
+return data_ptr;
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function finds (or more accurately tries to find) a single
@@ -353,6 +438,8 @@ default:
 	  }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // Given that an apparent horizon has been found, this function computes