diff options
Diffstat (limited to 'src/gr/expansion.cc')
-rw-r--r-- | src/gr/expansion.cc | 109 |
1 files changed, 70 insertions, 39 deletions
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc index 490ff99..1946cfa 100644 --- a/src/gr/expansion.cc +++ b/src/gr/expansion.cc @@ -59,20 +59,25 @@ void setup_xyz_posns(patch_system& ps, bool print_msg_flag); bool interpolate_geometry(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag); -void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag); +void convert_conformal_to_physical(patch_system& ps, + bool print_msg_flag); void Schwarzschild_EF_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, bool print_msg_flag); -bool h_is_finite(patch_system& ps, bool print_msg_flag); +bool h_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); bool geometry_is_finite(patch_system& ps, - bool print_msg_flag); + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); -bool compute_Theta(patch_system& ps, - bool Jacobian_flag, +bool compute_Theta(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag); } @@ -122,6 +127,11 @@ bool compute_Theta(patch_system& ps, // interpolator call are done, the latter with the number of // interpolation points is set to 0 and all the output array // pointers set to NULL. +// initial_flag = true if this is the first evaluation of expansion() +// for this horizon, +// false otherwise; +// this is used (only) to select which elements of error_info +// are relevant // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. // print_msg_flag = true to print status messages, @@ -145,6 +155,7 @@ bool compute_Theta(patch_system& ps, bool expansion(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool Jacobian_flag /* = false */, bool print_msg_flag /* = false */, jtutil::norm<fp>* Theta_norms_ptr /* = NULL */) @@ -164,8 +175,9 @@ if (active_flag) // fill in values of all ghosted gridfns in ghost zones ps_ptr->synchronize(); - if (gi.check_that_h_is_finite - && !h_is_finite(*ps_ptr, print_msg_flag)) + if (gi.check_that_h_is_finite && !h_is_finite(*ps_ptr, + error_info, initial_flag, + print_msg_flag)) then return false; // *** ERROR RETURN *** // set up xyz positions of grid points @@ -181,15 +193,19 @@ case geometry__local_interp_from_Cactus_grid: // ps_ptr (non-NULL vs NULL) to choose a normal vs dummy computation if (!interpolate_geometry(ps_ptr, cgi, gi, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** if (active_flag && cgi.Cactus_conformal_metric) - then convert_conformal_to_physical(*ps_ptr, print_msg_flag); + then convert_conformal_to_physical(*ps_ptr, + print_msg_flag); break; case geometry__Schwarzschild_EF: if (active_flag) - then Schwarzschild_EF_geometry(*ps_ptr, gi, print_msg_flag); + then Schwarzschild_EF_geometry(*ps_ptr, + gi, + print_msg_flag); break; default: @@ -203,7 +219,9 @@ default: if (active_flag) then { if (gi.check_that_geometry_is_finite - && !geometry_is_finite(*ps_ptr, print_msg_flag)) + && !geometry_is_finite(*ps_ptr, + error_info, initial_flag, + print_msg_flag)) then return false; // *** ERROR RETURN *** // compute remaining gridfns --> $\Theta$ @@ -211,6 +229,7 @@ if (active_flag) // by algebraic ops and angular finite differencing if (!compute_Theta(*ps_ptr, Jacobian_flag, Theta_norms_ptr, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** } @@ -322,6 +341,7 @@ namespace { bool interpolate_geometry(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { const bool active_flag = (ps_ptr != NULL); @@ -651,19 +671,26 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE) const char axis = "xyz"[error_axis]; assert((error_direction == -1) || (error_direction == +1)); - const char end = (error_direction == -1) ? '-' : '+'; + const char direction = (error_direction == -1) ? '-' : '+'; if (print_msg_flag) then { - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, + const int warn_level + = initial_flag + ? error_info.warn_level__point_outside__initial + : error_info.warn_level__point_outside__subsequent; + CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" the trial-horizon-surface point" -" global_(x,y,z)=(%g,%g,%g)" -" is outside the grid (or too close to the boundary) in the %c%c direction!" +"interpolate_geometry():\n" +" the trial-horizon-surface point\n" +" global_(x,y,z)=(%g,%g,%g)\n" +" is outside the grid (or too close to the boundary)" + " in the %c%c direction!\n" +" (this may mean that driver::ghost_size is too small,\n" +" especially if any horizon crosses a symmetry boundary)" , global_x, global_y, global_z, - end, axis); + direction, axis); } return false; // *** ERROR RETURN *** } @@ -811,7 +838,7 @@ if (print_msg_flag) // or more NaNs or infinities. // // Results: -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE // This function returns true if all the h values are finite, false // otherwise (i.e. if it contains any NaNs or infinities). #else @@ -820,12 +847,14 @@ if (print_msg_flag) #endif // namespace { -bool h_is_finite(patch_system& ps, bool print_msg_flag) +bool h_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " checking that h is finite"); -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -843,7 +872,7 @@ if (print_msg_flag) const fp sigma = p.sigma_of_isigma(isigma); const fp drho = jtutil::degrees_of_radians(rho); const fp dsigma = jtutil::degrees_of_radians(sigma); - CCTK_VWarn(expansion_warning_levels::nonfinite, + CCTK_VWarn(error_info.warn_level__nonfinite_geometry, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " h=%g isn't finite!\n" @@ -859,7 +888,7 @@ if (print_msg_flag) } return true; // *** all values finite *** #else -CCTK_VWarn(expansion_warning_levels::skip_nonfinite, +CCTK_VWarn(error_info.warn_level__skipping_finite_check, __LINE__, __FILE__, CCTK_THORNSTRING, " no finite() fn ==> skipping is-h-finite check"); return true; // *** no check possible *** @@ -878,7 +907,7 @@ return true; // *** no check possible *** // or more NaNs or infinities. // // Results: -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE // This function returns true if all the geometry variables are finite, // false otherwise (i.e. if they contain any NaNs or infinities). #else @@ -887,12 +916,14 @@ return true; // *** no check possible *** #endif // namespace { -bool geometry_is_finite(patch_system& ps, bool print_msg_flag) +bool geometry_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " checking that geometry is finite"); -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -989,7 +1020,7 @@ if (print_msg_flag) const fp global_x = ps.global_x_of_local_x(local_x); const fp global_y = ps.global_y_of_local_y(local_y); const fp global_z = ps.global_z_of_local_z(local_z); - CCTK_VWarn(expansion_warning_levels::nonfinite, + CCTK_VWarn(error_info.warn_level__nonfinite_geometry, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " geometry isn't finite at %s patch\n" @@ -1041,7 +1072,7 @@ if (print_msg_flag) } return true; // *** no NaNs found *** #else -CCTK_VWarn(expansion_warning_levels::skip_nonfinite, +CCTK_VWarn(error_info.warn_level__skipping_finite_check, __LINE__, __FILE__, CCTK_THORNSTRING, " no finite() ==> skipping is-geometry-finite check"); return true; // *** no check possible *** @@ -1073,9 +1104,9 @@ return true; // *** no check possible *** // g_ij isn't positive definite). // namespace { -bool compute_Theta(patch_system& ps, - bool Jacobian_flag, +bool compute_Theta(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { if (print_msg_flag) @@ -1159,16 +1190,16 @@ if (print_msg_flag) if (Theta_D <= 0) then { - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "Theta_D <= 0 at %s patch rho=%g sigma=%g!", - p.name(), double(rho), double(sigma)); - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "(this probably means the interpolated g_ij"); - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - " isn't positive definite)"); +const int warn_level + = initial_flag ? error_info.warn_level__gij_not_positive_definite__initial + : error_info.warn_level__gij_not_positive_definite__subsequent; +CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" compute_Theta(): Theta_D = $g^{ij} s_i s_j$ = %g <= 0\n" +" at %s patch rho=%g sigma=%g!\n" +" (i.e. the interpolated g_ij isn't positive definite)", + double(Theta_D), + p.name(), double(rho), double(sigma)); return false; // *** ERROR RETURN *** } |