diff options
| author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-02-16 22:21:49 +0000 |
|---|---|---|
| committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2003-02-16 22:21:49 +0000 |
| commit | cc221a875a640142697df54a2bcff6638f7b1a88 (patch) | |
| tree | a1477d8a3dd8b96d0dea7a2af4546556166389aa | |
| parent | 04ed769ddfacb1d6f5faedb66d3ff8d512ede7d4 (diff) | |
* general tidying of multiprocessor code
* add current iteration number to what's sent to processor #0 for
printing at each iteration
* add a new verbosity level "algorithm debug" for what used to be
#ifdef DEBUG debugging code
* add new parameters & logic to set CCTK_Warn() levels for various
warnings (eg AH is outside the gid, metric isn't +ve definite, etc)
* restore the ability to output the horizon shape to a data file,
which I mistakenly deleted a few days ago :(
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@947 f88db872-0e4f-0410-b76b-b9085cfa78c5
35 files changed, 738 insertions, 274 deletions
diff --git a/par/Kerr-tiny.par b/par/Kerr-tiny.par index 3c91a90..697f486 100644 --- a/par/Kerr-tiny.par +++ b/par/Kerr-tiny.par @@ -11,10 +11,10 @@ cactus::cctk_itlast = 0 ActiveThorns = "PUGH" -driver::ghost_size = 2 +driver::ghost_size = 3 driver::global_nx = 31 driver::global_ny = 31 -driver::global_nz = 17 +driver::global_nz = 19 ActiveThorns = "CartGrid3D" grid::domain = "bitant" diff --git a/par/Kerr.par b/par/Kerr.par index c9e1077..59f24b7 100644 --- a/par/Kerr.par +++ b/par/Kerr.par @@ -8,10 +8,10 @@ cactus::cctk_itlast = 0 ActiveThorns = "PUGH" -driver::ghost_size = 2 +driver::ghost_size = 3 driver::global_nx = 31 driver::global_ny = 31 -driver::global_nz = 17 +driver::global_nz = 19 ActiveThorns = "CartGrid3D" grid::domain = "bitant" @@ -38,7 +38,6 @@ AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "false" AHFinderDirect::print_timing_stats = "true" -AHFinderDirect::final_Theta_update_if_Delta_h_converged = "true" AHFinderDirect::how_often_to_output_Theta = 1 AHFinderDirect::h_base_file_name = "Kerr.h" AHFinderDirect::Theta_base_file_name = "Kerr.Theta" @@ -7,8 +7,9 @@ up Kerr/Kerr-Schild initial data and finds the apparent horizon in it, with the coordinate origin and initial guess offset to make this a nontrivial test of the apparent horizon finder. -misner-init.par is another fairly simple parameter file, which sets up -Misner initial data and then finds the apparent horizons in it. +misner1.2-025-init.par sets up Misner initial data (mu=1.2, so there's +a common horizon as well as two separate horizons) and then finds the +apparent horizons in it. misner-run.par is a more complicated parameter file, which sets up the same Misner initial data, but now time-evolves it, finding the apparent diff --git a/par/misner-run.par b/par/misner-run.par index a5c0caf..1ddbf79 100644 --- a/par/misner-run.par +++ b/par/misner-run.par @@ -14,7 +14,7 @@ ######################################## -ActiveThorns = "BAM_Elliptic ADMBase ADMConstraints StaticConformal Ellbase CartGrid3D Time PUGH PUGHreduce PUGHslab IOAscii IOUtil IOBasic ADMCoupling ADMMacros ADMAnalysis SpaceMask CoordGauge ADM_BSSN Boundary PUGHInterp Maximal NaNChecker AHFinder LocalInterp AHFinderDirect" +ActiveThorns = "BAM_Elliptic ADMBase ADMConstraints StaticConformal Ellbase CartGrid3D Time PUGH PUGHReduce PUGHSlab IOAscii IOUtil IOBasic ADMCoupling ADMMacros ADMAnalysis SpaceMask CoordGauge ADM_BSSN Boundary PUGHInterp Maximal NaNChecker AHFinder LocalInterp AHFinderDirect" ######################################## @@ -181,19 +181,27 @@ AHFinderDirect::find_AHs_at_postinitial = "true" AHFinderDirect::find_AHs_at_poststep = "true" AHFinderDirect::verbose_level = "algorithm highlights" -AHFinderDirect::final_Theta_update_if_Delta_h_converged = "true" AHFinderDirect::h_base_file_name = "misner.h" AHFinderDirect::how_often_to_output_h = 10 # time steps -AHFinderDirect::N_horizons = 1 +AHFinderDirect::N_horizons = 2 + AHFinderDirect::origin_x[1] = 0.0 AHFinderDirect::origin_y[1] = 0.0 AHFinderDirect::origin_z[1] = 1.0 - AHFinderDirect::initial_guess_method[1] = "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 -######################################## +AHFinderDirect::origin_x[2] = 0.0 +AHFinderDirect::origin_y[2] = 0.0 +AHFinderDirect::origin_z[2] = 0.0 +AHFinderDirect::initial_guess_method[1] = "coordinate ellipsoid" +AHFinderDirect::initial_guess__coord_ellipsoid__x_center[2] = 0.0 +AHFinderDirect::initial_guess__coord_ellipsoid__y_center[2] = 0.0 +AHFinderDirect::initial_guess__coord_ellipsoid__z_center[2] = 0.0 +AHFinderDirect::initial_guess__coord_ellipsoid__x_radius[2] = 0.8 +AHFinderDirect::initial_guess__coord_ellipsoid__y_radius[2] = 0.8 +AHFinderDirect::initial_guess__coord_ellipsoid__z_radius[2] = 1.4 @@ -99,6 +99,10 @@ keyword verbose_level \ # lots of details tracing what the code is doing "algorithm details" :: \ "physics details + lots of messages about the AH-finding algorithm" + +# even more details tracing what the code is doing :) +"algorithm debug" :: \ + "physics details + lots and lots of messages about the AH-finding algorithm" } "algorithm highlights" # n.b. printing timing stats is independent of verbose_level @@ -260,8 +264,11 @@ keyword horizon_file_format \ "HDF5" :: "HDF5 surface format (alas not implemented yet)" } "ASCII (gnuplot)" +# n.b. this refers to the patch system (angular) interpatch ghost zones, +# *not* the Cactus interprocessor ghost zones boolean output_ghost_zones_for_h \ - "should we include the ghost zones in h data files?" + "should we include the patch system (angular) interpatch ghost zones \ + in h data files?" { } "false" @@ -470,7 +477,7 @@ keyword patch_system_type[5] "what type of patch system should we use?" int N_zones_per_right_angle[5] "sets angular resolution of patch systems" { 1:* :: "any integer >= 1; must be even for patch systems other than full-sphere" -} 12 +} 18 # # This parameter sets the width of the interpatch ghost zones in the @@ -814,6 +821,80 @@ real initial_guess__coord_ellipsoid__z_radius[5] "z radius of ellipsoid" ################################################################################ # +# ***** parameters controlling how we handle errors/warnings ***** +# + +# +# We currently "handle" the following errors/warnings by CCTK_VWarn(); +# these parameters give the warning levels (in the sense of the -W flag +# to Cactus, and the first argument to CCTK_VWarn()) for various conditions: +# + +# point outside (or too close to boundary of) Cactus grid, +# i.e. geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE +# ... warning level if error occurs on first Newton iteration, +# i.e. when evaluating expansion/Jacobian for initial guess +# ... this may mean that driver::ghost_size is too small, especially +# if any horizon crosses a symmetry boundary +# ==> visible message by default +int warn_level__point_outside__initial \ + "warning level for point outside (or too close to boundary of) Cactus grid \ + (error occurs on first Newton iteration)" +{ +-1:* :: "any valid Cactus warning level" +} 1 +# ... warning level if error occurs on a subsequent Newton iteration +# ... this probably "just" means there's no horizon in the grid +# ==> no visible message by default +int warn_level__point_outside__subsequent \ + "warning level for point outside (or too close to boundary of) Cactus grid \ + (error occurs on subsequent Newton iteration)" +{ +-1:* :: "any valid Cactus warning level" +} 2 + +# the Cactus configure process didn't find a finite() function, +# so we're skipping the "check that the geometry is finite" tests +# even though the user set check_that_geometry_is_finite = "true" +int warn_level__skipping_finite_check \ + "warning level if the user sets check_that_geometry_is_finite \ + but the Cactus configure process doesn't find a finite() function \ + so we have to skip the finite-geometry check" +{ +-1:* :: "any valid Cactus warning level" +} 3 + +# infinity and/or NaN in interpolated {g_ij, partial_k g_ij, K_ij} +int warn_level__nonfinite_geometry \ + "warning level if we find infinity and/or NaN in the interpolated geometry \ + values {g_ij, partial_k g_ij, K_ij}" +{ +-1:* :: "any valid Cactus warning level" +} 1 + +# interpolated g_{ij} isn't positive definite +# (usually this means we're too close to a singularity) +# ... warning level if error occurs on first Newton iteration, +# i.e. when evaluating expansion/Jacobian for initial guess +int warn_level__gij_not_positive_definite__initial \ + "warning level if the interpolated g_{ij} isn't positive definite \ + (usually this means we're too close to a singularity) \ + (error occurs on first Newton iteration)" +{ +-1:* :: "any valid Cactus warning level" +} 2 +# ... warning level if error occurs on a subsequent Newton iteration +int warn_level__gij_not_positive_definite__subsequent \ + "warning level if the interpolated g_{ij} isn't positive definite \ + (usually this means we're too close to a singularity) \ + (error occurs on subsequent Newton iteration)" +{ +-1:* :: "any valid Cactus warning level" +} 2 + +################################################################################ + +# # ***** parameters for the test driver "src/patch/test_patch_system.cc" ***** # # By default this test driver isn't compiled into the cactus executable, diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc index ca411f8..5ec1418 100644 --- a/src/driver/BH_diagnostics.cc +++ b/src/driver/BH_diagnostics.cc @@ -143,7 +143,7 @@ case patch_system::patch_system__plus_xyz_octant_rotating: integral_z = integral_one * ps.origin_z(); break; default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n" " (this should never happen!)" @@ -247,7 +247,7 @@ snprintf(file_name_buffer, N_file_name_buffer, hn, IO_info.BH_diagnostics_file_name_extension); FILE *fileptr = fopen(file_name_buffer, "w"); if (fileptr == NULL) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " setup_BH_diagnostics_output_file():\n" " can't open BH-diagnostics output file\n" diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc index 1af34a0..2cde6c5 100644 --- a/src/driver/Newton.cc +++ b/src/driver/Newton.cc @@ -41,8 +41,6 @@ using jtutil::error_exit; #include "horizon_sequence.hh" #include "driver.hh" -#define DEBUG - //****************************************************************************** // @@ -54,14 +52,13 @@ bool reduce_and_broadcast(cGH* GH, int N_procs, int N_active_procs, int my_proc, bool my_active_flag, bool we_need_more_iterations, - int hn, fp rms_norm, fp infinity_norm, - int hn_buffer[], - fp rms_norm_buffer[], - fp infinity_norm_buffer[]); -void print_Theta_norms(int N_active_procs, - const int hn_buffer[], - const fp rms_norm_buffer[], - const fp infinity_norm_buffer[]); + int hn, int iteration, fp rms_norm, fp infinity_norm, + int hn_buffer[], int iteration_buffer[], + fp rms_norm_buffer[], fp infinity_norm_buffer[]); +void print_Theta_norms + (int N_active_procs, + const int hn_buffer[], const int iteration_buffer[], + const fp rms_norm_buffer[], const fp infinity_norm_buffer[]); void Newton_step(patch_system& ps, fp max_allowable_Delta_h_over_h, const struct verbose_info& verbose_info); @@ -73,7 +70,8 @@ void Newton_step(patch_system& ps, // This function tries to finds each horizon assigned to this processor, // by solving Theta(h) = 0 via Newton's method. For each horizon found, // it computes the BH diagnostics, optionally prints them (via CCTK_VInfo()), -// and optionally writes them to a BH diagnostics output file. +// and optionally writes them to a BH diagnostics output file. It also +// optionally writes the horizon shape itself to a data file. // // This function must be called synchronously across all processors, // and it operates synchronously, returning only when every processor @@ -96,6 +94,7 @@ void Newton(cGH* GH, const struct solver_info& solver_info, const struct IO_info& IO_info, const struct BH_diagnostics_info& BH_diagnostics_info, + const struct error_info& error_info, const struct verbose_info& verbose_info) { const bool my_active_flag = hs.has_genuine_horizons(); @@ -103,16 +102,28 @@ const int N_horizons = hs.N_horizons(); // buffers for reduce_and_broadcast() int* hn_buffer = NULL; +int* iteration_buffer = NULL; fp * rms_norm_buffer = NULL; fp * infinity_norm_buffer = NULL; if (hn_buffer == NULL) then { // allocate on first call hn_buffer = new int[N_active_procs]; + iteration_buffer = new int[N_active_procs]; rms_norm_buffer = new fp [N_active_procs]; infinity_norm_buffer = new fp [N_active_procs]; } +// print out which horizons we're finding on this processor +if (hs.has_genuine_horizons()) + then CCTK_VInfo(CCTK_THORNSTRING, + "proc %d: searching for horizon%s %s", + my_proc, + (N_horizons > 1 ? "s" : ""), hs.sequence_string(",")); + else CCTK_VInfo(CCTK_THORNSTRING, + "proc %d: dummy horizon(s) only", + my_proc); + // // each pass through this loop finds a single horizon, // or does corresponding dummy-horizon calls @@ -158,11 +169,10 @@ if (hn_buffer == NULL) // for (int iteration = 1 ; ; ++iteration) { - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - "Newton_solve(): iteration %d", - iteration); - #endif + if (verbose_info.print_algorithm_debug) + then CCTK_VInfo(CCTK_THORNSTRING, + "beginning iteration %d (horizon_is_genuine=%d)", + iteration, int(horizon_is_genuine)); // @@ -173,36 +183,33 @@ if (hn_buffer == NULL) && solver_info.debugging_output_at_each_Newton_iteration) then output_gridfn(*ps_ptr, gfns::gfn__h, IO_info, IO_info.h_base_file_name, - hn, verbose_info.print_algorithm_details, + hn, verbose_info.print_algorithm_highlights, iteration); - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - " computing Theta: horizon_is_genuine=%d", - int(horizon_is_genuine)); - #endif jtutil::norm<fp> Theta_norms; const bool Theta_is_ok = expansion(ps_ptr, cgi, gi, + error_info, (iteration == 1), true, // yes, we want the Jacobian coeffs verbose_info.print_algorithm_details, &Theta_norms); - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - "Newton_solve(): Theta_is_ok=%d", - Theta_is_ok); - if (horizon_is_genuine && Theta_is_ok) - then CCTK_VInfo(CCTK_THORNSTRING, - " Theta rms-norm %.1e, infinity-norm %.1e", - double(Theta_norms.rms_norm()), - double(Theta_norms.infinity_norm())); - #endif + if (verbose_info.print_algorithm_debug) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " Newton_solve(): Theta_is_ok=%d", + Theta_is_ok); + if (horizon_is_genuine && Theta_is_ok) + then CCTK_VInfo(CCTK_THORNSTRING, + " Theta rms-norm %.1e, infinity-norm %.1e", + double(Theta_norms.rms_norm()), + double(Theta_norms.infinity_norm())); + } if (horizon_is_genuine && Theta_is_ok && solver_info.debugging_output_at_each_Newton_iteration) then output_gridfn(*ps_ptr, gfns::gfn__Theta, IO_info, IO_info.Theta_base_file_name, - hn, verbose_info.print_algorithm_details, + hn, verbose_info.print_algorithm_highlights, iteration); @@ -218,35 +225,13 @@ if (hn_buffer == NULL) .Theta_norm_for_convergence); if (horizon_is_genuine) then AH_info_ptr->found_flag = found_this_horizon; - if (found_this_horizon) - then { - compute_BH_diagnostics(*ps_ptr, BH_diagnostics_info, - verbose_info, - AH_info_ptr->BH_diagnostics); - // FIXME: see header comment - print_BH_diagnostics(AH_info_ptr->BH_diagnostics, - N_horizons, hn, - verbose_info); - if (IO_info.output_BH_diagnostics) - then { - if (AH_info_ptr->BH_diagnostics_fileptr == NULL) - then AH_info_ptr->BH_diagnostics_fileptr - = setup_BH_diagnostics_output_file(IO_info, - N_horizons, - hn); - output_BH_diagnostics_fn(AH_info_ptr->BH_diagnostics, - IO_info, - hn, AH_info_ptr - ->BH_diagnostics_fileptr); - } - } // does *this* horizon need more iterations? // i.e. has this horizon's Newton iteration not yet converged? const bool this_horizon_needs_more_iterations = horizon_is_genuine && Theta_is_ok && !found_this_horizon - && (iteration <= max_iterations); + && (iteration < max_iterations); // do we (this processor) need to do more iterations // on this or a following horizon? @@ -257,37 +242,73 @@ if (hn_buffer == NULL) if (verbose_info.print_algorithm_details) then { CCTK_VInfo(CCTK_THORNSTRING, - "flags: found_this_horizon=%d", + " flags: found_this_horizon=%d", int(found_this_horizon)); CCTK_VInfo(CCTK_THORNSTRING, - " this_horizon_needs_more_iterations=%d", + " this_horizon_needs_more_iterations=%d", int(this_horizon_needs_more_iterations)); CCTK_VInfo(CCTK_THORNSTRING, - " we_need_more_iterations=%d", + " we_need_more_iterations=%d", int(we_need_more_iterations)); } // does any processor need more iterations // on its current or a following horizon? const bool any_proc_needs_more_iterations - = reduce_and_broadcast - (GH, N_procs, N_active_procs, - my_proc, my_active_flag, - we_need_more_iterations, - hn, horizon_is_genuine ? Theta_norms.rms_norm() : 0.0, - horizon_is_genuine ? Theta_norms.infinity_norm() : 0.0, - hn_buffer, rms_norm_buffer, infinity_norm_buffer); + = reduce_and_broadcast(GH, N_procs, N_active_procs, + my_proc, my_active_flag, + we_need_more_iterations, + hn, iteration, + ((horizon_is_genuine && Theta_is_ok) + ? Theta_norms.rms_norm() : 0.0), + ((horizon_is_genuine && Theta_is_ok) + ? Theta_norms.infinity_norm() : 0.0), + hn_buffer, iteration_buffer, + rms_norm_buffer, infinity_norm_buffer); if (verbose_info.print_algorithm_details) then { CCTK_VInfo(CCTK_THORNSTRING, - " ==> any_proc_needs_more_iterations=%d", + " ==> any_proc_needs_more_iterations=%d", int(any_proc_needs_more_iterations)); } if (my_active_flag && verbose_info.print_algorithm_highlights) then print_Theta_norms(N_active_procs, - hn_buffer, rms_norm_buffer, - infinity_norm_buffer); + hn_buffer, iteration_buffer, + rms_norm_buffer, infinity_norm_buffer); + if (found_this_horizon) + then { + compute_BH_diagnostics(*ps_ptr, BH_diagnostics_info, + verbose_info, + AH_info_ptr->BH_diagnostics); + // FIXME: see header comment + print_BH_diagnostics(AH_info_ptr->BH_diagnostics, + N_horizons, hn, + verbose_info); + if (IO_info.output_BH_diagnostics) + then { + if (AH_info_ptr->BH_diagnostics_fileptr == NULL) + then AH_info_ptr->BH_diagnostics_fileptr + = setup_BH_diagnostics_output_file(IO_info, + N_horizons, + hn); + output_BH_diagnostics_fn(AH_info_ptr->BH_diagnostics, + IO_info, + hn, AH_info_ptr + ->BH_diagnostics_fileptr); + } + + if (IO_info.output_h) + then output_gridfn(*ps_ptr, gfns::gfn__h, + IO_info, IO_info.h_base_file_name, + hn, verbose_info + .print_algorithm_highlights); + if (IO_info.output_Theta) + then output_gridfn(*ps_ptr, gfns::gfn__Theta, + IO_info, IO_info.Theta_base_file_name, + hn, verbose_info + .print_algorithm_highlights); + } // @@ -303,10 +324,9 @@ if (hn_buffer == NULL) } if ((N_procs == 1) && !this_horizon_needs_more_iterations) then { - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - " single processor: skipping to next horizon"); - #endif + if (verbose_info.print_algorithm_debug) + then CCTK_VInfo(CCTK_THORNSTRING, + "==> [single processor] Skipping to next horizon"); break; // *** LOOP EXIT *** } @@ -315,17 +335,17 @@ if (hn_buffer == NULL) // compute the Jacobian matrix // *** this is a synchronous operation across all processors *** // - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - " computing Jacobian: genuine/dummy flag %d", - int(this_horizon_needs_more_iterations)); - #endif + if (verbose_info.print_algorithm_debug) + then CCTK_VInfo(CCTK_THORNSTRING, + " computing Jacobian: genuine/dummy flag %d", + int(this_horizon_needs_more_iterations)); const bool Jacobian_is_ok = expansion_Jacobian(this_horizon_needs_more_iterations ? ps_ptr : NULL, this_horizon_needs_more_iterations ? Jac_ptr : NULL, cgi, gi, Jacobian_info, + error_info, (iteration == 1), verbose_info.print_algorithm_details); @@ -335,10 +355,9 @@ if (hn_buffer == NULL) // if (! (this_horizon_needs_more_iterations && Jacobian_is_ok)) then { - #ifdef DEBUG - CCTK_VInfo(CCTK_THORNSTRING, - " skipping to next horizon"); - #endif + if (verbose_info.print_algorithm_debug) + then CCTK_VInfo(CCTK_THORNSTRING, + " skipping to next horizon"); break; // *** LOOP EXIT *** } @@ -401,7 +420,8 @@ assert( false ); // This function does an inclusive-or--reduction of the "we need more // iterations (either on this or on another genuine horizon" flags across // all active processors. It also broadcasts which horizon we're working -// on, and the current error norms, to processor #0. +// on, the current iteration number, and the current error norms, to +// processor #0. // // The reduction and broadcast are combined in this function to allow // doing them all with a single reduction operation (and thus only a @@ -415,10 +435,11 @@ assert( false ); // we_need_more_iterations = "we need more iterations (either on this // or on another genuine horizon" flag to be // inclusive-or--reduced -// hn,{rms,infinity}_norm = horizon number and norms: if this is an active -// processor these are broadcast to processor #0, -// otherwise these are ignored -// {hn,{rms,infinity}_norm}_buffer[N_active_procs] +// hn,iteration,{rms,infinity}_norm = the values to be broadcast: +// if this is an active processor these +// are broadcast to processor #0, +// otherwise these are ignored +// {hn,iteration,{rms,infinity}_norm}_buffer[N_active_procs] // = buffers for horizon numbers and norms: if this is processor #0 // these are set to the horizon numbers and norms for all active // processors, otherwise these are set to undefined (garbage) values @@ -432,10 +453,9 @@ bool reduce_and_broadcast(cGH* GH, int N_procs, int N_active_procs, int my_proc, bool my_active_flag, bool we_need_more_iterations, - int hn, fp rms_norm, fp infinity_norm, - int hn_buffer[], - fp rms_norm_buffer[], - fp infinity_norm_buffer[]) + int hn, int iteration, fp rms_norm, fp infinity_norm, + int hn_buffer[], int iteration_buffer[], + fp rms_norm_buffer[], fp infinity_norm_buffer[]) { assert( my_proc >= 0 ); assert( my_proc < N_procs ); @@ -446,10 +466,12 @@ assert( my_proc < N_procs ); // doing a sum-reduction of the buffers. // // Alas, to do everything in a single operation, all the values have to -// be of a single datatype in a single array, so we convert hn to -// CCTK_REAL, and encode the "we need more iterations" in its sign. +// be of a single datatype in a single array, so we convert hn and +// iteration to CCTK_REAL, and encode the "we need more iterations" in +// hn's sign. // -const int N_reduce_buffer = 3*N_active_procs; +const int N_vars = 4; +const int N_reduce_buffer = N_vars*N_active_procs; static CCTK_REAL* reduce_buffer = NULL; if (reduce_buffer == NULL) then { @@ -467,9 +489,10 @@ if (reduce_buffer == NULL) if (my_active_flag) then { assert( my_proc < N_active_procs ); - reduce_buffer[3*my_proc + 0] = we_need_more_iterations ? hn : -hn; - reduce_buffer[3*my_proc + 1] = rms_norm; - reduce_buffer[3*my_proc + 2] = infinity_norm; + reduce_buffer[N_vars*my_proc + 0] = we_need_more_iterations ? hn : -hn; + reduce_buffer[N_vars*my_proc + 1] = iteration; + reduce_buffer[N_vars*my_proc + 2] = rms_norm; + reduce_buffer[N_vars*my_proc + 3] = infinity_norm; } // @@ -479,7 +502,7 @@ if (my_active_flag) // this name is appropriate for PUGHReduce, caveat user for other drivers :) const int reduction_handle = CCTK_ReductionArrayHandle("sum"); if (reduction_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "reduce_and_broadcast(): can't get sum-reduction handle!"); /*NOTREACHED*/ @@ -492,7 +515,7 @@ const int status N_reduce_buffer, CCTK_VARIABLE_REAL); if (status < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "reduce_and_broadcast(): error status %d from reduction!", status); /*NOTREACHED*/ @@ -502,12 +525,13 @@ if (status < 0) bool any_proc_needs_more_iterations = false; for (int proc = 0 ; proc < N_active_procs ; ++proc) { - const int i_proc = 3*proc; + const int i_proc = N_vars*proc; const bool proc_needs_more_iterations = reduce_buffer[i_proc + 0] > 0.0; any_proc_needs_more_iterations |= proc_needs_more_iterations; hn_buffer[proc] = (int) fabs(reduce_buffer[i_proc + 0]); - rms_norm_buffer[proc] = reduce_buffer[i_proc + 1]; - infinity_norm_buffer[proc] = reduce_buffer[i_proc + 2]; + iteration_buffer[proc] = (int) reduce_buffer[i_proc + 1]; + rms_norm_buffer[proc] = reduce_buffer[i_proc + 2]; + infinity_norm_buffer[proc] = reduce_buffer[i_proc + 3]; } return any_proc_needs_more_iterations; @@ -523,20 +547,22 @@ return any_proc_needs_more_iterations; // N_active_procs = The number of active processors. // hn_buffer[N_active_procs] = Gives the horizon each active processor // is currently working on. +// iteration_buffer[N_active_procs] = Gives the iteration number of each +// active processor on its current horizon. // {rms,infinity}_norm_buffer[N_active_procs] = Gives the Theta norms for // each active processor. // namespace { -void print_Theta_norms(int N_active_procs, - const int hn_buffer[], - const fp rms_norm_buffer[], - const fp infinity_norm_buffer[]) +void print_Theta_norms + (int N_active_procs, + const int hn_buffer[], const int iteration_buffer[], + const fp rms_norm_buffer[], const fp infinity_norm_buffer[]) { for (int proc = 0 ; proc < N_active_procs ; ++proc) { CCTK_VInfo(CCTK_THORNSTRING, - " proc %d/horizon %d: ||Theta|| rms=%.1e, infty=%.1e", - proc, int(hn_buffer[proc]), + " proc %d/horizon %d/it %d |Theta| rms=%.1e inf=%.1e", + proc, hn_buffer[proc], iteration_buffer[proc], double(rms_norm_buffer[proc]), double(infinity_norm_buffer[proc])); } diff --git a/src/driver/driver.hh b/src/driver/driver.hh index e5ad9aa..225130a 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -28,7 +28,8 @@ enum verbose_level verbose_level__physics_highlights, verbose_level__physics_details, verbose_level__algorithm_highlights, - verbose_level__algorithm_details // no comma + verbose_level__algorithm_details, + verbose_level__algorithm_debug // no comma }; // @@ -122,6 +123,9 @@ struct IO_info bool output_initial_guess; int how_often_to_output_h, how_often_to_output_Theta; + // based on the above, do we want to output things now (this time step)? + bool output_h, output_Theta; + bool output_ghost_zones_for_h; const char* ASCII_gnuplot_file_name_extension; const char* HDF5_file_name_extension; @@ -154,6 +158,7 @@ struct verbose_info bool print_physics_details; bool print_algorithm_highlights; bool print_algorithm_details; + bool print_algorithm_debug; }; //****************************************************************************** @@ -207,6 +212,7 @@ struct AH_info struct state { enum method method; + struct error_info error_info; struct verbose_info verbose_info; int timer_handle; int N_procs; // total number of processors @@ -275,6 +281,7 @@ void Newton(cGH* GH, const struct solver_info& solver_info, const struct IO_info& IO_info, const struct BH_diagnostics_info& BH_diagnostics_info, + const struct error_info& error_info, const struct verbose_info& verbose_info); // BH_diagnostics.cc diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index fd35afd..0324cf1 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -71,7 +71,7 @@ void do_evaluate_expansions(int my_proc, int N_horizons, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct IO_info& IO_info, - bool output_h, bool output_Theta, + const struct error_info& error_info, const struct verbose_info& verbose_info, int timer_handle); void do_test_expansion_Jacobians(int my_proc, int N_horizons, @@ -81,6 +81,7 @@ void do_test_expansion_Jacobians(int my_proc, int N_horizons, struct Jacobian_info& Jac_info, bool test_all_Jacobian_compute_methods, const struct IO_info& IO_info, + const struct error_info& error_info, const struct verbose_info& verbose_info, int timer_handle); @@ -117,6 +118,7 @@ const bool active_flag = hs.has_genuine_horizons(); const struct geometry_info& gi = state.gi; struct Jacobian_info& Jac_info = state.Jac_info; struct IO_info& IO_info = state.IO_info; +const struct error_info& error_info = state.error_info; const struct verbose_info& verbose_info = state.verbose_info; // what are the semantics of the Cactus gxx variables? (these may @@ -125,17 +127,17 @@ if (CCTK_Equals(metric_type, "physical")) then cgi.Cactus_conformal_metric = false; else if (CCTK_Equals(metric_type, "static conformal")) then cgi.Cactus_conformal_metric = (conformal_state > 0); -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!", metric_type); /*NOTREACHED*/ // get the Cactus time step and decide if we want to output h and/or Theta now IO_info.time_iteration = cctk_iteration; IO_info.time = cctk_time; -const bool output_h +IO_info.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_Theta +IO_info.output_Theta = (IO_info.how_often_to_output_Theta > 0) && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0); @@ -176,9 +178,9 @@ switch (state.method) case method__evaluate_expansions: do_evaluate_expansions(my_proc, state.N_horizons, *state.my_hs, state.my_AH_info, - cgi, gi, - IO_info, output_h, output_Theta, - state.verbose_info, state.timer_handle); + cgi, gi, IO_info, + error_info, verbose_info, + state.timer_handle); break; case method__test_expansion_Jacobians: @@ -186,7 +188,8 @@ case method__test_expansion_Jacobians: state.my_AH_info, cgi, gi, Jac_info, (test_all_Jacobian_compute_methods != 0), - IO_info, verbose_info, state.timer_handle); + IO_info, error_info, verbose_info, + state.timer_handle); break; case method__find_horizons: @@ -197,13 +200,13 @@ case method__find_horizons: *state.my_hs, state.my_AH_info, cgi, gi, Jac_info, state.solver_info, IO_info, state.BH_diagnostics_info, - verbose_info); + error_info, verbose_info); if (state.timer_handle >= 0) then CCTK_TimerStopI(state.timer_handle); break; default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " find_horizons(): unknown method=(int)%d!\n" " (this should never happen!)" @@ -213,7 +216,8 @@ default: if (state.timer_handle >= 0) then { - printf("timer stats for computation:\n"); + CCTK_VInfo(CCTK_THORNSTRING, + "timer stats for computation:"); CCTK_TimerPrintDataI(state.timer_handle, -1); } } @@ -279,7 +283,7 @@ const CCTK_REAL* data_ptr = static_cast<const fp*>( ); if (check_for_NULL && (data_ptr == NULL)) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " Cactus_gridfn_data_ptr(): got unexpected NULL data pointer\n" " for Cactus geometry gridfn!\n" @@ -317,7 +321,7 @@ void do_evaluate_expansions(int my_proc, int N_horizons, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct IO_info& IO_info, - bool output_h, bool output_Theta, + const struct error_info& error_info, const struct verbose_info& verbose_info, int timer_handle) { @@ -341,13 +345,14 @@ if (active_flag) jtutil::norm<fp> Theta_norms; const bool Theta_ok = expansion(&ps, cgi, gi, + error_info, true,// initial eval false, // no Jacobian coeffs true, // yes, print msgs &Theta_norms); if (timer_handle >= 0) then CCTK_TimerStopI(timer_handle); - if (output_h) + if (IO_info.output_h) then output_gridfn(ps, gfns::gfn__h, IO_info, IO_info.h_base_file_name, hn, verbose_info.print_algorithm_details); @@ -357,7 +362,7 @@ if (active_flag) CCTK_VInfo(CCTK_THORNSTRING, " Theta(h) rms-norm %.2e, infinity-norm %.2e", Theta_norms.rms_norm(), Theta_norms.infinity_norm()); - if (output_Theta) + if (IO_info.output_Theta) then output_gridfn(ps, gfns::gfn__Theta, IO_info, IO_info .Theta_base_file_name, @@ -369,8 +374,9 @@ if (active_flag) else { for (int i = 0 ; i < N_horizons ; ++i) { - expansion(NULL, // dummy computation - cgi, gi); + expansion(NULL, // dummy computation + cgi, gi, + error_info, true); // initial evaluation } } } @@ -410,6 +416,7 @@ void do_test_expansion_Jacobians(int my_proc, int N_horizons, struct Jacobian_info& Jac_info, bool test_all_Jacobian_compute_methods, const struct IO_info& IO_info, + const struct error_info& error_info, const struct verbose_info& verbose_info, int timer_handle) { @@ -427,10 +434,12 @@ patch_system* ps_ptr = active_flag ? AH_info_ptr->ps_ptr : NULL; // Jacobian* Jac_NP_ptr = active_flag ? AH_info_ptr->Jac_ptr : NULL; expansion(ps_ptr, - cgi, gi); + cgi, gi, + error_info, true); // initial evaluation Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation; expansion_Jacobian(ps_ptr, Jac_NP_ptr, cgi, gi, Jac_info, + error_info, true, // initial evaluation print_msg_flag); Jacobian* Jac_SD_FDdr_ptr = NULL; @@ -444,10 +453,12 @@ if (test_all_Jacobian_compute_methods) : NULL; expansion(ps_ptr, cgi, gi, + error_info, true, // initial evaluation true); // compute SD Jacobian coeffs Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr; expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr, cgi, gi, Jac_info, + error_info, true, // initial evaluation print_msg_flag); } diff --git a/src/driver/horizon_sequence.cc b/src/driver/horizon_sequence.cc index ca47489..fa88423 100644 --- a/src/driver/horizon_sequence.cc +++ b/src/driver/horizon_sequence.cc @@ -4,12 +4,17 @@ // // horizon_sequence::horizon_sequence // horizon_sequence::~horizon_sequence +// horizon_sequence::sequence_string // horizon_sequence::append_hn // horizon_sequence::next_posn // #include <stdio.h> #include <assert.h> + +#include "../include/stdc.h" +#include "../jtutil/util.hh" + #include "horizon_sequence.hh" //****************************************************************************** @@ -38,6 +43,35 @@ delete[] my_hn_; //****************************************************************************** // +// This function returns a (pointer to a) C-style string showing all +// the horizon numbers in the sequence, separated by the C-style string +// sep , eg. "2,3,5". +// +// Note the result points into a private static buffer. +// +char* horizon_sequence::sequence_string(const char sep[]) const +{ +const int N_hn_buffer = 10; +char hn_buffer[N_hn_buffer]; + +const int N_buffer = 100; +static char buffer[N_buffer]; + +buffer[0] = '\0'; + for (int pos = 0 ; pos < my_N_horizons_ ; ++pos) + { + if (pos > 0) + then AHFinderDirect_Strlcat(buffer, sep, N_buffer); + snprintf(hn_buffer, N_hn_buffer, "%d", my_hn_[pos]); + AHFinderDirect_Strlcat(buffer, hn_buffer, N_buffer); + } + +return buffer; +} + +//****************************************************************************** + +// // This function appends hn to the sequence. It returns the new value // of my_N_horizons(). // diff --git a/src/driver/horizon_sequence.hh b/src/driver/horizon_sequence.hh index 398f687..ca363a1 100644 --- a/src/driver/horizon_sequence.hh +++ b/src/driver/horizon_sequence.hh @@ -28,6 +28,11 @@ public: // are there any genuine horizons in the sequence (for this processor)? bool has_genuine_horizons() const { return my_N_horizons_ > 0; } + // C-style string showing all horizon numbers, + // separated by the C-style string sep, eg. "2,3,5" + // ... result points into a private static buffer + char* sequence_string(const char sep[]) const; + // is the current horizon in the sequence dummy/genuine? bool is_dummy () const { return posn_is_dummy (posn_); } bool is_genuine() const { return posn_is_genuine(posn_); } diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc index bc9ab89..7e00390 100644 --- a/src/driver/initial_guess.cc +++ b/src/driver/initial_guess.cc @@ -132,7 +132,7 @@ case initial_guess__coord_ellipsoid: break; default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown initial guess method=(int)%d!", int(igi.method)); /*NOTREACHED*/ } @@ -157,7 +157,7 @@ else if (STRING_EQUAL(initial_guess_method_string, "coordinate sphere")) then return initial_guess__coord_sphere; else if (STRING_EQUAL(initial_guess_method_string, "coordinate ellipsoid")) then return initial_guess__coord_ellipsoid; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " decode_initial_guess_method():\n" " unknown initial_guess_method_string=\"%s\"!", @@ -176,7 +176,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // and the horizon is worked out on page 13.2 of my AHFinderDirect notes. // // Arguments: -// [xyz]_posn = The position of the Kerr black hole. +// [xyz]_posn = The global-coordinate position of the Kerr black hole. // (m,a) = Describe the Kerr black hole. Note that my convention has // a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a). // Kerr_Schild_flag = false to use Kerr coordinates, @@ -215,9 +215,9 @@ if (verbose_info.print_algorithm_details) " setting coordinate %s", Kerr_Schild_flag ? "ellipsoid" : "sphere"); setup_coord_ellipsoid(ps, - x_posn, y_posn, z_posn, - xy_radius, xy_radius, z_radius, - verbose_info.print_algorithm_details); + x_posn, y_posn, z_posn, + xy_radius, xy_radius, z_radius, + verbose_info.print_algorithm_details); } } @@ -227,7 +227,7 @@ setup_coord_ellipsoid(ps, // This function sets up an ellipsoid initial guess, using the formulas // in "ellipsoid.maple" and the Maple-generated C code in "ellipsoid.c": // -// ellipsoid has center (A,B,C), radius (a,b,c) +// ellipsoid has global-coordinates center (A,B,C), radius (a,b,c) // angular coordinate system has center (U,V,W) // // direction cosines wrt angular coordinate center are (xcos,ycos,zcos) @@ -292,7 +292,8 @@ if (print_msg_flag) then r = r_plus; else if ((r_plus < 0.0) && (r_minus > 0.0)) then r = r_minus; - else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + else CCTK_VWarn(FATAL_ERROR, + __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " setup_coord_ellipsoid():\n" " expected exactly one r>0 solution to quadratic, got 0 or 2!\n" diff --git a/src/driver/io.cc b/src/driver/io.cc index 2bddaa9..3db7edf 100644 --- a/src/driver/io.cc +++ b/src/driver/io.cc @@ -88,11 +88,11 @@ case horizon_file_format__ASCII_gnuplot: break; case horizon_file_format__HDF5: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "input_gridfn(): HDF5 data files not implemented yet!"); /*NOTREACHED*/ default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " input_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" @@ -126,7 +126,7 @@ case horizon_file_format__ASCII_gnuplot: case gfns::gfn__h: if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " writing h to \"%s\"", file_name); + "writing h to \"%s\"", file_name); ps.print_ghosted_gridfn_with_xyz (unknown_gfn, true, gfns::gfn__h, @@ -137,19 +137,19 @@ case horizon_file_format__ASCII_gnuplot: case gfns::gfn__Theta: if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " writing Theta to \"%s\"", file_name); + "writing Theta to \"%s\"", 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); + "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 gfn=%d to \"%s\"", + "writing gfn=%d to \"%s\"", unknown_gfn, file_name); ps.print_gridfn(unknown_gfn, file_name); break; @@ -157,11 +157,11 @@ case horizon_file_format__ASCII_gnuplot: break; case horizon_file_format__HDF5: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "output_gridfn(): HDF5 data files not implemented yet!"); /*NOTREACHED*/ default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " output_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" @@ -206,7 +206,7 @@ if (print_msg_flag) FILE *fileptr = fopen(file_name, "w"); if (fileptr == NULL) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "output_Jacobians(): can't open output file \"%s\"!", file_name); /*NOTREACHED*/ @@ -334,7 +334,7 @@ case horizon_file_format__HDF5: file_name_extension = IO_info.HDF5_file_name_extension; break; default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " io_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" diff --git a/src/driver/makefile b/src/driver/makefile index bbc3443..852a6ed 100644 --- a/src/driver/makefile +++ b/src/driver/makefile @@ -1,5 +1,5 @@ # Makefile for driver code -# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/driver/makefile,v 1.3 2003-02-15 18:43:08 jthorn Exp $ +# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/driver/makefile,v 1.4 2003-02-16 22:21:48 jthorn Exp $ # # Targets: @@ -21,8 +21,10 @@ test : test_horizon_sequence .PHONY : test_horizon_sequence test_horizon_sequence : - $(CXX) $(CXXFLAGS) -o $@ test_horizon_sequence.cc horizon_sequence.cc + $(CXX) $(CXXFLAGS) -o $@ \ + test_horizon_sequence.cc horizon_sequence.cc \ + ../jtutil/miscstr.c .PHONY : clean clean : - -rm *core test_horizon_sequence + -rm test_horizon_sequence diff --git a/src/driver/setup.cc b/src/driver/setup.cc index df6e76e..d0b823e 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -105,6 +105,19 @@ CCTK_VInfo(CCTK_THORNSTRING, // state.method = decode_method(method); +state.error_info.warn_level__point_outside__initial + = warn_level__point_outside__initial; +state.error_info.warn_level__point_outside__subsequent + = warn_level__point_outside__subsequent; +state.error_info.warn_level__skipping_finite_check + = warn_level__skipping_finite_check; +state.error_info.warn_level__nonfinite_geometry + = warn_level__nonfinite_geometry; +state.error_info.warn_level__gij_not_positive_definite__initial + = warn_level__gij_not_positive_definite__initial; +state.error_info.warn_level__gij_not_positive_definite__subsequent + = warn_level__gij_not_positive_definite__subsequent; + struct verbose_info& verbose_info = state.verbose_info; verbose_info.verbose_level = decode_verbose_level(verbose_level); verbose_info.print_physics_highlights @@ -115,6 +128,8 @@ verbose_info.print_algorithm_highlights = (state.verbose_info.verbose_level >= verbose_level__algorithm_highlights); verbose_info.print_algorithm_details = (state.verbose_info.verbose_level >= verbose_level__algorithm_details); +verbose_info.print_algorithm_debug + = (state.verbose_info.verbose_level >= verbose_level__algorithm_debug); state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1; @@ -140,7 +155,7 @@ cgi.Cactus_conformal_metric = false; // dummy value, may change later cgi.coord_system_handle = CCTK_CoordSystemHandle(coordinate_system_name); if (cgi.coord_system_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_setup(): can't find Cactus coordinate system \"%s\"!", coordinate_system_name); /*NOTREACHED*/ cgi.g_dd_11_varindex = CCTK_VarIndex("ADMBase::gxx"); @@ -195,12 +210,12 @@ gi.geometry_method = decode_geometry_method(geometry_method); // parameters for geometry_method = interpolation methods gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name); if (gi.operator_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_setup(): couldn't find interpolator \"%s\"!", geometry_interpolator_name); /*NOTREACHED*/ gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars); if (gi.param_table_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_setup(): bad geometry-interpolator parameter(s) \"%s\"!", geometry_interpolator_pars); /*NOTREACHED*/ @@ -249,6 +264,8 @@ IO_info.horizon_file_format = decode_horizon_file_format(horizon_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_Theta = how_often_to_output_Theta; +IO_info.output_h = false; // dummy value +IO_info.output_Theta = false; // dummy value IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0); IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension; IO_info.HDF5_file_name_extension = HDF5_file_name_extension; @@ -308,13 +325,13 @@ if (verbose_info.print_algorithm_highlights) then CCTK_VInfo(CCTK_THORNSTRING, " setting up interpatch interpolator"); const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name); if (interp_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_setup(): couldn't find interpolator \"%s\"!", interpatch_interpolator_name); /*NOTREACHED*/ const int interp_param_table_handle = Util_TableCreateFromString(interpatch_interpolator_pars); if (interp_param_table_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "AHFinderDirect_setup(): bad interpatch-interpolator parameter(s) \"%s\"!", interpatch_interpolator_pars); /*NOTREACHED*/ @@ -441,7 +458,7 @@ else if (STRING_EQUAL(method_string, "test expansion Jacobians")) then return method__test_expansion_Jacobians; else if (STRING_EQUAL(method_string, "find horizons")) then return method__find_horizons; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_method(): unknown method_string=\"%s\"!", method_string); /*NOTREACHED*/ } @@ -465,7 +482,9 @@ else if (STRING_EQUAL(verbose_level_string, "algorithm highlights")) then return verbose_level__algorithm_highlights; else if (STRING_EQUAL(verbose_level_string, "algorithm details")) then return verbose_level__algorithm_details; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else if (STRING_EQUAL(verbose_level_string, "algorithm debug")) + then return verbose_level__algorithm_debug; +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_verbose_level(): unknown verbose_level_string=\"%s\"!", verbose_level_string); /*NOTREACHED*/ } @@ -485,7 +504,7 @@ if (STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)")) then return horizon_file_format__ASCII_gnuplot; else if (STRING_EQUAL(horizon_file_format_string, "HDF5")) then return horizon_file_format__HDF5; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!", horizon_file_format_string); /*NOTREACHED*/ } diff --git a/src/driver/test_horizon_sequence.cc b/src/driver/test_horizon_sequence.cc index 04e0978..d59a0f4 100644 --- a/src/driver/test_horizon_sequence.cc +++ b/src/driver/test_horizon_sequence.cc @@ -2,7 +2,10 @@ // $Header$ #include <stdio.h> +#include <string.h> #include <assert.h> + +#include "../include/stdc.h" #include "horizon_sequence.hh" // prototypes for functions local to this file @@ -34,6 +37,7 @@ horizon_sequence hs0(0); assert( hs0.N_horizons() == 0 ); assert( hs0.my_N_horizons() == 0 ); assert( ! hs0.has_genuine_horizons() ); +assert( STRING_EQUAL(hs0.sequence_string(","), "") ); assert( hs0.is_dummy() ); assert( ! hs0.is_genuine() ); assert( ! hs0.is_next_genuine() ); @@ -69,6 +73,7 @@ assert( hs1.append_hn(42) == 1 ); assert( hs1.N_horizons() == 1 ); assert( hs1.my_N_horizons() == 1 ); assert( hs1.has_genuine_horizons() ); +assert( STRING_EQUAL(hs1.sequence_string(","), "42") ); assert( hs1.is_genuine() ); assert( ! hs1.is_next_genuine() ); assert( hs1.dummy_number() == 0 ); @@ -117,6 +122,7 @@ assert( hs2.append_hn(69) == 2 ); assert( hs2.N_horizons() == 2 ); assert( hs2.my_N_horizons() == 2 ); assert( hs2.has_genuine_horizons() ); +assert( STRING_EQUAL(hs2.sequence_string(","), "42,69") ); assert( hs2.is_genuine() ); assert( hs2.is_next_genuine() ); assert( hs2.dummy_number() == 0 ); @@ -182,6 +188,8 @@ void test345() assert( hs.is_next_genuine() ); assert( hs.dummy_number() == 0 ); + assert( STRING_EQUAL(hs.sequence_string(","), "42,69,105") ); + for (int i = 1 ; i <= 4 ; ++i) { printf(" try %d...\n", i); diff --git a/src/elliptic/dense_Jacobian.cc b/src/elliptic/dense_Jacobian.cc index 1dd6fbb..3e4e582 100644 --- a/src/elliptic/dense_Jacobian.cc +++ b/src/elliptic/dense_Jacobian.cc @@ -133,7 +133,7 @@ dense_Jacobian::dense_Jacobian(patch_system& ps, { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " dense Jacobian matrix (N_rows_=%d)", + " dense Jacobian matrix (%d rows)", N_rows_); } #endif // HAVE_DENSE_JACOBIAN diff --git a/src/elliptic/row_sparse_Jacobian.cc b/src/elliptic/row_sparse_Jacobian.cc index 875ab7d..4927000 100644 --- a/src/elliptic/row_sparse_Jacobian.cc +++ b/src/elliptic/row_sparse_Jacobian.cc @@ -146,7 +146,7 @@ row_sparse_Jacobian::row_sparse_Jacobian(patch_system& ps, { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " row sparse matrix Jacobian (N_rows_=%d)", + " row sparse matrix Jacobian (%d rows)", N_rows_); zero_matrix(); @@ -459,7 +459,7 @@ integer istatus, ierror; #endif if (ierror != 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "***** row_sparse_Jacobian__ILUCG::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n" " error return istatus=%d from [sd]ilucg() routine!" , diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc index 490ff99..1946cfa 100644 --- a/src/gr/expansion.cc +++ b/src/gr/expansion.cc @@ -59,20 +59,25 @@ void setup_xyz_posns(patch_system& ps, bool print_msg_flag); bool interpolate_geometry(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag); -void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag); +void convert_conformal_to_physical(patch_system& ps, + bool print_msg_flag); void Schwarzschild_EF_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, bool print_msg_flag); -bool h_is_finite(patch_system& ps, bool print_msg_flag); +bool h_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); bool geometry_is_finite(patch_system& ps, - bool print_msg_flag); + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); -bool compute_Theta(patch_system& ps, - bool Jacobian_flag, +bool compute_Theta(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag); } @@ -122,6 +127,11 @@ bool compute_Theta(patch_system& ps, // interpolator call are done, the latter with the number of // interpolation points is set to 0 and all the output array // pointers set to NULL. +// initial_flag = true if this is the first evaluation of expansion() +// for this horizon, +// false otherwise; +// this is used (only) to select which elements of error_info +// are relevant // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. // print_msg_flag = true to print status messages, @@ -145,6 +155,7 @@ bool compute_Theta(patch_system& ps, bool expansion(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool Jacobian_flag /* = false */, bool print_msg_flag /* = false */, jtutil::norm<fp>* Theta_norms_ptr /* = NULL */) @@ -164,8 +175,9 @@ if (active_flag) // fill in values of all ghosted gridfns in ghost zones ps_ptr->synchronize(); - if (gi.check_that_h_is_finite - && !h_is_finite(*ps_ptr, print_msg_flag)) + if (gi.check_that_h_is_finite && !h_is_finite(*ps_ptr, + error_info, initial_flag, + print_msg_flag)) then return false; // *** ERROR RETURN *** // set up xyz positions of grid points @@ -181,15 +193,19 @@ case geometry__local_interp_from_Cactus_grid: // ps_ptr (non-NULL vs NULL) to choose a normal vs dummy computation if (!interpolate_geometry(ps_ptr, cgi, gi, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** if (active_flag && cgi.Cactus_conformal_metric) - then convert_conformal_to_physical(*ps_ptr, print_msg_flag); + then convert_conformal_to_physical(*ps_ptr, + print_msg_flag); break; case geometry__Schwarzschild_EF: if (active_flag) - then Schwarzschild_EF_geometry(*ps_ptr, gi, print_msg_flag); + then Schwarzschild_EF_geometry(*ps_ptr, + gi, + print_msg_flag); break; default: @@ -203,7 +219,9 @@ default: if (active_flag) then { if (gi.check_that_geometry_is_finite - && !geometry_is_finite(*ps_ptr, print_msg_flag)) + && !geometry_is_finite(*ps_ptr, + error_info, initial_flag, + print_msg_flag)) then return false; // *** ERROR RETURN *** // compute remaining gridfns --> $\Theta$ @@ -211,6 +229,7 @@ if (active_flag) // by algebraic ops and angular finite differencing if (!compute_Theta(*ps_ptr, Jacobian_flag, Theta_norms_ptr, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** } @@ -322,6 +341,7 @@ namespace { bool interpolate_geometry(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { const bool active_flag = (ps_ptr != NULL); @@ -651,19 +671,26 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE) const char axis = "xyz"[error_axis]; assert((error_direction == -1) || (error_direction == +1)); - const char end = (error_direction == -1) ? '-' : '+'; + const char direction = (error_direction == -1) ? '-' : '+'; if (print_msg_flag) then { - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, + const int warn_level + = initial_flag + ? error_info.warn_level__point_outside__initial + : error_info.warn_level__point_outside__subsequent; + CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" the trial-horizon-surface point" -" global_(x,y,z)=(%g,%g,%g)" -" is outside the grid (or too close to the boundary) in the %c%c direction!" +"interpolate_geometry():\n" +" the trial-horizon-surface point\n" +" global_(x,y,z)=(%g,%g,%g)\n" +" is outside the grid (or too close to the boundary)" + " in the %c%c direction!\n" +" (this may mean that driver::ghost_size is too small,\n" +" especially if any horizon crosses a symmetry boundary)" , global_x, global_y, global_z, - end, axis); + direction, axis); } return false; // *** ERROR RETURN *** } @@ -811,7 +838,7 @@ if (print_msg_flag) // or more NaNs or infinities. // // Results: -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE // This function returns true if all the h values are finite, false // otherwise (i.e. if it contains any NaNs or infinities). #else @@ -820,12 +847,14 @@ if (print_msg_flag) #endif // namespace { -bool h_is_finite(patch_system& ps, bool print_msg_flag) +bool h_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " checking that h is finite"); -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -843,7 +872,7 @@ if (print_msg_flag) const fp sigma = p.sigma_of_isigma(isigma); const fp drho = jtutil::degrees_of_radians(rho); const fp dsigma = jtutil::degrees_of_radians(sigma); - CCTK_VWarn(expansion_warning_levels::nonfinite, + CCTK_VWarn(error_info.warn_level__nonfinite_geometry, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " h=%g isn't finite!\n" @@ -859,7 +888,7 @@ if (print_msg_flag) } return true; // *** all values finite *** #else -CCTK_VWarn(expansion_warning_levels::skip_nonfinite, +CCTK_VWarn(error_info.warn_level__skipping_finite_check, __LINE__, __FILE__, CCTK_THORNSTRING, " no finite() fn ==> skipping is-h-finite check"); return true; // *** no check possible *** @@ -878,7 +907,7 @@ return true; // *** no check possible *** // or more NaNs or infinities. // // Results: -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE // This function returns true if all the geometry variables are finite, // false otherwise (i.e. if they contain any NaNs or infinities). #else @@ -887,12 +916,14 @@ return true; // *** no check possible *** #endif // namespace { -bool geometry_is_finite(patch_system& ps, bool print_msg_flag) +bool geometry_is_finite(patch_system& ps, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " checking that geometry is finite"); -#ifdef Theta_AVE_FINITE +#ifdef HAVE_FINITE for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -989,7 +1020,7 @@ if (print_msg_flag) const fp global_x = ps.global_x_of_local_x(local_x); const fp global_y = ps.global_y_of_local_y(local_y); const fp global_z = ps.global_z_of_local_z(local_z); - CCTK_VWarn(expansion_warning_levels::nonfinite, + CCTK_VWarn(error_info.warn_level__nonfinite_geometry, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " geometry isn't finite at %s patch\n" @@ -1041,7 +1072,7 @@ if (print_msg_flag) } return true; // *** no NaNs found *** #else -CCTK_VWarn(expansion_warning_levels::skip_nonfinite, +CCTK_VWarn(error_info.warn_level__skipping_finite_check, __LINE__, __FILE__, CCTK_THORNSTRING, " no finite() ==> skipping is-geometry-finite check"); return true; // *** no check possible *** @@ -1073,9 +1104,9 @@ return true; // *** no check possible *** // g_ij isn't positive definite). // namespace { -bool compute_Theta(patch_system& ps, - bool Jacobian_flag, +bool compute_Theta(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { if (print_msg_flag) @@ -1159,16 +1190,16 @@ if (print_msg_flag) if (Theta_D <= 0) then { - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "Theta_D <= 0 at %s patch rho=%g sigma=%g!", - p.name(), double(rho), double(sigma)); - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - "(this probably means the interpolated g_ij"); - CCTK_VWarn(expansion_warning_levels::failure, - __LINE__, __FILE__, CCTK_THORNSTRING, - " isn't positive definite)"); +const int warn_level + = initial_flag ? error_info.warn_level__gij_not_positive_definite__initial + : error_info.warn_level__gij_not_positive_definite__subsequent; +CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" compute_Theta(): Theta_D = $g^{ij} s_i s_j$ = %g <= 0\n" +" at %s patch rho=%g sigma=%g!\n" +" (i.e. the interpolated g_ij isn't positive definite)", + double(Theta_D), + p.name(), double(rho), double(sigma)); return false; // *** ERROR RETURN *** } diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc index 805ed97..ef2f7e8 100644 --- a/src/gr/expansion_Jacobian.cc +++ b/src/gr/expansion_Jacobian.cc @@ -48,12 +48,13 @@ using jtutil::error_exit; // namespace { -bool expansion_Jacobian_NP(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& ii, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); +bool expansion_Jacobian_NP + (patch_system& ps, Jacobian& Jac, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); void expansion_Jacobian_partial_SD(patch_system& ps, Jacobian& Jac, @@ -67,12 +68,13 @@ void add_ghost_zone_Jacobian(const patch_system& ps, const patch& xp, const ghost_zone& xmgz, int x_II, int xm_irho, int xm_isigma); -bool expansion_Jacobian_dr_FD(patch_system* ps_ptr, - Jacobian* Jac_ptr, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); +bool expansion_Jacobian_dr_FD + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag); } //****************************************************************************** @@ -108,6 +110,7 @@ bool expansion_Jacobian(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag /* = false */) { const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL); @@ -119,24 +122,23 @@ case Jacobian__numerical_perturbation: then { if (! expansion_Jacobian_NP(*ps_ptr, *Jac_ptr, cgi, gi, Jacobian_info, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** break; } - else CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" expansion_Jacobian():\n" + else error_exit(ERROR_EXIT, +"***** expansion_Jacobian():\n" " dummy computation isn't supported for\n" -" Jacobian_compute_method = \"numerical perturbation\"!");/*NOTREACHED*/ +" Jacobian_compute_method = \"numerical perturbation\"!\n"); + /*NOTREACHED*/ case Jacobian__symbolic_diff: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" expansion_Jacobian():\n" + error_exit(ERROR_EXIT, +"***** expansion_Jacobian():\n" " Jacobian_compute_method == \"symbolic differentiation\"\n" " isn't implemented (yet); reverting to\n" -" \"symbolic differentiation with finite diff d/dr\""); - // fall through +" \"symbolic differentiation with finite diff d/dr\"\n");/*NOTREACHED*/ case Jacobian__symbolic_diff_with_FD_dr: if (active_flag) @@ -147,6 +149,7 @@ case Jacobian__symbolic_diff_with_FD_dr: // to choose a normal vs dummy computation if (! expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, cgi, gi, Jacobian_info, + error_info, initial_flag, print_msg_flag)) then return false; // *** ERROR RETURN *** break; @@ -202,12 +205,13 @@ return true; // *** NORMAL RETURN *** // for details of the various ways the computation may fail. // namespace { -bool expansion_Jacobian_NP(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& ii, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) +bool expansion_Jacobian_NP + (patch_system& ps, Jacobian& Jac, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, @@ -234,7 +238,9 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma); yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon; - if (! expansion(&ps, cgi, ii)) + if (! expansion(&ps, + cgi, gi, + error_info, initial_flag)) then return false; // *** ERROR RETURN *** for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) @@ -523,12 +529,13 @@ patch& yp = ye.my_patch(); // FIXME: excision isn't implemented yet :( // namespace { -bool expansion_Jacobian_dr_FD(patch_system* ps_ptr, - Jacobian* Jac_ptr, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) +bool expansion_Jacobian_dr_FD + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, + bool print_msg_flag) { const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL); if (print_msg_flag) @@ -544,7 +551,9 @@ if (active_flag) ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h); } -if (! expansion(ps_ptr, cgi, gi)) +if (! expansion(ps_ptr, + cgi, gi, + error_info, initial_flag)) then return false; // *** ERROR RETURN *** if (active_flag) diff --git a/src/gr/gr.hh b/src/gr/gr.hh index 3014cd7..58275f5 100644 --- a/src/gr/gr.hh +++ b/src/gr/gr.hh @@ -42,9 +42,9 @@ struct cactus_grid_info cGH *GH; // --> Cactus grid hierarchy // this describes the semantics of the Cactus gij gridfns: - // true ==> the Cactus gij are conformal values as per + // true ==> the Cactus g_ij are conformal values as per // CactusEinstein/StaticConformal and the psi gridfn - // false ==> the Cactus gij are the physical metric + // false ==> the Cactus g_ij are the physical metric bool Cactus_conformal_metric; // @@ -136,6 +136,40 @@ struct Jacobian_info fp perturbation_amplitude; }; +// +// This struct holds information on what to do if various error/warning +// conditions occur. +// +struct error_info + { + // at present, we "handle" all these errors/warnings by CCTK_VWarn(); + // these parameters give the warning levels for various conditions: + + // point outside (or too close to boundary of) Cactus grid, + // i.e. geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE + // ... warning level if error occurs on first Newton iteration, + // i.e. when evaluating expansion/Jacobian for initial guess + int warn_level__point_outside__initial; + // ... warning level if error occurs on a subsequent Newton iteration + int warn_level__point_outside__subsequent; + + // the Cactus configure process didn't find a finite() function, + // so we're skipping the "check that the geometry is finite" tests + // even though the user set check_that_geometry_is_finite = "true" + int warn_level__skipping_finite_check; + + // infinity or NaN in interpolated {g_ij, partial_k g_ij, K_ij} + int warn_level__nonfinite_geometry; + + // interpolated g_ij isn't positive definite + // (usually this means we're too close to a singularity) + // ... warning level if error occurs on first Newton iteration, + // i.e. when evaluating expansion/Jacobian for initial guess + int warn_level__gij_not_positive_definite__initial; + // ... warning level if error occurs on a subsequent Newton iteration + int warn_level__gij_not_positive_definite__subsequent; + }; + //****************************************************************************** // @@ -148,18 +182,10 @@ struct Jacobian_info // ... returns false for failure (eg horizon outside grid); // if (print_msg_flag) then also reports CCTK_VWarn() at the // appropriate warning level -namespace expansion_warning_levels - { - enum { - failure = 2, // point outside grid etc - nonfinite = 2, // nonfinite gridfn detected - skip_nonfinite = 3 // no finite() function available - // ==> skipping nonfinite check - }; - } bool expansion(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, + const struct error_info& error_info, bool initial_flag, bool Jacobian_flag = false, bool print_msg_flag = false, jtutil::norm<fp>* H_norms_ptr = NULL); @@ -173,6 +199,7 @@ bool expansion_Jacobian(patch_system* ps_ptr, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, + const struct error_info& error_info, bool initial_flag, bool print_msg_flag = false); // Schwarzschild_EF.cc diff --git a/src/gr/misc.cc b/src/gr/misc.cc index 8a0ced4..846ee5d 100644 --- a/src/gr/misc.cc +++ b/src/gr/misc.cc @@ -56,7 +56,7 @@ else if (STRING_EQUAL(geometry_method_string, else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF")) then return geometry__Schwarzschild_EF; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_geometry_method(): unknown geometry_method_string=\"%s\"!", geometry_method_string); /*NOTREACHED*/ } @@ -78,7 +78,7 @@ case geometry__local_interp_from_Cactus_grid: case geometry__Schwarzschild_EF: return false; default: - CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "geometry_method_is_interp(): unknown geometry_method=(int)%d!", geometry_method); /*NOTREACHED*/ } @@ -103,7 +103,7 @@ else if (STRING_EQUAL(Jacobian_compute_method_string, else if (STRING_EQUAL(Jacobian_compute_method_string, "symbolic differentiation")) then return Jacobian__symbolic_diff; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " decode_Jacobian_compute_method():\n" " unknown Jacobian_compute_method_string=\"%s\"!", diff --git a/src/include/config.h b/src/include/config.h index a5334e9..e761015 100644 --- a/src/include/config.h +++ b/src/include/config.h @@ -32,6 +32,13 @@ typedef CCTK_INT integer; * definitions for C++, C, and Fortran */ +/* CCTK_VWarn() "warning level" for fatal errors (terminate execution) */ +/* note low-level software uses error_exit() instead of CCTK_VWarn(), */ +/* with exit codes defined in "config.h" */ +#define FATAL_ERROR (-1) + +/**************************************/ + /* * FIXME: this assumes fp == C 'double' * CCTK_REAL_PRECISION_{4,8,16} are helpful, but not quite enough) @@ -76,7 +83,7 @@ typedef CCTK_INT integer; /* store as row-oriented sparse matrix, solve with ILUCG */ #define HAVE_ROW_SPARSE_JACOBIAN__ILUCG -/**************************************/ +/******************************************************************************/ /* * The following #ifdefs and #defines shouldn't need changing; they diff --git a/src/include/stdc.h b/src/include/stdc.h index 972a870..872c012 100644 --- a/src/include/stdc.h +++ b/src/include/stdc.h @@ -12,8 +12,8 @@ /******************************************************************************/ /* - * I use this empty macro in all if statements -- I like the symmetry - * between the two branches of an if statement. + * I use this macro in all if statements -- I like the symmetry between + * the two branches of an if statement. * * if (foo) * then bar; diff --git a/src/jtutil/array.cc b/src/jtutil/array.cc index 5eb841a..c9d4176 100644 --- a/src/jtutil/array.cc +++ b/src/jtutil/array.cc @@ -17,7 +17,8 @@ // #include <assert.h> -#include <stddef.h> // for NULL +#include <stddef.h> // NULL +#include <stdlib.h> // size_t // we want to instantiate templates with CCTK_* types #ifdef STANDALONE_TEST diff --git a/src/jtutil/fuzzy.cc b/src/jtutil/fuzzy.cc index 34cffa2..d3e4314 100644 --- a/src/jtutil/fuzzy.cc +++ b/src/jtutil/fuzzy.cc @@ -10,6 +10,8 @@ // ***** template instantiations and specializations ***** // +#include <stdlib.h> + #include "stdc.h" #include "util.hh" diff --git a/src/jtutil/make.code.defn b/src/jtutil/make.code.defn index a3ecb79..29bf72e 100644 --- a/src/jtutil/make.code.defn +++ b/src/jtutil/make.code.defn @@ -2,7 +2,8 @@ # $Header$ # Source files in this directory -SRCS = array.cc cpm_map.cc fuzzy.cc linear_map.cc miscfp.cc norm.cc round.cc \ +SRCS = array.cc cpm_map.cc fuzzy.cc linear_map.cc norm.cc round.cc \ + miscstr.c miscfp.cc \ error_exit.cc # Subdirectories containing source files diff --git a/src/jtutil/makefile b/src/jtutil/makefile index ca3d738..abab2fb 100644 --- a/src/jtutil/makefile +++ b/src/jtutil/makefile @@ -1,5 +1,5 @@ # Makefile for standalone test drivers in this directory -# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/jtutil/makefile,v 1.9 2002-10-11 17:35:11 jthorn Exp $ +# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/jtutil/makefile,v 1.10 2003-02-16 22:21:49 jthorn Exp $ # # CC, CXX = C and C++ compilers. Defaults are gcc and g++ if # variables aren't set from command line or environment. @@ -26,7 +26,8 @@ ALL_TESTS := test_array test_array2 \ test_cpm_map test_linear_map \ test_fuzzy test_round \ test_modulo test_norm \ - test_error_exit + test_error_exit \ + test_strlcat ################################################################################ @@ -48,6 +49,8 @@ test_norm : test_norm.o norm.o \ fuzzy.o round.o -lm test_error_exit : test_error_exit.o error_exit.o +test_strlcat : test_strlcat.o miscstr.o + .PHONY : clean clean : -rm -f *.o diff --git a/src/jtutil/miscfp.cc b/src/jtutil/miscfp.cc index fb59352..0015644 100644 --- a/src/jtutil/miscfp.cc +++ b/src/jtutil/miscfp.cc @@ -8,6 +8,7 @@ // #include <math.h> +#include <stdlib.h> #include "stdc.h" #include "util.hh" diff --git a/src/jtutil/miscstr.c b/src/jtutil/miscstr.c new file mode 100644 index 0000000..6553b55 --- /dev/null +++ b/src/jtutil/miscstr.c @@ -0,0 +1,67 @@ +/* miscstr.c -- misc string routines */ +/* $Header$ */ + +#include <stdlib.h> +#include <string.h> + +#include "../include/stdc.h" + +/* + * exceptionally, this (C++) header file is safe to #include from C, + * and contains a prototype for AHFinderDirect_Strlcat() . + */ +#include "util.hh" + +/******************************************************************************/ + +/*@@ + @routine AHFinderDirect_Strlcat + @date 16.Feb.2003 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc This function implements the strcat() function + described in + http://www.openbsd.org/papers/strlcpy-paper.ps + + The strlcat(3) function appends the null-terminated string + src to the end of dst. It will append at most + size - strlen(dst) - 1 + characters, and always null-terminates the result. + (Hence this function never overflows the destination buffer.) + + strlcat(3) is a replacement for strncat(3). In comparison + to strncat(3), strlcat(3) is safer and easier to use: + it guarantees null termination of the destination buffer, + and the size parameter is easy to specify without danger + of off-by-one errors. + @enddesc + + @var dst + @vdesc A non-null pointer to the destination buffer. + @vtype char* dst + @endvar + + @var src + @vdesc A non-null pointer to the source string. + @vtype const char* dst + @endvar + + @var dst_size + @vdesc The size of the destination buffer. + @vtype size_t dst_size + @endvar + + @returntype size_t + @returndesc This function returns the length of the string it + tries to create, i.e. strlen(src) + strlen(dst() . + @endreturndesc + @@*/ +size_t AHFinderDirect_Strlcat(char* dst, const char* src, size_t dst_size) +{ +const size_t src_len = strlen(src); +const size_t dst_len = strlen(dst); +const int dst_remaining = dst_size - dst_len - 1; +if (dst_remaining > 0) + then strncat(dst, src, dst_remaining); + +return src_len + dst_len; +} diff --git a/src/jtutil/norm.cc b/src/jtutil/norm.cc index 23cac5c..edd3fd1 100644 --- a/src/jtutil/norm.cc +++ b/src/jtutil/norm.cc @@ -7,6 +7,7 @@ #include <math.h> #include <assert.h> +#include <stdlib.h> #include "util.hh" diff --git a/src/jtutil/round.cc b/src/jtutil/round.cc index b990ee0..6559944 100644 --- a/src/jtutil/round.cc +++ b/src/jtutil/round.cc @@ -9,6 +9,8 @@ // ***** template instantiations ***** // +#include <stdlib.h> + #include "stdc.h" #include "util.hh" diff --git a/src/jtutil/test_strlcat.c b/src/jtutil/test_strlcat.c new file mode 100644 index 0000000..0c69f57 --- /dev/null +++ b/src/jtutil/test_strlcat.c @@ -0,0 +1,96 @@ +/* test_strlcat -- test driver for Util_Strlcpy() */ +/* $Header$ */ + +#include <string.h> +#include <stdio.h> + +#include "../include/stdc.h" + +/* + * exceptionally, this (C++) header file is safe to #include from C, + * and contains a prototype for AHFinderDirect_Strlcat() . + */ +#include "util.hh" + +/******************************************************************************/ + +/* prototypes */ +size_t tryit(size_t dst_size, const char* src); +void nprint(int n_print, const char* buf); + +/* global data structures */ +static char buffer[100]; + +/******************************************************************************/ + +/* + * This program is a test driver for Util_Strlcpy() . + */ + +int main(void) +{ +size_t n; + +n = tryit(15, "world"); +printf("bufsize=15: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(11, "world"); +printf("bufsize=11: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(10, "world"); +printf("bufsize=10: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(9, "world"); +printf("bufsize=9: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(6, "world"); +printf("bufsize=6: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(5, "world"); +printf("bufsize=5: result=%d buffer=", (int) n); +nprint(20, buffer); + +n = tryit(4, "world"); +printf("bufsize=4: result=%d buffer=", (int) n); +nprint(20, buffer); + +return 0; +} + +/******************************************************************************/ + +size_t tryit(size_t dst_size, const char* src) +{ +const char hello[] = "hello\0xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"; +memcpy(buffer, hello, sizeof(hello)); +return AHFinderDirect_Strlcat(buffer, src, dst_size); +} + +/******************************************************************************/ + +/* print n_print characters of buf[], with visible indication of '\0' */ +void nprint(int n_print, const char* buf) +{ +int i; +int i_null = -1; + +printf("\""); + for (i = 0 ; i < n_print ; ++i) + { + if (buf[i] == '\0') + then { + if (i_null == -1) + then i_null = i; + printf("\\0"); + } + else printf("%c", buf[i]); + } +if (i_null >= 0) + then printf(" [null at i=%d]", i_null); +printf("\"\n"); +} diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh index df3daeb..d87330d 100644 --- a/src/jtutil/util.hh +++ b/src/jtutil/util.hh @@ -3,9 +3,22 @@ // // prerequisites: +// <stdlib.h> or <string.h> or <stdio.h> // size_t // "stdc.h" // +// +// misc string functions (in C, not C++) +// +#ifdef __cplusplus +extern "C" +#endif +size_t AHFinderDirect_Strlcat(char* dst, const char* src, size_t dst_size); + +// +// the rest of this file is C++ only +// +#ifdef __cplusplus namespace jtutil { @@ -185,3 +198,4 @@ public: //****************************************************************************** } // namespace jtutil +#endif /* __cplusplus */ diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc index 2f65378..3b3450b 100644 --- a/src/patch/test_patch_system.cc +++ b/src/patch/test_patch_system.cc @@ -150,13 +150,13 @@ DECLARE_CCTK_PARAMETERS CCTK_VInfo(CCTK_THORNSTRING, "setting up interpatch interpolator"); const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name); if (interp_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "couldn't find interpolator \"%s\"!", interpatch_interpolator_name); /*NOTREACHED*/ const int interp_par_table_handle = Util_TableCreateFromString(interpatch_interpolator_pars); if (interp_par_table_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "bad interpatch-interpolator parameter(s) \"%s\"!", interpatch_interpolator_pars); /*NOTREACHED*/ @@ -232,7 +232,7 @@ else if (STRING_EQUAL(which_test, "derivatives")) ps.print_gridfn(nominal_error_gfn, "nominal_error.dat"); } -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +else CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown which_test=\"%s\"!", which_test); /*NOTREACHED*/ @@ -316,7 +316,7 @@ setup_sym_fn_xyz(ps, NP_test_gfn, false); FILE *fileptr = fopen(Jacobian_file_name, "w"); if (fileptr == NULL) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "test_ghost_zone_Jacobians(): can't open output plot file \"%s\"!", Jacobian_file_name); /*NOTREACHED*/ fprintf(fileptr, "# column 1 = x patch number\n"); |