add new parameter

real min_horizon_size_for_mask \
  "only set mask if min horizon radius >= this number of local grid spacings"

docs will follow shortly...


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

5 files changed, 168 insertions, 6 deletions

diff --git a/par/bl-different-mask-min-radius=0.2.par b/par/bl-different-mask-min-radius=0.2.par
new file mode 100644
index 0000000..cee4170
--- /dev/null
+++ b/par/bl-different-mask-min-radius=0.2.par
@@ -0,0 +1,112 @@
+# This parameter file sets up Brill-Lindquist initial data and finds the
+# apparent horizons in it with AHFinderDirect.
+
+########################################
+ 
+#
+# ***** grid *****
+#
+
+ActiveThorns = "CoordBase CartGrid3D PUGH"
+
+Driver::global_nx   = 52
+Driver::global_ny   = 52
+Driver::global_nz   = 102
+PUGH::ghost_size    = 2
+
+Grid::domain        = "full"
+Grid::type          = "byspacing"
+Grid::dxyz          = 0.04
+Grid::avoid_origin  = "true"
+
+########################################
+
+#
+# physics
+#
+
+ActiveThorns = "ADMBase ADMCoupling CoordGauge StaticConformal SpaceMask"
+ActiveThorns = "IDAnalyticBH"
+
+ADMBase::initial_data              = "bl_bh"
+IDAnalyticBH::bl_nbh = 2
+
+IDAnalyticBH::bl_m_1 = 1.0
+IDAnalyticBH::bl_x0_1 = 0.0
+IDAnalyticBH::bl_y0_1 = 0.0
+IDAnalyticBH::bl_z0_1 = 1.0
+
+IDAnalyticBH::bl_m_2 = 0.5
+IDAnalyticBH::bl_x0_2 = 0.0
+IDAnalyticBH::bl_y0_2 = 0.0
+IDAnalyticBH::bl_z0_2 = -1.0
+
+SpaceMask::use_mask                = "true"
+
+########################################
+
+#
+# numerical methods
+#
+
+Cactus::cctk_itlast = 0
+ADMBase::metric_type            = "static conformal"
+StaticConformal::conformal_storage = "factor+derivs+2nd derivs"
+ADMBase::evolution_method       = "none"
+ADMBase::initial_lapse          = "Cadez"
+ADMBase::lapse_evolution_method = "static"
+ADMBase::initial_shift          = "zero"
+ADMBase::shift_evolution_method = "static"
+
+########################################
+
+ActiveThorns = "PUGHSlab IOBasic IOUtil IOASCII"
+IO::parfile_write = "no"
+IO::out_dir = "bl-different-mask-min-radius=0.2"
+IOASCII::out2D_every = 1
+IOASCII::out2D_style = "gnuplot f(x,y)"
+IOASCII::out2D_vars = "SpaceMask::emask"
+
+########################################
+
+#
+# AHFinderDirect
+#
+
+ActiveThorns = "always AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+
+AHFinderDirect::print_timing_stats = "true"
+AHFinderDirect::verbose_level = "algorithm details"
+
+AHFinderDirect::h_base_file_name              = "bl-different-mask-min-radius=0.2.h"
+AHFinderDirect::BH_diagnostics_base_file_name = "bl-different-mask-min-radius=0.2.BH_diagnostics"
+
+AHFinderDirect::min_horizon_size_for_mask = 5.0 # grid points
+AHFinderDirect::set_mask = "true"
+AHFinderDirect::set_old_style_mask = "true"
+# inside = horizon - 5 grid points
+# buffer = from there out to horizon
+AHFinderDirect::mask_radius_multiplier = 1.0
+AHFInderDirect::mask_radius_offset     = -5.0
+AHFInderDirect::mask_buffer_thickness  =  5.0
+AHFinderDirect::mask_is_noshrink = "true"
+
+AHFinderDirect::N_horizons = 2
+
+AHFinderDirect::origin_x[1] = 0.0
+AHFinderDirect::origin_y[1] = 0.0
+AHFinderDirect::origin_z[1] = +1.0
+AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = 0.0
+AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.0
+AHFinderDirect::initial_guess__coord_sphere__z_center[1] = +1.0
+AHFinderDirect::initial_guess__coord_sphere__radius[1] = 0.5
+
+AHFinderDirect::origin_x[2] = 0.0
+AHFinderDirect::origin_y[2] = 0.0
+AHFinderDirect::origin_z[2] = -1.0
+AHFinderDirect::initial_guess_method[2] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[2] = 0.0
+AHFinderDirect::initial_guess__coord_sphere__y_center[2] = 0.0
+AHFinderDirect::initial_guess__coord_sphere__z_center[2] = -1.0
+AHFinderDirect::initial_guess__coord_sphere__radius[2] = 0.25
diff --git a/param.ccl b/param.ccl
index 4c4e2e8..3c2a842 100644
--- a/param.ccl
+++ b/param.ccl
@@ -528,6 +528,26 @@ boolean mask_is_noshrink \
 } "true"
 
 #
+# If we're using mesh refinement (eg Carpet) it's possible to find a
+# horizon on a fnie grid, and then set a mask on a coarser grid where
+# the masked-off region is only a few grid points across.  This then
+# confuses various other thorns (eg LegoExcision) which require excised
+# regions to be "not too small".
+#
+# So...  This parameter may be used to specify that the mask not be set
+# for a given horizon (just as if that horizon wasn't found) on any grid
+# where the ratio
+#	BH_diagnostics.min_radius / max(Delta_x,Delta_y,Delta_z)
+# (where  Delta_[xyz]  are the *LOCAL* grid spacings on the current grid)
+# is too small (in particular, < this parameter).  
+#
+real min_horizon_size_for_mask \
+  "only set mask if min horizon radius >= this number of local grid spacings"
+{
+0.0:* :: "any real number >= 0.0"
+} 0.0
+
+#
 # This thorn supports two types of mask grid functions:
 # * the "old" style is that defined by CactusEinstein/AHFinder:
 #   the mask is a CCTK_REAL grid function, which is set to specified
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 04665b9..0ef2126 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -195,6 +195,7 @@ struct	mask_info
 	CCTK_REAL radius_multiplier, radius_offset;
 	CCTK_REAL buffer_thickness;
 	bool mask_is_noshrink;
+	CCTK_REAL min_horizon_size_for_mask;
 	bool set_old_style_mask, set_new_style_mask;
 	struct	old_style_mask_info
 		{
diff --git a/src/driver/mask.cc b/src/driver/mask.cc
index 129e035..bfcb9fd 100644
--- a/src/driver/mask.cc
+++ b/src/driver/mask.cc
@@ -79,12 +79,16 @@ struct	mask_grid_info
 	cGH *GH;			// --> Cactus grid hierarchy
 
 	// Cactus coordinate system
+	// (for the *current* grid if we are doing mesh refinement)
 	fp proc_coord_origin[N_GRID_DIMS];	// global (x,y,z) of
 						// *this processor's*
 						// (i,j,k) = (0,0,0)
 						// grid point
 	fp coord_delta[N_GRID_DIMS];		// (x,y,z) grid spacing
 
+	// maximum of x,y,z grid spacings
+	fp max_coord_delta;
+
 	// geometric mean of x,y,z grid spacings,
 	// on the *base* grid if we are doing mesh refinement
 	// (we need the base-grid semantics to make excision consistent
@@ -200,6 +204,11 @@ mgi.GH = cctkGH;
 mgi.coord_delta[X_AXIS] = CCTK_DELTA_SPACE(X_AXIS);
 mgi.coord_delta[Y_AXIS] = CCTK_DELTA_SPACE(Y_AXIS);
 mgi.coord_delta[Z_AXIS] = CCTK_DELTA_SPACE(Z_AXIS);
+mgi.max_coord_delta = jtutil::max(mgi.coord_delta[X_AXIS],
+		      jtutil::max(mgi.coord_delta[Y_AXIS],
+				  mgi.coord_delta[Z_AXIS]));
+
+// Cactus grid spacings on the *base* grid
 const fp base_grid_delta_product =   cctk_delta_space[X_AXIS]
 				   * cctk_delta_space[Y_AXIS]
 				   * cctk_delta_space[Z_AXIS];
@@ -391,22 +400,41 @@ set_mask_gridfn_to_outside_value(mgi,
 //
 // loop over each horizon's xyz bounding box
 // (intersected with (this processor's chunk of) the grid)
-// to set the mask accurately
+// to set the mask accurately (or skip this horizon if it's too small)
 //
 	for (int hn = 1 ; hn <= N_horizons ; ++hn)
 	{
 	const struct AH_data& AH_data = *AH_data_array[hn];
 	if (! AH_data.found_flag)
-	   then continue;			// *** LOOP CONTROL ***
+	   then continue;				// *** LOOP CONTROL ***
+
+	const patch_system& ps                      = *AH_data.ps_ptr;
+	const struct BH_diagnostics& BH_diagnostics =  AH_data.BH_diagnostics;
+
+	//
+	// skip this horizon if it's too small
+	//
+	if (BH_diagnostics.min_radius
+	    < mask_info.min_horizon_size_for_mask*mgi.max_coord_delta)
+	   then {
+		if (verbose_info.print_algorithm_details)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+   "   min_radius=%g < %g grid points ==> skipping mask for horizon %d",
+				   BH_diagnostics.min_radius,
+				   mask_info.min_horizon_size_for_mask,
+				   hn);
+		continue;				// *** LOOP CONTROL ***
+		}
+
+	//
+	// get to here ==> normal mask processing for this horizon
+	//
 
 	if (verbose_info.print_algorithm_details)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   setting mask grid function to \"buffer\"/\"inside\" for horizon %d",
+   "   setting mask grid function to \"buffer\"/\"inside\" for horizon %d",
 			   hn);
 
-	const patch_system& ps                      = *AH_data.ps_ptr;
-	const struct BH_diagnostics& BH_diagnostics =  AH_data.BH_diagnostics;
-
 
 	// horizon bounding box, rounded "out" to the next grid point
 	const int AH_min_i = mgi.ijk_floor_of_global_xyz(X_AXIS, BH_diagnostics.min_x);
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 42a5bfa..795010b 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -306,6 +306,7 @@ if (mask_info.set_mask)
 	mask_info.radius_offset      = mask_radius_offset;
 	mask_info.buffer_thickness   = mask_buffer_thickness;
 	mask_info.mask_is_noshrink   = mask_is_noshrink;
+	mask_info.min_horizon_size_for_mask = min_horizon_size_for_mask;
 	mask_info.set_old_style_mask = (set_old_style_mask != 0);
 	mask_info.set_new_style_mask = (set_new_style_mask != 0);
 	if (mask_info.set_old_style_mask)