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

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

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

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

//
// ***** prototypes for functions local to this file
//
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,
		  struct solver_info& solver_info,
		  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;

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

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

// if Cactus is using a conformal metric,
//    we need to re-fetch the data pointer for the conformal factor,
//    since the data array may have moved in memory since the last
//    time we were called due to storage being switched off/on
if (state.cgi.Cactus_conformal_metric)
   then state.cgi.psi_data = Cactus_gridfn_data_ptr(cctkGH,
						    "StaticConformal::psi");

state.IO_info.time_iteration = cctk_iteration;
	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;

	//
	// The first time we (try to) find a given horizon, our initial
	// guess is likely to be rather inaccurate, so we may need a
	// larger number of iterations.  But if we've found this horizon
	// before, then we have its previous position as an initial guess,
	// so we shouldn't need as many iterations.
	//
	state.solver_info.max_Newton_iterations
		= AH_info.AH_found ? max_Newton_iterations__subsequent
				   : max_Newton_iterations__initial;

	AH_info.AH_found
	   = find_horizon(state.method,
			  verbose_info, state.timer_handle,
			  state.IO_info, state.Jac_info, state.solver_info,
			  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 finds (or more accurately tries to find) a single
// apparent horizon.
//
// 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)
// 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,
		  struct solver_info& solver_info,
		  struct cactus_grid_info& cgi, struct geometry_info& gi,
		  patch_system& ps, Jacobian* Jac_ptr,
		  int hn, int N_horizons)
{
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 iv 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 (IO_info.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, IO_info,
			       hn, verbose_info);
	if (timer_handle >= 0)
	   then CCTK_TimerStopI(timer_handle);
	if (! status)
	   then return false;				// *** ERROR RETURN ***

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