add support for outputting te change in the horizon position from

one Newton iteration to the next, Delta_h, at each Newton iteration
(for help debugging convergence problems)

also move parameters for "should we output {h,H,Delta_h} from
struct IO_info to struct solver_info

also make the instructions on setting environment variables in
make.configuration.defn even more verbose...


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

1 files changed, 31 insertions, 23 deletions

diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 7df1fb3..7d38b34 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -61,7 +61,7 @@ bool find_horizon(enum method method,
 		  const struct verbose_info& verbose_info, int timer_handle,
 		  struct IO_info& IO_info,
 		  struct Jacobian_info& Jac_info,
-		  struct solver_info& solver_info,
+		  const struct solver_info& solver_info, bool initial_find_flag,
 		  struct cactus_grid_info& cgi, struct geometry_info& gi,
 		  patch_system& ps, Jacobian* Jac_ptr,
 		  int hn, int N_horizons);
@@ -83,7 +83,10 @@ extern "C"
 {
 DECLARE_CCTK_ARGUMENTS
 DECLARE_CCTK_PARAMETERS
+
 const struct verbose_info& verbose_info = state.verbose_info;
+      struct IO_info&      IO_info      = state.IO_info;
+const struct solver_info&  solver_info  = state.solver_info;
 
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
@@ -105,24 +108,13 @@ if (state.cgi.Cactus_conformal_metric)
    then state.cgi.psi_data = Cactus_gridfn_data_ptr(cctkGH,
 						    "StaticConformal::psi");
 
-state.IO_info.time_iteration = cctk_iteration;
+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];
 	patch_system& ps = *AH_info.ps_ptr;
 
 	//
-	// The first time we (try to) find a given horizon, our initial
-	// guess is likely to be rather inaccurate, so we may need a
-	// larger number of iterations.  But if we've found this horizon
-	// before, then we have its previous position as an initial guess,
-	// so we shouldn't need as many iterations.
-	//
-	state.solver_info.max_Newton_iterations
-		= AH_info.AH_found ? max_Newton_iterations__subsequent
-				   : max_Newton_iterations__initial;
-
-	//
 	// If this is our first attempt to find this horizon, or
 	//    if we've tried to find it before but we failed on our
 	//    immediately previous attempt, then we need to (re)set
@@ -132,16 +124,25 @@ state.IO_info.time_iteration = cctk_iteration;
 	//    then we can just reuse the previous horizon position as
 	//    our initial guess for this time around.
 	//
-	if (! AH_info.AH_found)
-	   then setup_initial_guess(ps,
+	const bool initial_find_flag = ! AH_info.AH_found;
+	if (initial_find_flag)
+	   then {
+		setup_initial_guess(ps,
 				    AH_info.initial_guess_info,
-				    state.IO_info,
+				    IO_info,
 				    hn, verbose_info);
+		if (solver_info.output_initial_guess)
+		   then output_gridfn(ps, gfns::gfn__h,
+				      IO_info, IO_info.h_base_file_name,
+				      hn, verbose_info
+					  .print_algorithm_highlights);
+		}
 
 	AH_info.AH_found
 	   = find_horizon(state.method,
 			  verbose_info, state.timer_handle,
-			  state.IO_info, state.Jac_info, state.solver_info,
+			  IO_info, state.Jac_info,
+			  solver_info, initial_find_flag,
 			  state.cgi, state.gi,
 			  ps, AH_info.Jac_ptr,
 			  hn, state.N_horizons);
@@ -184,6 +185,12 @@ if (state.timer_handle >= 0)
 // timer_handle = a valid Cactus timer handle if we want to time the
 //		  apparent horizon process, or -ve to skip this
 //		  (we only time the computation, not the file I/O)
+// initial_find_flag = true ==> This is the first attempt to find this
+//				horizon, or this is a subsequent attempt
+//				and the immediately previous attempt failed.
+//		       false ==> This isn't the first attempt to find this
+//				 horizon, and we found it successfully on
+//				 the immediately previous attempt.
 // Jac_ptr = may be NULL if no Jacobian is needed (depending on  method)
 // hn = the horizon number (used only in naming output files)
 // N_horizons = the total number of horizon(s) being searched for number
@@ -201,7 +208,7 @@ bool find_horizon(enum method method,
 		  const struct verbose_info& verbose_info, int timer_handle,
 		  struct IO_info& IO_info,
 		  struct Jacobian_info& Jac_info,
-		  struct solver_info& solver_info,
+		  const struct solver_info& solver_info, bool initial_find_flag,
 		  struct cactus_grid_info& cgi, struct geometry_info& gi,
 		  patch_system& ps, Jacobian* Jac_ptr,
 		  int hn, int N_horizons)
@@ -226,7 +233,7 @@ case method__horizon_function:
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "   H(h) rms-norm %.2e, infinity-norm %.2e",
 			   H_norms.rms_norm(), H_norms.infinity_norm());
-	if (IO_info.output_H)
+	if (solver_info.output_H)
 	   then output_gridfn(ps, gfns::gfn__H,
 			      IO_info, IO_info.H_base_file_name,
 			      hn, true);
@@ -297,18 +304,19 @@ case method__Newton_solve:
 	const bool status
 		= Newton_solve(ps, Jac,
 			       cgi, gi,
-			       Jac_info, solver_info, IO_info,
-			       hn, verbose_info);
+			       Jac_info,
+			       solver_info, initial_find_flag,
+			       IO_info, hn, verbose_info);
 	if (timer_handle >= 0)
 	   then CCTK_TimerStopI(timer_handle);
 	if (! status)
 	   then return false;				// *** ERROR RETURN ***
 
-	if (IO_info.output_h)
+	if (solver_info.output_h)
 	   then output_gridfn(ps, gfns::gfn__h,
 			      IO_info, IO_info.h_base_file_name,
 			      hn, verbose_info.print_algorithm_details);
-	if (IO_info.output_H)
+	if (solver_info.output_H)
 	   then output_gridfn(ps, gfns::gfn__H,
 			      IO_info, IO_info.H_base_file_name,
 			      hn, verbose_info.print_algorithm_details);