// 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
// Cactus_gridfn_data_ptr - get Cactus data pointer for CCTK_REAL gridfn
///
/// decode_method - decode the method parameter
/// decode_verbose_level - decode the verbose_level parameter
///
/// choose_patch_system_type - choose patch system type
///

#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"

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

//
// ***** 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 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 all our
// persistent data structures.  (These are stored in  struct state .)
//
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);

struct verbose_info& verbose_info = state.verbose_info;
verbose_info.verbose_level = decode_verbose_level(verbose_level);
verbose_info.print_physics_highlights
   = (state.verbose_info.verbose_level >= verbose_level__physics_highlights);
verbose_info.print_physics_details
   = (state.verbose_info.verbose_level >= verbose_level__physics_details);
verbose_info.print_algorithm_highlights
   = (state.verbose_info.verbose_level >= verbose_level__algorithm_highlights);
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);

struct IO_info& IO_info = state.IO_info;
IO_info.file_format = decode_file_format(file_format);
IO_info.output_ghost_zones_for_h  = (output_ghost_zones_for_h != 0);
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.Delta_h_base_file_name  = Delta_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.output_initial_guess = (output_initial_guess != 0);
solver_info.output_h              = (output_h != 0);
solver_info.output_H              = (output_H_of_h != 0);
solver_info.debugging_output_at_each_Newton_iteration
			= (debugging_output_at_each_Newton_iteration != 0);
solver_info.max_Newton_iterations__initial
			= max_Newton_iterations__initial;
solver_info.max_Newton_iterations__subsequent
			= max_Newton_iterations__subsequent;
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];

// dummy values -- real ones filled in in AHFinderDirect_find_horizons()
cgi.Cactus_conformal_metric = false;
cgi.psi_data = NULL;

cgi.g_dd_11_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxx");
cgi.g_dd_12_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxy");
cgi.g_dd_13_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gxz");
cgi.g_dd_22_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gyy");
cgi.g_dd_23_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gyz");
cgi.g_dd_33_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::gzz");
cgi.K_dd_11_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxx");
cgi.K_dd_12_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxy");
cgi.K_dd_13_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kxz");
cgi.K_dd_22_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kyy");
cgi.K_dd_23_data = Cactus_gridfn_data_ptr(cctkGH, "ADMBase::kyz");
cgi.K_dd_33_data = Cactus_gridfn_data_ptr(cctkGH, "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);

	// store the initial guess parameters
	AH_ptr->initial_guess_info.method
		= decode_initial_guess_method(initial_guess_method);
	// ... Kerr/Kerr
	AH_ptr->initial_guess_info.Kerr_Kerr_info.x_posn
		= initial_guess__Kerr_Kerr__x_posn[hn];
	AH_ptr->initial_guess_info.Kerr_Kerr_info.y_posn
		= initial_guess__Kerr_Kerr__y_posn[hn];
	AH_ptr->initial_guess_info.Kerr_Kerr_info.z_posn
		= initial_guess__Kerr_Kerr__z_posn[hn];
	AH_ptr->initial_guess_info.Kerr_Kerr_info.mass
		= initial_guess__Kerr_Kerr__mass[hn];
	AH_ptr->initial_guess_info.Kerr_Kerr_info.spin
		= initial_guess__Kerr_Kerr__spin[hn];
	// ... Kerr/Kerr-Schild
	AH_ptr->initial_guess_info.Kerr_KerrSchild_info.x_posn
		= initial_guess__Kerr_KerrSchild__x_posn[hn];
	AH_ptr->initial_guess_info.Kerr_KerrSchild_info.y_posn
		= initial_guess__Kerr_KerrSchild__y_posn[hn];
	AH_ptr->initial_guess_info.Kerr_KerrSchild_info.z_posn
		= initial_guess__Kerr_KerrSchild__z_posn[hn];
	AH_ptr->initial_guess_info.Kerr_KerrSchild_info.mass
		= initial_guess__Kerr_KerrSchild__mass[hn];
	AH_ptr->initial_guess_info.Kerr_KerrSchild_info.spin
		= initial_guess__Kerr_KerrSchild__spin[hn];
	// ... sphere
	AH_ptr->initial_guess_info.sphere_info.x_center
		= initial_guess__sphere__x_center[hn];
	AH_ptr->initial_guess_info.sphere_info.y_center
		= initial_guess__sphere__y_center[hn];
	AH_ptr->initial_guess_info.sphere_info.z_center
		= initial_guess__sphere__z_center[hn];
	AH_ptr->initial_guess_info.sphere_info.radius
		= initial_guess__sphere__radius[hn];
	// ... ellipsoid
	AH_ptr->initial_guess_info.ellipsoid_info.x_center
		= initial_guess__ellipsoid__x_center[hn];
	AH_ptr->initial_guess_info.ellipsoid_info.y_center
		= initial_guess__ellipsoid__y_center[hn];
	AH_ptr->initial_guess_info.ellipsoid_info.z_center
		= initial_guess__ellipsoid__z_center[hn];
	AH_ptr->initial_guess_info.ellipsoid_info.x_radius
		= initial_guess__ellipsoid__x_radius[hn];
	AH_ptr->initial_guess_info.ellipsoid_info.y_radius
		= initial_guess__ellipsoid__y_radius[hn];
	AH_ptr->initial_guess_info.ellipsoid_info.z_radius
		= initial_guess__ellipsoid__z_radius[hn];

	// initialize the BH diagnostics
	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 gets the Cactus data pointer for a gridfn, and checks
// to make sure this is non-NULL.
//
const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, const char gridfn_name[])
{
const int time_level = 0;

//
// CCTK_VarDataPtr() 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 *const data_ptr
	= static_cast<const fp*>(
	     const_cast<const void *>(
		CCTK_VarDataPtr(GH, time_level, gridfn_name)
				     )
				);

if (data_ptr == NULL)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"***** Cactus_gridfn_data_ptr(): can't find gridfn \"%s\"!",
		   gridfn_name);				/*NOTREACHED*/

return data_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 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_VInfo(CCTK_THORNSTRING,
		   "      bitant mode but patch system origin_z=%g != 0\n",
			   double(origin_z));
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      ==> using \"full sphere\" patch system");
		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;
	}
}
	  }