path: root/src/driver/setup.cc

// setup.cc -- top level driver to setup our persistent data structures
// $Header$
//
// <<<prototypes for functions local to this file>>>
// <<<access to persistent data>>>
//
// AHFinderDirect_driver - top-level driver to setup persistent data structures
///
/// decode_method - decode the method parameter
/// decode_verbose_level - decode the verbose_level parameter
/// decode_file_format - decode the file_format parameter
/// decode_Jacobian_method - decode the Jacobian_method parameter
/// decode_geometry_method - decode the geometry_method parameter
///
/// choose_patch_system_type - choose patch system type
///

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

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

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

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

//
// ***** prototypes for functions local to this file *****
//

namespace {
enum method
  decode_method(const char method_string[]);
enum verbose_level
  decode_verbose_level(const char verbose_level_string[]);
enum file_format
  decode_file_format(const char file_format_string[]);
enum Jacobian_method
  decode_Jacobian_method(const char Jacobian_method_string[]);
enum geometry_method
  decode_geometry_method(const char geometry_method_string[]);

enum patch_system::patch_system_type
  choose_patch_system_type(const char grid_domain[],
			   fp origin_x, fp origin_y, fp origin_z);
	  };

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

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

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

//
// This function is called by the Cactus scheduler to set up the
// following data structures:
//	state.gi
//	state.cgi
//
extern "C"
  void AHFinderDirect_setup(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS

CCTK_VInfo(CCTK_THORNSTRING,
	   "initializing AHFinderDirect data structures");

state.N_horizons = N_horizons;
CCTK_VInfo(CCTK_THORNSTRING,
           "             to search for %d horizon%s",
	   state.N_horizons, ((state.N_horizons == 1) ? "" : "s"));


//
// decode/copy parameters into structures
//

state.method = decode_method(method);
state.verbose_info.verbose_level = decode_verbose_level(verbose_level);
state.verbose_info.print_physics_highlights
   = (state.verbose_info.verbose_level >= verbose_level__physics_highlights);
state.verbose_info.print_physics_details
   = (state.verbose_info.verbose_level >= verbose_level__physics_details);
state.verbose_info.print_algorithm_highlights
   = (state.verbose_info.verbose_level >= verbose_level__algorithm_highlights);
state.verbose_info.print_algorithm_details
   = (state.verbose_info.verbose_level >= verbose_level__algorithm_details);
state.timer_handle = (print_timing_stats != 0)
		     ? CCTK_TimerCreateI()
		     : -1;
state.surface_integral_method
   = patch::decode_integration_method(surface_integral_method);

const struct verbose_info& verbose_info = state.verbose_info;

struct IO_info& IO_info = state.IO_info;
IO_info.output_initial_guess     = (output_initial_guess != 0);
IO_info.output_h_and_H_at_each_Newton_iteration
			      = (output_h_and_H_at_each_Newton_iteration != 0);
IO_info.output_h = (output_h != 0);
IO_info.output_H = (output_H_of_h != 0);
IO_info.file_format = decode_file_format(file_format);
IO_info.ASCII_file_name_extension = ASCII_file_name_extension;
IO_info.HDF5_file_name_extension  = HDF5_file_name_extension;
IO_info.h_base_file_name        = h_base_file_name;
IO_info.H_base_file_name        = H_of_h_base_file_name;
IO_info.Jacobian_base_file_name = Jacobian_base_file_name;
IO_info.time_iteration = 0;

struct Jacobian_info& Jac_info = state.Jac_info;
Jac_info.Jacobian_method = decode_Jacobian_method(Jacobian_method);
Jac_info.Jacobian_storage_method
	= decode_Jacobian_type(Jacobian_storage_method);
Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude;

struct solver_info& solver_info = state.solver_info;
solver_info.max_Newton_iterations = max_Newton_iterations;
solver_info.max_Delta_h_over_h           = max_Delta_h_over_h;
solver_info.H_norm_for_convergence       = H_norm_for_convergence;
solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
solver_info.final_H_update_if_exit_x_H_small
	= (final_H_update_if_exit_x_H_small != 0);

//
// set up the Cactus grid info
//
if (verbose_info.print_algorithm_highlights)
   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
struct cactus_grid_info& cgi = state.cgi;
cgi.GH = cctkGH;
cgi.coord_origin[0] = cctk_origin_space[0];
cgi.coord_origin[1] = cctk_origin_space[1];
cgi.coord_origin[2] = cctk_origin_space[2];
cgi.coord_delta[0] = cctk_delta_space[0];
cgi.coord_delta[1] = cctk_delta_space[1];
cgi.coord_delta[2] = cctk_delta_space[2];
cgi.gridfn_dims[0] = cctk_lsh[0];
cgi.gridfn_dims[1] = cctk_lsh[1];
cgi.gridfn_dims[2] = cctk_lsh[2];
// CCTK_VarDataPtr() returns a  void * , but we need a  const fp *;
// since  static_cast<>  won't change const-ness, we need a 2-stage cast
// for this:
#define DATA_PTR_CAST(void_ptr_)	\
	static_cast<const fp*>( const_cast<const void *>(void_ptr_) )
cgi.g_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxx"));
cgi.g_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy"));
cgi.g_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz"));
cgi.g_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy"));
cgi.g_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz"));
cgi.g_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz"));
cgi.K_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx"));
cgi.K_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy"));
cgi.K_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz"));
cgi.K_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy"));
cgi.K_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz"));
cgi.K_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));


//
// set up the geometry parameters and the geometry interpolator
//
struct geometry_info& gi = state.gi;
gi.geometry_method = decode_geometry_method(geometry_method);

// parameters for geometry_method = "interpolate from Cactus grid"
if (verbose_info.print_algorithm_highlights)
   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
if (gi.operator_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_setup(): couldn't find interpolator \"%s\"!",
		   geometry_interpolator_name);		/*NOTREACHED*/
gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
if (gi.param_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_setup(): bad geometry-interpolator parameter(s) \"%s\"!",
		   geometry_interpolator_pars);		/*NOTREACHED*/

// parameters for geometry_method = "Schwarzschild/EF"
gi.geometry__Schwarzschild_EF__mass     = geometry__Schwarzschild_EF__mass;
gi.geometry__Schwarzschild_EF__x_posn   = geometry__Schwarzschild_EF__x_posn;
gi.geometry__Schwarzschild_EF__y_posn   = geometry__Schwarzschild_EF__y_posn;
gi.geometry__Schwarzschild_EF__z_posn   = geometry__Schwarzschild_EF__z_posn;
gi.geometry__Schwarzschild_EF__epsilon  = geometry__Schwarzschild_EF__epsilon;
gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz;


//
// create AH-specific info for each AH
//

// set up the interpatch interpolator
if (verbose_info.print_algorithm_highlights)
   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up interpatch interpolator");
const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name);
if (interp_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_setup(): couldn't find interpolator \"%s\"!",
		   interpatch_interpolator_name);		/*NOTREACHED*/
const int interp_param_table_handle
	= Util_TableCreateFromString(interpatch_interpolator_pars);
if (interp_param_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_setup(): bad interpatch-interpolator parameter(s) \"%s\"!",
		   interpatch_interpolator_pars);		/*NOTREACHED*/

const bool ps_type_from_Cactus_grid_sym
   = STRING_EQUAL(patch_system_type, "match Cactus grid symmetry");

// we don't use the [0] array element
state.AH_info_ptrs.push_back(NULL);

	for (int hn = 1 ; hn <= N_horizons ; ++hn)
	{
	struct AH_info* AH_ptr = new AH_info;

	// decide what type of patch system this one should be
	const enum patch_system::patch_system_type ps_type
	   = ps_type_from_Cactus_grid_sym
	     ? // choose a patch system type based on grid::domain
	       // ... since we inherit from grid, we can see its
	       //     parameters, so grid::domain is just "domain" here
	       choose_patch_system_type(domain, origin_x[hn],
						origin_y[hn],
						origin_z[hn])
	     : patch_system::type_of_name(patch_system_type);
	if (verbose_info.print_algorithm_highlights)
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "   for horizon #%d/%d: %s patch system",
			   hn, int(N_horizons),
			   patch_system::name_of_type(ps_type));
		CCTK_VInfo(CCTK_THORNSTRING,
			   "                     at origin (%g,%g,%g)",
			   origin_x[hn], origin_y[hn], origin_z[hn]);
		}

	// create the patch system
	patch_system& ps
	   = *new
	      patch_system(origin_x[hn], origin_y[hn], origin_z[hn],
			   ps_type,
			   N_ghost_points, N_overlap_points, delta_drho_dsigma,
			   gfns::nominal_min_gfn, gfns::nominal_max_gfn,
			   gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
			   interp_handle, interp_param_table_handle,
			   verbose_info.print_algorithm_details);
	AH_ptr->ps_ptr = &ps;
	ps.set_gridfn_to_constant(1.0, gfns::gfn__one);

	// create the Jacobian matrix
	AH_ptr->Jac_ptr = (state.method == method__horizon_function)
			  ? NULL	// no Jacobian used for this case
			  : &new_Jacobian(ps, Jac_info.Jacobian_storage_method);

	AH_ptr->AH_found = false;
	AH_ptr->centroid_x = 0.0;
	AH_ptr->centroid_y = 0.0;
	AH_ptr->centroid_z = 0.0;
	AH_ptr->area = 0.0;
	AH_ptr->mass = 0.0;

	state.AH_info_ptrs.push_back(AH_ptr);
	}
}

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

//
// This function decodes the  method  parameter (string) into an
// internal enum for future use.
//
namespace {
enum method
  decode_method(const char method_string[])
{
if	(STRING_EQUAL(method_string, "horizon function"))
   then return method__horizon_function;
else if (STRING_EQUAL(method_string, "Jacobian test"))
   then return method__Jacobian_test;
else if (STRING_EQUAL(method_string, "Jacobian test (NP only)"))
   then return method__Jacobian_test_NP_only;
else if (STRING_EQUAL(method_string, "Newton solve"))
   then return method__Newton_solve;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "decode_method(): unknown method_string=\"%s\"!",
		   method_string);				/*NOTREACHED*/
}
	  }

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

//
// This function decodes the  verbose_level  parameter (string) into an
// internal enum for future use.
//
namespace {
enum verbose_level
  decode_verbose_level(const char verbose_level_string[])
{
if	(STRING_EQUAL(verbose_level_string, "physics highlights"))
   then return verbose_level__physics_highlights;
else if (STRING_EQUAL(verbose_level_string, "physics details"))
   then return verbose_level__physics_details;
else if (STRING_EQUAL(verbose_level_string, "algorithm highlights"))
   then return verbose_level__algorithm_highlights;
else if (STRING_EQUAL(verbose_level_string, "algorithm details"))
   then return verbose_level__algorithm_details;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"decode_verbose_level(): unknown verbose_level_string=\"%s\"!",
		   verbose_level_string);			/*NOTREACHED*/
}
	  }

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

//
// This function decodes the  file_format  parameter (string) into an
// internal enum for future use.
//
namespace {
enum file_format
  decode_file_format(const char file_format_string[])
{
if	(STRING_EQUAL(file_format_string, "ASCII"))
   then return file_format__ASCII;
else if (STRING_EQUAL(file_format_string, "HDF5"))
   then return file_format__HDF5;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"decode_file_format(): unknown file_format_string=\"%s\"!",
		   file_format_string);				/*NOTREACHED*/
}
	  }

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

//
// This function decodes the  Jacobian_method  parameter (string) into
// an internal enum for future use.
//
namespace {
enum Jacobian_method
  decode_Jacobian_method(const char Jacobian_method_string[])
{
if	(STRING_EQUAL(Jacobian_method_string,
		      "numerical perturbation"))
   then return Jacobian_method__numerical_perturb;
else if (STRING_EQUAL(Jacobian_method_string,
		      "symbolic differentiation with finite diff d/dr"))
   then return Jacobian_method__symbolic_diff_with_FD_dr;
else if (STRING_EQUAL(Jacobian_method_string,
		      "symbolic differentiation"))
   then return Jacobian_method__symbolic_diff;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"decode_Jacobian_method(): unknown Jacobian_method_string=\"%s\"!",
		   Jacobian_method_string);			/*NOTREACHED*/
}
	  }

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

//
// This function decodes the  geometry_method  parameter (string) into
// an internal enum for future use.
//
namespace {
enum geometry_method
  decode_geometry_method(const char geometry_method_string[])
{
if	(STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid"))
   then return geometry__interpolate_from_Cactus_grid;
else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF"))
   then return geometry__Schwarzschild_EF;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"decode_geometry_method(): unknown geometry_method_string=\"%s\"!",
		   geometry_method_string);			/*NOTREACHED*/
}
	  }

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

//
// This function chooses a patch system type based on the Cactus grid
// symmetry and the patch system's origin position.
//
// Arguments:
// grid_domain = The value of the Cactus grid::domain parameter.
// origin_[xyz] = The origin position for this patch system.
//
namespace {
enum patch_system::patch_system_type
  choose_patch_system_type(const char grid_domain[],
			   fp origin_x, fp origin_y, fp origin_z)
{
if	(STRING_EQUAL(grid_domain, "full"))
   then return patch_system::full_sphere_patch_system;

else if (STRING_EQUAL(grid_domain, "bitant"))
   then {
	if (origin_z == 0.0)
	   then return patch_system::plus_z_hemisphere_patch_system;
	   else {
		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   choose_patch_system_type()\n"
"        bitant mode but patch system origin_z=%g != 0\n"
"        ==> using \"full sphere\" patch system"
,
			   double(origin_z));
		return patch_system::full_sphere_patch_system;
		}
	}

else	{
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   choose_patch_system_type()\n"
"        grid::domain = \"%s\" not supported!\n"
"        ==> will try using \"full sphere\" patch system\n",
"            but this probably won't work :("
,
		   grid_domain);
	return patch_system::full_sphere_patch_system;
	}
}
	  }