add auto-determine of patch system type from Cactus grid symmetry

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

2 files changed, 112 insertions, 16 deletions

diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 62e3adc..86ca722 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -4,13 +4,15 @@
 // <<<prototypes for functions local to this file>>>
 // <<<access to persistent data>>>
 //
-// AHFinderDirect_driver - top-level driver
+// AHFinderDirect_driver - top-level driver to setup persistent data structures
 ///
 /// decode_method - decode the method parameter
 /// decode_verbose_level - decode the verbose_level parameter
 /// decode_Jacobian_method - decode the Jacobian_method parameter
 /// decode_geometry_method - decode the geometry_method parameter
 ///
+/// choose_patch_system_type - choose patch system type
+///
 
 #include <stdio.h>
 #include <assert.h>
@@ -44,7 +46,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
-#include "AHFinderDirect.hh"
+#include "driver.hh"
 
 //******************************************************************************
 
@@ -61,6 +63,10 @@ enum Jacobian_method
   decode_Jacobian_method(const char Jacobian_method_string[]);
 enum geometry_method
   decode_geometry_method(const char geometry_method_string[]);
+
+enum patch_system::patch_system_type
+  choose_patch_system_type(const char grid_domain[],
+			   fp origin_x, fp origin_y, fp origin_z);
 	  };
 
 //******************************************************************************
@@ -98,10 +104,22 @@ CCTK_VInfo(CCTK_THORNSTRING,
 //
 
 state.method = decode_method(method);
-state.verbose_level = decode_verbose_level(verbose_level);
+state.verbose_info.verbose_level = decode_verbose_level(verbose_level);
+state.verbose_info.print_physics_highlights
+   = (state.verbose_info.verbose_level >= verbose_level__physics_highlights);
+state.verbose_info.print_physics_details
+   = (state.verbose_info.verbose_level >= verbose_level__physics_details);
+state.verbose_info.print_algorithm_highlights
+   = (state.verbose_info.verbose_level >= verbose_level__algorithm_highlights);
+state.verbose_info.print_algorithm_details
+   = (state.verbose_info.verbose_level >= verbose_level__algorithm_details);
 state.timer_handle = (print_timing_stats != 0)
 		     ? CCTK_TimerCreateI()
 		     : -1;
+state.surface_integral_method
+   = patch::decode_integration_method(surface_integral_method);
+
+const struct verbose_info& verbose_info = state.verbose_info;
 
 struct IO_info& IO_info = state.IO_info;
 IO_info.h_base_file_name        = h_base_file_name;
@@ -125,11 +143,10 @@ solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
 solver_info.final_H_update_if_exit_x_H_small
 	= (final_H_update_if_exit_x_H_small != 0);
 
-
 //
 // set up the Cactus grid info
 //
-if (state.verbose_level >= verbose_level__algorithm_highlights)
+if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
 struct cactus_grid_info& cgi = state.cgi;
 cgi.GH = cctkGH;
@@ -168,7 +185,7 @@ struct geometry_info& gi = state.gi;
 gi.geometry_method = decode_geometry_method(geometry_method);
 
 // parameters for geometry_method = "interpolate from Cactus grid"
-if (state.verbose_level >= verbose_level__algorithm_highlights)
+if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
 gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
 if (gi.operator_handle < 0)
@@ -195,7 +212,7 @@ gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz;
 //
 
 // set up the interpatch interpolator
-if (state.verbose_level >= verbose_level__algorithm_highlights)
+if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up interpatch interpolator");
 const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name);
 if (interp_handle < 0)
@@ -209,20 +226,49 @@ if (interp_param_table_handle < 0)
 "AHFinderDirect_setup(): bad interpatch-interpolator parameter(s) \"%s\"!",
 		   interpatch_interpolator_pars);		/*NOTREACHED*/
 
-	for (int hn = 0 ; hn < N_horizons ; ++hn)
+const bool ps_type_from_Cactus_grid_sym
+   = STRING_EQUAL(patch_system_type, "match Cactus grid symmetry");
+
+// we don't use the [0] array element
+state.AH_info_ptrs.push_back(NULL);
+
+	for (int hn = 1 ; hn <= N_horizons ; ++hn)
 	{
 	struct AH_info* AH_ptr = new AH_info;
 
+	// decide what type of patch system this one should be
+	const enum patch_system::patch_system_type ps_type
+	   = ps_type_from_Cactus_grid_sym
+	     ? // choose a patch system type based on grid::domain
+	       // ... since we inherit from grid, we can see its
+	       //     parameters, so grid::domain is just "domain" here
+	       choose_patch_system_type(domain, origin_x[hn],
+						origin_y[hn],
+						origin_z[hn])
+	     : patch_system::type_of_name(patch_system_type);
+	if (verbose_info.print_algorithm_highlights)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "   for horizon #%d/%d: %s patch system",
+			   hn, int(N_horizons),
+			   patch_system::name_of_type(ps_type));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "                     at origin (%g,%g,%g)",
+			   origin_x[hn], origin_y[hn], origin_z[hn]);
+		}
+
 	// create the patch system
 	patch_system& ps
-	   = *new patch_system
-		     (origin_x[hn], origin_y[hn], origin_z[hn],
-		      patch_system::type_of_name(patch_system_type),
-		      N_ghost_points, N_overlap_points, delta_drho_dsigma,
-		      gfns::nominal_min_gfn, gfns::nominal_max_gfn,
-		      gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
-		      interp_handle, interp_param_table_handle);
+	   = *new
+	      patch_system(origin_x[hn], origin_y[hn], origin_z[hn],
+			   ps_type,
+			   N_ghost_points, N_overlap_points, delta_drho_dsigma,
+			   gfns::nominal_min_gfn, gfns::nominal_max_gfn,
+			   gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
+			   interp_handle, interp_param_table_handle,
+			   verbose_info.print_algorithm_details);
 	AH_ptr->ps_ptr = &ps;
+	ps.set_gridfn_to_constant(1.0, gfns::gfn__one);
 
 	// create the Jacobian matrix
 	AH_ptr->Jac_ptr = (state.method == method__horizon_function)
@@ -334,3 +380,53 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   geometry_method_string);			/*NOTREACHED*/
 }
 	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function chooses a patch system type based on the Cactus grid
+// symmetry and the patch system's origin position.
+//
+// Arguments:
+// grid_domain = The value of the Cactus grid::domain parameter.
+// origin_[xyz] = The origin position for this patch system.
+//
+namespace {
+enum patch_system::patch_system_type
+  choose_patch_system_type(const char grid_domain[],
+			   fp origin_x, fp origin_y, fp origin_z)
+{
+if	(STRING_EQUAL(grid_domain, "full"))
+   then return patch_system::full_sphere_patch_system;
+
+else if (STRING_EQUAL(grid_domain, "bitant"))
+   then {
+	if (origin_z == 0.0)
+	   then return patch_system::plus_z_hemisphere_patch_system;
+	   else {
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   choose_patch_system_type()\n"
+"        bitant mode but patch system origin_z=%g != 0\n"
+"        ==> using \"full sphere\" patch system"
+,
+			   double(origin_z));
+		return patch_system::full_sphere_patch_system;
+		}
+	}
+
+else	{
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   choose_patch_system_type()\n"
+"        grid::domain = \"%s\" not supported!\n"
+"        ==> will try using \"full sphere\" patch system\n",
+"            but this probably won't work :("
+,
+		   grid_domain);
+	return patch_system::full_sphere_patch_system;
+	}
+}
+	  }
diff --git a/src/driver/state.cc b/src/driver/state.cc
index 52d62c8..b1f78c7 100644
--- a/src/driver/state.cc
+++ b/src/driver/state.cc
@@ -33,7 +33,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
-#include "AHFinderDirect.hh"
+#include "driver.hh"
 
 //******************************************************************************