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
///
/// 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 "stl_vector.hh"

#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 find_horizon(enum method method,
		  const struct verbose_info& verbose_info, int timer_handle,
		  struct IO_info& IO_info,
		  struct Jacobian_info& Jac_info,
		  const struct solver_info& solver_info, bool initial_find_flag,
		  struct cactus_grid_info& cgi, struct geometry_info& gi,
		  patch_system& ps, Jacobian* Jac_ptr,
		  int hn, int N_horizons);

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

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

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

IO_info.time_iteration = cctk_iteration;

// what are the semantics of the Cactus gxx variables?
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*/

//
// we need to re-fetch the Cactus data pointers at least each time step,
// because they change each time Cactus rotates the time levels
//
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_ptrs[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, verbose_info);
		if (solver_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);
		}

	AH_info.AH_found
	   = find_horizon(state.method,
			  verbose_info, state.timer_handle,
			  IO_info, state.Jac_info,
			  solver_info, initial_find_flag,
			  state.cgi, state.gi,
			  ps, AH_info.Jac_ptr,
			  hn, state.N_horizons);

	if (AH_info.AH_found)
	   then {
		BH_diagnostics(state.surface_integral_method,
			       verbose_info,
			       AH_info);

		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));
			CCTK_VInfo(CCTK_THORNSTRING,
				   "estimated mass %.10g",
				   double(AH_info.mass));
			}
		}
	}

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 finds (or more accurately tries to find) a single
// apparent horizon.
//
// Note that if we decide to output h, we output it *after* any H(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)
// initial_find_flag = true ==> This is the first attempt to find this
//				horizon, or this is a subsequent attempt
//				and the immediately previous attempt failed.
//		       false ==> This isn't the first attempt to find this
//				 horizon, and we found it successfully on
//				 the immediately previous attempt.
// Jac_ptr = may be NULL if no Jacobian is needed (depending on  method)
// 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 computation succeeds; false if it fails.
// If  method  specifies finding the (an) apparent horizon, then "success"
// means finding an h satisfying H(h) = 0 to within the error tolerances.
// Otherwise, "success" means successfully evaluating the horizon function
// and/or its Jacobian, as appropriate.
//
namespace {
bool find_horizon(enum method method,
		  const struct verbose_info& verbose_info, int timer_handle,
		  struct IO_info& IO_info,
		  struct Jacobian_info& Jac_info,
		  const struct solver_info& solver_info, bool initial_find_flag,
		  struct cactus_grid_info& cgi, struct geometry_info& gi,
		  patch_system& ps, Jacobian* Jac_ptr,
		  int hn, int N_horizons)
{
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_H
   = (IO_info.how_often_to_output_H > 0)
     && ((IO_info.time_iteration % IO_info.how_often_to_output_H) == 0);

switch	(method)
	{
// just evaluate the horizon function
case method__horizon_function:
	  {
	jtutil::norm<fp> H_norms;

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

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

	if (output_h)
	   then output_gridfn(ps, gfns::gfn__h,
			      IO_info, IO_info.h_base_file_name,
			      hn, verbose_info.print_algorithm_details);
	if (output_H)
	   then output_gridfn(ps, gfns::gfn__H,
			      IO_info, IO_info.H_base_file_name,
			      hn, true);
	return true;					// *** NORMAL RETURN ***
	  }

// just compute/print the NP Jacobian
case method__Jacobian_test_NP_only:
	  {
	Jacobian& Jac_NP = *Jac_ptr;
	if (! horizon_function(ps, cgi, gi, true))
	   then return false;				// *** ERROR RETURN ***
	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
	if (!  horizon_Jacobian(ps, Jac_NP,
				cgi, gi, Jac_info,
				true))
	   then return false;				// *** ERROR RETURN ***

	print_Jacobians(ps,
			& Jac_NP, NULL,
			IO_info, IO_info.Jacobian_base_file_name,
			hn, true);
	return true;					// *** NORMAL RETURN ***
	  }

// compute/print the Jacobian by all possible methods
case method__Jacobian_test:
	  {
	Jacobian& Jac_NP = *Jac_ptr;
	if (!  horizon_function(ps, cgi, gi, true))
	   then return false;				// *** ERROR RETURN ***
	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
	if (! horizon_Jacobian(ps, Jac_NP,
			       cgi, gi, Jac_info,
			       true))
	   then return false;				// *** ERROR RETURN ***

	// symbolic differentiation with finite diff d/dr
	Jacobian& Jac_SD_FDdr
		= new_Jacobian(ps, Jac_info.Jacobian_storage_method);
	if (! horizon_function(ps, cgi, gi, true))
	   then return false;				// *** ERROR RETURN ***
	Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr;
	if (! horizon_Jacobian(ps, Jac_SD_FDdr,
			       cgi, gi, Jac_info,
			       true))
	   then return false;				// *** ERROR RETURN ***

	print_Jacobians(ps,
			& Jac_NP, & Jac_SD_FDdr,
			IO_info, IO_info.Jacobian_base_file_name,
			hn, true);
	return true;					// *** NORMAL RETURN ***
	  }

// find the apparent horizon via the Newton solver
case method__Newton_solve:
	  {
	Jacobian& Jac = *Jac_ptr;

	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, hn, verbose_info);
	if (timer_handle >= 0)
	   then CCTK_TimerStopI(timer_handle);
	if (! status)
	   then return false;				// *** ERROR RETURN ***

	if (output_h)
	   then output_gridfn(ps, gfns::gfn__h,
			      IO_info, IO_info.h_base_file_name,
			      hn, verbose_info.print_algorithm_details);
	if (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 ***
	  }

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

/*NOTREACHED*/
}
	  }

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

//
// Given that an apparent horizon has been found, this function computes
// various BH diagnostics.
//
// Inputs (gridfns)
// h		# ghosted
// one		# nominal
// global_[xyz]	# nominal
//
namespace {
void BH_diagnostics(enum patch::integration_method surface_integral_method,
		    const struct verbose_info& verbose_info,
		    struct AH_info& AH_info)
{
const patch_system& ps = * AH_info.ps_ptr;

//
// compute raw surface integrals
//
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
//
switch	(ps.type())
	{
case patch_system::full_sphere_patch_system:
	break;
case patch_system::plus_z_hemisphere_patch_system:
	integral_one *= 2.0;
	integral_x *= 2.0;
	integral_y *= 2.0;
	integral_z = integral_one * ps.origin_z();
	break;
case patch_system::plus_xy_quadrant_patch_system:
	integral_one *= 4.0;
	integral_x = integral_one * ps.origin_x();
	integral_y = integral_one * ps.origin_y();
	integral_z *= 4.0;
	break;
case patch_system::plus_xz_quadrant_patch_system:
	integral_one *= 4.0;
	integral_x = integral_one * ps.origin_x();
	integral_y *= 4.0;
	integral_z = integral_one * ps.origin_z();
	break;
case patch_system::plus_xyz_octant_patch_system:
	integral_one *= 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"
"   BH_diagnostics(): unknown ps.type()=(int)%d!\n"
"                     (this should never happen!)"
,
		   int(ps.type()));				/*NOTREACHED*/
	}

AH_info.area = integral_one;
AH_info.mass = sqrt(AH_info.area / (16.0*PI));

AH_info.centroid_x = integral_x / integral_one;
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);
}
	  }