diff options
| -rw-r--r-- | doc/documentation.tex | 13 | ||||
| -rw-r--r-- | par/misner-run.par | 10 | ||||
| -rw-r--r-- | param.ccl | 24 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/Kerr.par | 5 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/try-5.par | 5 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/try-7.5-debug.par | 4 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/try-7.5-horizon.par | 55 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/try-7.5.par | 4 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/misc/try-9.par | 4 | ||||
| -rw-r--r-- | src/driver/driver.hh | 18 | ||||
| -rw-r--r-- | src/driver/find_horizons.cc | 441 | ||||
| -rw-r--r-- | src/driver/setup.cc | 21 |
12 files changed, 371 insertions, 233 deletions
diff --git a/doc/documentation.tex b/doc/documentation.tex index 68fcbe2..5236c65 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -787,12 +787,13 @@ $\bar{\bf x}^i \equiv \int\! {\bf x}^i \, dS \big/ \int\! dS$ is computed about as accurately, but is quite coordinate-dependent and gives only an approximate estimate of the black hole position. -The apparent horizon mass is currently defined (very crudely) as -$m \equiv \sqrt{A/16\pi}$. Even for Kerr spacetime this formula -deviates significantly from the true black hole mass if the spin -is large. (For example, for spin~$0.6$, $0.8$, and $0.999$ this -gives $m/m_\text{true} = 0.949$, $0.894$, and $0.723$.) A better -technique would be to use the ``isolated horizons'' formalism of +At present the only mass estimate available for an apparent horizon +is the irreducible mass $m_\text{irreducible} \equiv \sqrt{A/16\pi}$. +Note that even for Kerr spacetime, this formula deviates significantly +from the ADM black hole mass if the spin is large. (For example, +for spin~$0.6$, $0.8$, and $0.999$, gives +$m_\text{irreducible}/m_\text{ADM} = 0.949$, $0.894$, and $0.723$.) +It would be better to (also) use the ``isolated horizons'' formalism of \cite{AEIDevelopment_AHFinderDirect_Dreyer-etal-2002-isolated-horizons}; at some point this thorn may be enhanced to do this. diff --git a/par/misner-run.par b/par/misner-run.par index dbaacfb..77a3910 100644 --- a/par/misner-run.par +++ b/par/misner-run.par @@ -191,10 +191,10 @@ AHFinderDirect::origin_y[1] = 0.0 AHFinderDirect::origin_z[1] = 1.0 AHFinderDirect::delta_drho_dsigma = 7.5 -AHFinderDirect::initial_guess_method = "sphere" -AHFinderDirect::initial_guess__sphere__x_center[1] = 0.0 -AHFinderDirect::initial_guess__sphere__y_center[1] = 0.0 -AHFinderDirect::initial_guess__sphere__z_center[1] = 1.0 -AHFinderDirect::initial_guess__sphere__radius[1] = 0.3 +AHFinderDirect::initial_guess_method = "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.3 ######################################## @@ -14,10 +14,11 @@ USES KEYWORD metric_type # all remaining parameters are private to this thorn private: + ################################################################################ # -# overall parameters for the apparent horizon finding algorithm itself +# overall parameters for the apparent horizon finding algorithm # boolean find_AHs_at_postinitial \ @@ -30,18 +31,15 @@ boolean find_AHs_at_poststep \ { } "true" -keyword method "top-level method used to find the apparent horizon" +keyword method "what should this thorn do for each apparent horizon?" { # these options are mostly for testing/debugging -"horizon function" :: "evaluate the LHS function H(h)" -"Jacobian test" :: \ - "compute/print the J[H(h)] Jacobian matrix by all possible methods" -"Jacobian test (NP only)" :: \ - "compute/print the J[H(h)] Jacobian matrix by numerical perturbation only" +"horizon function" :: "evaluate the LHS function H(h)" +"test Jacobian" :: "compute/print the J[H(h)] Jacobian matrix" # this is for normal apparent horizon finding -"Newton solve" :: "find the horizon via Newton's method" -} "Newton solve" +"find horizon" :: "find the apparent horizon" +} "find horizon" # # We support searching for up to N_horizons distinct apparent horizons @@ -144,6 +142,14 @@ real Jacobian_perturbation_amplitude \ (0.0:* :: "any real number > 0" } 1.0e-6 +# if AHFinderDirect::method = "test Jacobian", should we test all +# known methods for computing the Jacobian, or just the numerical perturbation +# method (the latter may be useful of some other methods are broken) +boolean test_all_Jacobian_methods \ + "should we test all Jacobian computation methods, or just NP?" +{ +} "true" + ################################################################################ # diff --git a/run/test-ahfinderdirect/misc/Kerr.par b/run/test-ahfinderdirect/misc/Kerr.par index bcee119..a727c02 100644 --- a/run/test-ahfinderdirect/misc/Kerr.par +++ b/run/test-ahfinderdirect/misc/Kerr.par @@ -34,11 +34,6 @@ AHFinderDirect::N_ghost_points = 2 AHFinderDirect::N_overlap_points = 1 AHFinderDirect::delta_drho_dsigma = 5.0 -AHFinderDirect::geometry_interpolator_name = "generalized polynomial interpolation" -AHFinderDirect::geometry_interpolator_pars = "order=3" -AHFinderDirect::interpatch_interpolator_name = "generalized polynomial interpolation" -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 diff --git a/run/test-ahfinderdirect/misc/try-5.par b/run/test-ahfinderdirect/misc/try-5.par index fb9f964..20b070c 100644 --- a/run/test-ahfinderdirect/misc/try-5.par +++ b/run/test-ahfinderdirect/misc/try-5.par @@ -31,8 +31,7 @@ Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "false" -AHFinderDirect::verbose_level = "algorithm details" -AHFinderDirect::print_timing_stats = "true" +AHFinderDirect::verbose_level = "algorithm highlights" AHFinderDirect::final_H_update_if_exit_x_H_small = "true" AHFinderDirect::h_base_file_name = "try-5.h" @@ -44,8 +43,6 @@ AHFinderDirect::origin_y[1] = 0.7 AHFinderDirect::origin_z[1] = 0.0 AHFinderDirect::delta_drho_dsigma = 5.0 -AHFinderDirect::geometry_interpolator_pars = "order=3 out_of_range_tolerance={-1.0 -1.0 -1.0 -1.0 -1.0 -1.0}" - AHFinderDirect::initial_guess_method = "coordinate sphere" AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2 AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.3 diff --git a/run/test-ahfinderdirect/misc/try-7.5-debug.par b/run/test-ahfinderdirect/misc/try-7.5-debug.par index 8b4d8f2..e337cbc 100644 --- a/run/test-ahfinderdirect/misc/try-7.5-debug.par +++ b/run/test-ahfinderdirect/misc/try-7.5-debug.par @@ -31,7 +31,7 @@ Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "false" -AHFinderDirect::verbose_level = "algorithm details" +AHFinderDirect::verbose_level = "algorithm highlights" AHFinderDirect::print_timing_stats = "true" AHFinderDirect::final_H_update_if_exit_x_H_small = "true" @@ -46,8 +46,6 @@ AHFinderDirect::origin_y[1] = 0.7 AHFinderDirect::origin_z[1] = 0.0 AHFinderDirect::delta_drho_dsigma = 7.5 -AHFinderDirect::geometry_interpolator_pars = "order=3 out_of_range_tolerance={-1.0 -1.0 -1.0 -1.0 -1.0 -1.0}" - AHFinderDirect::initial_guess_method = "coordinate sphere" AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2 AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.3 diff --git a/run/test-ahfinderdirect/misc/try-7.5-horizon.par b/run/test-ahfinderdirect/misc/try-7.5-horizon.par new file mode 100644 index 0000000..2de1ce3 --- /dev/null +++ b/run/test-ahfinderdirect/misc/try-7.5-horizon.par @@ -0,0 +1,55 @@ +# parameter file for patch system test +ActiveThorns = "CartGrid3D LocalInterp PUGH ADMBase ADMCoupling StaticConformal CoordGauge Exact AHFinderDirect" + +# flesh +cactus::cctk_itlast = 0 + +# PUGH +driver::ghost_size = 2 +driver::global_nx = 31 +driver::global_ny = 31 +driver::global_nz = 17 + +# CartGrid3D +grid::domain = "bitant" +grid::avoid_origin = "false" +grid::type = "byspacing" +grid::dxyz = 0.2 + +# ADMBase +ADMBase::initial_lapse = "exact" +ADMBase::initial_shift = "exact" +ADMBase::initial_data = "exact" +ADMBase::lapse_evolution_method = "static" +ADMBase::shift_evolution_method = "static" +ADMBase::metric_type = "physical" + +# Exact +Exact::exact_model = "Kerr/Kerr-Schild" +Exact::Kerr_KerrSchild__mass = 1.0 +Exact::Kerr_KerrSchild__spin = 0.6 + +AHFinderDirect::method = "horizon function" +AHFinderDirect::find_AHs_at_postinitial = "true" +AHFinderDirect::find_AHs_at_poststep = "false" +AHFinderDirect::verbose_level = "algorithm highlights" +AHFinderDirect::print_timing_stats = "true" + +AHFinderDirect::output_initial_guess = "true" +AHFinderDirect::how_often_to_output_H_of_h = 1 +AHFinderDirect::final_H_update_if_exit_x_H_small = "true" +AHFinderDirect::h_base_file_name = "try-7.5-horizon.h" +AHFinderDirect::H_of_h_base_file_name = "try-7.5-horizon.H" + +AHFinderDirect::N_horizons = 1 +AHFinderDirect::origin_x[1] = 0.5 +AHFinderDirect::origin_y[1] = 0.7 +AHFinderDirect::origin_z[1] = 0.0 +AHFinderDirect::delta_drho_dsigma = 7.5 + +AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" +AHFinderDirect::initial_guess__Kerr_KerrSchild__x_posn[1] = 0.0 +AHFinderDirect::initial_guess__Kerr_KerrSchild__y_posn[1] = 0.0 +AHFinderDirect::initial_guess__Kerr_KerrSchild__z_posn[1] = 0.0 +AHFinderDirect::initial_guess__Kerr_KerrSchild__mass[1] = 1.0 +AHFinderDirect::initial_guess__Kerr_KerrSchild__spin[1] = 0.6 diff --git a/run/test-ahfinderdirect/misc/try-7.5.par b/run/test-ahfinderdirect/misc/try-7.5.par index 5e80c22..eadb5b7 100644 --- a/run/test-ahfinderdirect/misc/try-7.5.par +++ b/run/test-ahfinderdirect/misc/try-7.5.par @@ -31,7 +31,7 @@ Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "false" -AHFinderDirect::verbose_level = "algorithm details" +AHFinderDirect::verbose_level = "algorithm highlights" AHFinderDirect::print_timing_stats = "true" AHFinderDirect::final_H_update_if_exit_x_H_small = "true" @@ -44,8 +44,6 @@ AHFinderDirect::origin_y[1] = 0.7 AHFinderDirect::origin_z[1] = 0.0 AHFinderDirect::delta_drho_dsigma = 7.5 -AHFinderDirect::geometry_interpolator_pars = "order=3 out_of_range_tolerance={-1.0 -1.0 -1.0 -1.0 -1.0 -1.0}" - AHFinderDirect::initial_guess_method = "coordinate sphere" AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2 AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.3 diff --git a/run/test-ahfinderdirect/misc/try-9.par b/run/test-ahfinderdirect/misc/try-9.par index f8d38f9..4df4812 100644 --- a/run/test-ahfinderdirect/misc/try-9.par +++ b/run/test-ahfinderdirect/misc/try-9.par @@ -31,7 +31,7 @@ Exact::Kerr_KerrSchild__spin = 0.6 AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "false" -AHFinderDirect::verbose_level = "algorithm details" +AHFinderDirect::verbose_level = "algorithm highlights" AHFinderDirect::print_timing_stats = "true" AHFinderDirect::final_H_update_if_exit_x_H_small = "true" @@ -44,8 +44,6 @@ AHFinderDirect::origin_y[1] = 0.7 AHFinderDirect::origin_z[1] = 0.0 AHFinderDirect::delta_drho_dsigma = 9.0 -AHFinderDirect::geometry_interpolator_pars = "order=3 out_of_range_tolerance={-1.0 -1.0 -1.0 -1.0 -1.0 -1.0}" - AHFinderDirect::initial_guess_method = "coordinate sphere" AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2 AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.3 diff --git a/src/driver/driver.hh b/src/driver/driver.hh index 547ec26..867577e 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -10,9 +10,8 @@ enum method { method__horizon_function, - method__Jacobian_test, - method__Jacobian_test_NP_only, - method__Newton_solve // no comma + method__test_Jacobian, + method__find_horizon // no comma }; // @@ -91,7 +90,6 @@ struct initial_guess_info // struct solver_info { - bool output_initial_guess; bool debugging_output_at_each_Newton_iteration; int max_Newton_iterations__initial, max_Newton_iterations__subsequent; @@ -102,11 +100,20 @@ struct solver_info }; // +// This struct holds info for computing black hole diagnostics. +// +struct diagnostics_info + { + enum patch::integration_method surface_integral_method; + }; + +// // This struct holds info for I/O // struct IO_info { enum file_format file_format; + bool output_initial_guess; int how_often_to_output_h, how_often_to_output_H; bool output_ghost_zones_for_h; @@ -165,7 +172,6 @@ struct state enum method method; struct verbose_info verbose_info; int timer_handle; - enum patch::integration_method surface_integral_method; struct IO_info IO_info; struct Jacobian_info Jac_info; @@ -173,6 +179,8 @@ struct state struct cactus_grid_info cgi; struct geometry_info gi; + struct diagnostics_info diagnostics_info; + int N_horizons; // this vector is of size N_horizons+1, diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index 4781372..9506634 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -9,6 +9,9 @@ /// Cactus_gridfn_data_ptr - get a single data pointer from a variable index /// /// find_horizon - find a horizon +/// do_horizon_function +/// do_test_Jacobian +/// do_find_horizon /// /// BH_diagnostics - compute BH diagnostics for a horizon /// surface_integral - compute surface integral of a gridfn over the horizon @@ -67,16 +70,30 @@ const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex, const char gridfn_name[], bool check_for_NULL = true); -bool find_horizon(enum method method, - const struct verbose_info& verbose_info, int timer_handle, - struct IO_info& IO_info, - struct Jacobian_info& Jac_info, - const struct solver_info& solver_info, bool initial_find_flag, - struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, - int hn, int N_horizons); - -void BH_diagnostics(enum patch::integration_method surface_integral_method, +bool do_horizon_function + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, bool output_h, bool output_H, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, + int hn, int N_horizons); +bool do_test_Jacobian + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, + bool test_all_Jacobian_methods, + struct Jacobian_info& Jac_info, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, Jacobian* Jac_ptr, + int hn, int N_horizons); +bool do_find_horizon + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, bool output_h, bool output_H, + const struct solver_info& solver_info, bool initial_find_flag, + const struct Jacobian_info& Jac_info, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, Jacobian* Jac_ptr, + int hn, int N_horizons); + +void BH_diagnostics(const struct diagnostics_info& diagnostics_info, const struct verbose_info& verbose_info, struct AH_info& AH_info); fp surface_integral(const patch_system& ps, int src_gfn, @@ -102,7 +119,14 @@ const struct solver_info& solver_info = state.solver_info; if (state.timer_handle >= 0) then CCTK_TimerResetI(state.timer_handle); +// get the Cactus time step and decide if we want to output [hH] now IO_info.time_iteration = cctk_iteration; +const bool output_h + = (IO_info.how_often_to_output_h > 0) + && ((IO_info.time_iteration % IO_info.how_often_to_output_h) == 0); +const bool output_H + = (IO_info.how_often_to_output_H > 0) + && ((IO_info.time_iteration % IO_info.how_often_to_output_H) == 0); // what are the semantics of the Cactus gxx variables? if (CCTK_Equals(metric_type, "physical")) @@ -141,40 +165,70 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi); AH_info.initial_guess_info, IO_info, hn, verbose_info); - if (solver_info.output_initial_guess) + if (IO_info.output_initial_guess) then output_gridfn(ps, gfns::gfn__h, IO_info, IO_info.h_base_file_name, hn, verbose_info .print_algorithm_highlights); } - AH_info.AH_found - = find_horizon(state.method, - verbose_info, state.timer_handle, - IO_info, state.Jac_info, - solver_info, initial_find_flag, - state.cgi, state.gi, - ps, AH_info.Jac_ptr, - hn, state.N_horizons); - - if (AH_info.AH_found) - then { - BH_diagnostics(state.surface_integral_method, - verbose_info, - AH_info); - - if (verbose_info.print_physics_details) + AH_info.AH_found = false; + switch (state.method) + { + case method__horizon_function: + do_horizon_function(verbose_info, state.timer_handle, + IO_info, output_h, output_H, + state.cgi, state.gi, + ps, + hn, state.N_horizons); + break; + case method__test_Jacobian: + do_test_Jacobian(verbose_info, state.timer_handle, + IO_info, + test_all_Jacobian_methods, + state.Jac_info, state.cgi, state.gi, + ps, AH_info.Jac_ptr, + hn, state.N_horizons); + break; + case method__find_horizon: + AH_info.AH_found + = do_find_horizon(verbose_info, state.timer_handle, + IO_info, output_h, output_H, + solver_info, initial_find_flag, + state.Jac_info, state.cgi, state.gi, + ps, AH_info.Jac_ptr, + hn, state.N_horizons); + if (AH_info.AH_found) then { - CCTK_VInfo(CCTK_THORNSTRING, - "AH found: A=%.10g at (%f,%f,%f)", - double(AH_info.area), - double(AH_info.centroid_x), - double(AH_info.centroid_y), - double(AH_info.centroid_z)); - CCTK_VInfo(CCTK_THORNSTRING, - "estimated mass %.10g", - double(AH_info.mass)); + BH_diagnostics(state.diagnostics_info, + verbose_info, + AH_info); + if (verbose_info.print_physics_details) + then { + CCTK_VInfo(CCTK_THORNSTRING, + "AH found: A=%.10g at (%f,%f,%f)", + double(AH_info.area), + double(AH_info.centroid_x), + double(AH_info.centroid_y), + double(AH_info.centroid_z)); + CCTK_VInfo(CCTK_THORNSTRING, + "m_irreducible=%.10g", + double(AH_info.mass)); + } } + else { + if (verbose_info.print_physics_details) + then CCTK_VInfo(CCTK_THORNSTRING, + "no apparent horizon found"); + } + break; + default: + CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "\n" + " find_horizons(): unknown method=(int)%d!\n" + " (this should never happen!)" + , + int(method)); /*NOTREACHED*/ } } @@ -263,8 +317,8 @@ return data_ptr; //****************************************************************************** // -// This function finds (or more accurately tries to find) a single -// apparent horizon. +// This function implements AHFinderDirect::method == "horizon function", +// that is, it evaluates the H(h) function for a single apparent horizon. // // Note that if we decide to output h, we output it *after* any H(h) // evaluation or horizon finding has been done, to ensure that all the @@ -274,166 +328,191 @@ return data_ptr; // timer_handle = a valid Cactus timer handle if we want to time the // apparent horizon process, or -ve to skip this // (we only time the computation, not the file I/O) -// initial_find_flag = true ==> This is the first attempt to find this -// horizon, or this is a subsequent attempt -// and the immediately previous attempt failed. -// false ==> This isn't the first attempt to find this -// horizon, and we found it successfully on -// the immediately previous attempt. -// Jac_ptr = may be NULL if no Jacobian is needed (depending on method) +// output_[hH] = flags to control whether or not we should output [hH] // hn = the horizon number (used only in naming output files) // N_horizons = the total number of horizon(s) being searched for number // (used only in formatting info messages) // // Results: -// This function returns true if the computation succeeds; false if it fails. -// If method specifies finding the (an) apparent horizon, then "success" -// means finding an h satisfying H(h) = 0 to within the error tolerances. -// Otherwise, "success" means successfully evaluating the horizon function -// and/or its Jacobian, as appropriate. +// This function returns true if the evaluation was successful, false +// otherwise. // namespace { -bool find_horizon(enum method method, - const struct verbose_info& verbose_info, int timer_handle, - struct IO_info& IO_info, - struct Jacobian_info& Jac_info, - const struct solver_info& solver_info, bool initial_find_flag, - struct cactus_grid_info& cgi, struct geometry_info& gi, - patch_system& ps, Jacobian* Jac_ptr, - int hn, int N_horizons) +bool do_horizon_function + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, bool output_h, bool output_H, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, + int hn, int N_horizons) { -const bool output_h - = (IO_info.how_often_to_output_h > 0) - && ((IO_info.time_iteration % IO_info.how_often_to_output_h) == 0); -const bool output_H - = (IO_info.how_often_to_output_H > 0) - && ((IO_info.time_iteration % IO_info.how_often_to_output_H) == 0); - -switch (method) - { -// just evaluate the horizon function -case method__horizon_function: - { - jtutil::norm<fp> H_norms; - - if (timer_handle >= 0) - then CCTK_TimerStartI(timer_handle); - const bool status - = horizon_function(ps, cgi, gi, false, true, &H_norms); - if (timer_handle >= 0) - then CCTK_TimerStopI(timer_handle); - if (!status) - then return false; // *** ERROR RETURN *** - - if (H_norms.is_nonempty()) // might be empty if H(h) eval failed - then CCTK_VInfo(CCTK_THORNSTRING, - " H(h) rms-norm %.2e, infinity-norm %.2e", - H_norms.rms_norm(), H_norms.infinity_norm()); - - if (output_h) - then output_gridfn(ps, gfns::gfn__h, - IO_info, IO_info.h_base_file_name, - hn, verbose_info.print_algorithm_details); - if (output_H) - then output_gridfn(ps, gfns::gfn__H, - IO_info, IO_info.H_base_file_name, - hn, true); - return true; // *** NORMAL RETURN *** - } - -// just compute/print the NP Jacobian -case method__Jacobian_test_NP_only: - { - Jacobian& Jac_NP = *Jac_ptr; - if (! horizon_function(ps, cgi, gi, true)) - then return false; // *** ERROR RETURN *** - Jac_info.Jacobian_method = Jacobian_method__numerical_perturb; - if (! horizon_Jacobian(ps, Jac_NP, - cgi, gi, Jac_info, - true)) - then return false; // *** ERROR RETURN *** - - print_Jacobians(ps, - & Jac_NP, NULL, - IO_info, IO_info.Jacobian_base_file_name, - hn, true); - return true; // *** NORMAL RETURN *** +jtutil::norm<fp> H_norms; + +if (timer_handle >= 0) + then CCTK_TimerStartI(timer_handle); +const bool status = horizon_function(ps, cgi, gi, false, true, &H_norms); +if (timer_handle >= 0) + then CCTK_TimerStopI(timer_handle); +if (!status) + then return false; // *** ERROR RETURN *** + +if (H_norms.is_nonempty()) // might be empty if H(h) eval failed + then CCTK_VInfo(CCTK_THORNSTRING, + " H(h) rms-norm %.2e, infinity-norm %.2e", + H_norms.rms_norm(), H_norms.infinity_norm()); + +if (output_h) + then output_gridfn(ps, gfns::gfn__h, + IO_info, IO_info.h_base_file_name, + hn, verbose_info.print_algorithm_details); +if (output_H) + then output_gridfn(ps, gfns::gfn__H, + IO_info, IO_info.H_base_file_name, + hn, verbose_info.print_algorithm_details); + +return true; // *** NORMAL RETURN *** +} } -// compute/print the Jacobian by all possible methods -case method__Jacobian_test: - { - Jacobian& Jac_NP = *Jac_ptr; - if (! horizon_function(ps, cgi, gi, true)) - then return false; // *** ERROR RETURN *** - Jac_info.Jacobian_method = Jacobian_method__numerical_perturb; - if (! horizon_Jacobian(ps, Jac_NP, - cgi, gi, Jac_info, - true)) - then return false; // *** ERROR RETURN *** +//****************************************************************************** +// +// This function implements AHFinderDirect::method == "test Jacobian", +// that is, it computes and prints the Jacobian matrix J[H(h)] for a +// single apparent horizon. The computation may optionally be done +// in several different ways, in which case all the resulting Jacobian +// matrices are printed, as are their differences. Alternatively, only +// the numerical perturbation computation may be done/printed. +// +// Arguments: +// timer_handle = a valid Cactus timer handle if we want to time the +// apparent horizon process, or -ve to skip this +// (we only time the computation, not the file I/O) +// test_all_Jacobian_methods = true ==> Test all known methods of computing +// the Jacobian matrix, and print all +// the resulting Jacobian matrices +// and their differences. +// false ==> Just test/print the numerical +// perturbation calculation. (This +// may be useful if one or more of +// the other methods is broken.) +// hn = the horizon number (used only in naming output files) +// N_horizons = the total number of horizon(s) being searched for number +// (used only in formatting info messages) +// +// Results: +// This function returns true if the Jacobian computation was successful, +// false otherwise. +// +namespace { +bool do_test_Jacobian + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, + bool test_all_Jacobian_methods, + struct Jacobian_info& Jac_info, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, Jacobian* Jac_ptr, + int hn, int N_horizons) +{ +// numerical perturbation +Jacobian* Jac_NP_ptr = Jac_ptr; +if (! horizon_function(ps, cgi, gi, true)) + then return false; // *** ERROR RETURN *** +Jac_info.Jacobian_method = Jacobian_method__numerical_perturb; +if (! horizon_Jacobian(ps, *Jac_NP_ptr, + cgi, gi, Jac_info, + true)) + then return false; // *** ERROR RETURN *** + +Jacobian* Jac_SD_FDdr_ptr = NULL; +if (test_all_Jacobian_methods) + then { // symbolic differentiation with finite diff d/dr - Jacobian& Jac_SD_FDdr - = new_Jacobian(ps, Jac_info.Jacobian_storage_method); + Jac_SD_FDdr_ptr = & new_Jacobian(ps, Jac_info.Jacobian_storage_method); if (! horizon_function(ps, cgi, gi, true)) then return false; // *** ERROR RETURN *** Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr; - if (! horizon_Jacobian(ps, Jac_SD_FDdr, + if (! horizon_Jacobian(ps, *Jac_SD_FDdr_ptr, cgi, gi, Jac_info, true)) then return false; // *** ERROR RETURN *** + } - print_Jacobians(ps, - & Jac_NP, & Jac_SD_FDdr, - IO_info, IO_info.Jacobian_base_file_name, - hn, true); - return true; // *** NORMAL RETURN *** - } +print_Jacobians(ps, + Jac_NP_ptr, Jac_SD_FDdr_ptr, + IO_info, IO_info.Jacobian_base_file_name, + hn, true); -// find the apparent horizon via the Newton solver -case method__Newton_solve: - { - Jacobian& Jac = *Jac_ptr; - - if (verbose_info.print_algorithm_highlights) - then CCTK_VInfo(CCTK_THORNSTRING, - "searching for horizon #%d/%d", - hn, N_horizons); - - if (timer_handle >= 0) - then CCTK_TimerStartI(timer_handle); - const bool status - = Newton_solve(ps, Jac, - cgi, gi, - Jac_info, - solver_info, initial_find_flag, - IO_info, hn, verbose_info); - if (timer_handle >= 0) - then CCTK_TimerStopI(timer_handle); - if (! status) - then return false; // *** ERROR RETURN *** - - if (output_h) - then output_gridfn(ps, gfns::gfn__h, - IO_info, IO_info.h_base_file_name, - hn, verbose_info.print_algorithm_details); - if (output_H) - then output_gridfn(ps, gfns::gfn__H, - IO_info, IO_info.H_base_file_name, - hn, verbose_info.print_algorithm_details); - return true; // *** NORMAL RETURN *** +return true; // *** NORMAL RETURN *** +} } -default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" find_horizons(): unknown method=(int)%d!\n" -" (this should never happen!)" -, - int(method)); /*NOTREACHED*/ - } +//****************************************************************************** -/*NOTREACHED*/ +// +// This function implements AHFinderDirect::method == "find horizon", +// that is, it finds the (an) apparent horizon. +// +// Arguments: +// timer_handle = a valid Cactus timer handle if we want to time the +// apparent horizon process, or -ve to skip this +// (we only time the computation, not the file I/O) +// initial_find_flag = true ==> This is the first time we've tried to find +// this horizon, or we've tried before but +// failed on our most recent previous attempt. +// Thus we should use "initial" parameters +// for the Newton iteration. +// flase ==> We've tried to find this horizon before, +// and we succeeded on our most recent +// previous attempt. Thus we should use +// "subsequent" parameters for the Newton +// iteration. +// hn = the horizon number (used only in naming output files) +// N_horizons = the total number of horizon(s) being searched for number +// (used only in formatting info messages) +// +// Results: +// This function returns true if the apparent horizon was found +// successfully, false otherwise. +// +namespace { +bool do_find_horizon + (const struct verbose_info& verbose_info, int timer_handle, + struct IO_info& IO_info, bool output_h, bool output_H, + const struct solver_info& solver_info, bool initial_find_flag, + const struct Jacobian_info& Jac_info, + struct cactus_grid_info& cgi, struct geometry_info& gi, + patch_system& ps, Jacobian* Jac_ptr, + int hn, int N_horizons) +{ +Jacobian& Jac = *Jac_ptr; + +if (verbose_info.print_algorithm_highlights) + then CCTK_VInfo(CCTK_THORNSTRING, + "searching for horizon #%d of %d", + hn, N_horizons); + +if (timer_handle >= 0) + then CCTK_TimerStartI(timer_handle); +const bool status = Newton_solve(ps, Jac, + cgi, gi, + Jac_info, + solver_info, initial_find_flag, + IO_info, hn, verbose_info); +if (timer_handle >= 0) + then CCTK_TimerStopI(timer_handle); +if (! status) + then return false; // *** ERROR RETURN *** + +if (output_h) + then output_gridfn(ps, gfns::gfn__h, + IO_info, IO_info.h_base_file_name, + hn, verbose_info.print_algorithm_details); +if (output_H) + then output_gridfn(ps, gfns::gfn__H, + IO_info, IO_info.H_base_file_name, + hn, verbose_info.print_algorithm_details); + +return true; // *** NORMAL RETURN *** } } @@ -443,7 +522,7 @@ default: // // Given that an apparent horizon has been found, this function computes -// various BH diagnostics. +// various black hole diagnostics. // // Inputs (gridfns) // h # ghosted @@ -451,7 +530,7 @@ default: // global_[xyz] # nominal // namespace { -void BH_diagnostics(enum patch::integration_method surface_integral_method, +void BH_diagnostics(const struct diagnostics_info& diagnostics_info, const struct verbose_info& verbose_info, struct AH_info& AH_info) { @@ -461,13 +540,13 @@ const patch_system& ps = * AH_info.ps_ptr; // compute raw surface integrals // fp integral_one = surface_integral(ps, gfns::gfn__one, - surface_integral_method); + diagnostics_info.surface_integral_method); fp integral_x = surface_integral(ps, gfns::gfn__global_x, - surface_integral_method); + diagnostics_info.surface_integral_method); fp integral_y = surface_integral(ps, gfns::gfn__global_y, - surface_integral_method); + diagnostics_info.surface_integral_method); fp integral_z = surface_integral(ps, gfns::gfn__global_z, - surface_integral_method); + diagnostics_info.surface_integral_method); // // adjust integrals to take into account patch system not covering full sphere diff --git a/src/driver/setup.cc b/src/driver/setup.cc index d2a59b1..34cc0cd 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -111,11 +111,10 @@ verbose_info.print_algorithm_details = (state.verbose_info.verbose_level >= verbose_level__algorithm_details); state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1; -state.surface_integral_method - = patch::decode_integration_method(surface_integral_method); struct IO_info& IO_info = state.IO_info; IO_info.file_format = decode_file_format(file_format); +IO_info.output_initial_guess = (output_initial_guess != 0); IO_info.how_often_to_output_h = how_often_to_output_h; IO_info.how_often_to_output_H = how_often_to_output_H_of_h; IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0); @@ -134,7 +133,6 @@ Jac_info.Jacobian_storage_method Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude; struct solver_info& solver_info = state.solver_info; -solver_info.output_initial_guess = (output_initial_guess != 0); solver_info.debugging_output_at_each_Newton_iteration = (debugging_output_at_each_Newton_iteration != 0); solver_info.max_Newton_iterations__initial @@ -213,6 +211,13 @@ gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0); // +// other misc setup +// +state.diagnostics_info.surface_integral_method + = patch::decode_integration_method(surface_integral_method); + + +// // create AH-specific info for each AH // @@ -354,12 +359,10 @@ enum method { if (STRING_EQUAL(method_string, "horizon function")) then return method__horizon_function; -else if (STRING_EQUAL(method_string, "Jacobian test")) - then return method__Jacobian_test; -else if (STRING_EQUAL(method_string, "Jacobian test (NP only)")) - then return method__Jacobian_test_NP_only; -else if (STRING_EQUAL(method_string, "Newton solve")) - then return method__Newton_solve; +else if (STRING_EQUAL(method_string, "test Jacobian")) + then return method__test_Jacobian; +else if (STRING_EQUAL(method_string, "find horizon")) + then return method__find_horizon; else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_method(): unknown method_string=\"%s\"!", method_string); /*NOTREACHED*/ |