diff options
| -rw-r--r-- | param.ccl | 101 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset/large-Jacobian.par | 4 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset/large-ahf.par | 1 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset/large.par | 5 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset/small-Jacobian.par | 4 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset/small.par | 4 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Schw/large-ahf.par | 1 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Schw/offset.par | 10 | ||||
| -rw-r--r-- | src/gr/AHFinderDirect.hh | 46 | ||||
| -rw-r--r-- | src/gr/Newton.cc | 69 | ||||
| -rw-r--r-- | src/gr/Schwarzschild_EF.cc | 33 | ||||
| -rw-r--r-- | src/gr/driver.cc | 47 | ||||
| -rw-r--r-- | src/gr/gr_gfas.minc | 9 | ||||
| -rw-r--r-- | src/gr/horizon_Jacobian.cc | 2 | ||||
| -rw-r--r-- | src/gr/horizon_function.cc | 20 | ||||
| -rw-r--r-- | src/gr/setup_gr_gfas.maple | 7 | ||||
| -rw-r--r-- | src/include/config.hh | 9 | ||||
| -rw-r--r-- | src/patch/patch_system.cc | 2 |
18 files changed, 267 insertions, 107 deletions
@@ -64,6 +64,12 @@ int max_Newton_iterations "maximum number of Newton iterations before giving up" (0:* :: "any positive integer" } 10 +# if this is used, the file names are the usual ones with ".it%d" appended +boolean output_h_and_H_at_each_Newton_iteration \ + "should we output h and H at each Newton iteration (for debugging)?" +{ +} "false" + # # we declare convergence if *either* of the following two criteria are met # @@ -78,6 +84,12 @@ real Delta_h_norm_for_convergence \ (0.0:* :: "any positive real number" } 1.0e-10 +boolean final_H_update_if_exit_x_H_small \ + "should we do a final H(h) update after a h += Delta_h update which is\ + so small it meets the Delta_h_norm_for_convergence convergence criterion?" +{ +} "false" + ################################################################################ # @@ -149,47 +161,81 @@ real delta_drho_dsigma "angular grid spacing of patches, in degrees" ################################################################################ # -# parameters for (optionally) hard-wiring the geometry to -# Schwarzschild spacetime in Eddington-Finkelstein coordinates -# (useful for testing smoothness requirements of interpolation) +# parameters for how we compute the slice's geometry +# (gij, Kij, partial_k gij) # -boolean hardwire_Schwarzschild_EF \ - "should we hard-wire Schwarzschild/EF geometry, \ - bypassing the usual $g_{ij}$/$K_{ij}$ interpolation from the Cactus grid?" +keyword geometry_method "how do we compute the slice's geometry?" { -} "false" +"interpolate from Cactus grid" :: "interpolate gij and Kij from Cactus grid" +"Schwarzschild/EF" :: \ + "hard-wire to Schwarzschild spacetime / Eddington-Finkelstein slice" +} "interpolate from Cactus grid" -real hardwire_Schwarzschild_EF__mass "mass of Schwarzschild BH" +######################################## + +# +# parameters for geometry_method = "interpolate from Cactus grid" +# + +string geometry_interpolator_name \ + "name under which the geometry interpolation operator is registered in Cactus" +{ +.* :: "any string" +} "generalized polynomial interpolation" + +string geometry_interpolator_pars \ + "parameters for the geometry interpolator" +{ +.* :: "any string acceptable to Util_TableSetFromString()" +} "" +######################################## + +# +# parameters for geometry_method = "Schwarzschild/EF" +# + +real geometry__Schwarzschild_EF__mass "mass of Schwarzschild BH" { (0.0:* :: "BH mass = any real number > 0" } 1.0 -real hardwire_Schwarzschild_EF__x_posn "x coordinate of Schwarzschild BH" +real geometry__Schwarzschild_EF__x_posn "x coordinate of Schwarzschild BH" { *:* :: "any real number" } 0.0 -real hardwire_Schwarzschild_EF__y_posn "y coordinate of Schwarzschild BH" +real geometry__Schwarzschild_EF__y_posn "y coordinate of Schwarzschild BH" { *:* :: "any real number" } 0.0 -real hardwire_Schwarzschild_EF__z_posn "z coordinate of Schwarzschild BH" +real geometry__Schwarzschild_EF__z_posn "z coordinate of Schwarzschild BH" { *:* :: "any real number" } 0.0 -real hardwire_Schwarzschild_EF__epsilon \ +# some of the formulas have 0/0 limits on the z axis; this parameter controls +# where we switch from the generic formulas to the L'Hopital's-rule z axis +# limits +# - don't set this parameter too small or roundoff errors will be excessive +# - don't set this parameter too large or finite differencing errors will +# be excessive +# in practice the default value should be fine +# n.b. this is used for centered finite differencing, unlike the Jacobian +real geometry__Schwarzschild_EF__epsilon \ "threshold for sin^2 theta = (x^2+y^2)/r^2 below which we use z axis limits" { -(0.0:* :: "this should be a little bit above the floating-point roundoff level" +(0.0:* :: "this should be somewhat above the floating-point roundoff level" } 1.0e-12 -# notes on this parameter: -# - don't set it too small or roundoff errors will become large -# - don't set it too large or finite differencing errors will become large +# we compute partial_k g_ij by numerical finite differencing of the exact +# analytical g_ij values; this parameter sets the "grid spacing" for this +# - don't set this parameter too small or roundoff errors will be excessive +# - don't set this parameter too large or finite differencing errors will +# be excessive # in practice the default value should be fine -# n.b. this is used for centered finite differencing, unlike the Jacobian -real hardwire_Schwarzschild_EF__Delta_xyz \ +# ... n.b. this finite differencing is *centered*, unlike that in the +# Jacobian computation +real geometry__Schwarzschild_EF__Delta_xyz \ "finite diff pseuo-grid spacing for computing $\partial_k g_{ij}$" { (0.0:* :: "any real number > 0" @@ -198,25 +244,6 @@ real hardwire_Schwarzschild_EF__Delta_xyz \ ################################################################################ # -# parameters for the interpolator used to interpolate the "geometry" -# $g_{ij}$ and $K_{ij}$ to the apparent horizon surface position -# - -string geometry_interpolator_name \ - "name under which the geometry interpolation operator is registered in Cactus" -{ -.* :: "any string" -} "generalized polynomial interpolation" - -string geometry_interpolator_pars \ - "parameters for the geometry interpolator" -{ -.* :: "any string acceptable to Util_TableSetFromString()" -} "" - -################################################################################ - -# # parameters for the interpatch interpolator # string interpatch_interpolator_name \ diff --git a/run/test-ahfinderdirect/Kerr-offset/large-Jacobian.par b/run/test-ahfinderdirect/Kerr-offset/large-Jacobian.par index 63f7c93..9fa37de 100644 --- a/run/test-ahfinderdirect/Kerr-offset/large-Jacobian.par +++ b/run/test-ahfinderdirect/Kerr-offset/large-Jacobian.par @@ -44,5 +44,5 @@ AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpola AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" -AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 -AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::initial_guess__Kerr__mass = 1.0 +AHFinderDirect::initial_guess__Kerr__spin = 0.6 diff --git a/run/test-ahfinderdirect/Kerr-offset/large-ahf.par b/run/test-ahfinderdirect/Kerr-offset/large-ahf.par index c28ebc5..fb5b390 100644 --- a/run/test-ahfinderdirect/Kerr-offset/large-ahf.par +++ b/run/test-ahfinderdirect/Kerr-offset/large-ahf.par @@ -27,6 +27,7 @@ Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::method = "Newton solve" +AHFinderDirect::final_H_update_if_exit_x_H_small = "true" AHFinderDirect::h_file_name = "large-ahf.h.dat" AHFinderDirect::H_of_h_file_name = "large-ahf.H.dat" diff --git a/run/test-ahfinderdirect/Kerr-offset/large.par b/run/test-ahfinderdirect/Kerr-offset/large.par index 6220e09..bef93fc 100644 --- a/run/test-ahfinderdirect/Kerr-offset/large.par +++ b/run/test-ahfinderdirect/Kerr-offset/large.par @@ -27,6 +27,7 @@ Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::method = "horizon function" +AHFinderDirect::final_H_update_if_exit_x_H_small = "true" AHFinderDirect::h_file_name = "large.h.dat" AHFinderDirect::H_of_h_file_name = "large.H.dat" @@ -44,5 +45,5 @@ AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpola AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" -AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 -AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::initial_guess__Kerr__mass = 1.0 +AHFinderDirect::initial_guess__Kerr__spin = 0.6 diff --git a/run/test-ahfinderdirect/Kerr-offset/small-Jacobian.par b/run/test-ahfinderdirect/Kerr-offset/small-Jacobian.par index 953fb07..e7b270c 100644 --- a/run/test-ahfinderdirect/Kerr-offset/small-Jacobian.par +++ b/run/test-ahfinderdirect/Kerr-offset/small-Jacobian.par @@ -44,5 +44,5 @@ AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpola AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" -AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 -AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::initial_guess__Kerr__mass = 1.0 +AHFinderDirect::initial_guess__Kerr__spin = 0.6 diff --git a/run/test-ahfinderdirect/Kerr-offset/small.par b/run/test-ahfinderdirect/Kerr-offset/small.par index 3b7f8f4..6759004 100644 --- a/run/test-ahfinderdirect/Kerr-offset/small.par +++ b/run/test-ahfinderdirect/Kerr-offset/small.par @@ -44,5 +44,5 @@ AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpola AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" -AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 -AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::initial_guess__Kerr__mass = 1.0 +AHFinderDirect::initial_guess__Kerr__spin = 0.6 diff --git a/run/test-ahfinderdirect/Schw/large-ahf.par b/run/test-ahfinderdirect/Schw/large-ahf.par index 536d893..c7daeb7 100644 --- a/run/test-ahfinderdirect/Schw/large-ahf.par +++ b/run/test-ahfinderdirect/Schw/large-ahf.par @@ -26,6 +26,7 @@ Exact::exact_model = "Schwarzschild/EF" Exact::Schwarzschild_EF__mass = 1.0 AHFinderDirect::method = "Newton solve" +AHFinderDirect::final_H_update_if_exit_x_H_small = "true" AHFinderDirect::h_file_name = "large-ahf.h.dat" AHFinderDirect::H_of_h_file_name = "large-ahf.H.dat" diff --git a/run/test-ahfinderdirect/Schw/offset.par b/run/test-ahfinderdirect/Schw/offset.par index 8b32933..c16dcfc 100644 --- a/run/test-ahfinderdirect/Schw/offset.par +++ b/run/test-ahfinderdirect/Schw/offset.par @@ -23,6 +23,12 @@ ADMBase::metric_type = "physical" # Exact Exact::exact_model = "Schwarzschild/EF" +Exact::Schwarzschild_EF__mass = 1.0 + +AHFinderDirect::method = "horizon function" +AHFinderDirect::final_H_update_if_exit_x_H_small = "true" +AHFinderDirect::h_file_name = "offset.h.dat" +AHFinderDirect::H_of_h_file_name = "offset.H.dat" AHFinderDirect::origin_x = 0.5 AHFinderDirect::origin_y = 0.7 @@ -38,4 +44,6 @@ AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpola AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "ellipsoid" -AHFinderDirect::method = "horizon" +AHFinderDirect::initial_guess__ellipsoid__x_radius = 2.0 +AHFinderDirect::initial_guess__ellipsoid__y_radius = 2.0 +AHFinderDirect::initial_guess__ellipsoid__z_radius = 2.0 diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index 403ea39..22e713c 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -20,6 +20,15 @@ enum Jacobian_method symbolic_differentiation // no comma }; +// +// this enum holds the (a) decoded geometry_method parameter, +// i.e. it specifies how we compute the (a) geometry +// +enum geometry_method + { + geometry__interpolate_from_Cactus_grid, + geometry__Schwarzschild_EF // no comma + }; // // This structure holds information for computing the spacetime geometry. @@ -29,31 +38,22 @@ enum Jacobian_method // struct geometry_info { - // - // parameters for hard-wiring Schwarzschild/EF geometry - // - bool hardwire_Schwarzschild_EF; // should we hard-wire the - // Schwarzschild/EF geometry? - fp hardwire_Schwarzschild_EF__x_posn; // x posn of Schwarzschild BH - fp hardwire_Schwarzschild_EF__y_posn; // y posn of Schwarzschild BH - fp hardwire_Schwarzschild_EF__z_posn; // z posn of Schwarzschild BH - fp hardwire_Schwarzschild_EF__mass; // mass of Schwarzschild BH - fp hardwire_Schwarzschild_EF__epsilon; // threshold for sin^2 theta - // = (x^2+y^2)/r^2 below which - // we use z axis limits - fp Delta_xyz; // pseudo-grid spacing for finite differencing - // computation of $\partial_k g_{ij}$ - - // - // parameters for normal interpolation from Cactus grid - // + enum geometry_method geometry_method; + + // parameters for geometry_method = interpolate_from_Cactus_grid int operator_handle; // Cactus handle to interpolation op int param_table_handle; // Cactus handle to parameter table // for the interpolator - // this doesn't really belong in this structure (it doesn't - // have any logical connection to the geometry calculations), - // but it's convenient to put it here anyway... + // parameters for geometry_method = Schwarzschild_EF + fp geometry__Schwarzschild_EF__x_posn; // x posn of Schwarzschild BH + fp geometry__Schwarzschild_EF__y_posn; // y posn of Schwarzschild BH + fp geometry__Schwarzschild_EF__z_posn; // z posn of Schwarzschild BH + fp geometry__Schwarzschild_EF__mass; // mass of Schwarzschild BH + fp geometry__Schwarzschild_EF__epsilon; // use z axis limits if + // (x^2+y^2)/r^2 <= this + fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD + // computation of partial_k g_ij }; // @@ -108,8 +108,12 @@ struct solver_info { // stuff for Newton_solve() int max_Newton_iterations; + bool output_h_and_H_at_each_Newton_iteration; + const char *h_file_name; + const char *H_of_h_file_name; fp H_norm_for_convergence; fp Delta_h_norm_for_convergence; + bool final_H_update_if_exit_x_H_small; }; //****************************************************************************** diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc index d82c752..dd030d7 100644 --- a/src/gr/Newton.cc +++ b/src/gr/Newton.cc @@ -72,6 +72,33 @@ bool Newton_solve(patch_system& ps, " H rms-norm=%.2e, infinity-norm=%.2e", H_norms.rms_norm(), H_norms.infinity_norm()); + if (solver_info.output_h_and_H_at_each_Newton_iteration) + then { + const char N_buffer = 100; + char file_name_buffer[100]; + + snprintf(file_name_buffer, N_buffer, + "%s.it%d", + solver_info.h_file_name, iteration); + CCTK_VInfo(CCTK_THORNSTRING, + " writing h to \"%s\"", + file_name_buffer); + ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h, + true, gfns::gfn__h, + file_name_buffer, + false); // no ghost zones + + snprintf(file_name_buffer, N_buffer, + "%s.it%d", + solver_info.H_of_h_file_name, iteration); + CCTK_VInfo(CCTK_THORNSTRING, + " writing H(h) to \"%s\"", + file_name_buffer); + ps.print_gridfn_with_xyz(gfns::gfn__H, + true, gfns::gfn__h, + file_name_buffer); + } + if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence) then { CCTK_VInfo(CCTK_THORNSTRING, @@ -80,14 +107,27 @@ bool Newton_solve(patch_system& ps, return true; // *** NORMAL RETURN *** } - // take a Newton step + // compute the Newton step horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac); CCTK_VInfo(CCTK_THORNSTRING, - "solving linear system for %d unknowns", Jac.NN()); - Jac.solve_linear_system(gfns::gfn__H, // rhs gridfn - gfns::gfn__Delta_h); // soln gridfn - CCTK_VInfo(CCTK_THORNSTRING, - "done solving linear system"); + "solving linear system for %d unknowns Delta_h", + Jac.NN()); + const fp rcond = Jac.solve_linear_system(gfns::gfn__H, + gfns::gfn__Delta_h); + if (rcond == 0.0) + then { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Newton_solve: Jacobian matrix is numerically singular!"); + return false; // *** ERROR RETURN *** + } + if (rcond > 0.0) + then CCTK_VInfo(CCTK_THORNSTRING, + "done solving linear system (condition number %.2e)", + double(rcond)); + else CCTK_VInfo(CCTK_THORNSTRING, + "done solving linear system"); + + // take the Newton step jtutil::norm<fp> Delta_h_norms; for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { @@ -107,7 +147,7 @@ bool Newton_solve(patch_system& ps, } } CCTK_VInfo(CCTK_THORNSTRING, - "moved h by Delta_h (rhs-norm=%.2e, infinity-norm=%.2e)", + "h += Delta_h (rms-norm=%.2e, infinity-norm=%.2e)", Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm()); if (Delta_h_norms.infinity_norm() <= solver_info @@ -116,12 +156,25 @@ bool Newton_solve(patch_system& ps, CCTK_VInfo(CCTK_THORNSTRING, "==> finished (||Delta_h|| < %g)", double(solver_info.Delta_h_norm_for_convergence)); + if (solver_info.final_H_update_if_exit_x_H_small) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " doing final H(h) evaluation"); + horizon_function(ps, cgi, gi, false, true); + } return true; // *** NORMAL RETURN *** } } + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Newton_solve: no convergence in %d iterations!\n", + "Newton_solve: no convergence in %d iterations!", solver_info.max_Newton_iterations); +if (solver_info.final_H_update_if_exit_x_H_small) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " doing final H(h) evaluation"); + horizon_function(ps, cgi, gi, false, true); + } return false; // *** ERROR RETURN *** } diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc index f0e18d4..1bce005 100644 --- a/src/gr/Schwarzschild_EF.cc +++ b/src/gr/Schwarzschild_EF.cc @@ -107,11 +107,12 @@ void Schwarzschild_EF_geometry(patch_system& ps, const struct geometry_info& gi, bool msg_flag) { -const fp mass = gi.hardwire_Schwarzschild_EF__mass; -const fp epsilon = gi.hardwire_Schwarzschild_EF__epsilon; -const fp x_posn = gi.hardwire_Schwarzschild_EF__x_posn; -const fp y_posn = gi.hardwire_Schwarzschild_EF__y_posn; -const fp z_posn = gi.hardwire_Schwarzschild_EF__z_posn; +const fp mass = gi.geometry__Schwarzschild_EF__mass; +const fp x_posn = gi.geometry__Schwarzschild_EF__x_posn; +const fp y_posn = gi.geometry__Schwarzschild_EF__y_posn; +const fp z_posn = gi.geometry__Schwarzschild_EF__z_posn; +const fp epsilon = gi.geometry__Schwarzschild_EF__epsilon; +const fp Delta_xyz = gi.geometry__Schwarzschild_EF__Delta_xyz; if (msg_flag) then { @@ -165,16 +166,16 @@ if (msg_flag) // compute partial_x g_ij by finite differencing in xyz Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel+gi.Delta_xyz, y_rel, z_rel, + x_rel+Delta_xyz, y_rel, z_rel, g11p, g12p, g13p, g22p, g23p, g33p); Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel-gi.Delta_xyz, y_rel, z_rel, + x_rel-Delta_xyz, y_rel, z_rel, g11m, g12m, g13m, g22m, g23m, g33m); - const fp fx = 1.0 / (2.0*gi.Delta_xyz); + const fp fx = 1.0 / (2.0*Delta_xyz); p.gridfn(gfns::gfn__partial_d_g_dd_111,irho,isigma) = fx * (g11p-g11m); p.gridfn(gfns::gfn__partial_d_g_dd_112,irho,isigma) = fx * (g12p-g12m); p.gridfn(gfns::gfn__partial_d_g_dd_113,irho,isigma) = fx * (g13p-g13m); @@ -184,16 +185,16 @@ if (msg_flag) // compute partial_y g_ij by finite differencing in xyz Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel, y_rel+gi.Delta_xyz, z_rel, + x_rel, y_rel+Delta_xyz, z_rel, g11p, g12p, g13p, g22p, g23p, g33p); Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel, y_rel-gi.Delta_xyz, z_rel, + x_rel, y_rel-Delta_xyz, z_rel, g11m, g12m, g13m, g22m, g23m, g33m); - const fp fy = 1.0 / (2.0*gi.Delta_xyz); + const fp fy = 1.0 / (2.0*Delta_xyz); p.gridfn(gfns::gfn__partial_d_g_dd_211,irho,isigma) = fy * (g11p-g11m); p.gridfn(gfns::gfn__partial_d_g_dd_212,irho,isigma) = fy * (g12p-g12m); p.gridfn(gfns::gfn__partial_d_g_dd_213,irho,isigma) = fy * (g13p-g13m); @@ -203,16 +204,16 @@ if (msg_flag) // compute partial_x g_ij by finite differencing in xyz Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel, y_rel, z_rel+gi.Delta_xyz, + x_rel, y_rel, z_rel+Delta_xyz, g11p, g12p, g13p, g22p, g23p, g33p); Schwarzschild_EF_gij_xyz(mass, epsilon, - x_rel, y_rel, z_rel-gi.Delta_xyz, + x_rel, y_rel, z_rel-Delta_xyz, g11m, g12m, g13m, g22m, g23m, g33m); - const fp fz = 1.0 / (2.0*gi.Delta_xyz); + const fp fz = 1.0 / (2.0*Delta_xyz); p.gridfn(gfns::gfn__partial_d_g_dd_311,irho,isigma) = fz * (g11p-g11m); p.gridfn(gfns::gfn__partial_d_g_dd_312,irho,isigma) = fz * (g12p-g12m); p.gridfn(gfns::gfn__partial_d_g_dd_313,irho,isigma) = fz * (g13p-g13m); @@ -336,7 +337,7 @@ K_theta_theta = 2.0 * temp; // Arguments: // epsilon = Tolerance parameter for deciding when to switch from // generic expressions to z-axis limits: we switch iff -// sin^2 theta = (x^2+y^2)/r^2 is < epsilon +// sin^2 theta = (x^2+y^2)/r^2 is <= epsilon // namespace { void xform_from_rthetaphi_to_xyz(fp epsilon, @@ -351,7 +352,7 @@ const fp y2 = y*y; const fp x2py2 = x2 + y2; const fp z2 = z*z; const fp r2 = x2 + y2 + z2; const fp r4 = r2*r2; -const bool z_flag = x2py2/r2 < epsilon; +const bool z_flag = x2py2/r2 <= epsilon; T_xx = (x2/r2)*T_rr + (z_flag ? 1.0/r2 : (x2*z2/r2 + y2) / (r2*x2py2)) *T_theta_theta; diff --git a/src/gr/driver.cc b/src/gr/driver.cc index a9b9a13..eb4d629 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -4,6 +4,7 @@ // <<<prototypes for functions local to this file>>> // AHFinderDirect_driver - top-level driver /// decode_Jacobian_method - decode the Jacobian_method parameter +/// decode_geometry_method - decode the geometry_method parameter /// /// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords) /// setup_ellipsoid - setup up a coordiante ellipsoid in h @@ -52,6 +53,9 @@ using jtutil::error_exit; namespace { enum Jacobian_method decode_Jacobian_method(const char Jacobian_method_string[]); +enum geometry_method + decode_geometry_method(const char geometry_method_string[]); + void setup_Kerr_horizon(patch_system& ps, fp x_posn, fp y_posn, fp z_posn, fp m, fp a, @@ -83,14 +87,9 @@ CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures"); // set up the geometry parameters and the geometry interpolator // struct geometry_info gi; -gi.hardwire_Schwarzschild_EF = (hardwire_Schwarzschild_EF != 0); -gi.hardwire_Schwarzschild_EF__x_posn = hardwire_Schwarzschild_EF__x_posn; -gi.hardwire_Schwarzschild_EF__y_posn = hardwire_Schwarzschild_EF__y_posn; -gi.hardwire_Schwarzschild_EF__z_posn = hardwire_Schwarzschild_EF__z_posn; -gi.hardwire_Schwarzschild_EF__mass = hardwire_Schwarzschild_EF__mass; -gi.hardwire_Schwarzschild_EF__epsilon = hardwire_Schwarzschild_EF__epsilon; -gi.Delta_xyz = hardwire_Schwarzschild_EF__Delta_xyz; +gi.geometry_method = decode_geometry_method(geometry_method); +// parameters for geometry_method = "interpolate from Cactus grid" CCTK_VInfo(CCTK_THORNSTRING, " setting up geometry interpolator"); gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name); if (gi.operator_handle < 0) @@ -103,6 +102,14 @@ if (gi.param_table_handle < 0) "bad geometry-interpolator parameter(s) \"%s\"!", geometry_interpolator_pars); /*NOTREACHED*/ +// parameters for geometry_method = "Schwarzschild/EF" +gi.geometry__Schwarzschild_EF__mass = geometry__Schwarzschild_EF__mass; +gi.geometry__Schwarzschild_EF__x_posn = geometry__Schwarzschild_EF__x_posn; +gi.geometry__Schwarzschild_EF__y_posn = geometry__Schwarzschild_EF__y_posn; +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; + // // set up the interpatch interpolator @@ -227,8 +234,14 @@ Jacobian_info.perturbation_amplitude = Jacobian_perturbation_amplitude; struct solver_info solver_info; solver_info.max_Newton_iterations = max_Newton_iterations; +solver_info.output_h_and_H_at_each_Newton_iteration + = (output_h_and_H_at_each_Newton_iteration != 0); +solver_info.h_file_name = h_file_name; +solver_info.H_of_h_file_name = H_of_h_file_name; solver_info.H_norm_for_convergence = H_norm_for_convergence; solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence; +solver_info.final_H_update_if_exit_x_H_small + = (final_H_update_if_exit_x_H_small != 0); // // find the apparent horizon @@ -329,6 +342,26 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, } //****************************************************************************** + +// +// This function decodes the geometry_method parameter (string) into +// an internal enum for future use. +// +namespace { +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, + "unknown geometry_method_string=\"%s\"!", + geometry_method_string); /*NOTREACHED*/ +} + } + +//****************************************************************************** //****************************************************************************** //****************************************************************************** diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc index e47146e..5a7ebc1 100644 --- a/src/gr/gr_gfas.minc +++ b/src/gr/gr_gfas.minc @@ -21,12 +21,17 @@ partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, -H, H__fnd, - HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, +H, H__fnd, + +partial_d_HA, partial_d_HA__fnd, +partial_d_HB, partial_d_HB__fnd, +partial_d_HC, partial_d_HC__fnd, +partial_d_HD, partial_d_HD__fnd, +partial_d_H, partial_d_H__fnd, partial_HA_wrt_partial_d_h, partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h, partial_HB_wrt_partial_d_h__fnd, diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 0b9e0e0..f0b2af1 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -219,7 +219,7 @@ void horizon_Jacobian_SD(patch_system& ps, CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian (symbolic differentiation)"); ps.compute_synchronize_Jacobian(); -Jac.zero(); +Jac.zero_matrix(); // diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 02310ce..41eaf23 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -123,10 +123,22 @@ ps.synchronize(); // set up xyz positions of grid points setup_xyz_posns(ps, msg_flag); -// compute g_ij, K_ij, and partial_k g_ij -if (gi.hardwire_Schwarzschild_EF) - then Schwarzschild_EF_geometry(ps, gi, msg_flag); - else interpolate_geometry(ps, cgi, gi, msg_flag); +// compute the "geometry" g_ij, K_ij, and partial_k g_ij +switch (gi.geometry_method) + { +case geometry__interpolate_from_Cactus_grid: + interpolate_geometry(ps, cgi, gi, msg_flag); + break; +case geometry__Schwarzschild_EF: + Schwarzschild_EF_geometry(ps, gi, msg_flag); + break; +default: + error_exit(PANIC_EXIT, +"***** horizon_function(): unknown gi.geometry_method=(int)%d!\n" +" (this should never happen!)\n" +, + int(gi.geometry_method)); /*NOTREACHED*/ + } // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple index dc079cb..e833338 100644 --- a/src/gr/setup_gr_gfas.maple +++ b/src/gr/setup_gr_gfas.maple @@ -54,6 +54,13 @@ make_gfa('HD', {inert,fnd}, [], none); # LHS of apparent horizon equation make_gfa('H', {inert,fnd}, [], none); +# 1st derivatives of H[ABCD] and of H +make_gfa('partial_d_HA', {inert,fnd}, [1..N], none); +make_gfa('partial_d_HB', {inert,fnd}, [1..N], none); +make_gfa('partial_d_HC', {inert,fnd}, [1..N], none); +make_gfa('partial_d_HD', {inert,fnd}, [1..N], none); +make_gfa('partial_d_H', {inert,fnd}, [1..N], none); + # Jacobian coefficients for H[ABCD] # these are the partial derivatives of H[ABCD] # ... wrt Diff(h,x_rs[x]) diff --git a/src/include/config.hh b/src/include/config.hh index 6cf3693..2729a62 100644 --- a/src/include/config.hh +++ b/src/include/config.hh @@ -15,7 +15,14 @@ typedef CCTK_INT integer; // (CCTK_REAL_PRECISION_{4,8,16} are helpful, but not quite enough) #undef FP_IS_FLOAT #define FP_IS_DOUBLE -#define FP_SCANF_FORMAT "%lf" + +#if defined(FP_IS_FLOAT) + #define FP_SCANF_FORMAT "%f" +#elif defined(FP_IS_DOUBLE) + #define FP_SCANF_FORMAT "%lf" +#else + #error "don't know fp datatype!" +#endif // // The angular finite differencing in our multipatch system can be diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index 2fb1f17..a824e4a 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -1356,7 +1356,7 @@ if (output_fp == NULL) const fp global_y = global_y_of_local_y(local_y); const fp global_z = global_z_of_local_z(local_z); fprintf(output_fp, - "\t%#g\t%#g\t%#g", + "\t%#.10g\t%#.10g\t%#.10g", global_x, global_y, global_z); } fprintf(output_fp, "\n"); |