diff options
Diffstat (limited to 'src/driver/find_horizons.cc')
| -rw-r--r-- | src/driver/find_horizons.cc | 441 |
1 files changed, 260 insertions, 181 deletions
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 |