diff options
-rw-r--r-- | src/gr/horizon_function.cc | 208 |
1 files changed, 120 insertions, 88 deletions
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 88ebfc7..0821eaf 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -563,10 +563,10 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) then error_exit(ERROR_EXIT, "***** interpolate_geometry():\n" " point out of range when interpolating geometry info from 3-D grid!\n" -" ==> the trial horizon surface is (at least partially)\n" -" outside the grid and/or in an excised region!\n" -" (unable to get info about which point is out of range:\n" -" maybe an interpolator problem?)\n"); /*NOTREACHED*/ +" (one or more points on the trial horizon surface\n" +" are outside the grid or too close to the boundary),\n" +" but we can't get info about the out-of-range point!\n" +" ==> maybe an interpolator problem?\n"); /*NOTREACHED*/ assert(out_of_range_pt >= 0); assert(out_of_range_pt < ps.N_grid_points()); @@ -588,17 +588,12 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) then { CCTK_VWarn(horizon_function_warning_levels::failure, __LINE__, __FILE__, CCTK_THORNSTRING, - "the trial-horizon-surface point"); - CCTK_VWarn(horizon_function_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - " (%g,%g,%g)", - global_x, global_y, global_z); - CCTK_VWarn(horizon_function_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "is outside the grid, or too close to the boundary"); - CCTK_VWarn(horizon_function_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "in the %c%c direction!", +"\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!" + , + global_x, global_y, global_z, end, axis); } return false; // *** ERROR RETURN *** @@ -778,12 +773,17 @@ if (print_msg_flag) then { const fp rho = p.rho_of_irho(irho); 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(horizon_function_warning_levels::nonfinite, __LINE__, __FILE__, CCTK_THORNSTRING, - "h=%g isn't finite at %s patch rho=%g sigma=%g!", +"\n" +" h=%g isn't finite!\n" +" %s patch (rho,sigma)=(%g,%g) (drho,dsigma)=(%g,%g)\n" + , double(h), - p.name(), - double(rho), double(sigma)); + p.name(), double(rho), double(sigma), + double(drho), double(dsigma)); return false; // *** found a NaN *** } } @@ -841,56 +841,58 @@ if (print_msg_flag) 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 (!finite(g_dd_11) || !finite(g_dd_12) || !finite(g_dd_13) - || !finite(g_dd_22) || !finite(g_dd_23) - || !finite(g_dd_33)) - then { - const fp rho = p.rho_of_irho(irho); - const fp sigma = p.sigma_of_isigma(isigma); - CCTK_VWarn(horizon_function_warning_levels::nonfinite, - __LINE__, __FILE__, CCTK_THORNSTRING, - "g_ij isn't finite at %s patch rho=%g sigma=%g!", - p.name(), double(rho), double(sigma)); - return false; // *** 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 ( !finite(partial_d_g_dd_111) + 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); + + 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 ( !finite(g_dd_11) || !finite(g_dd_12) || !finite(g_dd_13) + || !finite(g_dd_22) || !finite(g_dd_23) + || !finite(g_dd_33) + || !finite(K_dd_11) || !finite(K_dd_12) || !finite(K_dd_13) + || !finite(K_dd_22) || !finite(K_dd_23) + || !finite(K_dd_33) + || !finite(partial_d_g_dd_111) || !finite(partial_d_g_dd_112) || !finite(partial_d_g_dd_113) || !finite(partial_d_g_dd_122) @@ -909,31 +911,61 @@ if (print_msg_flag) || !finite(partial_d_g_dd_323) || !finite(partial_d_g_dd_333) ) then { + const fp h = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); const fp rho = p.rho_of_irho(irho); const fp sigma = p.sigma_of_isigma(isigma); + const fp drho = jtutil::degrees_of_radians(rho); + const fp dsigma = jtutil::degrees_of_radians(sigma); + fp local_x, local_y, local_z; + p.xyz_of_r_rho_sigma(h,rho,sigma, local_x,local_y,local_z); + 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(horizon_function_warning_levels::nonfinite, __LINE__, __FILE__, CCTK_THORNSTRING, - "partial_k g_ij isn't finite at %s patch rho=%g sigma=%g!", - p.name(), double(rho), double(sigma)); - return false; // *** 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 (!finite(K_dd_11) || !finite(K_dd_12) || !finite(K_dd_13) - || !finite(K_dd_22) || !finite(K_dd_23) - || !finite(K_dd_33)) - then { - const fp rho = p.rho_of_irho(irho); - const fp sigma = p.sigma_of_isigma(isigma); - CCTK_VWarn(horizon_function_warning_levels::nonfinite, - __LINE__, __FILE__, CCTK_THORNSTRING, - "K_ij isn't finite at %s patch rho=%g sigma=%g!", - p.name(), double(rho), double(sigma)); +"\n" +" geometry isn't finite at %s patch\n" +" h=%g (rho,sigma)=(%g,%g) (drho,dsigma)=(%g,%g)\n" +" local_(x,y,z)=(%g,%g,%g)\n" +" global_(x,y,z)=(%g,%g,%g)\n" +" g_dd_11=%g _12=%g _13=%g\n" +" _22=%g _23=%g _33=%g\n" +" K_dd_11=%g _12=%g _13=%g\n" +" _22=%g _23=%g _33=%g\n" +" partial_d_g_dd_111=%g _112=%g _113=%g\n" +" _122=%g _123=%g _133=%g\n" +" partial_d_g_dd_211=%g _212=%g _213=%g\n" +" _222=%g _223=%g _233=%g\n" +" partial_d_g_dd_311=%g _312=%g _313=%g\n" +" _322=%g _323=%g _333=%g\n" + , + p.name(), + double(h), double(rho), double(sigma), + double(drho), double(dsigma), + double(local_x), double(local_y), double(local_z), + double(global_x), double(global_y), double(global_z), + double(g_dd_11), double(g_dd_12), double(g_dd_13), + double(g_dd_22), double(g_dd_23), double(g_dd_33), + double(K_dd_11), double(K_dd_12), double(K_dd_13), + double(K_dd_22), double(K_dd_23), double(K_dd_33), + double(partial_d_g_dd_111), + double(partial_d_g_dd_112), + double(partial_d_g_dd_113), + double(partial_d_g_dd_122), + double(partial_d_g_dd_123), + double(partial_d_g_dd_133), + double(partial_d_g_dd_211), + double(partial_d_g_dd_212), + double(partial_d_g_dd_213), + double(partial_d_g_dd_222), + double(partial_d_g_dd_223), + double(partial_d_g_dd_233), + double(partial_d_g_dd_311), + double(partial_d_g_dd_312), + double(partial_d_g_dd_313), + double(partial_d_g_dd_322), + double(partial_d_g_dd_323), + double(partial_d_g_dd_333)); return false; // *** found a NaN *** } } |