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

#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 "horizon_sequence.hh"
#include "BH_diagnostics.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);

void do_evaluate_expansions(int my_proc, int N_horizons,
			    horizon_sequence& hs,
			       struct AH_info* const my_AH_info[],
			    const struct cactus_grid_info& cgi,
			    const struct geometry_info& gi,
			    const struct IO_info& IO_info,
			    const struct error_info& error_info,
			    const struct verbose_info& verbose_info,
			    int timer_handle);
void do_test_expansion_Jacobians(int my_proc, int N_horizons,
				 struct AH_info* const my_AH_info[],
				 const struct cactus_grid_info& cgi,
				 const struct geometry_info& gi,
				       struct Jacobian_info& Jac_info,
				 bool test_all_Jacobian_compute_methods,
				 const struct IO_info& IO_info,
				 const struct error_info& error_info,
				 const struct verbose_info& verbose_info,
				 int timer_handle);
	  }

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

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

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

const int my_proc = state.my_proc;
horizon_sequence& hs = *state.my_hs;
const bool active_flag = hs.has_genuine_horizons();

      struct cactus_grid_info&          cgi =  state.cgi;
const struct    geometry_info&           gi =  state.gi;
      struct    Jacobian_info&     Jac_info =  state.Jac_info;
      struct          IO_info&      IO_info =  state.IO_info;
const struct       error_info&   error_info =  state.error_info;
const struct     verbose_info& verbose_info =  state.verbose_info;

// 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 cgi.Cactus_conformal_metric = false;
else if (CCTK_Equals(metric_type, "static conformal"))
   then cgi.Cactus_conformal_metric = (conformal_state > 0);
else	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
		   metric_type);				/*NOTREACHED*/

// get the Cactus time step and decide if we want to output h and/or Theta now
IO_info.time_iteration = cctk_iteration;
IO_info.time           = cctk_time;
IO_info.output_h
   = (IO_info.how_often_to_output_h > 0)
     && ((IO_info.time_iteration % IO_info.how_often_to_output_h) == 0);
IO_info.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 (gi.geometry_method == geometry__local_interp_from_Cactus_grid)
   then setup_Cactus_gridfn_data_ptrs(cctkGH, cgi);

// set initial guess for any (genuine) horizons that need it,
// i.e. for any (genuine) horizons where we didn't find the horizon previously
	for (int hn = hs.init_hn() ; hs.is_genuine() ; hn = hs.next_hn())
	{
	assert( state.my_AH_info[hn] != NULL );
	struct AH_info& AH_info = *state.my_AH_info[hn];
	if (AH_info.found_flag)
	   then AH_info.initial_find_flag = false;
	   else {
		patch_system& ps = *AH_info.ps_ptr;
		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);
		AH_info.initial_find_flag = true;
		}
	}

//
// now the main horizon finding (or other computation)
//
switch	(state.method)
	{
case method__evaluate_expansions:
	do_evaluate_expansions(my_proc, state.N_horizons,
			       *state.my_hs, state.my_AH_info,
			       cgi, gi, IO_info,
			       error_info, verbose_info,
			       state.timer_handle);
	break;

case method__test_expansion_Jacobians:
	do_test_expansion_Jacobians(my_proc, state.N_horizons,
				    state.my_AH_info,
				    cgi, gi, Jac_info, 
				    (test_all_Jacobian_compute_methods != 0),
				    IO_info, error_info, verbose_info,
				    state.timer_handle);
	break;

case method__find_horizons:
	if (state.timer_handle >= 0)
	   then CCTK_TimerStartI(state.timer_handle);
	Newton(cctkGH,
	       state.N_procs, state.N_active_procs, my_proc,
	       *state.my_hs, state.my_AH_info,
	       cgi, gi, Jac_info, state.solver_info,
	       IO_info, state.BH_diagnostics_info,
	       error_info, verbose_info);
	if (state.timer_handle >= 0)
	   then CCTK_TimerStopI(state.timer_handle);
	break;

default:
	CCTK_VWarn(FATAL_ERROR, __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 {
	CCTK_VInfo(CCTK_THORNSTRING,
		   "timer stats for computation:");
	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(FATAL_ERROR, __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":
// On processor #0 it evaluates the Theta(h) function for each apparent
// horizon (and does any I/O desired); on other processors it does N_horizons
// dummy evaluations on horizon #0.
//
// 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)
//
namespace {
void do_evaluate_expansions(int my_proc, int N_horizons,
			    horizon_sequence& hs,
			       struct AH_info* const my_AH_info[],
			    const struct cactus_grid_info& cgi,
			    const struct geometry_info& gi,
			    const struct IO_info& IO_info,
			    const struct error_info& error_info,
			    const struct verbose_info& verbose_info,
			    int timer_handle)
{
const bool active_flag = (my_proc == 0);

if (active_flag)
   then {
	assert( hs.N_horizons() == N_horizons );
	assert( hs.my_N_horizons() == N_horizons );

		for (int hn = hs.init_hn() ;
		     hs.is_genuine() ;
		     hn = hs.next_hn())
		{
		assert( my_AH_info[hn] != NULL );
		struct AH_info& AH_info = *my_AH_info[hn];
		patch_system& ps = *AH_info.ps_ptr;

		if (timer_handle >= 0)
		   then CCTK_TimerStartI(timer_handle);
		jtutil::norm<fp> Theta_norms;
		const bool Theta_ok = expansion(&ps,
						cgi, gi,
						error_info, true,// initial eval
						false,	// no Jacobian coeffs
						true,	// yes, print msgs
						&Theta_norms);
		if (timer_handle >= 0)
		   then CCTK_TimerStopI(timer_handle);

		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 (Theta_ok)
		   then {
			CCTK_VInfo(CCTK_THORNSTRING,
			   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
			if (IO_info.output_Theta)
			   then output_gridfn(ps, gfns::gfn__Theta,
					      IO_info, IO_info
						       .Theta_base_file_name,
					      hn, verbose_info
						  .print_algorithm_details);
			}
		}
	}
   else {
		for (int i = 0 ; i < N_horizons ; ++i)
		{
		expansion(NULL,			// dummy computation
			  cgi, gi,
			  error_info, true);	// initial evaluation
		}
	}
}
	  }

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

//
// This function implements
//  AHFinderDirect::method == "test expansion Jacobians":
// On processor #0 it computes and prints the Jacobian matrix J[Theta(h)]
// function for horizon #1; on other processors it does dummy Jacobian
// computations.
//
// The Jacobian 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_compute_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.) 
//
namespace {
void do_test_expansion_Jacobians(int my_proc, int N_horizons,
				 struct AH_info* const my_AH_info[],
				 const struct cactus_grid_info& cgi,
				 const struct geometry_info& gi,
				       struct Jacobian_info& Jac_info,
				 bool test_all_Jacobian_compute_methods,
				 const struct IO_info& IO_info,
				 const struct error_info& error_info,
				 const struct verbose_info& verbose_info,
				 int timer_handle)
{
const bool active_flag = (my_proc == 0);
assert(N_horizons >= 1);

const bool print_msg_flag = true;
const int hn = 1;

struct AH_info* AH_info_ptr = active_flag ? my_AH_info[hn]      : NULL;
patch_system*        ps_ptr = active_flag ? AH_info_ptr->ps_ptr : NULL;

//
// numerical-perturbation Jacobian
//
Jacobian* Jac_NP_ptr = active_flag ? AH_info_ptr->Jac_ptr : NULL;
expansion(ps_ptr,
	  cgi, gi,
	  error_info, true);		// initial evaluation
Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
expansion_Jacobian(ps_ptr, Jac_NP_ptr,
		   cgi, gi, Jac_info,
		   error_info, true,	// initial evaluation
		   print_msg_flag);

Jacobian* Jac_SD_FDdr_ptr = NULL;
if (test_all_Jacobian_compute_methods)
   then {
	// symbolic differentiation with finite diff d/dr
	Jac_SD_FDdr_ptr = active_flag
			  ? new_Jacobian(Jac_info.Jacobian_store_solve_method,
					 *ps_ptr,
					 verbose_info.print_algorithm_details)
			  : NULL;
	expansion(ps_ptr,
		  cgi, gi,
		  error_info, true,	// initial evaluation
		  true);		// compute SD Jacobian coeffs
	Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
			   cgi, gi, Jac_info,
			   error_info, true,	// initial evaluation
			   print_msg_flag);
	}

if (active_flag)
   then output_Jacobians(*ps_ptr,
			 Jac_NP_ptr, Jac_SD_FDdr_ptr,
			 IO_info, IO_info.Jacobian_base_file_name,
			 hn, print_msg_flag);
}
	  }