diff options
Diffstat (limited to 'src/driver/find_horizons.cc')
| -rw-r--r-- | src/driver/find_horizons.cc | 84 |
1 files changed, 54 insertions, 30 deletions
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc index 3f3279e..0739a3c 100644 --- a/src/driver/find_horizons.cc +++ b/src/driver/find_horizons.cc @@ -4,6 +4,9 @@ // <<<access to persistent data>>> // <<<prototypes for functions local to this file>>> // AHFinderDirect_find_horizons - top-level driver to find apparent horizons +// find_horizon - find a horizon +// BH_diagnostics - compute BH diagnostics for a horizon +// surface_integral - compute surface integral of a gridfn over the horizon // #include <stdio.h> @@ -24,14 +27,14 @@ #include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../util/coords.hh" -#include "../util/grid.hh" -#include "../util/fd_grid.hh" -#include "../util/patch.hh" -#include "../util/patch_edge.hh" -#include "../util/patch_interp.hh" -#include "../util/ghost_zone.hh" -#include "../util/patch_system.hh" +#include "../patch/coords.hh" +#include "../patch/grid.hh" +#include "../patch/fd_grid.hh" +#include "../patch/patch.hh" +#include "../patch/patch_edge.hh" +#include "../patch/patch_interp.hh" +#include "../patch/ghost_zone.hh" +#include "../patch/patch_system.hh" #include "../elliptic/Jacobian.hh" @@ -64,6 +67,8 @@ bool find_horizon(enum method method, void BH_diagnostics(enum patch::integration_method surface_integral_method, const struct verbose_info& verbose_info, struct AH_info& AH_info); +fp surface_integral(const patch_system& ps, int src_gfn, + enum patch::integration_method method); }; //****************************************************************************** @@ -105,7 +110,7 @@ if (state.timer_handle >= 0) if (verbose_info.print_physics_details) then CCTK_VInfo(CCTK_THORNSTRING, - "AH found: A=%g m=%g at (%g,%g,%g)", + "AH found: A=%.10g m=%.10g at (%g,%g,%g)", double(AH_info.area), double(AH_info.mass), double(AH_info.centroid_x), double(AH_info.centroid_y), @@ -172,9 +177,10 @@ case method__horizon_function: then CCTK_VInfo(CCTK_THORNSTRING, " H(h) rms-norm %.2e, infinity-norm %.2e", H_norms.rms_norm(), H_norms.infinity_norm()); - output_gridfn(ps, gfns::gfn__H, - IO_info, IO_info.H_base_file_name, - hn, true); + if (IO_info.output_H) + then output_gridfn(ps, gfns::gfn__H, + IO_info, IO_info.H_base_file_name, + hn, true); return true; // *** NORMAL RETURN *** } @@ -249,12 +255,14 @@ case method__Newton_solve: if (! status) then return false; // *** ERROR RETURN *** - output_gridfn(ps, gfns::gfn__h, - IO_info, IO_info.h_base_file_name, - hn, verbose_info.print_algorithm_details); - output_gridfn(ps, gfns::gfn__H, - IO_info, IO_info.H_base_file_name, - hn, verbose_info.print_algorithm_details); + if (IO_info.output_h) + then output_gridfn(ps, gfns::gfn__h, + IO_info, IO_info.h_base_file_name, + hn, verbose_info.print_algorithm_details); + if (IO_info.output_H) + then output_gridfn(ps, gfns::gfn__H, + IO_info, IO_info.H_base_file_name, + hn, verbose_info.print_algorithm_details); return true; // *** NORMAL RETURN *** } @@ -292,18 +300,14 @@ const patch_system& ps = * AH_info.ps_ptr; // // compute raw surface integrals // -fp integral_one = ps.integrate_gridfn(gfns::gfn__one, - true, // area weighting - surface_integral_method); -fp integral_x = ps.integrate_gridfn(gfns::gfn__global_x, - true, // area weighting - surface_integral_method); -fp integral_y = ps.integrate_gridfn(gfns::gfn__global_y, - true, // area weighting - surface_integral_method); -fp integral_z = ps.integrate_gridfn(gfns::gfn__global_z, - true, // area weighting - surface_integral_method); +fp integral_one = surface_integral(ps, gfns::gfn__one, + surface_integral_method); +fp integral_x = surface_integral(ps, gfns::gfn__global_x, + surface_integral_method); +fp integral_y = surface_integral(ps, gfns::gfn__global_y, + surface_integral_method); +fp integral_z = surface_integral(ps, gfns::gfn__global_z, + surface_integral_method); // // adjust integrals to take into account patch system not covering full sphere @@ -353,3 +357,23 @@ AH_info.centroid_y = integral_y / integral_one; AH_info.centroid_z = integral_z / integral_one; } } + +//****************************************************************************** + +// +// This function computes the surface integral of a gridfn over the +// horizon. +// +namespace { +fp surface_integral(const patch_system& ps, int src_gfn, + enum patch::integration_method method) +{ +return ps.integrate_gridfn + (src_gfn, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method); +} + } |