diff options
| -rw-r--r-- | interface.ccl | 6 | ||||
| -rw-r--r-- | param.ccl | 127 | ||||
| -rw-r--r-- | src/driver/Newton.cc | 2 | ||||
| -rw-r--r-- | src/driver/driver.hh | 70 | ||||
| -rw-r--r-- | src/driver/find_horizons.cc | 98 | ||||
| -rw-r--r-- | src/driver/initial_guess.cc | 2 | ||||
| -rw-r--r-- | src/driver/io.cc | 114 | ||||
| -rw-r--r-- | src/driver/setup.cc | 22 |
8 files changed, 324 insertions, 117 deletions
diff --git a/interface.ccl b/interface.ccl index df8b800..8e2dfe6 100644 --- a/interface.ccl +++ b/interface.ccl @@ -16,7 +16,8 @@ protected: # all the remaining diagnostics are arrays subscripted by the # (1-origin) horizon number hn ; subscript 0 is unused # -int BH_diagnostics__int TYPE=array DIM=1 SIZE=N_horizons+1 DISTRIB=constant +int BH_diagnostics__int TYPE=array DIM=1 SIZE=N_horizons+1 \ + DISTRIB=constant TIMELEVELS=1 { # this is a Boolean flag: 0=false, 1=true AH_found # was this AH found the last time we searched for it? @@ -27,7 +28,8 @@ AH_found # was this AH found the last time we searched for it? # these diagnostics are only defined for those apparent horizons # where AH_found[hn] is true # -real BH_diagnostics__real TYPE=array DIM=1 SIZE=N_horizons+1 DISTRIB=constant +real BH_diagnostics__real TYPE=array DIM=1 SIZE=N_horizons+1 \ + DISTRIB=constant TIMELEVELS=1 { # FIXME: it would be nice to gather these together in some way centroid_x, centroid_y, centroid_z # centroid position @@ -3,6 +3,10 @@ ################################################################################ +# +# ***** interface to the rest of Cactus ***** +# + # we may need to look at grid::domain to choose our patch system symmetries shares: grid USES KEYWORD domain @@ -18,7 +22,7 @@ private: ################################################################################ # -# overall parameters for the apparent horizon finding algorithm +# ***** overall parameters ***** # boolean find_AHs_at_postinitial \ @@ -94,7 +98,7 @@ boolean print_timing_stats \ ################################################################################ # -# parameters for the Jacobian matrix +# ***** parameters for the Jacobian matrix ***** # keyword Jacobian_method "how do we compute the Jacobian matrix?" @@ -153,7 +157,7 @@ boolean test_all_Jacobian_methods \ ################################################################################ # -# parameters for the Newton's-method solution of H(h) = 0 +# ***** parameters for the Newton's-method solution of H(h) = 0 ***** # # @@ -215,22 +219,14 @@ boolean final_H_update_if_exit_x_H_small \ ################################################################################ # -# I/O parameters +# ***** I/O parameters ***** # -# this is mainly useful for debugging purposes -boolean output_initial_guess \ - "should we output the initial guess back to the h data file?" -{ -} "false" +######################################## -# for debugging convergence failures, 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" +# +# parameters for horizon-shape and other similar data files +# # the next two parameters control how often we write full-sized output # files (format controlled by horizon_file_format ); we still *find* @@ -270,11 +266,11 @@ boolean output_ghost_zones_for_h \ string ASCII_gnuplot_file_name_extension \ "extension for ASCII (gnuplot) data files" { -.* :: "any string" +.+ :: "any nonempty string" } "gp" string HDF5_file_name_extension "extension for HDF5 data files" { -.* :: "any string" +.+ :: "any nonempty string" } "hdf5" # @@ -282,7 +278,7 @@ string HDF5_file_name_extension "extension for HDF5 data files" # file names being given by a printf() format "%s.t%d.ah%d[.it%d].%s", # where # - the first %s is the base file name, -# - the first %d is the global Cactus time iteration number, +# - the first %d is the global Cactus time iteration number cctk_iteration, # - the second %d is the apparent horizon number, ande # - the optional third %d is the horizon finder iteration number # - the second %s is the file name extension as set by the @@ -293,18 +289,85 @@ string h_base_file_name \ "base file name for horizon shape h input/output file(s)" { .+ :: "any nonempty string" -} "h.dat" +} "h" string H_of_h_base_file_name "base file name for H(h) output file(s)" { .+ :: "any nonempty string" -} "H.dat" +} "H" string Delta_h_base_file_name \ "base file name for horizon-shape-update Delta_h output file(s)" { .+ :: "any nonempty string" -} "Delta_h.dat" +} "Delta_h" + +######################################## + +# +# parameters for BH diagnostics +# + +# +# The file format is currently hard-wired to a simple ASCII format: +# - there is one file per AH +# - after some header comments (starting with '#'), +# there is one line per successful-AH-finding +# (there is *no* line if we fail to find an AH) +# - each line contains the tab-separated fields +# cctk_iteration +# cctk_time +# centroid_x +# centroid_y +# centroid_z +# area +# m_irreducible +# +boolean output_BH_diagnostics \ + "should we output BH diagnostics to a data file for each AH found?" +{ +} "true" + +# +# These file names are actually just "base" file names, with the full +# file names being given by a printf() format "%s.ah%d.%s", +# where +# - the first %s is the base file name, +# - the %d is the global Cactus time iteration number cctk_iteration, +# - the second %s is the file name extension as set by the +# BH_diagnostics_file_name_extension parameter +# +string BH_diagnostics_base_file_name \ + "base file name for BH diagnostics output file(s)" +{ +.+ :: "any nonempty string" +} "BH_diagnostics" + +string BH_diagnostics_file_name_extension \ + "extension for BH diagnostics data files" +{ +.+ :: "any nonempty string" +} "gp" + +######################################## + +# +# parameters mainly for debugging +# + +# this is mainly useful for debugging purposes +boolean output_initial_guess \ + "should we output the initial guess back to the h data file?" +{ +} "false" + +# for debugging convergence failures, 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" string Jacobian_base_file_name "base file name for Jacobian output file(s)" { @@ -314,9 +377,8 @@ string Jacobian_base_file_name "base file name for Jacobian output file(s)" ################################################################################ # -# parameters to define the patch systems +# ***** parameters to define the patch system(s) ***** # -private: # # For each apparent horizon, you need to set these parameters to the @@ -420,7 +482,7 @@ real delta_drho_dsigma "angular grid spacing of patches, in degrees" ################################################################################ # -# parameters for how we compute surface integrals over the horizon +# ***** parameters for computing surface integrals over the horizon ***** # # ... N is the number of grid zones in a patch, in either direction @@ -445,8 +507,8 @@ keyword surface_integral_method \ ################################################################################ # -# parameters for how we compute the slice's geometry -# (gij, Kij, partial_k gij) +# ***** parameters for how we compute the slice's geometry ***** +# ***** (gij, Kij, partial_k gij) ***** # keyword geometry_method "how do we compute the slice's geometry?" @@ -569,7 +631,9 @@ boolean check_that_geometry_is_finite \ ################################################################################ # -# parameters for the interpatch interpolator +# ***** parameters for the interpatch interpolator ***** +# + # # This 1D interpolator is used to interpolate the h function between # angular patches. Because any given patch boundary only interpolates @@ -595,7 +659,10 @@ string interpatch_interpolator_pars \ ################################################################################ # -# parameters specifying the initial guess for the apparent horizon shape +# ***** parameters for the initial guess for the apparent horizon shape ***** +# + +# # (Note that if at any time we fail to find the (an) apparent horizon, # then we reset our trial horizon surface to this initial guess before # next attempting to find this horizon.) @@ -714,7 +781,7 @@ real initial_guess__coord_ellipsoid__z_radius[5] "z radius of ellipsoid" ################################################################################ # -# parameters for the test driver "src/patch/test_patch_system.cc" +# ***** parameters for the test driver "src/patch/test_patch_system.cc" ***** # # By default this test driver isn't compiled into the cactus executable, # and these parameters are ignored. To compile this test driver into diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc index 47bfd64..26391cb 100644 --- a/src/driver/Newton.cc +++ b/src/driver/Newton.cc @@ -64,7 +64,7 @@ bool Newton_solve(patch_system& ps, const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, bool initial_find_flag, - struct IO_info& IO_info, + const struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info) { const int max_Newton_iterations diff --git a/src/driver/driver.hh b/src/driver/driver.hh index 1e128cf..1811ce6 100644 --- a/src/driver/driver.hh +++ b/src/driver/driver.hh @@ -1,6 +1,11 @@ // driver.hh -- header file for driver code // $Header$ +// +// prerequisites: +// <stdio.h> +// + //****************************************************************************** // @@ -40,17 +45,6 @@ enum initial_guess_method initial_guess__coord_ellipsoid // no comma }; -// -// this enum holds the decoded horizon_file_format parameter, i.e. -// it specifies what format of input/output file(s) we should use -// for h and H (and other angular grid functions) -// -enum horizon_file_format - { - horizon_file_format__ASCII_gnuplot, - horizon_file_format__HDF5 // no comma - }; - //****************************************************************************** // @@ -103,11 +97,24 @@ struct solver_info // // This struct holds info for computing black hole diagnostics. // -struct diagnostics_info +struct BH_diagnostics_info { enum patch::integration_method surface_integral_method; }; +//****************************************************************************** + +// +// this enum holds the decoded horizon_file_format parameter, i.e. +// it specifies what format of input/output file(s) we should use +// for h and H (and other angular grid functions) +// +enum horizon_file_format + { + horizon_file_format__ASCII_gnuplot, + horizon_file_format__HDF5 // no comma + }; + // // This struct holds info for I/O // @@ -125,9 +132,14 @@ struct IO_info const char* Delta_h_base_file_name; const char* Jacobian_base_file_name; + bool output_BH_diagnostics; + const char* BH_diagnostics_base_file_name; + const char* BH_diagnostics_file_name_extension; + // this is used to choose file names int time_iteration; // the Cactus time interation number // (cctk_iteration) + fp time; // the Cactus time coordinate (cctk_time) }; // @@ -149,6 +161,18 @@ struct verbose_info //****************************************************************************** // +// (A single copy of) this struct holds all of our black hole diagnostics +// for a single apparent horizon. These diagnostics are only meaningful +// if the apparent horizon has indeed been found. +// +struct BH_diagnostics + { + fp centroid_x, centroid_y, centroid_z; + fp area; + fp m_irreducible; + }; + +// // (A single copy of) this struct holds all of our information about // a single apparent horizon. // @@ -160,9 +184,8 @@ struct AH_info struct initial_guess_info initial_guess_info; bool AH_found; - fp centroid_x, centroid_y, centroid_z; - fp area; - fp m_irreducible; + struct BH_diagnostics BH_diagnostics; + FILE *BH_diagnostics_fileptr; }; // @@ -181,7 +204,7 @@ struct state struct cactus_grid_info cgi; struct geometry_info gi; - struct diagnostics_info diagnostics_info; + struct BH_diagnostics_info BH_diagnostics_info; int N_horizons; @@ -210,7 +233,7 @@ extern "C" // initial_guess.cc void setup_initial_guess(patch_system& ps, const struct initial_guess_info& igi, - struct IO_info& IO_info, + const struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info); enum initial_guess_method decode_initial_guess_method(const char initial_guess_method_string[]); @@ -224,19 +247,24 @@ bool Newton_solve(patch_system& ps, const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, bool initial_find_flag, - struct IO_info& IO_info, + const struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info); // io.cc void input_gridfn(patch_system& ps, int unknown_gfn, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration = 0); void output_gridfn(patch_system& ps, int unknown_gfn, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration = 0); void print_Jacobians(const patch_system& ps, const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration = 0); +FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info, + int hn, int N_horizons); +void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics, + const struct IO_info& IO_info, + int hn, FILE* fileptr); enum horizon_file_format decode_horizon_file_format(const char horizon_file_format_string[]); diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index 32b30b5..bf21f0d 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -13,7 +13,7 @@ /// do_test_Jacobian /// do_find_horizon /// -/// BH_diagnostics - compute BH diagnostics for a horizon +/// compute_BH_diagnostics - compute BH diagnostics for a horizon /// surface_integral - compute surface integral of a gridfn over the horizon // @@ -93,9 +93,11 @@ bool do_find_horizon patch_system& ps, Jacobian* Jac_ptr, int hn, int N_horizons); -void BH_diagnostics(const struct diagnostics_info& diagnostics_info, - const struct verbose_info& verbose_info, - struct AH_info& AH_info); +void compute_BH_diagnostics + (const patch_system& ps, + const struct BH_diagnostics_info& BH_diagnostics_info, + const struct verbose_info& verbose_info, + struct BH_diagnostics& BH_diagnostics); fp surface_integral(const patch_system& ps, int src_gfn, enum patch::integration_method method); } @@ -131,6 +133,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // get the Cactus time step and decide if we want to output [hH] now IO_info.time_iteration = cctk_iteration; +IO_info.time = cctk_time; const bool output_h = (IO_info.how_often_to_output_h > 0) && ((IO_info.time_iteration % IO_info.how_often_to_output_h) == 0); @@ -200,41 +203,56 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi); ps, AH_info.Jac_ptr, hn, state.N_horizons); - // compute BH diagnostics? +#ifdef NOT_YET + // store found flag in Cactus array variable + AH_found[hn] = AH_info.AH_found; +#endif + if (AH_info.AH_found) then { - BH_diagnostics(state.diagnostics_info, - verbose_info, - AH_info); + // compute BH diagnostics + compute_BH_diagnostics(ps, + state.BH_diagnostics_info, + verbose_info, + AH_info.BH_diagnostics); + const struct BH_diagnostics& BH_diagnostics + = AH_info.BH_diagnostics; + + // print BH diagnostics? if (verbose_info.print_physics_details) then { CCTK_VInfo(CCTK_THORNSTRING, "AH found: A=%.10g at (%f,%f,%f)", - double(AH_info.area), - double(AH_info.centroid_x), - double(AH_info.centroid_y), - double(AH_info.centroid_z)); + double(BH_diagnostics.area), + double(BH_diagnostics.centroid_x), + double(BH_diagnostics.centroid_y), + double(BH_diagnostics.centroid_z)); CCTK_VInfo(CCTK_THORNSTRING, "m_irreducible=%.10g", - double(AH_info.m_irreducible)); + double(BH_diagnostics.m_irreducible)); } + +#ifdef NOT_YET + // store BH diagnostics in Cactus array variables + centroid_x[hn] = BH_diagnostics.centroid_x; + centroid_y[hn] = BH_diagnostics.centroid_y; + centroid_z[hn] = BH_diagnostics.centroid_z; + area[hn] = BH_diagnostics.area; + m_irreducible[hn] = BH_diagnostics.m_irreducible; +#endif + + // output BH diagnostics? + if (IO_info.output_BH_diagnostics) + then output_BH_diagnostics_fn + (BH_diagnostics, + IO_info, hn, + AH_info.BH_diagnostics_fileptr); } else { if (verbose_info.print_physics_details) then CCTK_VInfo(CCTK_THORNSTRING, "no apparent horizon found"); } - - // store results in Cactus array variables - AH_found[hn] = AH_info.AH_found; - if (AH_info.AH_found) - then { - centroid_x[hn] = AH_info.centroid_x; - centroid_y[hn] = AH_info.centroid_y; - centroid_z[hn] = AH_info.centroid_z; - area[hn] = AH_info.area; - m_irreducible[hn] = AH_info.m_irreducible; - } break; default: CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, @@ -544,23 +562,23 @@ return true; // *** NORMAL RETURN *** // global_[xyz] # nominal // namespace { -void BH_diagnostics(const struct diagnostics_info& diagnostics_info, - const struct verbose_info& verbose_info, - struct AH_info& AH_info) +void compute_BH_diagnostics + (const patch_system& ps, + const struct BH_diagnostics_info& BH_diagnostics_info, + const struct verbose_info& verbose_info, + struct BH_diagnostics& BH_diagnostics) { -const patch_system& ps = * AH_info.ps_ptr; - // // compute raw surface integrals // fp integral_one = surface_integral(ps, gfns::gfn__one, - diagnostics_info.surface_integral_method); + BH_diagnostics_info.surface_integral_method); fp integral_x = surface_integral(ps, gfns::gfn__global_x, - diagnostics_info.surface_integral_method); + BH_diagnostics_info.surface_integral_method); fp integral_y = surface_integral(ps, gfns::gfn__global_y, - diagnostics_info.surface_integral_method); + BH_diagnostics_info.surface_integral_method); fp integral_z = surface_integral(ps, gfns::gfn__global_z, - diagnostics_info.surface_integral_method); + BH_diagnostics_info.surface_integral_method); // // adjust integrals to take into account patch system not covering full sphere @@ -596,18 +614,18 @@ case patch_system::plus_xyz_octant_patch_system: default: CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" BH_diagnostics(): unknown ps.type()=(int)%d!\n" -" (this should never happen!)" +" compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n" +" (this should never happen!)" , int(ps.type())); /*NOTREACHED*/ } -AH_info.centroid_x = integral_x / integral_one; -AH_info.centroid_y = integral_y / integral_one; -AH_info.centroid_z = integral_z / integral_one; +BH_diagnostics.centroid_x = integral_x / integral_one; +BH_diagnostics.centroid_y = integral_y / integral_one; +BH_diagnostics.centroid_z = integral_z / integral_one; -AH_info.area = integral_one; -AH_info.m_irreducible = sqrt(AH_info.area / (16.0*PI)); +BH_diagnostics.area = integral_one; +BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI)); } } diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc index 9c92ca9..c72590d 100644 --- a/src/driver/initial_guess.cc +++ b/src/driver/initial_guess.cc @@ -70,7 +70,7 @@ void setup_coord_ellipsoid(patch_system& ps, // void setup_initial_guess(patch_system& ps, const struct initial_guess_info& igi, - struct IO_info& IO_info, + const struct IO_info& IO_info, int hn, const struct verbose_info& verbose_info) { if (verbose_info.print_algorithm_highlights) diff --git a/src/driver/io.cc b/src/driver/io.cc index 5b788f7..5b4aaaa 100644 --- a/src/driver/io.cc +++ b/src/driver/io.cc @@ -2,11 +2,14 @@ // $Header$ // // <<<prototypes for functions local to this file>>> -// input_gridfn -// output_gridfn -// output_Jacobians +// input_gridfn - read an angular grid function from a data file +// output_gridfn - write an angular grid function to a data file +// output_Jacobians - write a Jacobian matrix or matrices to a data file // decode_horizon_file_format - decode the horizon_file_format parameter // +// setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file +// output_BH_diagnostics_fn - append BH diagnostics to file +// /// io_file_name // @@ -52,7 +55,8 @@ using jtutil::error_exit; // namespace { -const char* io_file_name(struct IO_info& IO_info, const char base_file_name[], +const char* io_file_name(const struct IO_info& IO_info, + const char base_file_name[], int hn, int AHF_iteration = 0); } @@ -65,7 +69,7 @@ const char* io_file_name(struct IO_info& IO_info, const char base_file_name[], // present in the data file. // void input_gridfn(patch_system& ps, int unknown_gfn, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration /* = 0 */) { const char* file_name = io_file_name(IO_info, base_file_name, @@ -98,7 +102,7 @@ default: " input_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" , - int(IO_info.horizon_file_format)); /*NOTREACHED*/ + int(IO_info.horizon_file_format)); /*NOTREACHED*/ } } @@ -113,7 +117,7 @@ default: // FIXME: if the gridfn is not h, we assume that it's nominal-grid. // void output_gridfn(patch_system& ps, int unknown_gfn, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration /* = 0 */) { const char* file_name = io_file_name(IO_info, base_file_name, @@ -167,7 +171,7 @@ default: " output_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" , - int(IO_info.horizon_file_format)); /*NOTREACHED*/ + int(IO_info.horizon_file_format)); /*NOTREACHED*/ } } @@ -185,7 +189,7 @@ default: // void print_Jacobians(const patch_system& ps, const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr, - struct IO_info& IO_info, const char base_file_name[], + const struct IO_info& IO_info, const char base_file_name[], int hn, bool print_msg_flag, int AHF_iteration /* = 0 */) { if (Jac_NP == NULL) @@ -309,7 +313,86 @@ else if (STRING_EQUAL(horizon_file_format_string, "HDF5")) then return horizon_file_format__HDF5; else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!", - horizon_file_format_string); /*NOTREACHED*/ + horizon_file_format_string); /*NOTREACHED*/ +} + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function creates a BH-diagnostics output file, writes a suitable +// header comment, and returns a stdio file pointer which can be used to +// append data to the file. +// +FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info, + int hn, int N_horizons) +{ +const int N_file_name_buffer = 200; +char file_name_buffer[N_file_name_buffer]; + +snprintf(file_name_buffer, N_file_name_buffer, + "%s.ah%d.%s", + IO_info.BH_diagnostics_base_file_name, + 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, +"\n" +" setup_BH_diagnostics_output_file():\n" +" can't open BH-diagnostics output file\n" +" \"%s\"!" + , + file_name_buffer); /*NOTREACHED*/ + +fprintf(fileptr, "# apparent horizon %d of %d\n", hn, N_horizons); +fprintf(fileptr, "#\n"); +fprintf(fileptr, "# column 1 = cctk_iteration\n"); +fprintf(fileptr, "# column 2 = cctk_time\n"); +fprintf(fileptr, "# column 3 = centroid_x\n"); +fprintf(fileptr, "# column 4 = centroid_y\n"); +fprintf(fileptr, "# column 5 = centroid_z\n"); +fprintf(fileptr, "# column 6 = area\n"); +fprintf(fileptr, "# column 7 = m_irreducible\n"); +fflush(fileptr); + +return fileptr; +} + +//****************************************************************************** + +// +// This function appends a BH-diagnostics line to a data file. It +// attempts to ensure that the newly-written line is flushed to disk, +// so the output file can be examined while the Cactus run is still in +// progress. +// +// (The "_fn" in the function name is to distinguish it from the Cactus +// parameter output_BH_diagnostics .) +// +// Arguments: +// BH_diagnostics = The BH diagnostics to be written +// fileptr = The stdio file pointer to append to +// +void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics, + const struct IO_info& IO_info, + int hn, FILE* fileptr) +{ +assert(fileptr != NULL); + +// cctk_iteration +// == cctk_time +// == ==== centroid_[xyz] +// == ==== ========== area m_irreducible +// == ==== ========== ===== ===== +fprintf(fileptr, "%d\t%.3f\t%f\t%f\t%f\t%.10g\t%.10g\n", + IO_info.time_iteration, double(IO_info.time), + BH_diagnostics.centroid_x, + BH_diagnostics.centroid_y, + BH_diagnostics.centroid_z, + BH_diagnostics.area, BH_diagnostics.m_irreducible); + +fflush(fileptr); } //****************************************************************************** @@ -317,13 +400,11 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, //****************************************************************************** // -// This function encapsulates our file-naming conventions. +// This function encapsulates our file-naming conventions for "horizon" +// output files (those used for h, H, and other angular grid functions). // // Arguments: -// IO_info = note this is *not* const since we may keep things like -// file positions etc here // base_file_name[] = from the parameter file -// time_iteration = the Cactus time iteration number (cctk_iteration) // hn = the horizon number // AHF_iteration = the apparent horizon finder's internal iteration // number (>= 1) if this is an intermediate iterate, @@ -335,7 +416,8 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // result points into a private static buffer; the usual caveats apply. // namespace { -const char* io_file_name(struct IO_info& IO_info, const char base_file_name[], +const char* io_file_name(const struct IO_info& IO_info, + const char base_file_name[], int hn, int AHF_iteration /* = 0 */) { const int N_file_name_buffer = 200; @@ -356,7 +438,7 @@ default: " io_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n" " (this should never happen!)" , - int(IO_info.horizon_file_format)); /*NOTREACHED*/ + int(IO_info.horizon_file_format)); /*NOTREACHED*/ } if (AHF_iteration == 0) diff --git a/src/driver/setup.cc b/src/driver/setup.cc index d5d51b3..125b716 100644 --- a/src/driver/setup.cc +++ b/src/driver/setup.cc @@ -124,7 +124,11 @@ 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.output_BH_diagnostics = (output_BH_diagnostics != 0); +IO_info.BH_diagnostics_base_file_name = BH_diagnostics_base_file_name; +IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension; IO_info.time_iteration = 0; +IO_info.time = 0.0; struct Jacobian_info& Jac_info = state.Jac_info; Jac_info.Jacobian_method = decode_Jacobian_method(Jacobian_method); @@ -213,7 +217,7 @@ gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0); // // other misc setup // -state.diagnostics_info.surface_integral_method +state.BH_diagnostics_info.surface_integral_method = patch::decode_integration_method(surface_integral_method); @@ -335,11 +339,17 @@ state.AH_info_ptrs.push_back(NULL); // initialize the BH diagnostics AH_ptr->AH_found = false; - AH_ptr->area = 0.0; - AH_ptr->m_irreducible = 0.0; - AH_ptr->centroid_x = 0.0; - AH_ptr->centroid_y = 0.0; - AH_ptr->centroid_z = 0.0; + AH_ptr->BH_diagnostics.centroid_x = 0.0; + AH_ptr->BH_diagnostics.centroid_y = 0.0; + AH_ptr->BH_diagnostics.centroid_z = 0.0; + AH_ptr->BH_diagnostics.area = 0.0; + AH_ptr->BH_diagnostics.m_irreducible = 0.0; + + // set up the BH-diagnostics output files? + AH_ptr->BH_diagnostics_fileptr + = output_BH_diagnostics + ? setup_BH_diagnostics_output_file(IO_info, hn, N_horizons) + : NULL; state.AH_info_ptrs.push_back(AH_ptr); } |