add new parameter and code to allow computing surfaces of constant expansion

(the case expansion = 0 is the AH, others should be surfaces nested inside
or outside the AH)


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

8 files changed, 55 insertions, 34 deletions

diff --git a/param.ccl b/param.ccl
index 415358c..21cd349 100644
--- a/param.ccl
+++ b/param.ccl
@@ -105,6 +105,15 @@ boolean print_timing_stats \
 {
 } "false"
 
+# By default, we find apparent horizons by solving the equation
+#  Theta(h) = 0
+# This parameter allows the RHS to be set to any specified constant,
+# to find a surface of constant expansion.
+real surface_expansion[11] "search for a surface with this (constant) expansion"
+{
+*:* :: "any real number"
+} 0.0
+
 ################################################################################
 
 #
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 94eb9cf..7dddda7 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -145,8 +145,12 @@ if (hs.has_genuine_horizons())
 
     struct AH_data* AH_data_ptr
 	                  = horizon_is_genuine ? AH_data_array[hn]   : NULL;
-    patch_system*  ps_ptr = horizon_is_genuine ? AH_data_ptr->ps_ptr : NULL;
-    Jacobian*     Jac_ptr = horizon_is_genuine ? AH_data_ptr->Jac_ptr: NULL;
+    patch_system* const  ps_ptr = horizon_is_genuine
+				  ? AH_data_ptr->ps_ptr : NULL;
+    Jacobian*     const Jac_ptr = horizon_is_genuine
+				  ? AH_data_ptr->Jac_ptr: NULL;
+    const fp add_to_expansion = horizon_is_genuine
+				? -AH_data_ptr->surface_expansion : 0.0;
 
     const int max_iterations
        = horizon_is_genuine
@@ -180,7 +184,7 @@ if (hs.has_genuine_horizons())
 
 	jtutil::norm<fp> Theta_norms;
 	const enum expansion_status
-	  expansion_status = expansion(ps_ptr,
+	  expansion_status = expansion(ps_ptr, add_to_expansion,
 				       cgi, gi,
 				       error_info, (iteration == 1),
 				       true,	// yes, we want Jacobian coeffs
@@ -396,6 +400,7 @@ if (hs.has_genuine_horizons())
 	     = expansion_Jacobian
 		 (this_horizon_needs_more_iterations ? ps_ptr : NULL,
 		  this_horizon_needs_more_iterations ? Jac_ptr : NULL,
+		  add_to_expansion,
 		  cgi, gi, Jacobian_info,
 		  error_info, (iteration == 1),
 		  verbose_info.print_algorithm_details);
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 7a3c734..5826fa0 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -262,6 +262,7 @@ struct	AH_data
 	//
 	patch_system* ps_ptr;
 	Jacobian* Jac_ptr;
+	fp surface_expansion;
 
 	// are we finding this horizon "from scratch"?
 	// ... true if this is the first time we've tried to find it,
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index f3f3803..5f891fd 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -316,7 +316,7 @@ if (active_flag)
 		if (timer_handle >= 0)
 		   then CCTK_TimerStartI(timer_handle);
 		jtutil::norm<fp> Theta_norms;
-		const bool Theta_ok = expansion(&ps,
+		const bool Theta_ok = expansion(&ps, -AH_data.surface_expansion,
 						cgi, gi,
 						error_info, true,// initial eval
 						false,	// no Jacobian coeffs
@@ -347,7 +347,7 @@ if (active_flag)
    else {
 		for (int i = 0 ; i < N_horizons ; ++i)
 		{
-		expansion(NULL,			// dummy computation
+		expansion(NULL, 0.0,		// dummy computation
 			  cgi, gi,
 			  error_info, true);	// initial evaluation
 		}
@@ -399,18 +399,19 @@ assert(N_horizons >= 1);
 const bool print_msg_flag = true;
 const int hn = 1;
 
-struct AH_data* AH_data_ptr = active_flag ? AH_data_array[hn]      : NULL;
-patch_system*        ps_ptr = active_flag ? AH_data_ptr->ps_ptr : NULL;
+struct AH_data* const AH_data_ptr = active_flag ? AH_data_array[hn]   : NULL;
+patch_system*   const      ps_ptr = active_flag ? AH_data_ptr->ps_ptr : NULL;
+const fp add_to_expansion = active_flag ? -AH_data_ptr->surface_expansion : 0.0;
 
 //
 // numerical-perturbation Jacobian
 //
 Jacobian* Jac_NP_ptr = active_flag ? AH_data_ptr->Jac_ptr : NULL;
-expansion(ps_ptr,
+expansion(ps_ptr, add_to_expansion,
 	  cgi, gi,
 	  error_info, true);		// initial evaluation
 Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
-expansion_Jacobian(ps_ptr, Jac_NP_ptr,
+expansion_Jacobian(ps_ptr, Jac_NP_ptr, add_to_expansion,
 		   cgi, gi, Jac_info,
 		   error_info, true,	// initial evaluation
 		   print_msg_flag);
@@ -424,12 +425,12 @@ if (test_all_Jacobian_compute_methods)
 					 *ps_ptr,
 					 verbose_info.print_algorithm_details)
 			  : NULL;
-	expansion(ps_ptr,
+	expansion(ps_ptr, add_to_expansion,
 		  cgi, gi,
 		  error_info, true,	// initial evaluation
 		  true);		// compute SD Jacobian coeffs
 	Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
-	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
+	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr, add_to_expansion,
 			   cgi, gi, Jac_info,
 			   error_info, true,	// initial evaluation
 			   print_msg_flag);
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 84cb548..f7137e1 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -492,6 +492,8 @@ if (strlen(surface_interpolator_name) > 0)
 					 verbose_info.print_algorithm_details)
 			  : NULL;
 
+	AH_data.surface_expansion = surface_expansion[hn];
+
 	AH_data.initial_find_flag = genuine_flag;
 
 	if (genuine_flag)
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index 51af9ed..f26cac1 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -76,8 +76,8 @@ bool geometry_is_finite(patch_system& ps,
 			const struct error_info& error_info, bool initial_flag,
 			bool print_msg_flag);
 
-bool compute_Theta(patch_system& ps, bool Jacobian_flag,
-		   jtutil::norm<fp>* Theta_norms_ptr,
+bool compute_Theta(patch_system& ps, fp add_to_expansion,
+		   bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr,
 		   const struct error_info& error_info, bool initial_flag,
 		   bool print_msg_flag);
 	  }
@@ -128,6 +128,8 @@ bool compute_Theta(patch_system& ps, bool Jacobian_flag,
 //	      interpolator call are done, the latter with the number of
 //	      interpolation points is set to 0 and all the output array
 //	      pointers set to NULL.
+// add_to_expansion = A real number which is added to the computed expansion
+//		      at each grid point.
 // initial_flag = true if this is the first evaluation of  expansion()
 //		       for this horizon,
 //		  false otherwise;
@@ -147,7 +149,7 @@ bool compute_Theta(patch_system& ps, bool Jacobian_flag,
 // succeeded or failed, and if the latter, what caused the failure.
 //
 enum expansion_status
-  expansion(patch_system* ps_ptr,
+  expansion(patch_system* ps_ptr, fp add_to_expansion,
 	    const struct cactus_grid_info& cgi,
 	    const struct geometry_info& gi,
 	    const struct error_info& error_info, bool initial_flag,
@@ -229,7 +231,7 @@ if (active_flag)
 	// compute remaining gridfns --> $\Theta$
 	// and optionally also the Jacobian coefficients
 	// by algebraic ops and angular finite differencing
-	if (!compute_Theta(*ps_ptr,
+	if (!compute_Theta(*ps_ptr, add_to_expansion,
 			   Jacobian_flag, Theta_norms_ptr,
 			   error_info, initial_flag,
 			   print_msg_flag))
@@ -1150,8 +1152,8 @@ return true;					// *** no check possible ***
 // g_ij isn't positive definite).
 //
 namespace {
-bool compute_Theta(patch_system& ps, bool Jacobian_flag,
-		   jtutil::norm<fp>* Theta_norms_ptr,
+bool compute_Theta(patch_system& ps, fp add_to_expansion,
+		   bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr,
 		   const struct error_info& error_info, bool initial_flag,
 		   bool print_msg_flag)
 {
@@ -1254,7 +1256,8 @@ CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING,
 		Theta = + Theta_A/(Theta_D*sqrt_Theta_D)
 			+ Theta_B/sqrt_Theta_D
 			+ Theta_C/Theta_D
-			- K;
+			- K
+			+ add_to_expansion;
 
 		// update running norms of Theta(h) function
 		if (Theta_norms_ptr != NULL)
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index 7d102f3..d360076 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -50,15 +50,14 @@ using jtutil::error_exit;
 namespace {
 enum expansion_status
   expansion_Jacobian_NP
-	(patch_system& ps, Jacobian& Jac,
+	(patch_system& ps, Jacobian& Jac, fp add_to_expansion,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
 	 const struct Jacobian_info& Jacobian_info,
 	 const struct error_info& error_info, bool initial_flag,
 	 bool print_msg_flag);
 
-void expansion_Jacobian_partial_SD(patch_system& ps,
-				   Jacobian& Jac,
+void expansion_Jacobian_partial_SD(patch_system& ps, Jacobian& Jac,
 				   const struct cactus_grid_info& cgi,
 				   const struct geometry_info& gi,
 				   const struct Jacobian_info& Jacobian_info,
@@ -71,7 +70,7 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
 			     int xm_irho, int xm_isigma);
 enum expansion_status
   expansion_Jacobian_dr_FD
-	(patch_system* ps_ptr, Jacobian* Jac_ptr,
+	(patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
 	 const struct Jacobian_info& Jacobian_info,
@@ -99,14 +98,15 @@ enum expansion_status
 // Arguments:
 // ps_ptr --> The patch system, or == NULL to do (only) a dummy computation.
 // Jac_ptr --> The Jacobian, or == NULL to do (only) a dummy computation.
+// add_to_expansion = A real number to add to the expansion.
 //
 // Results:
 // This function returns a status code indicating whether the computation
 // succeeded or failed, and if the latter, what caused the failure.
 //
 enum expansion_status
-  expansion_Jacobian(patch_system* ps_ptr,
-		     Jacobian* Jac_ptr,
+  expansion_Jacobian(patch_system* ps_ptr, Jacobian* Jac_ptr,
+		     fp add_to_expansion,
 		     const struct cactus_grid_info& cgi,
 		     const struct geometry_info& gi,
 		     const struct Jacobian_info& Jacobian_info,
@@ -122,6 +122,7 @@ case Jacobian__numerical_perturbation:
 	if (active_flag)
 	   then {
 		status = expansion_Jacobian_NP(*ps_ptr, *Jac_ptr,
+					       add_to_expansion,
 					       cgi, gi, Jacobian_info,
 					       error_info, initial_flag,
 					       print_msg_flag);
@@ -149,7 +150,7 @@ case Jacobian__symbolic_diff_with_FD_dr:
 	// this function looks at ps_ptr and Jac_ptr (non-NULL vs NULL)
 	// to choose a normal vs dummy computation
 	  {
-	status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr,
+	status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, add_to_expansion,
 					  cgi, gi, Jacobian_info,
 					  error_info, initial_flag,
 					  print_msg_flag);
@@ -210,7 +211,7 @@ return expansion_success;				// *** NORMAL RETURN ***
 namespace {
 enum expansion_status
   expansion_Jacobian_NP
-	(patch_system& ps, Jacobian& Jac,
+	(patch_system& ps, Jacobian& Jac, fp add_to_expansion,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
 	 const struct Jacobian_info& Jacobian_info,
@@ -243,7 +244,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 	const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
 	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
 	const
-	  enum expansion_status status = expansion(&ps,
+	  enum expansion_status status = expansion(&ps, add_to_expansion,
 						   cgi, gi,
 						   error_info, initial_flag);
 	if (status != expansion_success)
@@ -299,8 +300,7 @@ return expansion_success;				// *** NORMAL RETURN ***
 //	The Jacobian matrix is stored in the Jacobian object Jac.
 //
 namespace {
-void expansion_Jacobian_partial_SD(patch_system& ps,
-				   Jacobian& Jac,
+void expansion_Jacobian_partial_SD(patch_system& ps, Jacobian& Jac,
 				   const struct cactus_grid_info& cgi,
 				   const struct geometry_info& gi,
 				   const struct Jacobian_info& Jacobian_info,
@@ -532,7 +532,7 @@ patch& yp = ye.my_patch();
 namespace {
 enum expansion_status
   expansion_Jacobian_dr_FD
-	(patch_system* ps_ptr, Jacobian* Jac_ptr,
+	(patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion,
 	 const struct cactus_grid_info& cgi,
 	 const struct geometry_info& gi,
 	 const struct Jacobian_info& Jacobian_info,
@@ -554,7 +554,7 @@ if (active_flag)
 	ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
 	}
 const
-  enum expansion_status status = expansion(ps_ptr,
+  enum expansion_status status = expansion(ps_ptr, add_to_expansion,
 					   cgi, gi,
 					   error_info, initial_flag);
 if (status != expansion_success)
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index bcc8702..11ab58e 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -233,7 +233,7 @@ struct	error_info
 
 // expansion.cc
 enum expansion_status
-  expansion(patch_system* ps_ptr,
+  expansion(patch_system* ps_ptr, fp add_to_expansion,
 	    const struct cactus_grid_info& cgi,
 	    const struct geometry_info& gi,
 	    const struct error_info& error_info, bool initial_flag,
@@ -243,8 +243,8 @@ enum expansion_status
 
 // expansion_Jacobian.cc
 enum expansion_status
-  expansion_Jacobian(patch_system* ps_ptr,
-		     Jacobian* Jac_ptr,
+  expansion_Jacobian(patch_system* ps_ptr, Jacobian* Jac_ptr,
+		     fp add_to_expansion,
 		     const struct cactus_grid_info& cgi,
 		     const struct geometry_info& gi,
 		     const struct Jacobian_info& Jacobian_info,