diff options
| -rw-r--r-- | README | 4 | ||||
| -rw-r--r-- | param.ccl | 18 | ||||
| -rw-r--r-- | src/driver/Newton.cc | 27 | ||||
| -rw-r--r-- | src/driver/driver.hh | 11 | ||||
| -rw-r--r-- | src/driver/find_horizons.cc | 54 | ||||
| -rw-r--r-- | src/driver/initial_guess.cc | 6 | ||||
| -rw-r--r-- | src/driver/io.cc | 15 | ||||
| -rw-r--r-- | src/driver/setup.cc | 33 | ||||
| -rw-r--r-- | src/make.configuration.defn | 5 |
9 files changed, 108 insertions, 65 deletions
@@ -96,8 +96,8 @@ g2c: Note that all of these settings are of shell environment variables, using the syntax (eg) - $ export LAPACK_DIR=/usr/bin # bash - % setenv LAPACK_DIR /usr/bin # csh,tcsh + % export LAPACK_DIR=/usr/lib # bash + % setenv LAPACK_DIR /usr/lib # csh,tcsh Alas, setting these in your ~/.cactus/config file doesn't work. :( @@ -214,9 +214,11 @@ boolean output_initial_guess \ { } "false" -# 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)?" +# for debugging failure to converge, we can optionally output +# h, H, and delta_h at each Newton iteration +# (the file names are the usual ones with ".it%d" appended) +boolean debugging_output_at_each_Newton_iteration \ + "should we output {h, H, delta_h} at each Newton iteration?" { } "false" @@ -264,14 +266,20 @@ string HDF5_file_name_extension "extension for HDF5 data files" string h_base_file_name \ "base file name for horizon shape h input/output file(s)" { -.* :: "any string" +.+ :: "any nonempty string" } "h.dat" string H_of_h_base_file_name "base file name for H(h) output file(s)" { -.* :: "any string" +.+ :: "any nonempty string" } "H.dat" +string Delta_h_base_file_name \ + "base file name for horizon-shape-update Delta_h output file(s)" +{ +.+ :: "any nonempty string" +} "Delta_h.dat" + string Jacobian_base_file_name "base file name for Jacobian output file(s)" { .+ :: "any valid file name" diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc index aa31876..f3efb56 100644 --- a/src/driver/Newton.cc +++ b/src/driver/Newton.cc @@ -44,6 +44,12 @@ using jtutil::error_exit; // This function solves H(h) = 0 for h via Newton's method. // // Arguments: +// 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. // hn = the horizon number (used only in naming output files) // // Results: @@ -57,18 +63,23 @@ bool Newton_solve(patch_system& ps, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, + bool initial_find_flag, struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info) { +const int max_Newton_iterations + = initial_find_flag ? solver_info.max_Newton_iterations__initial + : solver_info.max_Newton_iterations__subsequent; for (int iteration = 1 ; - iteration <= solver_info.max_Newton_iterations ; + iteration <= max_Newton_iterations ; ++iteration) { if (verbose_info.print_algorithm_details) then CCTK_VInfo(CCTK_THORNSTRING, - "Newton iteration %d", iteration); + "Newton iteration %d of %d", + iteration, max_Newton_iterations); - if (IO_info.output_h_and_H_at_each_Newton_iteration) + if (solver_info.debugging_output_at_each_Newton_iteration) then output_gridfn(ps, gfns::gfn__h, IO_info, IO_info.h_base_file_name, hn, verbose_info.print_algorithm_details, @@ -87,7 +98,7 @@ bool Newton_solve(patch_system& ps, iteration, H_norms.rms_norm(), H_norms.infinity_norm()); - if (IO_info.output_h_and_H_at_each_Newton_iteration) + if (solver_info.debugging_output_at_each_Newton_iteration) then output_gridfn(ps, gfns::gfn__H, IO_info, IO_info.H_base_file_name, hn, verbose_info.print_algorithm_details, @@ -129,6 +140,12 @@ bool Newton_solve(patch_system& ps, "done solving linear system"); } + if (solver_info.debugging_output_at_each_Newton_iteration) + then output_gridfn(ps, gfns::gfn__Delta_h, + IO_info, IO_info.Delta_h_base_file_name, + hn, verbose_info.print_algorithm_details, + iteration); + // if the Newton step is too large, scale it down jtutil::norm<fp> h_norms; ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms); @@ -199,7 +216,7 @@ bool Newton_solve(patch_system& ps, CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "Newton_solve: no convergence in %d iterations!", - solver_info.max_Newton_iterations); + max_Newton_iterations); if (solver_info.final_H_update_if_exit_x_H_small) then { if (verbose_info.print_algorithm_details) diff --git a/src/driver/driver.hh b/src/driver/driver.hh index 84134da..d5f9ea7 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -91,7 +91,11 @@ struct initial_guess_info // struct solver_info { - int max_Newton_iterations; + bool output_initial_guess; + bool output_h, output_H; + bool debugging_output_at_each_Newton_iteration; + int max_Newton_iterations__initial, + max_Newton_iterations__subsequent; fp max_Delta_h_over_h; fp H_norm_for_convergence; fp Delta_h_norm_for_convergence; @@ -103,15 +107,13 @@ struct solver_info // struct IO_info { - bool output_initial_guess; - bool output_h_and_H_at_each_Newton_iteration; - bool output_h, output_H; enum file_format file_format; bool output_ghost_zones_for_h; const char* ASCII_file_name_extension; const char* HDF5_file_name_extension; const char* h_base_file_name; const char* H_base_file_name; + const char* Delta_h_base_file_name; const char* Jacobian_base_file_name; // this is used to choose file names @@ -212,6 +214,7 @@ bool Newton_solve(patch_system& ps, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, + bool initial_find_flag, struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info); diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index 7df1fb3..7d38b34 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -61,7 +61,7 @@ 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, - struct solver_info& solver_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); @@ -83,7 +83,10 @@ extern "C" { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + const struct verbose_info& verbose_info = state.verbose_info; + struct IO_info& IO_info = state.IO_info; +const struct solver_info& solver_info = state.solver_info; if (state.timer_handle >= 0) then CCTK_TimerResetI(state.timer_handle); @@ -105,24 +108,13 @@ if (state.cgi.Cactus_conformal_metric) then state.cgi.psi_data = Cactus_gridfn_data_ptr(cctkGH, "StaticConformal::psi"); -state.IO_info.time_iteration = cctk_iteration; +IO_info.time_iteration = cctk_iteration; for (int hn = 1 ; hn <= state.N_horizons ; ++hn) { struct AH_info& AH_info = * state.AH_info_ptrs[hn]; patch_system& ps = *AH_info.ps_ptr; // - // The first time we (try to) find a given horizon, our initial - // guess is likely to be rather inaccurate, so we may need a - // larger number of iterations. But if we've found this horizon - // before, then we have its previous position as an initial guess, - // so we shouldn't need as many iterations. - // - state.solver_info.max_Newton_iterations - = AH_info.AH_found ? max_Newton_iterations__subsequent - : max_Newton_iterations__initial; - - // // If this is our first attempt to find this horizon, or // if we've tried to find it before but we failed on our // immediately previous attempt, then we need to (re)set @@ -132,16 +124,25 @@ state.IO_info.time_iteration = cctk_iteration; // then we can just reuse the previous horizon position as // our initial guess for this time around. // - if (! AH_info.AH_found) - then setup_initial_guess(ps, + const bool initial_find_flag = ! AH_info.AH_found; + if (initial_find_flag) + then { + setup_initial_guess(ps, AH_info.initial_guess_info, - state.IO_info, + IO_info, hn, verbose_info); + if (solver_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, - state.IO_info, state.Jac_info, state.solver_info, + IO_info, state.Jac_info, + solver_info, initial_find_flag, state.cgi, state.gi, ps, AH_info.Jac_ptr, hn, state.N_horizons); @@ -184,6 +185,12 @@ if (state.timer_handle >= 0) // 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) // hn = the horizon number (used only in naming output files) // N_horizons = the total number of horizon(s) being searched for number @@ -201,7 +208,7 @@ 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, - struct solver_info& solver_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) @@ -226,7 +233,7 @@ case method__horizon_function: then CCTK_VInfo(CCTK_THORNSTRING, " H(h) rms-norm %.2e, infinity-norm %.2e", H_norms.rms_norm(), H_norms.infinity_norm()); - if (IO_info.output_H) + if (solver_info.output_H) then output_gridfn(ps, gfns::gfn__H, IO_info, IO_info.H_base_file_name, hn, true); @@ -297,18 +304,19 @@ case method__Newton_solve: const bool status = Newton_solve(ps, Jac, cgi, gi, - Jac_info, solver_info, IO_info, - hn, verbose_info); + 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 (IO_info.output_h) + if (solver_info.output_h) then output_gridfn(ps, gfns::gfn__h, IO_info, IO_info.h_base_file_name, hn, verbose_info.print_algorithm_details); - if (IO_info.output_H) + if (solver_info.output_H) then output_gridfn(ps, gfns::gfn__H, IO_info, IO_info.H_base_file_name, hn, verbose_info.print_algorithm_details); diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc index 73de3e6..0c11742 100644 --- a/src/driver/initial_guess.cc +++ b/src/driver/initial_guess.cc @@ -134,12 +134,6 @@ default: "unknown initial guess method=(int)%d!", int(igi.method)); /*NOTREACHED*/ } - -// write initial guess back to the data file? -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); } //****************************************************************************** diff --git a/src/driver/io.cc b/src/driver/io.cc index 9a9e01d..af7d5f9 100644 --- a/src/driver/io.cc +++ b/src/driver/io.cc @@ -139,15 +139,20 @@ case file_format__ASCII: if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " writing H to \"%s\"", file_name); - ps.print_gridfn(unknown_gfn, - file_name); + ps.print_gridfn(unknown_gfn, file_name); + break; + case gfns::gfn__Delta_h: + if (print_msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " writing Delta_h to \"%s\"", file_name); + ps.print_gridfn(unknown_gfn, file_name); break; default: if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " writing \"%s\"", file_name); - ps.print_gridfn(unknown_gfn, - file_name); + " writing gfn=%d to \"%s\"", + unknown_gfn, file_name); + ps.print_gridfn(unknown_gfn, file_name); break; } break; diff --git a/src/driver/setup.cc b/src/driver/setup.cc index daf3963..f3e859d 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -98,35 +98,32 @@ CCTK_VInfo(CCTK_THORNSTRING, // state.method = decode_method(method); -state.verbose_info.verbose_level = decode_verbose_level(verbose_level); -state.verbose_info.print_physics_highlights + +struct verbose_info& verbose_info = state.verbose_info; +verbose_info.verbose_level = decode_verbose_level(verbose_level); +verbose_info.print_physics_highlights = (state.verbose_info.verbose_level >= verbose_level__physics_highlights); -state.verbose_info.print_physics_details +verbose_info.print_physics_details = (state.verbose_info.verbose_level >= verbose_level__physics_details); -state.verbose_info.print_algorithm_highlights +verbose_info.print_algorithm_highlights = (state.verbose_info.verbose_level >= verbose_level__algorithm_highlights); -state.verbose_info.print_algorithm_details +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); -const struct verbose_info& verbose_info = state.verbose_info; - struct IO_info& IO_info = state.IO_info; -IO_info.output_initial_guess = (output_initial_guess != 0); -IO_info.output_h_and_H_at_each_Newton_iteration - = (output_h_and_H_at_each_Newton_iteration != 0); -IO_info.output_h = (output_h != 0); -IO_info.output_H = (output_H_of_h != 0); IO_info.file_format = decode_file_format(file_format); IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0); IO_info.ASCII_file_name_extension = ASCII_file_name_extension; IO_info.HDF5_file_name_extension = HDF5_file_name_extension; IO_info.h_base_file_name = h_base_file_name; IO_info.H_base_file_name = H_of_h_base_file_name; +IO_info.Delta_h_base_file_name = Delta_h_base_file_name; IO_info.Jacobian_base_file_name = Jacobian_base_file_name; IO_info.time_iteration = 0; @@ -137,9 +134,15 @@ Jac_info.Jacobian_storage_method Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude; struct solver_info& solver_info = state.solver_info; -solver_info.max_Newton_iterations = 0; // dummy value; actual value - // will be filled in by - // AHFinderDirect_find_horizons() +solver_info.output_initial_guess = (output_initial_guess != 0); +solver_info.output_h = (output_h != 0); +solver_info.output_H = (output_H_of_h != 0); +solver_info.debugging_output_at_each_Newton_iteration + = (debugging_output_at_each_Newton_iteration != 0); +solver_info.max_Newton_iterations__initial + = max_Newton_iterations__initial; +solver_info.max_Newton_iterations__subsequent + = max_Newton_iterations__subsequent; solver_info.max_Delta_h_over_h = max_Delta_h_over_h; solver_info.H_norm_for_convergence = H_norm_for_convergence; solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence; diff --git a/src/make.configuration.defn b/src/make.configuration.defn index 932a37a..f5ab6fb 100644 --- a/src/make.configuration.defn +++ b/src/make.configuration.defn @@ -27,6 +27,11 @@ MissingLAPACK_DIR: @echo '*** LAPACK_EXTRA_LIBS = g2c' @echo '*** LAPACK_EXTRA_LIBDIRS = `g77 --print-file-name=libg2c.a | xargs dirname`' @echo '***' + @echo '*** Note that all of these settings are shell environment variables, e.g.,' + @echo '*** % export LAPACK_DIR=/usr/lib # bash' + @echo '*** % setenv LAPACK_DIR /usr/lib # csh, tcsh' + @echo '*** Alas, setting these in your ~/.cactus/config file doesn'\''t work. :(' + @echo '***' exit 2 endif |