add support for (mirrored, not rotating) quadrant and octant grids

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

4 files changed, 129 insertions, 37 deletions

diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 1811ce6..3db180d 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -234,7 +234,8 @@ extern "C"
 void setup_initial_guess(patch_system& ps,
 			 const struct initial_guess_info& igi,
 			 const struct IO_info& IO_info,
-			 int hn, const struct verbose_info& verbose_info);
+			 int hn, int N_horizons,
+			 const struct verbose_info& verbose_info);
 enum initial_guess_method
   decode_initial_guess_method(const char initial_guess_method_string[]);
 
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 446eae1..bfeedef 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -168,7 +168,7 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 		setup_initial_guess(ps,
 				    AH_info.initial_guess_info,
 				    IO_info,
-				    hn, verbose_info);
+				    hn, state.N_horizons, verbose_info);
 		if (IO_info.output_initial_guess)
 		   then output_gridfn(ps, gfns::gfn__h,
 				      IO_info, IO_info.h_base_file_name,
@@ -186,6 +186,7 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 				    ps,
 				    hn, state.N_horizons);
 		break;
+
 	case method__test_Jacobian:
 		do_test_Jacobian(verbose_info, state.timer_handle,
 				 IO_info,
@@ -194,6 +195,7 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 				 ps, AH_info.Jac_ptr,
 				 hn, state.N_horizons);
 		break;
+
 	case method__find_horizon:
 		AH_info.AH_found
 		   = do_find_horizon(verbose_info, state.timer_handle,
@@ -250,6 +252,7 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 					   "no apparent horizon found");
 			}
 		break;
+
 	default:
 		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 			   "\n"
@@ -581,27 +584,29 @@ fp integral_z = surface_integral(ps, gfns::gfn__global_z,
 //
 switch	(ps.type())
 	{
-case patch_system::full_sphere_patch_system:
+case patch_system::patch_system__full_sphere:
 	break;
-case patch_system::plus_z_hemisphere_patch_system:
+case patch_system::patch_system__plus_z_hemisphere:
 	integral_one *= 2.0;
 	integral_x *= 2.0;
 	integral_y *= 2.0;
 	integral_z = integral_one * ps.origin_z();
 	break;
-case patch_system::plus_xy_quadrant_patch_system:
+case patch_system::patch_system__plus_xy_quadrant_mirrored:
+case patch_system::patch_system__plus_xy_quadrant_rotating:
 	integral_one *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
 	integral_z *= 4.0;
 	break;
-case patch_system::plus_xz_quadrant_patch_system:
+case patch_system::patch_system__plus_xz_quadrant_rotating:
 	integral_one *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y *= 4.0;
 	integral_z = integral_one * ps.origin_z();
 	break;
-case patch_system::plus_xyz_octant_patch_system:
+case patch_system::patch_system__plus_xyz_octant_mirrored:
+case patch_system::patch_system__plus_xyz_octant_rotating:
 	integral_one *= 8.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
index c72590d..369aaeb 100644
--- a/src/driver/initial_guess.cc
+++ b/src/driver/initial_guess.cc
@@ -71,11 +71,13 @@ void setup_coord_ellipsoid(patch_system& ps,
 void setup_initial_guess(patch_system& ps,
 			 const struct initial_guess_info& igi,
 			 const struct IO_info& IO_info,
-			 int hn, const struct verbose_info& verbose_info)
+			 int hn, int N_horizons,
+			 const struct verbose_info& verbose_info)
 {
 if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "setting initial guess for horizon shape h");
+		   "setting initial guess for horizon %d of %d",
+		   hn, N_horizons);
 
 switch	(igi.method)
 	{
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 125b716..2de1959 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -62,6 +62,9 @@ enum verbose_level
 
 enum patch_system::patch_system_type
   choose_patch_system_type(const char grid_domain[],
+			   const char grid_bitant_plane[],
+			   const char grid_quadrant_direction[],
+			   const char grid_rotation_axis[],
 			   fp origin_x, fp origin_y, fp origin_z);
 	  }
 
@@ -85,11 +88,11 @@ DECLARE_CCTK_ARGUMENTS
 DECLARE_CCTK_PARAMETERS
 
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "initializing AHFinderDirect data structures");
+	   "setting up AHFinderDirect data structures");
 
 state.N_horizons = N_horizons;
 CCTK_VInfo(CCTK_THORNSTRING,
-           "             to search for %d horizon%s",
+           "           to search for %d horizon%s",
 	   state.N_horizons, ((state.N_horizons == 1) ? "" : "s"));
 
 
@@ -248,28 +251,29 @@ state.AH_info_ptrs.push_back(NULL);
 
 	for (int hn = 1 ; hn <= N_horizons ; ++hn)
 	{
+	if (verbose_info.print_algorithm_highlights)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   setting up data structures for horizon %d of %d",
+			   hn, int(N_horizons));
+
 	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])
+	     ? // choose a patch system type based on grid:: parameters
+	       // ... we inherit from grid, and we ask for some of its
+	       //     parameters in our param.ccl file; these appear as
+	       //     ordinary Cactus parameters here, so (eg)
+	       //     grid::domain is just "domain" here
+	       choose_patch_system_type(/* grid:: */ domain,
+					/* grid:: */ bitant_plane,
+					/* grid:: */ quadrant_direction,
+					/* grid:: */ rotation_axis,
+					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
@@ -280,15 +284,22 @@ state.AH_info_ptrs.push_back(NULL);
 			   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);
+			   true, verbose_info.print_algorithm_details);
 	AH_ptr->ps_ptr = &ps;
 	ps.set_gridfn_to_constant(1.0, gfns::gfn__one);
 
 	// create the Jacobian matrix
+	if (verbose_info.print_algorithm_details)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "      constructing Jacobian matrix");
 	AH_ptr->Jac_ptr = (state.method == method__horizon_function)
 			  ? NULL	// no Jacobian used for this case
 			  : &new_Jacobian(ps, Jac_info.Jacobian_storage_method);
 
+	if (verbose_info.print_algorithm_details)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "      setting initial guess parameters etc");
+
 	// store the initial guess parameters
 	AH_ptr->initial_guess_info.method
 		= decode_initial_guess_method(initial_guess_method);
@@ -347,7 +358,8 @@ state.AH_info_ptrs.push_back(NULL);
 
 	// set up the BH-diagnostics output files?
 	AH_ptr->BH_diagnostics_fileptr
-		= output_BH_diagnostics
+		= ( output_BH_diagnostics
+		    && (state.method == method__find_horizon) )
 		  ? setup_BH_diagnostics_output_file(IO_info, hn, N_horizons)
 		  : NULL;
 
@@ -412,28 +424,100 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 // symmetry and the patch system's origin position.
 //
 // Arguments:
-// grid_domain = The value of the Cactus grid::domain parameter.
+// grid_* = The values of the Cactus parameters
+//		grid::domain
+//		grid::bitant_plane
+//		grid::quadrant_direction
+//		grid::rotation_axis
 // origin_[xyz] = The origin position for this patch system.
 //
 namespace {
 enum patch_system::patch_system_type
   choose_patch_system_type(const char grid_domain[],
+			   const char grid_bitant_plane[],
+			   const char grid_quadrant_direction[],
+			   const char grid_rotation_axis[],
 			   fp origin_x, fp origin_y, fp origin_z)
 {
 if	(STRING_EQUAL(grid_domain, "full"))
-   then return patch_system::full_sphere_patch_system;
+   then return patch_system::patch_system__full_sphere;
 
-else if (STRING_EQUAL(grid_domain, "bitant"))
+else if ( STRING_EQUAL(grid_domain, "bitant")
+	  && STRING_EQUAL(grid_bitant_plane, "xy") )
    then {
 	if (origin_z == 0.0)
-	   then return patch_system::plus_z_hemisphere_patch_system;
+	   then return patch_system::patch_system__plus_z_hemisphere;
 	   else {
 		CCTK_VInfo(CCTK_THORNSTRING,
-		   "      bitant mode but patch system origin_z=%g != 0",
+			   "      Cactus has bitant mode about xy plane");
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "      but patch system origin_z=%g != 0",
 			   double(origin_z));
 		CCTK_VInfo(CCTK_THORNSTRING,
+			   "      ==> using \"full sphere\" patch system");
+		return patch_system::patch_system__full_sphere;
+		}
+	}
+
+else if ( STRING_EQUAL(grid_domain, "quadrant")
+	  && STRING_EQUAL(grid_quadrant_direction, "z") )
+   then {
+	if ((origin_x == 0.0) && (origin_y == 0.0))
+	   then return patch_system::patch_system__plus_xy_quadrant_mirrored;
+	   else {
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      Cactus has quadrant mode about z axis");
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      but patch system origin_(x,y)=(%g,%g) != (0,0)",
+		   double(origin_x), double(origin_y));
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      ==> using \"full sphere\" patch system");
+		return patch_system::patch_system__full_sphere;
+		}
+	}
+
+else if ( STRING_EQUAL(grid_domain, "quadrant_reflect_rotate")
+	  && STRING_EQUAL(grid_quadrant_direction, "z")
+	  && STRING_EQUAL(grid_rotation_axis, "z") )
+   then {
+	if ((origin_x == 0.0) && (origin_y == 0.0))
+	   then return patch_system::patch_system__plus_xy_quadrant_rotating;
+	   else {
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      Cactus has rotating quadrant mode about z axis");
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      but patch system origin_(x,y)=(%g,%g) != (0,0)",
+		   double(origin_x), double(origin_y));
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      ==> using \"full sphere\" patch system");
+		return patch_system::patch_system__full_sphere;
+		}
+	}
+
+else if (STRING_EQUAL(grid_domain, "octant"))
+   then {
+	if	((origin_x == 0.0) && (origin_y == 0.0) && (origin_z == 0.0))
+	   then return patch_system::patch_system__plus_xyz_octant_mirrored;
+	else if ((origin_x == 0.0) && (origin_y == 0.0))
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      Cactus has mirrored octant mode");
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      but patch system origin_z=%g != 0",
+		   double(origin_z));
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      ==> using \"+xy quadrant (mirrored)\" patch system");
+		return patch_system::patch_system__plus_xy_quadrant_mirrored;
+		}
+	else	{
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      Cactus has mirrored octant mode");
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      but patch system origin_(x,y,z)=(%g,%g,%g) != (0,0,0)",
+		   double(origin_x), double(origin_y), double(origin_z));
+		CCTK_VInfo(CCTK_THORNSTRING,
 		   "      ==> using \"full sphere\" patch system");
-		return patch_system::full_sphere_patch_system;
+		return patch_system::patch_system__full_sphere;
 		}
 	}
 
@@ -443,10 +527,10 @@ else	{
 "   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 :("
+"            but this may or may not work..."
 ,
 		   grid_domain);
-	return patch_system::full_sphere_patch_system;
+	return patch_system::patch_system__full_sphere;
 	}
 }
 	  }