// test_patch_system.cc -- test driver for patch_system::
// $Id$

//
// <<<misc global data>>>
// <<<prototypes>>>
//
// driver - Cactus interface
//
/// setup_sym_fn_xyz - set up symmetrized test function fn(global_[xyz])
/// setup_fn_rho_sigma - set up test function fn(rho,sigma)
/// finite_diff - compute linear combination of finite differences
/// analytic_derivs - compute linear combination of analytic derivatives
///
/// sym_fn_xyz - symmetrized test function fn(x,y,z) + fn(-y,x,z) + ...
/// fn_xyz - test function fn(x,y,z)
///
/// gridfn_minus - compute gridfn x - gridfn y --> gridfn z
/// ghosted_gridfn_minus - compute [ghosted] gridfn x - gridfn y --> gridfn z
///
/// fn_rho_sigma - test function fn(rho,sigma)
/// finite_diff_fn - finite differences of fn(rho,sigma)
/// analytic_deriv_fn - analytical derivs of fn(rho,sigma) (via Maple)
//

#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 "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;

#include "fp.hh"
#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "patch_interp.hh"
#include "ghost_zone.hh"
#include "patch_system.hh"

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

//
// misc global data
//

// which test are we going to do?
static const int which_deriv_fn          = 0x1;
static const int which_deriv_rho         = 0x2;
static const int which_deriv_sigma       = 0x4;
static const int which_deriv_rho_rho     = 0x8;
static const int which_deriv_rho_sigma   = 0x10;
static const int which_deriv_sigma_sigma = 0x20;

// when testing multiple derivatives, we combine them together
// with these "quasi-random" weights (chosen from digits of pi)
static const fp deriv_weight_fn          = 3.14,
		deriv_weight_rho         = 1.59,
		deriv_weight_sigma       = 2.65,
		deriv_weight_rho_rho     = 3.58,
		deriv_weight_rho_sigma   = 9.79,
		deriv_weight_sigma_sigma = 3.23;

// n.b. nominal gfns must all be < 0, ghosted > 0
// nominal gridfns
static const int FD_derivs_gfn = -1;
static const int analytic_derivs_gfn = -2;
static const int nominal_error_gfn = -3;
static const int nominal_min_gfn = -3;
static const int nominal_max_gfn = -1;

// ghosted gridfns
static const int test_fn_gfn = 1;
static const int test_fn_copy_gfn = 2;
static const int ghosted_error_gfn = 3;
static const int ghosted_min_gfn = 1;
static const int ghosted_max_gfn = 3;

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

//
// ***** prototypes *****
//

extern "C"
  void test_patch_system(CCTK_ARGUMENTS);

namespace {
void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones);
void setup_fn_rho_sigma(patch_system& ps, int ghosted_gfn);
void finite_diff(patch_system& ps,
		 int ghosted_gfn_src, int gfn_dst,
		 int which_derivs);
void analytic_derivs(patch_system& ps, int gfn_dst, int which_derivs);
void gridfn_minus(patch_system& ps, 
		  int gfn_x, int gfn_y, int gfn_dst);
void ghosted_gridfn_minus(patch_system& ps, 
			  int ghosted_gfn_x, int ghosted_gfn_y,
			  int ghosted_gfn_dst);

fp sym_fn_xyz(enum patch_system::patch_system_type type, fp x, fp y, fp z);
fp fn_xyz(fp x, fp y, fp z);

fp fn_rho_sigma(fp rho, fp sigma);
fp finite_diff_fn(const patch& p,
		  int ghosted_gfn_src, int which_derivs,
		  int irho, int isigma);
fp analytic_deriv_fn(fp rho, fp sigma, int which_derivs);
	  };

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

//
// This function is the Cactus interface for the test driver.
//
extern "C"
  void test_patch_system(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS

//
// set up the interpatch interpolator
//
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,
		   "couldn't find interpolator \"%s\"!",
		   interpatch_interpolator_name);		/*NOTREACHED*/
const int interp_par_table_handle
	= Util_TableCreateFromString(interpatch_interpolator_pars);
if (interp_par_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "bad interpatch-interpolator parameter(s) \"%s\"!",
		   interpatch_interpolator_pars);		/*NOTREACHED*/

//
// create the patch system
//
CCTK_VInfo(CCTK_THORNSTRING, "about to create patch system");
patch_system ps(origin_x, origin_y, origin_z,
		patch_system::type_of_name(patch_system_type),
		N_ghost_points, N_overlap_points, delta_drho_dsigma,
		nominal_min_gfn, nominal_max_gfn,
		ghosted_min_gfn, ghosted_max_gfn,
		interp_handle, interp_par_table_handle);
CCTK_VInfo(CCTK_THORNSTRING, "patch system created ok");

//
// do the actual tests
//
if      (STRING_EQUAL(which_test, "gridfn"))
   then {
	setup_sym_fn_xyz(ps, test_fn_gfn, true);
	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn");
	}
else if (STRING_EQUAL(which_test, "synchronize"))
   then {
	setup_sym_fn_xyz(ps, test_fn_gfn, false);
	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_init.dat");
	setup_sym_fn_xyz(ps, test_fn_copy_gfn, true);
	ps.print_ghosted_gridfn(test_fn_copy_gfn, "test_fn_copy.dat");
	ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn);
	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_sync.dat");
	ghosted_gridfn_minus(ps,
			     test_fn_gfn, test_fn_copy_gfn,
			     ghosted_error_gfn);
	ps.print_ghosted_gridfn(ghosted_error_gfn, "ghosted_error.dat");
	}
else if (STRING_EQUAL(which_test, "derivatives"))
   then {
	setup_fn_rho_sigma(ps, test_fn_gfn);
	ps.print_gridfn(nominal_error_gfn, "test_fn.dat");
	finite_diff(ps, test_fn_gfn, FD_derivs_gfn, which_derivs);
	ps.print_gridfn(FD_derivs_gfn, "FD_derivs.dat");
	analytic_derivs(ps, analytic_derivs_gfn, which_derivs);
	ps.print_gridfn(analytic_derivs_gfn, "analytic_derivs.dat");
	gridfn_minus(ps,
		     FD_derivs_gfn, analytic_derivs_gfn,
		     nominal_error_gfn);
	ps.print_gridfn(nominal_error_gfn, "nominal_error.dat");
	}
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown which_test=\"%s\"!",
		   which_test);					/*NOTREACHED*/

CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system");
}

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

//
// This function sets up the ghosted test function for the gridfn and
// synchronize tests, symmetrizing the test function to match the
// symmetry of the patch system.
//
// Arguments:
// ps = The patch system.
// ghosted_gfn = Specifies the gridfn to set up.
// want_ghost_zones = true ==> Set up on ghosted part of ghosted grid
//		      false ==> Set up on nominal part of ghosted grid
//
namespace {
void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones)
{
CCTK_VInfo(CCTK_THORNSTRING, "setting up ghosted test fn(x,y,z)");
CCTK_VInfo(CCTK_THORNSTRING, "   on %s part of ghosted grid",
	   want_ghost_zones ? "ghosted" : "nominal");

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

		for (int irho = p.effective_min_irho(want_ghost_zones) ;
		     irho <= p.effective_max_irho(want_ghost_zones) ;
		     ++irho)
		{
		for (int isigma = p.effective_min_isigma(want_ghost_zones) ;
		     isigma <= p.effective_max_isigma(want_ghost_zones) ;
		     ++isigma)
		{
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);

		fp local_x, local_y, local_z;
		p.xyz_of_r_rho_sigma(1.0, rho, sigma,
				     local_x, local_y, local_z);

		p.ghosted_gridfn(ghosted_gfn, irho,isigma)
			= sym_fn_xyz(ps.type(), local_x, local_y, local_z);
		}
		}
	}
}
	  }

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

//
// This function sets up the test function for the finite differencing
// tests.
//
// Arguments:
// ps = The patch system.
// ghosted_gfn = Specifies the gridfn to set up.
//
namespace {
void setup_fn_rho_sigma(patch_system& ps, int ghosted_gfn)
{
CCTK_VInfo(CCTK_THORNSTRING, "setting up ghosted test fn(rho,sigma)");

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

		for (int irho = p.ghosted_min_irho() ;
		     irho <= p.ghosted_max_irho() ;
		     ++irho)
		{
		for (int isigma = p.ghosted_min_isigma() ;
		     isigma <= p.ghosted_max_isigma() ;
		     ++isigma)
		{
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);

		p.ghosted_gridfn(ghosted_gfn, irho,isigma)
			= fn_rho_sigma(rho, sigma);
		}
		}
	}
}
	  }

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

//
// This function computes (on the nominal grid only) the specified
// linear combination of finite derivatives of a test function.
//
// Arguments:
// ps = The patch system.
// gfn_src = Specifies the gridfn to finite difference.
// gfn_dst = Specifies the gridfn in which to store the result.
// which = Specifies which finite derivatives to include in the test.
//
namespace {
void finite_diff(patch_system& ps,
		 int ghosted_gfn_src, int gfn_dst,
		 int which_derivs)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "finite differencing ghosted_gfn=%d --> gfn=%d (which_derivs=%d)",
	   ghosted_gfn_src, gfn_dst, which_derivs);

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

	for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
	{
	for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; ++isigma)
	{
	p.gridfn(gfn_dst, irho,isigma)
		= finite_diff_fn(p,
				 ghosted_gfn_src, which_derivs,
				 irho,isigma);
	}
	}

	}
}
	  }

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

//
// This function computes the specified linear combination of
// derivatives of the test function, everywhere on the nominal grid.
//
namespace {
void analytic_derivs(patch_system& ps,
		     int gfn_dst,
		     int which_derivs)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "computing analytic derivatives (which_derivs=%d)",
	   which_derivs);

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

	for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
	{
	for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; ++isigma)
	{
	const fp rho = p.rho_of_irho(irho);
	const fp sigma = p.sigma_of_isigma(isigma);
	p.gridfn(gfn_dst, irho,isigma)
		= analytic_deriv_fn(rho,sigma, which_derivs);
	}
	}

	}
}
	  }

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

//
// This function computes [nominal]  gridfn_x - gridfn_y --> gridfn_z .
//
namespace {
void gridfn_minus(patch_system& ps, 
		  int gfn_x, int gfn_y, int gfn_dst)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "[nominal] gfn=%d - gfn=%d --> gfn=%d",
	   gfn_x, gfn_y, gfn_dst);

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

	for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
	{
	for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; ++isigma)
	{
	p.gridfn(gfn_dst, irho,isigma)
		=   p.gridfn(gfn_x, irho,isigma)
		  - p.gridfn(gfn_y, irho,isigma);
	}
	}

	}
}
	  }

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

//
// This function computes [ghosted]  gridfn_x - gridfn_y --> gridfn_z .
//
namespace {
void ghosted_gridfn_minus(patch_system& ps, 
			  int ghosted_gfn_x, int ghosted_gfn_y,
			  int ghosted_gfn_dst)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "[ghosted] gfn=%d - gfn=%d --> gfn=%d",
	   ghosted_gfn_x, ghosted_gfn_y, ghosted_gfn_dst);

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

		for (int irho = p.ghosted_min_irho() ;
		     irho <= p.ghosted_max_irho() ;
		     ++irho)
		{
		for (int isigma = p.ghosted_min_isigma() ;
		     isigma <= p.ghosted_max_isigma() ;
		     ++isigma)
		{
		p.ghosted_gridfn(ghosted_gfn_dst, irho,isigma)
			=   p.ghosted_gridfn(ghosted_gfn_x, irho,isigma)
			  - p.ghosted_gridfn(ghosted_gfn_y, irho,isigma);
		}
		}

	}
}
	  }

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

//
// This function symmetrizes  fn_xyz()  (about the origin) to match
// the patch system's symmetries.
//
// To rotate f(x,y) by 90, 180, or 270 degrees:
//
//	      (-y,x) |
//	             |
//	             |
//	             |
//	             |        (x,y)
//	             |
//	-------------+-------------
//	             |
//	(-x,-y)      |
//	             |
//	             |
//	             |
//	             | (y,-x)
//
namespace {
fp sym_fn_xyz(enum patch_system::patch_system_type type, fp x, fp y, fp z)
{
switch	(type)
	{
case patch_system::full_sphere_patch_system:
	return fn_xyz(x,y,z);
	break;
case patch_system::plus_z_hemisphere_patch_system:
	return fn_xyz(x,y,+z) + fn_xyz(x,y,-z);
	break;
case patch_system::plus_xy_quadrant_patch_system:
	return   fn_xyz(+x,+y,z)
	       + fn_xyz(-y,+x,z)
	       + fn_xyz(-x,-y,z)
	       + fn_xyz(+y,-x,z);
	break;
case patch_system::plus_xz_quadrant_patch_system:
	return   fn_xyz(+x,+y,+z) + fn_xyz(-x,-y,+z)
	       + fn_xyz(+x,+y,-z) + fn_xyz(-x,-y,-z);
	break;
case patch_system::plus_xyz_octant_patch_system:
	return   fn_xyz(+x,+y,+z) + fn_xyz(+x,+y,-z)
	       + fn_xyz(-y,+x,+z) + fn_xyz(-y,+x,-z)
	       + fn_xyz(-x,-y,+z) + fn_xyz(-x,-y,-z)
	       + fn_xyz(+y,-x,+z) + fn_xyz(+y,-x,-z);
	break;
default:
	error_exit(PANIC_EXIT,
"***** sym_fn_xyz(): impossible type=(int)%d!\n",
		   int(type));					/*NOTREACHED*/
			}
}
	  };

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

//
// This is the underlying test function for our function and ghost-zone tests.
//
namespace {
fp fn_xyz(fp x, fp y, fp z)
{
return (x*(x+0.238) + 2.417*y*(y-0.917) + 1.38*z*(z-0.472))
       * tanh(jtutil::pow3(cos(z)));
}
	  };

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

//
// This is the underlying test function for our finite differencing
// tests.
//
namespace {
fp fn_rho_sigma(fp rho, fp sigma)
{
return exp(sin(1.38*rho)) * tanh(0.17+0.83*jtutil::pow2(sin(sigma)));
}
	  };

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

//
// This function computes the sum of our linear combination of various
// finite difference approximations to derivatives of fn_rho_sigma().
//
// Arguments:
// ghosted_gfn_src = Specifies the gridfn to finite difference.
//
namespace {
fp finite_diff_fn(const patch& p,
		  int ghosted_gfn_src, int which_derivs,
		  int irho, int isigma)
{
fp sum = 0.0;

if (which_derivs & which_deriv_fn)
   then sum += deriv_weight_fn
	       * p.ghosted_gridfn(ghosted_gfn_src, irho,isigma);

if (which_derivs & which_deriv_rho)
   then sum += deriv_weight_rho
	       * p.partial_rho(ghosted_gfn_src, irho,isigma);
if (which_derivs & which_deriv_sigma)
   then sum += deriv_weight_sigma
	       * p.partial_sigma(ghosted_gfn_src, irho,isigma);

if (which_derivs & which_deriv_rho_rho)
   then sum += deriv_weight_rho_rho
	       * p.partial_rho_rho(ghosted_gfn_src, irho,isigma);
if (which_derivs & which_deriv_rho_sigma)
   then sum += deriv_weight_rho_sigma
	       * p.partial_rho_sigma(ghosted_gfn_src, irho,isigma);
if (which_derivs & which_deriv_sigma_sigma)
   then sum += deriv_weight_sigma_sigma
	       * p.partial_sigma_sigma(ghosted_gfn_src, irho,isigma);

return sum;
}
	  };

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

//
// This function computes the sum of a specified linear combination of
// various (analytical) derivatives of fn_rho_sigma().
//
// The derivatives were machine-generated via Maple's codegen[C]() function:
// "deriv_patch_system.maple" is the Maple input
// "deriv_patch_system.out" is the Maple input; code here is cut-n-pasted
//			   from the Maple codegen[C]() output there
//
namespace {
fp analytic_deriv_fn(fp rho, fp sigma, int which_derivs)
{
fp sum = 0.0;

if (which_derivs & which_deriv_fn)
   then sum += deriv_weight_fn
	       * fn_rho_sigma(rho,sigma);

if (which_derivs & which_deriv_rho)
   then sum += deriv_weight_rho
	       * 0.138E1*cos(0.138E1*rho)*exp(sin(0.138E1*rho))
		 *tanh(0.17+0.83*pow(sin(sigma),2.0));
if (which_derivs & which_deriv_sigma)
   then sum += deriv_weight_sigma
	       * 0.166E1*exp(sin(0.138E1*rho))
		 *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
		 *sin(sigma)*cos(sigma);

if (which_derivs & which_deriv_rho_rho)
   then sum += deriv_weight_rho_rho
	       * (
		 -0.19044E1*sin(0.138E1*rho)*exp(sin(0.138E1*rho))
		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
		 +0.19044E1*pow(cos(0.138E1*rho),2.0)*exp(sin(0.138E1*rho))
		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
		 );
if (which_derivs & which_deriv_rho_sigma)
   then sum += deriv_weight_rho_sigma
	       * 0.22908E1*cos(0.138E1*rho)*exp(sin(0.138E1*rho))
		 *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
		 *sin(sigma)*cos(sigma);
if (which_derivs & which_deriv_sigma_sigma)
   then sum += deriv_weight_sigma_sigma
	       * (
		 -0.55112E1*exp(sin(0.138E1*rho))
		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
		  *pow(sin(sigma),2.0)*pow(cos(sigma),2.0)
		 +0.166E1*exp(sin(0.138E1*rho))
		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
		  *pow(cos(sigma),2.0)
		 -0.166E1*exp(sin(0.138E1*rho))
		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
		  *pow(sin(sigma),2.0)
		 );

return sum;
}
	  };