diff options
Diffstat (limited to 'src/driver/BH_diagnostics.cc')
| -rw-r--r-- | src/driver/BH_diagnostics.cc | 521 |
1 files changed, 415 insertions, 106 deletions
diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc index 125177c..dde7d5c 100644 --- a/src/driver/BH_diagnostics.cc +++ b/src/driver/BH_diagnostics.cc @@ -12,15 +12,22 @@ // print - print a line or two summarizing the diagnostics // setup_output_file - create/open output file, write header describing fields // output - write a (long) line of all the diagnostics +// store - copy the surface and the diagnostics into the SphericalSurface +// variables +// save - copy the surface and the diagnostics into the Cactus variables +// load - set the surface and the diagnostics from the Cactus variables // #include <stdio.h> #include <assert.h> #include <math.h> +#include <vector> + #include "util_Table.h" #include "cctk.h" #include "cctk_Arguments.h" +#include "cctk_Parameters.h" #include "config.h" #include "stdc.h" @@ -28,6 +35,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "../patch/coords.hh" #include "../patch/grid.hh" @@ -50,7 +58,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** @@ -67,17 +74,28 @@ extern struct state state; // This function initializes a struct BH_diagnostics to all zeros. // BH_diagnostics::BH_diagnostics() - : centroid_x(0.0), centroid_y(0.0), centroid_z(0.0), - quadrupole_xx(0.0), quadrupole_xy(0.0), quadrupole_xz(0.0), - quadrupole_yy(0.0), quadrupole_yz(0.0), - quadrupole_zz(0.0), - min_radius(0.0), max_radius(0.0), - mean_radius(0.0), + : origin_x(0.0), origin_y(0.0), origin_z(0.0), + centroid_x(0.0), centroid_y(0.0), centroid_z(0.0), + quadrupole_xx(0.0), quadrupole_xy(0.0), quadrupole_xz(0.0), + quadrupole_yy(0.0), quadrupole_yz(0.0), quadrupole_zz(0.0), + min_radius(0.0), max_radius(0.0), mean_radius(0.0), min_x(0.0), max_x(0.0), min_y(0.0), max_y(0.0), min_z(0.0), max_z(0.0), circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0), - area(0.0), irreducible_mass(0.0), areal_radius(0.0) // no comma + area(1.0), + expansion(0.0), + inner_expansion(0.0), + product_expansion(0.0), + mean_curvature(0.0), + area_gradient(0.0), + expansion_gradient(0.0), + inner_expansion_gradient(0.0), + product_expansion_gradient(0.0), + mean_curvature_gradient(0.0), + mean_curvature_minimum(0.0), + mean_curvature_maximum(0.0), + mean_curvature_integral(0.0) { } //****************************************************************************** @@ -90,6 +108,10 @@ BH_diagnostics::BH_diagnostics() void BH_diagnostics::copy_to_buffer(CCTK_REAL buffer[N_buffer]) const { +buffer[posn__origin_x] = this->origin_x; +buffer[posn__origin_y] = this->origin_y; +buffer[posn__origin_z] = this->origin_z; + buffer[posn__centroid_x] = centroid_x; buffer[posn__centroid_y] = centroid_y; buffer[posn__centroid_z] = centroid_z; @@ -112,13 +134,24 @@ buffer[posn__max_y] = max_y; buffer[posn__min_z] = min_z; buffer[posn__max_z] = max_z; -buffer[posn__circumference_xy] = circumference_xy; -buffer[posn__circumference_xz] = circumference_xz; -buffer[posn__circumference_yz] = circumference_yz; - -buffer[posn__area] = area; -buffer[posn__irreducible_mass] = irreducible_mass; -buffer[posn__areal_radius] = areal_radius; +buffer[posn__circumference_xy] = circumference_xy; +buffer[posn__circumference_xz] = circumference_xz; +buffer[posn__circumference_yz] = circumference_yz; +buffer[posn__area] = area; +buffer[posn__expansion] = expansion; +buffer[posn__inner_expansion] = inner_expansion; +buffer[posn__product_expansion] = product_expansion; +buffer[posn__mean_curvature] = mean_curvature; + +buffer[posn__area_gradient] = area_gradient; +buffer[posn__expansion_gradient] = expansion_gradient; +buffer[posn__inner_expansion_gradient] = inner_expansion_gradient; +buffer[posn__product_expansion_gradient] = product_expansion_gradient; +buffer[posn__mean_curvature_gradient] = mean_curvature_gradient; + +buffer[posn__mean_curvature_minimum] = mean_curvature_minimum; +buffer[posn__mean_curvature_maximum] = mean_curvature_maximum; +buffer[posn__mean_curvature_integral] = mean_curvature_integral; } //****************************************************************************** @@ -128,6 +161,10 @@ buffer[posn__areal_radius] = areal_radius; // void BH_diagnostics::copy_from_buffer(const CCTK_REAL buffer[N_buffer]) { +this->origin_x = buffer[posn__origin_x]; +this->origin_y = buffer[posn__origin_y]; +this->origin_z = buffer[posn__origin_z]; + centroid_x = buffer[posn__centroid_x]; centroid_y = buffer[posn__centroid_y]; centroid_z = buffer[posn__centroid_z]; @@ -150,13 +187,24 @@ max_y = buffer[posn__max_y]; min_z = buffer[posn__min_z]; max_z = buffer[posn__max_z]; -circumference_xy = buffer[posn__circumference_xy]; -circumference_xz = buffer[posn__circumference_xz]; -circumference_yz = buffer[posn__circumference_yz]; - - area = buffer[posn__area]; - irreducible_mass = buffer[posn__irreducible_mass]; - areal_radius = buffer[posn__areal_radius]; + circumference_xy = buffer[posn__circumference_xy]; + circumference_xz = buffer[posn__circumference_xz]; + circumference_yz = buffer[posn__circumference_yz]; + area = buffer[posn__area]; + expansion = buffer[posn__expansion]; + inner_expansion = buffer[posn__inner_expansion]; +product_expansion = buffer[posn__product_expansion]; + mean_curvature = buffer[posn__mean_curvature]; + + area_gradient = buffer[posn__area_gradient]; + expansion_gradient = buffer[posn__expansion_gradient]; + inner_expansion_gradient = buffer[posn__inner_expansion_gradient]; +product_expansion_gradient = buffer[posn__product_expansion_gradient]; + mean_curvature_gradient = buffer[posn__mean_curvature_gradient]; + + mean_curvature_minimum = buffer[posn__mean_curvature_minimum]; + mean_curvature_maximum = buffer[posn__mean_curvature_maximum]; + mean_curvature_integral = buffer[posn__mean_curvature_integral]; } //****************************************************************************** @@ -165,13 +213,12 @@ circumference_yz = buffer[posn__circumference_yz]; // // Given that an apparent horizon has been found, this function computes -// the black hole diagnostics. +// various black hole diagnostics. // // Inputs (gridfns) -// h # ghosted -// one # nominal -// global_[xyz] # nominal -// global_{xx,xy,xz,yy,yz,zz} # nominal +// h # ghosted +// one # nominal +// global_[xyz] # nominal // // Bugs: // The computation is rather inefficient -- we make many passes over the @@ -179,6 +226,16 @@ circumference_yz = buffer[posn__circumference_yz]; // void BH_diagnostics::compute (const patch_system& ps, + const fp the_area, + const fp mean_expansion, + const fp mean_inner_expansion, + const fp mean_product_expansion, + const fp mean_mean_curvature, + const fp the_area_gradient, + const fp mean_expansion_gradient, + const fp mean_inner_expansion_gradient, + const fp mean_product_expansion_gradient, + const fp mean_mean_curvature_gradient, const struct BH_diagnostics_info& BH_diagnostics_info) { // @@ -194,19 +251,20 @@ max_radius = h_norms.max_abs_value(); // xyz bounding box of horizon // -// compute bounding box of nominal grid (for stored part of the horizon only) +// compute bounding box of nominal grid +// ... this is only the stored part of the horizon if there are symmetries jtutil::norm<fp> x_norms; -jtutil::norm<fp> y_norms; -jtutil::norm<fp> z_norms; - ps.gridfn_norms(gfns::gfn__global_x, x_norms); -ps.gridfn_norms(gfns::gfn__global_y, y_norms); -ps.gridfn_norms(gfns::gfn__global_z, z_norms); - min_x = x_norms.min_value(); max_x = x_norms.max_value(); + +jtutil::norm<fp> y_norms; +ps.gridfn_norms(gfns::gfn__global_y, y_norms); min_y = y_norms.min_value(); max_y = y_norms.max_value(); + +jtutil::norm<fp> z_norms; +ps.gridfn_norms(gfns::gfn__global_z, z_norms); min_z = z_norms.min_value(); max_z = z_norms.max_value(); @@ -262,23 +320,31 @@ const fp integral_z = surface_integral(ps, gfns::gfn__global_z, false, true, true, BH_diagnostics_info.integral_method); const fp integral_xx = surface_integral(ps, - gfns::gfn__global_xx, true, true, true, - BH_diagnostics_info.integral_method); + gfns::gfn__global_xx, true, true, true, + BH_diagnostics_info.integral_method); const fp integral_xy = surface_integral(ps, - gfns::gfn__global_xy, true, false, false, - BH_diagnostics_info.integral_method); + gfns::gfn__global_xy, true, false, false, + BH_diagnostics_info.integral_method); const fp integral_xz = surface_integral(ps, - gfns::gfn__global_xz, false, true, false, - BH_diagnostics_info.integral_method); + gfns::gfn__global_xz, false, true, false, + BH_diagnostics_info.integral_method); const fp integral_yy = surface_integral(ps, - gfns::gfn__global_yy, true, true, true, - BH_diagnostics_info.integral_method); + gfns::gfn__global_yy, true, true, true, + BH_diagnostics_info.integral_method); const fp integral_yz = surface_integral(ps, - gfns::gfn__global_yz, false, false, true, - BH_diagnostics_info.integral_method); + gfns::gfn__global_yz, false, false, true, + BH_diagnostics_info.integral_method); const fp integral_zz = surface_integral(ps, - gfns::gfn__global_zz, true, true, true, - BH_diagnostics_info.integral_method); + gfns::gfn__global_zz, true, true, true, + BH_diagnostics_info.integral_method); + + +// +// originds +// +this->origin_x = ps.origin_x(); +this->origin_y = ps.origin_y(); +this->origin_z = ps.origin_z(); // @@ -290,32 +356,52 @@ centroid_z = integral_z / integral_one; // -// quadrupoles (taken about centroid position) +// quadrupoles // -quadrupole_xx = integral_xx / integral_one - centroid_x * centroid_x; -quadrupole_xy = integral_xy / integral_one - centroid_x * centroid_y; -quadrupole_xz = integral_xz / integral_one - centroid_x * centroid_z; -quadrupole_yy = integral_yy / integral_one - centroid_y * centroid_y; -quadrupole_yz = integral_yz / integral_one - centroid_y * centroid_z; -quadrupole_zz = integral_zz / integral_one - centroid_z * centroid_z; +quadrupole_xx = integral_xx / integral_one; +quadrupole_xy = integral_xy / integral_one; +quadrupole_xz = integral_xz / integral_one; +quadrupole_yy = integral_yy / integral_one; +quadrupole_yz = integral_yz / integral_one; +quadrupole_zz = integral_zz / integral_one; // -// mean radius of horizon +// area, mean radius, and mass // mean_radius = integral_h / integral_one; // -// surface area and quantities derived from it +// expansion +// +area = the_area; +expansion = mean_expansion; +inner_expansion = mean_inner_expansion; +product_expansion = mean_product_expansion; +mean_curvature = mean_mean_curvature; + +area_gradient = the_area_gradient; +expansion_gradient = mean_expansion_gradient; +inner_expansion_gradient = mean_inner_expansion_gradient; +product_expansion_gradient = mean_product_expansion_gradient; +mean_curvature_gradient = mean_mean_curvature_gradient; + +// +// minimum, maximum and the integral of the mean curvature // -area = integral_one; -irreducible_mass = sqrt(area / (16.0*PI)); -areal_radius = sqrt(area / ( 4.0*PI)); +jtutil::norm<fp> mean_curvature_norms; +ps.gridfn_norms(gfns::gfn__mean_curvature, mean_curvature_norms); +mean_curvature_minimum = mean_curvature_norms.min_value(); +mean_curvature_maximum = mean_curvature_norms.max_value(); +mean_curvature_integral = surface_integral(ps, + gfns::gfn__mean_curvature, + true, true, true, + BH_diagnostics_info.integral_method); // -// proper circumferences +// circumferences // circumference_xy = ps.circumference("xy", gfns::gfn__h, @@ -373,15 +459,16 @@ return ps.integrate_gridfn void BH_diagnostics::print(int N_horizons, int hn) const { +const fp m_irreducible = sqrt(area / (16*PI)); CCTK_VInfo(CCTK_THORNSTRING, "AH %d/%d: r=%g at (%f,%f,%f)", hn, N_horizons, double(mean_radius), double(centroid_x), double(centroid_y), double(centroid_z)); CCTK_VInfo(CCTK_THORNSTRING, - "AH %d/%d: area=%.10g irreducible_mass=%.10g", + "AH %d/%d: area=%.10g m_irreducible=%.10g", hn, N_horizons, - double(area), double(irreducible_mass)); + double(area), double(m_irreducible)); } //****************************************************************************** @@ -419,10 +506,13 @@ snprintf(file_name_buffer, IO_info::file_name_buffer_size, directory, IO_info.BH_diagnostics_base_file_name, hn, IO_info.BH_diagnostics_file_name_extension); -const char *const file_open_mode = IO_TruncateOutputFiles(state.cgi.GH) - ? "w" : "a"; +const char *openMode; +if (IO_TruncateOutputFiles(state.cgi.GH) == 1) + openMode = "w"; +else + openMode = "a"; -FILE *fileptr = fopen(file_name_buffer, file_open_mode); +FILE *fileptr = fopen(file_name_buffer, openMode); if (fileptr == NULL) then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" @@ -460,17 +550,20 @@ fprintf(fileptr, "# column 23 = yz-plane circumference\n"); fprintf(fileptr, "# column 24 = ratio of xz/xy-plane circumferences\n"); fprintf(fileptr, "# column 25 = ratio of yz/xy-plane circumferences\n"); fprintf(fileptr, "# column 26 = area\n"); -fprintf(fileptr, "# column 27 = irreducible mass\n"); +fprintf(fileptr, "# column 27 = m_irreducible\n"); fprintf(fileptr, "# column 28 = areal radius\n"); -fprintf(fileptr, "# column 29 = [not implemented yet] (outer) expansion Theta_(l)\n"); -fprintf(fileptr, "# column 30 = [not implemented yet] inner expansion Theta_(n)\n"); -fprintf(fileptr, "# column 31 = [not implemented yet] product of inner and outer expansions\n"); -fprintf(fileptr, "# column 32 = [not implemented yet] mean curvature\n"); -fprintf(fileptr, "# column 33 = [not implemented yet] d/d(coordinate radius) of area\n"); -fprintf(fileptr, "# column 34 = [not implemented yet] d/d(coordinate radius) of (outer) expansion Theta_(l)\n"); -fprintf(fileptr, "# column 35 = [not implemented yet] d/d(coordinate radius) of inner expansion Theta_(n)\n"); -fprintf(fileptr, "# column 36 = [not implemented yet] d/d(coordinate radius) of product of inner and outer expansions\n"); -fprintf(fileptr, "# column 37 = [not implemented yet] d/d(coordinate radius) of mean curvature\n"); +fprintf(fileptr, "# column 29 = expansion Theta_(l)\n"); +fprintf(fileptr, "# column 30 = inner expansion Theta_(n)\n"); +fprintf(fileptr, "# column 31 = product of the expansions\n"); +fprintf(fileptr, "# column 32 = mean curvature\n"); +fprintf(fileptr, "# column 33 = gradient of the areal radius\n"); +fprintf(fileptr, "# column 34 = gradient of the expansion Theta_(l)\n"); +fprintf(fileptr, "# column 35 = gradient of the inner expansion Theta_(n)\n"); +fprintf(fileptr, "# column 36 = gradient of the product of the expansions\n"); +fprintf(fileptr, "# column 37 = gradient of the mean curvature\n"); +fprintf(fileptr, "# column 38 = minimum of the mean curvature\n"); +fprintf(fileptr, "# column 39 = maximum of the mean curvature\n"); +fprintf(fileptr, "# column 40 = integral of the mean curvature\n"); fflush(fileptr); return fileptr; @@ -493,22 +586,25 @@ void BH_diagnostics::output(FILE*fileptr, const struct IO_info& IO_info) assert(fileptr != NULL); fprintf(fileptr, - // cctk_iteration min radius mean radius - // == cctk_time ====== max radius - // == == centroid_[xyz] ====== ====== - // == == ========== ====== ====== ====== - "%d\t%f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t", + // cctk_iteration min radius mean radius + // == cctk_time ====== max radius + // == ==== centroid_[xyz] ====== ====== + // == ==== ========== ====== ====== ====== + "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t", IO_info.time_iteration, double(IO_info.time), double(centroid_x), double(centroid_y), double(centroid_z), double(min_radius), double(max_radius), double(mean_radius)); fprintf(fileptr, // quadrupole_{xx,xy,xz,yy,yz,zz} - // ====== ====== ====== ====== ====== ====== + // ================================================ "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t", - double(quadrupole_xx), double(quadrupole_xy), double(quadrupole_xz), - double(quadrupole_yy), double(quadrupole_yz), - double(quadrupole_zz)); + double(quadrupole_xx - centroid_x*centroid_x), + double(quadrupole_xy - centroid_x*centroid_y), + double(quadrupole_xz - centroid_x*centroid_z), + double(quadrupole_yy - centroid_y*centroid_y), + double(quadrupole_yz - centroid_y*centroid_z), + double(quadrupole_zz - centroid_z*centroid_z)); fprintf(fileptr, // {min,max}_x {min,max}_y {min,max}_z @@ -530,35 +626,248 @@ fprintf(fileptr, double(circumference_xz / circumference_xy), double(circumference_yz / circumference_xy)); +const fp m_irreducible = sqrt(area / (16*PI));; +const fp areal_radius = sqrt(area / (4*PI)); +const fp areal_radius_gradient = sqrt(1 / (16*PI*area)) * area_gradient; fprintf(fileptr, - // area areal_radius - // ====== irreducible_mass - // ====== ====== ====== - "%#.10g\t%#.10g\t%#.10g\t", - double(area), double(irreducible_mass), double(areal_radius)); - -// these diagnostics aren't implemented yet :( -fprintf(fileptr, - // expansion product_expansion - // ====== inner_expansion - // ====== ====== ====== mean_curvature - // ====== ====== ====== ====== - "%#.10g\t%#.10g\t%#.10g\t%#.10g\t", - 0.0, 0.0, 0.0, 0.0); - -// these diagnostics aren't implemented yet :( -fprintf(fileptr, - // dr_area dr_inner_expansion - // ====== dr_expansion dr_product_expansion - // ====== ====== ====== ====== dr_mean_curvature - // ====== ====== ====== ====== ====== - "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n", - 0.0, 0.0, 0.0, 0.0, 0.0); + // area m_irre- areal expan- inner prod. mean areal expan- inner prod. mean mean mean mean + // ducible radius sion expan- of the curva- radius sion expan- of the curva- curva- curva- curva- + // sion expan- ture grad. grad. sion exp.s ture ture ture ture + // sions grad. grad. grad. min. max. integ. + // ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== + // + // ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== + "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n", + double(area), double(m_irreducible), double(areal_radius), + double(expansion), double(inner_expansion), double(product_expansion), double(mean_curvature), + double(areal_radius_gradient), + double(expansion_gradient), double(inner_expansion_gradient), double(product_expansion_gradient), double(mean_curvature_gradient), + double(mean_curvature_minimum), double(mean_curvature_maximum), double(mean_curvature_integral)); fflush(fileptr); } //****************************************************************************** + +// +// This function copies the BH-diagnostics into the SphericalSurface variables. +// +// Arguments: +// CCTK_ARGUMENTS +// +void BH_diagnostics::store(CCTK_ARGUMENTS, + const int horizon_number, const int surface_number) + const +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + assert (surface_number >= 0 && surface_number < nsurfaces); + + assert(state.AH_data_array[horizon_number] != NULL); + const struct AH_data& AH_data = *state.AH_data_array[horizon_number]; + + // tell the world if we don't look for this horizon any more + const int my_dont_find_after = dont_find_after_individual[horizon_number]; + const fp my_find_after_time = find_after_individual_time[horizon_number]; + const fp my_dont_find_after_time = dont_find_after_individual_time[horizon_number]; + if ((my_dont_find_after >= 0 and cctk_iteration > my_dont_find_after) or + (my_dont_find_after_time > my_find_after_time and cctk_time > my_dont_find_after_time)) + { + assert (! AH_data.search_flag); + sf_active[surface_number] = 0; + } + + // only try to copy AH info if we've found AHs at this time level + if (! AH_data.search_flag) { + sf_valid[surface_number] = 0; + return; + } + + // did we actually *find* this horizon? + if (! AH_data.found_flag) { + sf_valid[surface_number] = -1; + return; + } + + sf_active [surface_number] = 1; + sf_valid [surface_number] = 1; + sf_origin_x [surface_number] = this->origin_x; + sf_origin_y [surface_number] = this->origin_y; + sf_origin_z [surface_number] = this->origin_z; + sf_mean_radius [surface_number] = mean_radius; + sf_min_radius [surface_number] = min_radius; + sf_max_radius [surface_number] = max_radius; + sf_centroid_x [surface_number] = centroid_x; + sf_centroid_y [surface_number] = centroid_y; + sf_centroid_z [surface_number] = centroid_z; + sf_quadrupole_xx[surface_number] = quadrupole_xx - centroid_x*centroid_x; + sf_quadrupole_xy[surface_number] = quadrupole_xy - centroid_x*centroid_y; + sf_quadrupole_xz[surface_number] = quadrupole_xz - centroid_x*centroid_z; + sf_quadrupole_yy[surface_number] = quadrupole_yy - centroid_y*centroid_y; + sf_quadrupole_yz[surface_number] = quadrupole_yz - centroid_y*centroid_z; + sf_quadrupole_zz[surface_number] = quadrupole_zz - centroid_z*centroid_z; + sf_min_x [surface_number] = min_x; + sf_max_x [surface_number] = max_x; + sf_min_y [surface_number] = min_y; + sf_max_y [surface_number] = max_y; + sf_min_z [surface_number] = min_z; + sf_max_z [surface_number] = max_z; + sf_area [surface_number] = area; + + std::vector<CCTK_REAL> xx, yy, zz, rr; + xx.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]); + yy.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]); + zz.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]); + rr.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]); + size_t n; + n = 0; + for (int j=0; j<sf_nphi[surface_number]; ++j) { + for (int i=0; i<sf_ntheta[surface_number]; ++i) { + const CCTK_REAL theta + = sf_origin_theta[surface_number] + i * sf_delta_theta[surface_number]; + const CCTK_REAL phi + = sf_origin_phi[surface_number] + j * sf_delta_phi[surface_number]; + xx[n] = sf_origin_x[surface_number] + sin(theta) * cos(phi); + yy[n] = sf_origin_y[surface_number] + sin(theta) * sin(phi); + zz[n] = sf_origin_z[surface_number] + cos(theta); + ++n; + } + } + AHFinderDirect_radius_in_direction + (horizon_number, sf_ntheta[surface_number] * sf_nphi[surface_number], + &xx[0], &yy[0], &zz[0], &rr[0]); + n = 0; + for (int j=0; j<sf_nphi[surface_number]; ++j) { + for (int i=0; i<sf_ntheta[surface_number]; ++i) { + sf_radius[i + maxntheta * (j + maxnphi * surface_number)] = rr[n]; + ++n; + } + } +} + +//****************************************************************************** + +// +// This function copies the BH-diagnostics into the Cactus variables. +// +// Arguments: +// CCTK_ARGUMENTS +// +void BH_diagnostics::save(CCTK_ARGUMENTS, + const int horizon_number) + const +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + assert(state.AH_data_array[horizon_number] != NULL); + const struct verbose_info& verbose_info = state.verbose_info; + const struct AH_data& AH_data = *state.AH_data_array[horizon_number]; + patch_system& ps = *AH_data.ps_ptr; + + if (AH_data.search_flag && + AH_data.found_flag && + cctk_iteration > ah_centroid_iteration[horizon_number-1]) + { + ah_centroid_x_p[horizon_number-1] = ah_centroid_x[horizon_number-1]; + ah_centroid_y_p[horizon_number-1] = ah_centroid_y[horizon_number-1]; + ah_centroid_z_p[horizon_number-1] = ah_centroid_z[horizon_number-1]; + ah_centroid_t_p[horizon_number-1] = ah_centroid_t[horizon_number-1]; + ah_centroid_valid_p[horizon_number-1] = ah_centroid_valid[horizon_number-1]; + ah_centroid_iteration_p[horizon_number-1] = + ah_centroid_iteration[horizon_number-1]; + + ah_centroid_x[horizon_number-1] = centroid_x; + ah_centroid_y[horizon_number-1] = centroid_y; + ah_centroid_z[horizon_number-1] = centroid_z; + ah_centroid_t[horizon_number-1] = cctk_time; + ah_centroid_valid[horizon_number-1] = AH_data.found_flag; + ah_centroid_iteration[horizon_number-1] = cctk_iteration; + } + + ah_origin_x[horizon_number-1] = ps.origin_x(); + ah_origin_y[horizon_number-1] = ps.origin_y(); + ah_origin_z[horizon_number-1] = ps.origin_z(); + + for (int pn = 0; pn < ps.N_patches(); ++pn) { + assert (pn < 6); + patch& p = ps.ith_patch(pn); + for (int i = 0, irho = p.min_irho(); irho <= p.max_irho(); ++i, ++irho) { + assert (i <= max_N_zones_per_right_angle); + for (int j = 0, isigma = p.min_isigma(); isigma <= p.max_isigma(); ++j, ++isigma) { + assert (j <= max_N_zones_per_right_angle); + ah_radius[i + (max_N_zones_per_right_angle+1) * (j + (max_N_zones_per_right_angle+1) * (pn + 6 * (horizon_number-1)))] + = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + + ah_initial_find_flag [horizon_number-1] = AH_data.initial_find_flag; + ah_really_initial_find_flag[horizon_number-1] = AH_data.really_initial_find_flag; + ah_search_flag [horizon_number-1] = AH_data.search_flag; + ah_found_flag [horizon_number-1] = AH_data.found_flag; + if (verbose_info.print_algorithm_highlights) { + printf ("AHF BH_diagnostics::save[%d] initial_find_flag=%d\n", horizon_number, (int) AH_data.initial_find_flag); + printf ("AHF BH_diagnostics::save[%d] really_initial_find_flag=%d\n", horizon_number, (int) AH_data.really_initial_find_flag); + printf ("AHF BH_diagnostics::save[%d] search_flag=%d\n", horizon_number, (int) AH_data.search_flag); + printf ("AHF BH_diagnostics::save[%d] found_flag=%d\n", horizon_number, (int) AH_data.found_flag); + } +} + +//****************************************************************************** + +// +// This function copies the BH-diagnostics from the Cactus variables. +// +// Arguments: +// CCTK_ARGUMENTS +// +void BH_diagnostics::load(CCTK_ARGUMENTS, + const int horizon_number) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + assert(state.AH_data_array[horizon_number] != NULL); + const struct verbose_info& verbose_info = state.verbose_info; + struct AH_data& AH_data = *state.AH_data_array[horizon_number]; + patch_system& ps = *AH_data.ps_ptr; + + ps.origin_x(ah_origin_x[horizon_number-1]); + ps.origin_y(ah_origin_y[horizon_number-1]); + ps.origin_z(ah_origin_z[horizon_number-1]); + + for (int pn = 0; pn < ps.N_patches(); ++pn) { + assert (pn < 6); + patch& p = ps.ith_patch(pn); + for (int i = 0, irho = p.min_irho(); irho <= p.max_irho(); ++i, ++irho) { + assert (i <= max_N_zones_per_right_angle); + for (int j = 0, isigma = p.min_isigma(); isigma <= p.max_isigma(); ++j, ++isigma) { + assert (j <= max_N_zones_per_right_angle); + p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = + ah_radius[i + (max_N_zones_per_right_angle+1) * (j + (max_N_zones_per_right_angle+1) * (pn + 6 * (horizon_number-1)))]; + } + } + } + + // recover the full ghosted-grid horizon shape + // (we only save and load the nominal-grid shape) + ps.synchronize(); + + AH_data.initial_find_flag = ah_initial_find_flag [horizon_number-1]; + AH_data.really_initial_find_flag = ah_really_initial_find_flag[horizon_number-1]; + AH_data.search_flag = ah_search_flag [horizon_number-1]; + AH_data.found_flag = ah_found_flag [horizon_number-1]; + if (verbose_info.print_algorithm_highlights) { + printf ("AHF BH_diagnostics::load[%d] initial_find_flag=%d\n", horizon_number, (int) AH_data.initial_find_flag); + printf ("AHF BH_diagnostics::load[%d] really_initial_find_flag=%d\n", horizon_number, (int) AH_data.really_initial_find_flag); + printf ("AHF BH_diagnostics::load[%d] search_flag=%d\n", horizon_number, (int) AH_data.search_flag); + printf ("AHF BH_diagnostics::load[%d] found_flag=%d\n", horizon_number, (int) AH_data.found_flag); + } +} + +//****************************************************************************** //****************************************************************************** //****************************************************************************** |