// 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[],
			   const char grid_bitant_plane[],
			   const char grid_quadrant_direction[],
			   const char grid_rotation_axis[],
			   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,
	   "setting up 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;

struct IO_info& IO_info = state.IO_info;
IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
IO_info.output_initial_guess = (output_initial_guess != 0);
IO_info.how_often_to_output_h = how_often_to_output_h;
IO_info.how_often_to_output_H = how_often_to_output_H_of_h;
IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0);
IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_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.output_BH_diagnostics              = (output_BH_diagnostics != 0);
IO_info.BH_diagnostics_base_file_name      = BH_diagnostics_base_file_name;
IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension;
IO_info.time_iteration = 0;
IO_info.time           = 0.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.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];

cgi.g_dd_11_varindex = CCTK_VarIndex("ADMBase::gxx");
cgi.g_dd_12_varindex = CCTK_VarIndex("ADMBase::gxy");
cgi.g_dd_13_varindex = CCTK_VarIndex("ADMBase::gxz");
cgi.g_dd_22_varindex = CCTK_VarIndex("ADMBase::gyy");
cgi.g_dd_23_varindex = CCTK_VarIndex("ADMBase::gyz");
cgi.g_dd_33_varindex = CCTK_VarIndex("ADMBase::gzz");
cgi.K_dd_11_varindex = CCTK_VarIndex("ADMBase::kxx");
cgi.K_dd_12_varindex = CCTK_VarIndex("ADMBase::kxy");
cgi.K_dd_13_varindex = CCTK_VarIndex("ADMBase::kxz");
cgi.K_dd_22_varindex = CCTK_VarIndex("ADMBase::kyy");
cgi.K_dd_23_varindex = CCTK_VarIndex("ADMBase::kyz");
cgi.K_dd_33_varindex = CCTK_VarIndex("ADMBase::kzz");
cgi.psi_varindex     = CCTK_VarIndex("StaticConformal::psi");


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

gi.check_that_h_is_finite        = (check_that_h_is_finite        != 0);
gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0);


//
// other misc setup
//
state.BH_diagnostics_info.surface_integral_method
   = patch::decode_integration_method(surface_integral_method);


//
// 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)
	{
	if (verbose_info.print_algorithm_highlights)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "   setting up data structures for horizon %d of %d",
			   hn, int(N_horizons));

	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:: parameters
	       // ... we inherit from grid, and we ask for some of its
	       //     parameters in our param.ccl file; these appear as
	       //     ordinary Cactus parameters here, so (eg)
	       //     grid::domain is just "domain" here
	       choose_patch_system_type(/* grid:: */ domain,
					/* grid:: */ bitant_plane,
					/* grid:: */ quadrant_direction,
					/* grid:: */ rotation_axis,
					origin_x[hn],
					origin_y[hn],
					origin_z[hn])
	     : patch_system::type_of_name(patch_system_type);

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

	// create the Jacobian matrix
	if (verbose_info.print_algorithm_details)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "      constructing 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);

	if (verbose_info.print_algorithm_details)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "      setting initial guess parameters etc");

	// 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];
	// ... coordinate sphere
	AH_ptr->initial_guess_info.coord_sphere_info.x_center
		= initial_guess__coord_sphere__x_center[hn];
	AH_ptr->initial_guess_info.coord_sphere_info.y_center
		= initial_guess__coord_sphere__y_center[hn];
	AH_ptr->initial_guess_info.coord_sphere_info.z_center
		= initial_guess__coord_sphere__z_center[hn];
	AH_ptr->initial_guess_info.coord_sphere_info.radius
		= initial_guess__coord_sphere__radius[hn];
	// ... coordinate ellipsoid
	AH_ptr->initial_guess_info.coord_ellipsoid_info.x_center
		= initial_guess__coord_ellipsoid__x_center[hn];
	AH_ptr->initial_guess_info.coord_ellipsoid_info.y_center
		= initial_guess__coord_ellipsoid__y_center[hn];
	AH_ptr->initial_guess_info.coord_ellipsoid_info.z_center
		= initial_guess__coord_ellipsoid__z_center[hn];
	AH_ptr->initial_guess_info.coord_ellipsoid_info.x_radius
		= initial_guess__coord_ellipsoid__x_radius[hn];
	AH_ptr->initial_guess_info.coord_ellipsoid_info.y_radius
		= initial_guess__coord_ellipsoid__y_radius[hn];
	AH_ptr->initial_guess_info.coord_ellipsoid_info.z_radius
		= initial_guess__coord_ellipsoid__z_radius[hn];

	// initialize the BH diagnostics
	AH_ptr->AH_found = false;
	AH_ptr->BH_diagnostics.centroid_x = 0.0;
	AH_ptr->BH_diagnostics.centroid_y = 0.0;
	AH_ptr->BH_diagnostics.centroid_z = 0.0;
	AH_ptr->BH_diagnostics.area = 0.0;
	AH_ptr->BH_diagnostics.m_irreducible = 0.0;

	// set up the BH-diagnostics output files?
	AH_ptr->BH_diagnostics_fileptr
		= ( output_BH_diagnostics
		    && (state.method == method__find_horizon) )
		  ? setup_BH_diagnostics_output_file(IO_info, hn, N_horizons)
		  : NULL;

	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, "test Jacobian"))
   then return method__test_Jacobian;
else if (STRING_EQUAL(method_string, "find horizon"))
   then return method__find_horizon;
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_* = The values of the Cactus parameters
//		grid::domain
//		grid::bitant_plane
//		grid::quadrant_direction
//		grid::rotation_axis
// origin_[xyz] = The origin position for this patch system.
//
namespace {
enum patch_system::patch_system_type
  choose_patch_system_type(const char grid_domain[],
			   const char grid_bitant_plane[],
			   const char grid_quadrant_direction[],
			   const char grid_rotation_axis[],
			   fp origin_x, fp origin_y, fp origin_z)
{
if	(STRING_EQUAL(grid_domain, "full"))
   then return patch_system::patch_system__full_sphere;

else if ( STRING_EQUAL(grid_domain, "bitant")
	  && STRING_EQUAL(grid_bitant_plane, "xy") )
   then {
	if (origin_z == 0.0)
	   then return patch_system::patch_system__plus_z_hemisphere;
	   else {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "      Cactus has bitant mode about xy plane");
		CCTK_VInfo(CCTK_THORNSTRING,
			   "      but patch system origin_z=%g != 0",
			   double(origin_z));
		CCTK_VInfo(CCTK_THORNSTRING,
			   "      ==> using \"full sphere\" patch system");
		return patch_system::patch_system__full_sphere;
		}
	}

else if ( STRING_EQUAL(grid_domain, "quadrant")
	  && STRING_EQUAL(grid_quadrant_direction, "z") )
   then {
	if ((origin_x == 0.0) && (origin_y == 0.0))
	   then return patch_system::patch_system__plus_xy_quadrant_mirrored;
	   else {
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      Cactus has quadrant mode about z axis");
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      but patch system origin_(x,y)=(%g,%g) != (0,0)",
		   double(origin_x), double(origin_y));
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      ==> using \"full sphere\" patch system");
		return patch_system::patch_system__full_sphere;
		}
	}

else if ( STRING_EQUAL(grid_domain, "quadrant_reflect_rotate")
	  && STRING_EQUAL(grid_quadrant_direction, "z")
	  && STRING_EQUAL(grid_rotation_axis, "z") )
   then {
	if ((origin_x == 0.0) && (origin_y == 0.0))
	   then return patch_system::patch_system__plus_xy_quadrant_rotating;
	   else {
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      Cactus has rotating quadrant mode about z axis");
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      but patch system origin_(x,y)=(%g,%g) != (0,0)",
		   double(origin_x), double(origin_y));
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      ==> using \"full sphere\" patch system");
		return patch_system::patch_system__full_sphere;
		}
	}

else if (STRING_EQUAL(grid_domain, "octant"))
   then {
	if	((origin_x == 0.0) && (origin_y == 0.0) && (origin_z == 0.0))
	   then return patch_system::patch_system__plus_xyz_octant_mirrored;
	else if ((origin_x == 0.0) && (origin_y == 0.0))
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      Cactus has mirrored octant mode");
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      but patch system origin_z=%g != 0",
		   double(origin_z));
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      ==> using \"+xy quadrant (mirrored)\" patch system");
		return patch_system::patch_system__plus_xy_quadrant_mirrored;
		}
	else	{
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      Cactus has mirrored octant mode");
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      but patch system origin_(x,y,z)=(%g,%g,%g) != (0,0,0)",
		   double(origin_x), double(origin_y), double(origin_z));
		CCTK_VInfo(CCTK_THORNSTRING,
		   "      ==> using \"full sphere\" patch system");
		return patch_system::patch_system__full_sphere;
		}
	}

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 may or may not work..."
,
		   grid_domain);
	return patch_system::patch_system__full_sphere;
	}
}
	  }