diff options
| -rw-r--r-- | param.ccl | 14 | ||||
| -rw-r--r-- | src/driver/setup.cc | 3 | ||||
| -rw-r--r-- | src/gr/gr.hh | 4 | ||||
| -rw-r--r-- | src/gr/horizon_function.cc | 260 |
4 files changed, 259 insertions, 22 deletions
@@ -502,6 +502,20 @@ real geometry__Schwarzschild_EF__Delta_xyz \ (0.0:* :: "any real number > 0" } 1.0e-6 +######################################## + +# these are pretty cheap tests to catch assorted wierdness, so it's +# probably worth leaving them on unless you're trying to squeeze every +# last nanosecond... +boolean check_h_for_NaNs \ + "should we check the horizon shape function h for NaNs?" +{ +} "true" +boolean check_geometry_for_NaNs \ + "should we check the interpolated geometry variables for NaNs?" +{ +} "true" + ################################################################################ # diff --git a/src/driver/setup.cc b/src/driver/setup.cc index 34e180c..958676d 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -213,6 +213,9 @@ gi.geometry__Schwarzschild_EF__z_posn = geometry__Schwarzschild_EF__z_posn; gi.geometry__Schwarzschild_EF__epsilon = geometry__Schwarzschild_EF__epsilon; gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz; +gi.check_h_for_NaNs = (check_h_for_NaNs != 0); +gi.check_geometry_for_NaNs = (check_geometry_for_NaNs != 0); + // // create AH-specific info for each AH diff --git a/src/gr/gr.hh b/src/gr/gr.hh index 6c54ffe..9109e7b 100644 --- a/src/gr/gr.hh +++ b/src/gr/gr.hh @@ -99,6 +99,10 @@ struct geometry_info // (x^2+y^2)/r^2 <= this fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD // computation of partial_k g_ij + + // should we check various gridfns for NaNs? + bool check_h_for_NaNs; + bool check_geometry_for_NaNs; }; // diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index bc49a1b..92bc397 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -8,6 +8,10 @@ /// setup_xyz_posns - setup global xyz posns of grid points /// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid /// convert_conformal_to_physical - convert conformal gij to physical +/// +/// h_contains_NaNs - does the horizon shape function h contain any NaNs? +/// geometry_contains_NaNs - do the geometry variables contain any NaNs? +/// /// compute_H - compute H(h) given earlier setup /// @@ -46,6 +50,8 @@ using jtutil::pow4; #include "gr.hh" //****************************************************************************** +//****************************************************************************** +//****************************************************************************** // // ***** prototypes for functions local to this file ***** @@ -62,6 +68,10 @@ void Schwarzschild_EF_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, bool print_msg_flag); + +bool h_contains_NaNs(patch_system& ps, bool print_msg_flag); +bool geometry_contains_NaNs(patch_system& ps, bool print_msg_flag); + void compute_H(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* H_norms_ptr, @@ -69,23 +79,7 @@ void compute_H(patch_system& ps, } //****************************************************************************** - -// -// This function decodes the geometry_method parameter (string) into -// an internal enum for future use. -// -enum geometry_method - decode_geometry_method(const char geometry_method_string[]) -{ -if (STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid")) - then return geometry__interpolate_from_Cactus_grid; -else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF")) - then return geometry__Schwarzschild_EF; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, -"decode_geometry_method(): unknown geometry_method_string=\"%s\"!", - geometry_method_string); /*NOTREACHED*/ -} - +//****************************************************************************** //****************************************************************************** // @@ -132,9 +126,11 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // // Results: // This function returns true for a successful computation, or false -// if the computation failed because the geometry interpolation would -// need data outside the Cactus grid, or data from an excised region. -// FIXME: excision isn't implemented yet :( +// if the computation failed for any of the following reasons: +// - the geometry interpolation would need data outside the Cactus grid, +// or from an excised region +// FIXME: excision isn't implemented yet :( +// - NaNs are found in h or in the interpolate geometry variables // bool horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, @@ -149,6 +145,10 @@ if (print_msg_flag) // fill in values of all ghosted gridfns in ghost zones ps.synchronize(); +if (gi.check_h_for_NaNs + && h_contains_NaNs(ps, print_msg_flag)) + then return false; // *** ERROR RETURN *** + // set up xyz positions of grid points setup_xyz_posns(ps, print_msg_flag); @@ -156,7 +156,7 @@ setup_xyz_posns(ps, print_msg_flag); switch (gi.geometry_method) { case geometry__interpolate_from_Cactus_grid: - if (! interpolate_geometry(ps, cgi, gi, print_msg_flag)) + if (!interpolate_geometry(ps, cgi, gi, print_msg_flag)) then return false; // *** ERROR RETURN *** if (cgi.Cactus_conformal_metric) then convert_conformal_to_physical(ps, print_msg_flag); @@ -172,6 +172,10 @@ default: int(gi.geometry_method)); /*NOTREACHED*/ } +if (gi.check_geometry_for_NaNs + && geometry_contains_NaNs(ps, print_msg_flag)) + then return false; // *** ERROR RETURN *** + // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag); @@ -180,6 +184,24 @@ return true; // *** NORMAL RETURN *** } //****************************************************************************** + +// +// This function decodes the geometry_method parameter (string) into +// an internal enum for future use. +// +enum geometry_method + decode_geometry_method(const char geometry_method_string[]) +{ +if (STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid")) + then return geometry__interpolate_from_Cactus_grid; +else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF")) + then return geometry__Schwarzschild_EF; +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +"decode_geometry_method(): unknown geometry_method_string=\"%s\"!", + geometry_method_string); /*NOTREACHED*/ +} + +//****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -710,6 +732,197 @@ if (print_msg_flag) } //****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function whether the horizon shape function h (nominal grid only) +// contain any NaNs. +// +// Results: +// This function returns true if any NaNs are found, false otherwise. +// +namespace { +bool h_contains_NaNs(patch_system& ps, bool print_msg_flag) +{ +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " checking h for NaNs"); + +#ifdef HAVE_ISNAN + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + const fp h = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + if (isnan(h)) + then { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "h = NaN at %s patch rho=%g sigma=%g!", + p.name(), double(rho), double(sigma)); + return true; // *** found a NaN *** + } + } + } + } +return false; // *** no NaNs found *** +#else +CCTK_VWarn(3, __LINE__, __FILE__, CCTK_THORNSTRING, + " no isnan() ==> skipping NaN check of h"); +return false; // *** no check possible *** +#endif +} + } + +//****************************************************************************** + +// +// This function whether the geometry variables +// g_dd_ij +// partial_d_g_dd_kij +// K_dd_ij +// contain any NaNs. +// +// Results: +// This function returns true if any NaNs are found, false otherwise. +// +namespace { +bool geometry_contains_NaNs(patch_system& ps, bool print_msg_flag) +{ +if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " checking geometry for NaNs"); + +#ifdef HAVE_ISNAN + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + const fp g_dd_11 = p.gridfn(gfns::gfn__g_dd_11, irho,isigma); + const fp g_dd_12 = p.gridfn(gfns::gfn__g_dd_12, irho,isigma); + const fp g_dd_13 = p.gridfn(gfns::gfn__g_dd_13, irho,isigma); + const fp g_dd_22 = p.gridfn(gfns::gfn__g_dd_22, irho,isigma); + const fp g_dd_23 = p.gridfn(gfns::gfn__g_dd_23, irho,isigma); + const fp g_dd_33 = p.gridfn(gfns::gfn__g_dd_33, irho,isigma); + if (isnan(g_dd_11) || isnan(g_dd_12) || isnan(g_dd_13) + || isnan(g_dd_22) || isnan(g_dd_23) + || isnan(g_dd_33)) + then { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "NaN in g_ij at %s patch rho=%g sigma=%g!", + p.name(), double(rho), double(sigma)); + return true; // *** found a NaN *** + } + + const fp partial_d_g_dd_111 = p.gridfn(gfns::gfn__partial_d_g_dd_111, + irho,isigma); + const fp partial_d_g_dd_112 = p.gridfn(gfns::gfn__partial_d_g_dd_112, + irho,isigma); + const fp partial_d_g_dd_113 = p.gridfn(gfns::gfn__partial_d_g_dd_113, + irho,isigma); + const fp partial_d_g_dd_122 = p.gridfn(gfns::gfn__partial_d_g_dd_122, + irho,isigma); + const fp partial_d_g_dd_123 = p.gridfn(gfns::gfn__partial_d_g_dd_123, + irho,isigma); + const fp partial_d_g_dd_133 = p.gridfn(gfns::gfn__partial_d_g_dd_133, + irho,isigma); + const fp partial_d_g_dd_211 = p.gridfn(gfns::gfn__partial_d_g_dd_211, + irho,isigma); + const fp partial_d_g_dd_212 = p.gridfn(gfns::gfn__partial_d_g_dd_212, + irho,isigma); + const fp partial_d_g_dd_213 = p.gridfn(gfns::gfn__partial_d_g_dd_213, + irho,isigma); + const fp partial_d_g_dd_222 = p.gridfn(gfns::gfn__partial_d_g_dd_222, + irho,isigma); + const fp partial_d_g_dd_223 = p.gridfn(gfns::gfn__partial_d_g_dd_223, + irho,isigma); + const fp partial_d_g_dd_233 = p.gridfn(gfns::gfn__partial_d_g_dd_233, + irho,isigma); + const fp partial_d_g_dd_311 = p.gridfn(gfns::gfn__partial_d_g_dd_311, + irho,isigma); + const fp partial_d_g_dd_312 = p.gridfn(gfns::gfn__partial_d_g_dd_312, + irho,isigma); + const fp partial_d_g_dd_313 = p.gridfn(gfns::gfn__partial_d_g_dd_313, + irho,isigma); + const fp partial_d_g_dd_322 = p.gridfn(gfns::gfn__partial_d_g_dd_322, + irho,isigma); + const fp partial_d_g_dd_323 = p.gridfn(gfns::gfn__partial_d_g_dd_323, + irho,isigma); + const fp partial_d_g_dd_333 = p.gridfn(gfns::gfn__partial_d_g_dd_333, + irho,isigma); + if ( isnan(partial_d_g_dd_111) + || isnan(partial_d_g_dd_112) + || isnan(partial_d_g_dd_113) + || isnan(partial_d_g_dd_122) + || isnan(partial_d_g_dd_123) + || isnan(partial_d_g_dd_133) + || isnan(partial_d_g_dd_211) + || isnan(partial_d_g_dd_212) + || isnan(partial_d_g_dd_213) + || isnan(partial_d_g_dd_222) + || isnan(partial_d_g_dd_223) + || isnan(partial_d_g_dd_233) + || isnan(partial_d_g_dd_311) + || isnan(partial_d_g_dd_312) + || isnan(partial_d_g_dd_313) + || isnan(partial_d_g_dd_322) + || isnan(partial_d_g_dd_323) + || isnan(partial_d_g_dd_333) ) + then { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "NaN in partial_k g_ij at %s patch rho=%g sigma=%g!", + p.name(), double(rho), double(sigma)); + return true; // *** found a NaN *** + } + + const fp K_dd_11 = p.gridfn(gfns::gfn__K_dd_11, irho,isigma); + const fp K_dd_12 = p.gridfn(gfns::gfn__K_dd_12, irho,isigma); + const fp K_dd_13 = p.gridfn(gfns::gfn__K_dd_13, irho,isigma); + const fp K_dd_22 = p.gridfn(gfns::gfn__K_dd_22, irho,isigma); + const fp K_dd_23 = p.gridfn(gfns::gfn__K_dd_23, irho,isigma); + const fp K_dd_33 = p.gridfn(gfns::gfn__K_dd_33, irho,isigma); + if (isnan(K_dd_11) || isnan(K_dd_12) || isnan(K_dd_13) + || isnan(K_dd_22) || isnan(K_dd_23) + || isnan(K_dd_33)) + then { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "NaN in K_ij at %s patch rho=%g sigma=%g!", + p.name(), double(rho), double(sigma)); + return true; // *** found a NaN *** + } + } + } + } +return false; // *** no NaNs found *** +#else +CCTK_VWarn(3, __LINE__, __FILE__, CCTK_THORNSTRING, + " no isnan() ==> skipping NaN check of geometry"); +return false; // *** no check possible *** +#endif +} + } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** // // This function computes H(h), and optionally its Jacobian coefficients, @@ -717,6 +930,9 @@ if (print_msg_flag) // uses a mixture of algebraic operations and (rho,sigma) finite differencing. // The computation is done (entirely) on the nominal angular grid. // +// N.b. This function #includes "cg.hh", which defines "dangerous" macros +// which will stay in effect for the rest of this compilation unit! +// // Arguments: // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. @@ -778,7 +994,7 @@ if (print_msg_flag) // and so must be inside its own set of { } braces // - // gridfn #defins + // gridfn #defines #include "cg.hh" { |