path: root/src/driver/find_horizons.cc

// find_horizons.cc -- top level driver for finding apparent horizons
// $Header$
//
// <<<access to persistent data>>>
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_find_horizons - top-level driver to find apparent horizons
///
/// setup_Cactus_gridfn_data_ptrs - get all data pointers given variable indices
/// Cactus_gridfn_data_ptr - get a single data pointer from a variable index
///
/// find_horizon - find a horizon
/// do_evaluate_expansion
/// do_test_expansion_Jacobian
/// do_find_horizon
///
/// compute_BH_diagnostics - compute BH diagnostics for a horizon
/// surface_integral - compute surface integral of a gridfn over the horizon
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#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"
#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"

#include "../gr/gfns.hh"
#include "../gr/gr.hh"

#include "driver.hh"

//******************************************************************************

//
// ***** access to persistent data *****
//
extern struct state state;

//******************************************************************************

//
// ***** prototypes for functions local to this file
//
namespace {
void setup_Cactus_gridfn_data_ptrs(const cGH *GH, struct cactus_grid_info& cgi);
const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
					const char gridfn_name[],
					bool check_for_NULL = true);

bool do_evaluate_expansion
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info, bool output_h, bool output_Theta,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps,
	 bool active_flag, int hn, int N_horizons);
bool do_test_expansion_Jacobian
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info,
	 bool test_all_Jacobian_methods,
	 struct Jacobian_info& Jac_info,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps, Jacobian_matrix& Jac,
	 bool active_flag, int hn, int N_horizons);
bool do_find_horizon
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info, bool output_h, bool output_Theta,
	 const struct solver_info& solver_info, bool initial_find_flag,
	 const struct Jacobian_info& Jac_info,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps, Jacobian_matrix& Jac,
	 bool active_flag, int hn, int N_horizons);

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);
	  }

//******************************************************************************

//
// This function is called by the Cactus scheduler to find the apparent
// horizon or horizons in the current slice.
//
extern "C"
  void AHFinderDirect_find_horizons(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS

const struct verbose_info& verbose_info = state.verbose_info;
      struct IO_info&      IO_info      = state.IO_info;
const struct solver_info&  solver_info  = state.solver_info;

if (state.timer_handle >= 0)
   then CCTK_TimerResetI(state.timer_handle);

// what are the semantics of the Cactus gxx variables? (these may
// change from one call to another, so we have to re-check each time)
if      (CCTK_Equals(metric_type, "physical"))
   then state.cgi.Cactus_conformal_metric = false;
else if (CCTK_Equals(metric_type, "static conformal"))
   then state.cgi.Cactus_conformal_metric = (conformal_state > 0);
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
		   metric_type);				/*NOTREACHED*/

// only processor #0 is "active"
const bool active_flag = (state.my_proc == 0);

// 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);
const bool output_Theta
   = (IO_info.how_often_to_output_Theta > 0)
     && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0);

//
// if we're using them,
//    we need to re-fetch the Cactus data pointers at least each time step,
//    because they change each time Cactus rotates the time levels
//
if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
   then setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);

	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
	{
	struct AH_info& AH_info = state.AH_info_array[hn];
	patch_system& ps = *AH_info.ps_ptr;

	//
	// If this is our first attempt to find this horizon, or
	//    if we've tried to find it before but we failed on our
	//    immediately previous attempt, then we need to (re)set
	//    the initial guess.
	// Otherwise (i.e. if we've tried to find this horizon before,
	//    and we succeeded on our immediately previous attempt),
	//    then we can just reuse the previous horizon position as
	//    our initial guess for this time around.
	//
	const bool initial_find_flag = ! AH_info.AH_found;
	if (initial_find_flag)
	   then {
		setup_initial_guess(ps,
				    AH_info.initial_guess_info,
				    IO_info,
				    hn, state.N_horizons, verbose_info);
		if (active_flag && IO_info.output_initial_guess)
		   then output_gridfn(ps, gfns::gfn__h,
				      IO_info, IO_info.h_base_file_name,
				      hn, verbose_info
					  .print_algorithm_highlights);
		}

	// if we're going to need a Jacobian
	//    and it doesn't exist, create it.
	if ( (state.method != method__evaluate_expansion)
	     && (AH_info.Jac_ptr == NULL) )
	   then AH_info.Jac_ptr
		   = create_Jacobian_matrix(ps,
					    state.cgi, state.gi,
					    state.Jac_info,
					    verbose_info.print_algorithm_details);

	AH_info.AH_found = false;
	switch	(state.method)
		{
	case method__evaluate_expansion:
		do_evaluate_expansion(verbose_info, state.timer_handle,
				      IO_info, output_h, output_Theta,
				      state.cgi, state.gi,
				      ps,
				      active_flag, hn, state.N_horizons);
		break;

	case method__test_expansion_Jacobian:
		do_test_expansion_Jacobian(verbose_info, state.timer_handle,
					   IO_info,
					   test_all_Jacobian_methods,
					   state.Jac_info, state.cgi, state.gi,
					   ps, *AH_info.Jac_ptr,
					   active_flag, hn, state.N_horizons);
		break;

	case method__find_horizon:
		AH_info.AH_found
		   = do_find_horizon(verbose_info, state.timer_handle,
				     IO_info, output_h, output_Theta,
				     solver_info, initial_find_flag,
				     state.Jac_info, state.cgi, state.gi,
				     ps, *AH_info.Jac_ptr,
				     active_flag, hn, state.N_horizons);

		// store found flag in Cactus array variable
		AH_found[hn] = AH_info.AH_found;

		if (AH_info.AH_found)
		   then {
// 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 %d/%d: r=%g at (%f,%f,%f)",
		   hn, state.N_horizons,
		   double(BH_diagnostics.mean_radius),
		   double(BH_diagnostics.centroid_x),
		   double(BH_diagnostics.centroid_y),
		   double(BH_diagnostics.centroid_z));
	CCTK_VInfo(CCTK_THORNSTRING,
		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
		   hn, state.N_horizons,
		   double(BH_diagnostics.area),
		   double(BH_diagnostics.m_irreducible));
	}

// 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;

// output BH diagnostics?
if (active_flag && IO_info.output_BH_diagnostics)
   then {
	if (AH_info.BH_diagnostics_fileptr == NULL)
	   then AH_info.BH_diagnostics_fileptr
		   = setup_BH_diagnostics_output_file(IO_info,
						      hn, state.N_horizons);
	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");
			}
		break;

	default:
		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
			   "\n"
			   "   find_horizons(): unknown method=(int)%d!\n"
			   "                    (this should never happen!)"
			   ,
			   int(state.method));			/*NOTREACHED*/
		}
	}

if (state.timer_handle >= 0)
   then {
	printf("timer stats for computation:\n");
	CCTK_TimerPrintDataI(state.timer_handle, -1);
	}
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function sets up the geometry data pointers in a
//  struct cactus_grid_info .
//

namespace {
void setup_Cactus_gridfn_data_ptrs(const cGH *GH, struct cactus_grid_info& cgi)
{
cgi.g_dd_11_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_11_varindex, "g_11");
cgi.g_dd_12_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_12_varindex, "g_12");
cgi.g_dd_13_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_13_varindex, "g_13");
cgi.g_dd_22_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_22_varindex, "g_22");
cgi.g_dd_23_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_23_varindex, "g_23");
cgi.g_dd_33_data = Cactus_gridfn_data_ptr(GH, cgi.g_dd_33_varindex, "g_33");
cgi.K_dd_11_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_11_varindex, "K_11");
cgi.K_dd_12_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_12_varindex, "K_12");
cgi.K_dd_13_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_13_varindex, "K_13");
cgi.K_dd_22_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_22_varindex, "K_22");
cgi.K_dd_23_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_23_varindex, "K_23");
cgi.K_dd_33_data = Cactus_gridfn_data_ptr(GH, cgi.K_dd_33_varindex, "K_33");
cgi.psi_data     = Cactus_gridfn_data_ptr(GH, cgi.psi_varindex,     "psi",
					  cgi.Cactus_conformal_metric);
}
	  }

//******************************************************************************

//
// This function gets the Cactus data pointer for a single gridfn, and
// optionally checks to make sure this is non-NULL.
//
// Arguments:
// gridfn_name[] = The character-string name of the grid function;
//		   this is used only for formatting error messages.
// check_for_NULL = true ==> check to make sure the data pointer is non-NULL
//		    false ==> skip this check (presumably because a NULL
//			      pointer is ok)
//
namespace {
const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
					const char gridfn_name[],
					bool check_for_NULL /* = true */)
{
const int time_level = 0;

//
// CCTK_VarDataPtrI() returns a  void * , but we need a  const CCTK_REAL*;
// since  static_cast<>  won't change const-ness, we need a 2-stage cast
// for this:
//
const CCTK_REAL* data_ptr = static_cast<const fp*>(
			       const_cast<const void *>(
				  CCTK_VarDataPtrI(GH, time_level, varindex)
						       )
						  );

if (check_for_NULL && (data_ptr == NULL))
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   Cactus_gridfn_data_ptr(): got unexpected NULL data pointer\n"
"                             for Cactus geometry gridfn!\n"
"                             name=\"%s\" varindex=%d"
		   ,
		   gridfn_name, varindex);			/*NOTREACHED*/

return data_ptr;
}
	  }

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function implements  AHFinderDirect::method == "horizon function",
// that is, it evaluates the Theta(h) function for a single apparent horizon.
//
// Note that if we decide to output h, we output it *after* any Theta(h)
// evaluation or horizon finding has been done, to ensure that all the
// ghost zones are filled in in case we need to print them.
//
// Arguments:
// timer_handle = a valid Cactus timer handle if we want to time the
//		  apparent horizon process, or -ve to skip this
//		  (we only time the computation, not the file I/O)
// output_[hH] = flags to control whether or not we should output [hH]
// active_flag = Controls whether this processor is "active":
//		 true ==> Normal operation.
//		 false ==> We evaluate Theta(h) as usual, but don't write
//			   any output files.
// hn = the horizon number (used only in naming output files)
// N_horizons = the total number of horizon(s) being searched for number
//		(used only in formatting info messages)
//
// Results:
// This function returns true if the evaluation was successful, false
// otherwise.
//
namespace {
bool do_evaluate_expansion
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info, bool output_h, bool output_Theta,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps,
	 bool active_flag, int hn, int N_horizons)
{
jtutil::norm<fp> Theta_norms;

if (timer_handle >= 0)
   then CCTK_TimerStartI(timer_handle);
const bool status = expansion(ps, cgi, gi, false, true, &Theta_norms);
if (timer_handle >= 0)
   then CCTK_TimerStopI(timer_handle);
if (!status)
   then return false;					// *** ERROR RETURN ***

if (Theta_norms.is_nonempty())	// might be empty if Theta(h) eval failed
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
		   Theta_norms.rms_norm(), Theta_norms.infinity_norm());

if (active_flag && output_h)
   then output_gridfn(ps, gfns::gfn__h,
		      IO_info, IO_info.h_base_file_name,
		      hn, verbose_info.print_algorithm_details);
if (active_flag && output_Theta)
   then output_gridfn(ps, gfns::gfn__Theta,
		      IO_info, IO_info.Theta_base_file_name,
		      hn, verbose_info.print_algorithm_details);

return true;						// *** NORMAL RETURN ***
}
	  }

//******************************************************************************

//
// This function implements  AHFinderDirect::method == "test Jacobian",
// that is, it computes and prints the Jacobian matrix J[Theta(h)] for a
// single apparent horizon.  The computation may optionally be done
// in several different ways, in which case all the resulting Jacobian
// matrices are printed, as are their differences.  Alternatively, only
// the numerical perturbation computation may be done/printed.
//
// Arguments:
// timer_handle = a valid Cactus timer handle if we want to time the
//		  apparent horizon process, or -ve to skip this
//		  (we only time the computation, not the file I/O)
// test_all_Jacobian_methods = true ==> Test all known methods of computing
//					the Jacobian matrix, and print all
//					the resulting Jacobian matrices
//					and their differences.
//			       false ==> Just test/print the numerical
//					 perturbation calculation.  (This
//					 may be useful if one or more of
//					 the other methods is broken.) 
// active_flag = Controls whether this processor is "active":
//		 true ==> Normal operation.
//		 false ==> We compute the Jacobian(s) as usual, but don't
//			   write any output files.
// hn = the horizon number (used only in naming output files)
// N_horizons = the total number of horizon(s) being searched for number
//		(used only in formatting info messages)
//
// Results:
// This function returns true if the Jacobian computation was successful,
// false otherwise.
//
namespace {
bool do_test_expansion_Jacobian
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info,
	 bool test_all_Jacobian_methods,
	 struct Jacobian_info& Jac_info,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps, Jacobian_matrix& Jac,
	 bool active_flag, int hn, int N_horizons)
{
// numerical perturbation
Jacobian_matrix* Jac_NP_ptr = & Jac;
if (!  expansion(ps, cgi, gi, true))
   then return false;					// *** ERROR RETURN ***
Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
if (! expansion_Jacobian(ps, *Jac_NP_ptr,
			 cgi, gi, Jac_info,
			 true))	// print msgs
   then return false;					// *** ERROR RETURN ***

Jacobian_matrix* Jac_SD_FDdr_ptr = NULL;
if (test_all_Jacobian_methods)
   then {
	// symbolic differentiation with finite diff d/dr
	Jac_SD_FDdr_ptr
	   = create_Jacobian_matrix(ps,
				    cgi, gi,
				    Jac_info,
				    verbose_info.print_algorithm_details);
	if (! expansion(ps, cgi, gi, true))
	   then return false;				// *** ERROR RETURN ***
	Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr;
	if (! expansion_Jacobian(ps, *Jac_SD_FDdr_ptr,
				 cgi, gi, Jac_info,
				 true))	// print msgs
	   then return false;				// *** ERROR RETURN ***
	}

if (active_flag)
   then output_Jacobians(ps,
			 Jac_NP_ptr, Jac_SD_FDdr_ptr,
			 IO_info, IO_info.Jacobian_base_file_name,
			 hn, true);

return true;						// *** NORMAL RETURN ***
}
	  }

//******************************************************************************

//
// This function implements  AHFinderDirect::method == "find horizon",
// that is, it finds the (an) apparent horizon.
//
// Arguments:
// timer_handle = a valid Cactus timer handle if we want to time the
//		  apparent horizon process, or -ve to skip this
//		  (we only time the computation, not the file I/O)
// initial_find_flag = true ==> This is the first time we've tried to find
//				this horizon, or we've tried before but
//				failed on our most recent previous attempt.
//				Thus we should use "initial" parameters
//				for the Newton iteration.
//		       flase ==> We've tried to find this horizon before,
//				 and we succeeded on our most recent
//				 previous attempt.  Thus we should use
//				 "subsequent" parameters for the Newton
//				 iteration.
// active_flag = Controls whether this processor is "active":
//		 true ==> Normal operation.
//		 false ==> We find the horizon as usual, but don't write
//			   any output files.
// hn = the horizon number (used only in naming output files)
// N_horizons = the total number of horizon(s) being searched for number
//		(used only in formatting info messages)
//
// Results:
// This function returns true if the apparent horizon was found
// successfully, false otherwise.
//
namespace {
bool do_find_horizon
	(const struct verbose_info& verbose_info, int timer_handle,
	 struct IO_info& IO_info, bool output_h, bool output_Theta,
	 const struct solver_info& solver_info, bool initial_find_flag,
	 const struct Jacobian_info& Jac_info,
	 struct cactus_grid_info& cgi, struct geometry_info& gi,
	 patch_system& ps, Jacobian_matrix& Jac,
	 bool active_flag, int hn, int N_horizons)
{
if (verbose_info.print_algorithm_highlights)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "searching for horizon %d/%d",
		   hn, N_horizons);

if (timer_handle >= 0)
   then CCTK_TimerStartI(timer_handle);
const bool status = Newton_solve(ps, Jac,
				 cgi, gi, Jac_info,
				 solver_info, initial_find_flag, IO_info,
				 active_flag, hn, verbose_info);
if (timer_handle >= 0)
   then CCTK_TimerStopI(timer_handle);
if (! status)
   then return false;					// *** ERROR RETURN ***

if (active_flag && output_h)
   then output_gridfn(ps, gfns::gfn__h,
		      IO_info, IO_info.h_base_file_name,
		      hn, verbose_info.print_algorithm_details);
if (active_flag && output_Theta)
   then output_gridfn(ps, gfns::gfn__Theta,
		      IO_info, IO_info.Theta_base_file_name,
		      hn, verbose_info.print_algorithm_details);

return true;						// *** NORMAL RETURN ***
}
	  }

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// Given that an apparent horizon has been found, this function computes
// various black hole diagnostics.
//
// Inputs (gridfns)
// h		# ghosted
// one		# nominal
// global_[xyz]	# nominal
//
namespace {
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)
{
//
// compute raw surface integrals
//
fp integral_one = surface_integral(ps, gfns::gfn__one,
				   BH_diagnostics_info.integral_method);
fp integral_h = surface_integral(ps, gfns::gfn__h,
				 BH_diagnostics_info.integral_method);
fp integral_x = surface_integral(ps, gfns::gfn__global_x,
				 BH_diagnostics_info.integral_method);
fp integral_y = surface_integral(ps, gfns::gfn__global_y,
				 BH_diagnostics_info.integral_method);
fp integral_z = surface_integral(ps, gfns::gfn__global_z,
				 BH_diagnostics_info.integral_method);

//
// adjust integrals to take into account patch system not covering full sphere
//
switch	(ps.type())
	{
case patch_system::patch_system__full_sphere:
	break;
case patch_system::patch_system__plus_z_hemisphere:
	integral_one *= 2.0;
	integral_h *= 2.0;
	integral_x *= 2.0;
	integral_y *= 2.0;
	integral_z = integral_one * ps.origin_z();
	break;
case patch_system::patch_system__plus_xy_quadrant_mirrored:
case patch_system::patch_system__plus_xy_quadrant_rotating:
	integral_one *= 4.0;
	integral_h *= 4.0;
	integral_x = integral_one * ps.origin_x();
	integral_y = integral_one * ps.origin_y();
	integral_z *= 4.0;
	break;
case patch_system::patch_system__plus_xz_quadrant_rotating:
	integral_one *= 4.0;
	integral_h *= 4.0;
	integral_x = integral_one * ps.origin_x();
	integral_y *= 4.0;
	integral_z = integral_one * ps.origin_z();
	break;
case patch_system::patch_system__plus_xyz_octant_mirrored:
case patch_system::patch_system__plus_xyz_octant_rotating:
	integral_one *= 8.0;
	integral_h *= 8.0;
	integral_x = integral_one * ps.origin_x();
	integral_y = integral_one * ps.origin_y();
	integral_z = integral_one * ps.origin_z();
	break;
default:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n"
"                             (this should never happen!)"
,
		   int(ps.type()));				/*NOTREACHED*/
	}

BH_diagnostics.centroid_x = integral_x / integral_one;
BH_diagnostics.centroid_y = integral_y / integral_one;
BH_diagnostics.centroid_z = integral_z / integral_one;

BH_diagnostics.area = integral_one;
BH_diagnostics.circumference_xy = ps.xy_circumference
    (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,
     BH_diagnostics_info.integral_method);
BH_diagnostics.circumference_xz = ps.xz_circumference
    (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,
     BH_diagnostics_info.integral_method);
BH_diagnostics.circumference_yz = ps.yz_circumference
    (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,
     BH_diagnostics_info.integral_method);
BH_diagnostics.mean_radius = integral_h / integral_one;
BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
}
	  }

//******************************************************************************

//
// 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);
}
	  }