diff options
Diffstat (limited to 'src/driver/BH_diagnostics.cc')
| -rw-r--r-- | src/driver/BH_diagnostics.cc | 202 |
1 files changed, 158 insertions, 44 deletions
diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc index a641dde..163ef49 100644 --- a/src/driver/BH_diagnostics.cc +++ b/src/driver/BH_diagnostics.cc @@ -7,7 +7,7 @@ // BH_diagnostics::copy_from_buffer - copy buffer to diagnostics // // BH_diagnostics::compute - compute BH diagnostics after an AH has been found -/// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere +// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere // // print - print a line or two summarizing the diagnostics // setup_output_file - create/open output file, write header describing fields @@ -18,6 +18,8 @@ #include <assert.h> #include <math.h> +#include <vector> + #include "util_Table.h" #include "cctk.h" #include "cctk_Arguments.h" @@ -57,17 +59,27 @@ namespace AHFinderDirect //****************************************************************************** // +// ***** access to persistent data ***** +// +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), - min_radius(0.0), max_radius(0.0), mean_radius(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), - m_irreducible(0.0) + area(0.0), irreducible_mass(0.0), areal_radius(0.0) // no comma { } //****************************************************************************** @@ -84,6 +96,13 @@ buffer[posn__centroid_x] = centroid_x; buffer[posn__centroid_y] = centroid_y; buffer[posn__centroid_z] = centroid_z; +buffer[posn__quadrupole_xx] = quadrupole_xx; +buffer[posn__quadrupole_xy] = quadrupole_xy; +buffer[posn__quadrupole_xz] = quadrupole_xz; +buffer[posn__quadrupole_yy] = quadrupole_yy; +buffer[posn__quadrupole_xz] = quadrupole_yz; +buffer[posn__quadrupole_zz] = quadrupole_zz; + buffer[posn__min_radius] = min_radius; buffer[posn__max_radius] = max_radius; buffer[posn__mean_radius] = mean_radius; @@ -98,8 +117,10 @@ 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__m_irreducible] = m_irreducible; + +buffer[posn__area] = area; +buffer[posn__irreducible_mass] = irreducible_mass; +buffer[posn__areal_radius] = areal_radius; } //****************************************************************************** @@ -113,6 +134,13 @@ centroid_x = buffer[posn__centroid_x]; centroid_y = buffer[posn__centroid_y]; centroid_z = buffer[posn__centroid_z]; +quadrupole_xx = buffer[posn__quadrupole_xx]; +quadrupole_xy = buffer[posn__quadrupole_xy]; +quadrupole_xz = buffer[posn__quadrupole_xz]; +quadrupole_yy = buffer[posn__quadrupole_yy]; +quadrupole_yz = buffer[posn__quadrupole_yz]; +quadrupole_zz = buffer[posn__quadrupole_zz]; + min_radius = buffer[posn__min_radius]; max_radius = buffer[posn__max_radius]; mean_radius = buffer[posn__mean_radius]; @@ -127,8 +155,10 @@ 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]; - m_irreducible = buffer[posn__m_irreducible]; + irreducible_mass = buffer[posn__irreducible_mass]; + areal_radius = buffer[posn__areal_radius]; } //****************************************************************************** @@ -137,12 +167,13 @@ circumference_yz = buffer[posn__circumference_yz]; // // Given that an apparent horizon has been found, this function computes -// various black hole diagnostics. +// the black hole diagnostics. // // Inputs (gridfns) -// h # ghosted -// one # nominal -// global_[xyz] # nominal +// h # ghosted +// one # nominal +// global_[xyz] # nominal +// global_{xx,xy,xz,yy,yz,zz} # nominal // // Bugs: // The computation is rather inefficient -- we make many passes over the @@ -165,20 +196,19 @@ max_radius = h_norms.max_abs_value(); // xyz bounding box of horizon // -// compute bounding box of nominal grid -// ... this is only the stored part of the horizon if there are symmetries +// compute bounding box of nominal grid (for stored part of the horizon only) 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(); @@ -233,6 +263,24 @@ const fp integral_y = surface_integral(ps, 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); +const fp integral_xy = surface_integral(ps, + 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); +const fp integral_yy = surface_integral(ps, + 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); +const fp integral_zz = surface_integral(ps, + gfns::gfn__global_zz, true, true, true, + BH_diagnostics_info.integral_method); // @@ -244,15 +292,32 @@ centroid_z = integral_z / integral_one; // -// area, mean radius, and mass +// quadrupoles (taken about centroid position) +// +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; + + +// +// mean radius of horizon // -area = integral_one; mean_radius = integral_h / integral_one; -m_irreducible = sqrt(area / (16.0*PI)); // -// circumferences +// surface area and quantities derived from it +// +area = integral_one; +irreducible_mass = sqrt(area / (16.0*PI)); +areal_radius = sqrt(area / ( 4.0*PI)); + + +// +// proper circumferences // circumference_xy = ps.circumference("xy", gfns::gfn__h, @@ -310,15 +375,16 @@ return ps.integrate_gridfn void BH_diagnostics::print(int N_horizons, int hn) const { +const fp irreducible_mass = 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 m_irreducible=%.10g", + "AH %d/%d: area=%.10g irreducible_mass=%.10g", hn, N_horizons, - double(area), double(m_irreducible)); + double(area), double(irreducible_mass)); } //****************************************************************************** @@ -375,19 +441,35 @@ fprintf(fileptr, "# column 5 = centroid_z\n"); fprintf(fileptr, "# column 6 = min radius\n"); fprintf(fileptr, "# column 7 = max radius\n"); fprintf(fileptr, "# column 8 = mean radius\n"); -fprintf(fileptr, "# column 9 = min x\n"); -fprintf(fileptr, "# column 10 = max x\n"); -fprintf(fileptr, "# column 11 = min y\n"); -fprintf(fileptr, "# column 12 = max y\n"); -fprintf(fileptr, "# column 13 = min z\n"); -fprintf(fileptr, "# column 14 = max z\n"); -fprintf(fileptr, "# column 15 = xy-plane circumference\n"); -fprintf(fileptr, "# column 16 = xz-plane circumference\n"); -fprintf(fileptr, "# column 17 = yz-plane circumference\n"); -fprintf(fileptr, "# column 18 = ratio of xz/xy-plane circumferences\n"); -fprintf(fileptr, "# column 19 = ratio of yz/xy-plane circumferences\n"); -fprintf(fileptr, "# column 20 = area\n"); -fprintf(fileptr, "# column 21 = m_irreducible\n"); +fprintf(fileptr, "# column 9 = quadrupole_xx\n"); +fprintf(fileptr, "# column 10 = quadrupole_xy\n"); +fprintf(fileptr, "# column 11 = quadrupole_xz\n"); +fprintf(fileptr, "# column 12 = quadrupole_yy\n"); +fprintf(fileptr, "# column 13 = quadrupole_yz\n"); +fprintf(fileptr, "# column 14 = quadrupole_zz\n"); +fprintf(fileptr, "# column 15 = min x\n"); +fprintf(fileptr, "# column 16 = max x\n"); +fprintf(fileptr, "# column 17 = min y\n"); +fprintf(fileptr, "# column 18 = max y\n"); +fprintf(fileptr, "# column 19 = min z\n"); +fprintf(fileptr, "# column 20 = max z\n"); +fprintf(fileptr, "# column 21 = xy-plane circumference\n"); +fprintf(fileptr, "# column 22 = xz-plane circumference\n"); +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 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"); fflush(fileptr); return fileptr; @@ -420,6 +502,14 @@ fprintf(fileptr, 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)); + +fprintf(fileptr, // {min,max}_x {min,max}_y {min,max}_z // ============== ============== ============== "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t", @@ -428,17 +518,41 @@ fprintf(fileptr, double(min_z), double(max_z)); fprintf(fileptr, - // {xy,xz,yz}-plane xz/xy yz/xy area m_irreducible - // circumferences circumference ====== ====== + // {xy,xz,yz}-plane xz/xy yz/xy + // circumferences circumference // ratios - // ====================== ============== ====== ====== - "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n", + // ====================== ============== + "%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t", double(circumference_xy), double(circumference_xz), double(circumference_yz), double(circumference_xz / circumference_xy), - double(circumference_yz / circumference_xy), - double(area), double(m_irreducible)); + double(circumference_yz / circumference_xy)); + +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); fflush(fileptr); } |