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

//
// misc global data
// main
///
/// 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
/// gridfn_minus - compute gridfn x - gridfn y --> gridfn z
///
/// sym_fn_xyz - symmetrized test function fn(x,y,z) + fn(-y,x,z) + ...
/// fn_xyz - test function fn(x,y,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 "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
#include "jt/interpolate.hh"

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

using jtutil::error_exit;

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

//
// 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;
static const int which_deriv_all
   = which_deriv_fn
     | which_deriv_rho | which_deriv_sigma
     | which_deriv_rho_rho | which_deriv_rho_sigma | which_deriv_sigma_sigma;

// 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.1,
		deriv_weight_rho         = 4.1,
		deriv_weight_sigma       = 5.9,
		deriv_weight_rho_rho     = 2.6,
		deriv_weight_rho_sigma   = 5.3,
		deriv_weight_sigma_sigma = 5.8;

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

//
// function prototypes
//

namespace {
void setup_sym_fn_xyz(patch_system& ps, int gfn, bool want_ghost_zones);
void setup_fn_rho_sigma(patch_system& ps, int gfn, bool want_ghost_zones);
void finite_diff(patch_system& ps, int 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,
		  bool want_ghost_zones);

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 irho, int isigma,
		  int gfn_src, int which_derivs);
fp analytic_deriv_fn(fp rho, fp sigma, int which_derivs);
	  };

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

//
// main test driver -- see --help message below for usage
//
int main(int argc, const char *argv[])
{
const char *help_msg =
"\
Usage:\n\
  test_patch_system\n\
    { full-sphere | +z-hemisphere | +xy-quadrant | +xz-quadrant | octant }\n\
    N_ghost_points  N_overlap_points  delta_drho_dsigma\n\
    interpolator_order\n\
    { fn | ghost-zone\
 | deriv.{fn|rho|sigma|rho.rho|rho.sigma|sigma.sigma|all} }\n\
    [ fn | ghosted-fn | deriv-fd | deriv-analytic | error ]\n\
\n\
This program tests patch_system:: and its subsidiary classes.  After\n\
creating a patch system specified by the first 5 arguments, it does\n\
the test specified by the 6th argument:\n\
fn\n\
   Test gridfn storage, indexing, and coordinates:\n\
   * set up a test function on the ghosted grid with symmetry just matching\n\
     that of the patch system
   * [default] print the gridfn\n\
ghost-zone\n\
   Test extending a gridfn to ghost zones:\n\
   * set up test function on the nominal grid with symmetry just matching\n\
     that of the patch system\n\
   * call  patch_system::extend_scalar_gridfn_to_all_ghost_zones()\n\
   * set up the same test function on the ghosted grid\n\
   * compute error\n\
   * [default] print the error\n\
deriv.*\n\
   Test finite differencing:\n\
   * set up a test function (separately specified as fn(rho,sigma)\n\
     for each patch) on the ghosted grid\n\
   * compute specified linear combination of finite differences\n\
   * compute specified linear combination of true analytical derivatives\n\
   * compute error\n\
   * [default] print the error\n\
\n\
The optional last (7th) argument specifies which gridfn to print; for\n\
each test this has the default given above.\n\
";


// gridfn numbers
static const int gfn_min = -1;
static const int gfn_max = 2;
static const int gfn_fn = -1;			// common to all tests
static const int gfn_ghosted_fn = 0;		// for ghost-zone setup tests
static const int gfn_deriv_fd = 0;		// for finite diff tests
static const int gfn_deriv_analytic = 1;	// for finite diff tests
static const int gfn_error = 2;			// for finite diff tests


//
// ***** command line parsing *****
//

if ((argc == 2) && STRING_EQUAL(argv[1], "--help"))
   then {
	printf("%s", help_msg);
	return 0;						/*NOTREACHED*/
	}

if (! ((argc == 7) || (argc == 8)) )
   then error_exit(ERROR_EXIT, "%s", help_msg);			/*NOTREACHED*/

enum patch_system::patch_system_type type;
if	(STRING_EQUAL(argv[1], "full-sphere"))
   then type = patch_system::full_sphere_patch_system;
else if (STRING_EQUAL(argv[1], "+z-hemisphere"))
   then type = patch_system::plus_z_hemisphere_patch_system;
else if (STRING_EQUAL(argv[1], "+xy-quadrant"))
   then type = patch_system::plus_xy_quadrant_patch_system;
else if (STRING_EQUAL(argv[1], "+xz-quadrant"))
   then type = patch_system::plus_xz_quadrant_patch_system;
else if (STRING_EQUAL(argv[1], "octant"))
   then type = patch_system::octant_patch_system;
else	error_exit(ERROR_EXIT, "%s", help_msg);			/*NOTREACHED*/

int N_ghost_points, N_overlap_points;
fp delta_drho_dsigma;
int interpolator_order;
if (! (    (sscanf(argv[2], "%d", &N_ghost_points) == 1)
	&& (sscanf(argv[3], "%d", &N_overlap_points) == 1)
	&& (sscanf(argv[4], FP_SCANF_FORMAT, &delta_drho_dsigma) == 1)
	&& (sscanf(argv[5], "%d", &interpolator_order) == 1)    ) )
   then error_exit(ERROR_EXIT, "%s", help_msg);			/*NOTREACHED*/

enum {test_fn, test_ghost_zone, test_deriv} which_test;
int which_derivs;
if	(STRING_EQUAL(argv[6], "fn"))
   then which_test = test_fn;
else if (STRING_EQUAL(argv[6], "ghost-zone"))
   then which_test = test_ghost_zone;
else if (STRING_EQUAL(argv[6], "deriv.fn"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_fn;
	}
else if (STRING_EQUAL(argv[6], "deriv.rho"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_rho;
	}
else if (STRING_EQUAL(argv[6], "deriv.sigma"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_sigma;
	}
else if (STRING_EQUAL(argv[6], "deriv.rho.rho"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_rho_rho;
	}
else if (STRING_EQUAL(argv[6], "deriv.rho.sigma"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_rho_sigma;
	}
else if (STRING_EQUAL(argv[6], "deriv.sigma.sigma"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_sigma_sigma;
	}
else if (STRING_EQUAL(argv[6], "deriv.all"))
   then {
	which_test = test_deriv;
	which_derivs = which_deriv_all;
	}
else	error_exit(ERROR_EXIT, "%s", help_msg);		/*NOTREACHED*/

int gfn_to_print = (which_test == test_fn) ? gfn_fn : gfn_error;
if (argc >= 8)
   then {
	if	(STRING_EQUAL(argv[7], "fn"))
	   then gfn_to_print = gfn_fn;
	else if (STRING_EQUAL(argv[7], "ghosted-fn"))
	   then gfn_to_print = gfn_ghosted_fn;
	else if (STRING_EQUAL(argv[7], "deriv-fd"))
	   then gfn_to_print = gfn_deriv_fd;
	else if (STRING_EQUAL(argv[7], "deriv-analytic"))
	   then gfn_to_print = gfn_deriv_analytic;
	else if (STRING_EQUAL(argv[7], "error"))
	   then gfn_to_print = gfn_error;
	}

//
// ***** end of command line parsing *****
//


printf("##");
	for (int ap = 0 ; ap < argc ; ++ap)
	{
	printf(" %s", argv[ap]);
	}
printf("\n");

const fp origin_x = 0.314, origin_y = 0.159, origin_z = 0.265;

printf("## creating patch_system...\n");
patch_system ps(origin_x, origin_y, origin_z,
		type,
		N_ghost_points, N_overlap_points, delta_drho_dsigma,
		gfn_min, gfn_max,
		interpolator_order);

switch	(which_test)
	{
case test_fn:
	setup_sym_fn_xyz(ps, gfn_fn, true);
	ps.print_gridfn(gfn_to_print, true);
	break;
case test_ghost_zone:
	setup_sym_fn_xyz(ps, gfn_fn, false);
	ps.extend_scalar_gridfn_to_all_ghost_zones(gfn_fn);
	setup_sym_fn_xyz(ps, gfn_ghosted_fn, true);
	gridfn_minus(ps, gfn_fn, gfn_ghosted_fn, gfn_error, true);
	ps.print_gridfn(gfn_to_print, true);
	break;
case test_deriv:
	setup_fn_rho_sigma(ps, gfn_fn, true);	// fn(rho,sigma), ghost zones
	finite_diff(ps, gfn_fn, gfn_deriv_fd, which_derivs);
	analytic_derivs(ps, gfn_deriv_analytic, which_derivs);
	gridfn_minus(ps, gfn_deriv_fd, gfn_deriv_analytic, gfn_error, false);
	ps.print_gridfn(gfn_to_print, false);
	break;
default:
	error_exit(PANIC_EXIT,
"main(): impossible which_test=(int)%d!\n",
		   int(which_test));				/*NOTREACHED*/
	}
}

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

//
// This function sets up the test function for the function and
// ghost-zone tests, symmetrizing the test function to match the
// symmetry of the patch system.
//
// Arguments:
// ps = The patch system.
// gfn = Specifies the gridfn to set up.
// want_ghost_zones = true ==> Set up on ghosted grid
//		      false ==> Set up on nominal grid
//
namespace {
void setup_sym_fn_xyz(patch_system& ps, int gfn, bool want_ghost_zones)
{
printf("## setting up test fn(x,y,z) on %s grid...\n",
       (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.gridfn(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.
// gfn = Specifies the gridfn to set up.
// want_ghost_zones = true ==> Set up on ghosted grid
//		      false ==> Set up on nominal grid
//
namespace {
void setup_fn_rho_sigma(patch_system& ps, int gfn, bool want_ghost_zones)
{
printf("## setting up test fn(rho,sigma) on each patch of %s grid...\n",
       (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);

		p.gridfn(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 gfn_src, int gfn_dst,
		 int which_derivs)
{
printf(
"## finite differencing (gfn_src=%d, gfn_dst=%d, which_derivs=0x%02x)...\n",
       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, irho,isigma, gfn_src, which_derivs);
	}
	}

	}
}
	  }

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

//
// This function computes (on the nominal grid only) the specified
// linear combination of derivatives of the test function.
//
namespace {
void analytic_derivs(patch_system& ps,
		     int gfn_dst,
		     int which_derivs)
{
printf("## computing analytic derivatives(gfn_dst=%d, which_derivs=0x%02x...\n",
       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)
	{
	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  gridfn_x - gridfn_y --> gridfn_z
// on either the nominal or the full grid.
//
namespace {
void gridfn_minus(patch_system& ps, 
		  int gfn_x, int gfn_y, int gfn_dst,
		  bool want_ghost_zones)
{
printf("## gridfn_minus(gfn_x=%d, gfn_y=%d, gfn_dst=%d...\n",
       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.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)
		{
		p.gridfn(gfn_dst, irho,isigma)
			= p.gridfn(gfn_x, irho,isigma)
			  - p.gridfn(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::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(4*rho)) * tanh(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:
// gfn = Specifies the gridfn to set up.
//
namespace {
fp finite_diff_fn(const patch& p, int irho, int isigma,
		  int gfn_src, int which_derivs)
{
fp sum = 0.0;

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

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

if (which_derivs & which_deriv_rho_rho)
   then sum += deriv_weight_rho_rho
	       * p.partial_rho_rho(gfn_src, irho,isigma);
if (which_derivs & which_deriv_rho_sigma)
   then sum += deriv_weight_rho_sigma
	       * p.partial_rho_sigma(gfn_src, irho,isigma);
if (which_derivs & which_deriv_sigma_sigma)
   then sum += deriv_weight_sigma_sigma
	       * p.partial_sigma_sigma(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
	       * - cos(rho)*exp(sin(rho))*sinh(-1.0+pow(cos(sigma),2.0))
	           / cosh(-1.0+pow(cos(sigma),2.0));
if (which_derivs & which_deriv_sigma)
   then sum += deriv_weight_sigma
	       * 2.0*exp(sin(rho))*sin(sigma)*cos(sigma)
		 / pow(cosh(-1.0+pow(cos(sigma),2.0)),2.0);

if (which_derivs & which_deriv_rho_rho)
   then sum += deriv_weight_rho_rho
	       * exp(sin(rho))*sinh(-1.0+pow(cos(sigma),2.0))
		 * (sin(rho)-pow(cos(rho),2.0))
		 / cosh(-1.0+pow(cos(sigma),2.0));
if (which_derivs & which_deriv_rho_sigma)
   then sum += deriv_weight_rho_sigma
	       * 2.0*cos(rho)*exp(sin(rho))*sin(sigma)*cos(sigma)
		 / pow(cosh(-1.0+pow(cos(sigma),2.0)),2.0);
if (which_derivs & which_deriv_sigma_sigma)
   then sum += deriv_weight_sigma_sigma
	       * -2.0*exp(sin(rho))
		  * ( - 4.0*sinh(-1.0+pow(cos(sigma),2.0))*pow(cos(sigma),2.0)
		      + 4.0*sinh(-1.0+pow(cos(sigma),2.0))*pow(cos(sigma),4.0)
		      - 2.0*pow(cos(sigma),2.0)*cosh(-1.0+pow(cos(sigma),2.0))
		      + cosh(-1.0+pow(cos(sigma),2.0)) )
		  / pow(cosh(-1.0+pow(cos(sigma),2.0)),3.0);

return sum;
}
	  };