path: root/src/patch/patch.cc

// patch.cc -- describes a coordinate/grid patch
// $Id$

//
// patch::patch
// patch::~patch
// z_patch::z_patch
// x_patch::x_patch
// y_patch::y_patch
//
// patch::minmax_ang_ghost_zone
// patch::ghost_zone_on_edge
// patch::ghost_zone_ptr
// patch::patch_frontier_ptr
// patch::setup_mirror_symmetry_ghost_zone
// patch::setup_periodic_symmetry_ghost_zone
// patch::setup_interpatch_ghost_zone
// patch::set_ghost_zone
// patch::minmax_ang_ghost_zone
// patch::ghost_zone_on_edge
// patch::frontier_ptr_on_edge
// patch::edge_adjacent_to_patch
// patch::interpatch_ghost_zone_on_edge
// patch::assert_all_ghost_zones_fully_setup
//
// patch::print_gridfn
//
// patch_info::grid_array_pars
// patch_info::grid_pars
//

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

using jtutil::error_exit;

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

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

//
// This function constructs a  patch  object.
//
patch::patch(patch_system &my_patch_system_in, int patch_number_in,
	     const char *name_in, bool is_plus_in, char ctype_in,
	     local_coords::coords_set coords_set_rho_in,
	     local_coords::coords_set coords_set_sigma_in,
	     local_coords::coords_set coords_set_tau_in,
	     const grid_arrays::grid_array_pars& grid_array_pars_in,
	     const grid::grid_pars& grid_pars_in)

	: fd_grid(grid_array_pars_in, grid_pars_in),

	  my_patch_system_(my_patch_system_in), patch_number_(patch_number_in),
	  name_(name_in), is_plus_(is_plus_in), ctype_(ctype_in),

	  coords_set_rho_  (coords_set_rho_in  ),
	  coords_set_sigma_(coords_set_sigma_in),
	  coords_set_tau_  (coords_set_tau_in  ),

	  min_rho_patch_edge_(*new patch_edge(*this, side_is_min, side_is_rho)),
	  max_rho_patch_edge_(*new patch_edge(*this, side_is_max, side_is_rho)),
	  min_sigma_patch_edge_(*new patch_edge(*this,
						side_is_min, side_is_sigma)),
	  max_sigma_patch_edge_(*new patch_edge(*this,
						side_is_max, side_is_sigma)),

	  min_rho_ghost_zone_(NULL),
	  max_rho_ghost_zone_(NULL),
	  min_sigma_ghost_zone_(NULL),
	  max_sigma_ghost_zone_(NULL),

	  min_rho_patch_frontier_(NULL),
	  max_rho_patch_frontier_(NULL),
	  min_sigma_patch_frontier_(NULL),
	  max_sigma_patch_frontier_(NULL) // no comma

{ }

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

//
// This function destroys a  patch  object.
//
patch::~patch()
{
// no need to check for null pointers, since  delete NULL  is a silent no-op

delete max_sigma_ghost_zone_;
delete min_sigma_ghost_zone_;
delete max_rho_ghost_zone_;
delete min_rho_ghost_zone_;

delete max_sigma_patch_frontier_;
delete min_sigma_patch_frontier_;
delete max_rho_patch_frontier_;
delete min_rho_patch_frontier_;

delete & max_sigma_patch_edge_;
delete & min_sigma_patch_edge_;
delete & max_rho_patch_edge_;
delete & min_rho_patch_edge_;
}

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

//
// This function constructs a  z_patch  object.
//
z_patch::z_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'z',
		local_coords::set_mu, local_coords::set_nu,
		local_coords::set_phi,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function constructs an  x_patch  object.
//
x_patch::x_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'x',
		local_coords::set_nu, local_coords::set_phi,
		local_coords::set_mu,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function constructs a  y_patch  object.
//
y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'y',
		local_coords::set_mu, local_coords::set_phi,
		local_coords::set_nu,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function returns a reference to the specified ghost zone
// of this patch.
//
ghost_zone& patch::minmax_ang_ghost_zone(bool want_min, bool want_rho)
	const
{
return want_min ? (want_rho ? min_rho_ghost_zone()
			    : min_sigma_ghost_zone())
		: (want_rho ? max_rho_ghost_zone()
			    : max_sigma_ghost_zone());
}

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

//
// This function returns a reference to the specified ghost zone
// of this patch.
//
ghost_zone& patch::ghost_zone_on_edge(const patch_edge &edge)
const
{
assert(& edge.my_patch() == this);
return minmax_ang_ghost_zone(edge.is_min(), edge.is_rho());
}

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

//
// This function returns a reference to the specified ghost zone
// pointer data member of this patch, asserting that the current
// value of this pointer is NULL.  This is used (only) by
//	patch::setup_*_ghost_zone()
// to set the ghost zone pointers; the assert() ensures that it can't
// be used to overwrite a non-NULL pointer.
//
ghost_zone*& patch::ghost_zone_ptr(const patch_edge& edge)
{
assert(& edge.my_patch() == this);

ghost_zone*& gzp = edge.is_min() ? (edge.is_rho() ? min_rho_ghost_zone_
						  : min_sigma_ghost_zone_)
				 : (edge.is_rho() ? max_rho_ghost_zone_
						  : max_sigma_ghost_zone_);
assert(gzp == NULL);
return gzp;
}

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

//
// This function returns a reference to the specified patch frontier
// pointer data member of this patch, asserting that the current
// value of this pointer is NULL.  This is used (only) by the other
// patch's
//	interpatch_ghost_zone::setup_other_frontier()
// to set this patch's patch frontier pointers; the assert() ensures that
// it can't be used to overwrite a non-NULL pointer.
//
patch_frontier*& patch::patch_frontier_ptr(const patch_edge& edge)
{
assert(& edge.my_patch() == this);

patch_frontier*& pfp
	= edge.is_min() ? (edge.is_rho() ? min_rho_patch_frontier_
					 : min_sigma_patch_frontier_)
			: (edge.is_rho() ? max_rho_patch_frontier_
					 : max_sigma_patch_frontier_);
assert(pfp == NULL);
return pfp;
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then sets it up as a mirror-symmetry
// ghost zone.  It returns a reference to the newly-set-up symmetry
// ghost zone object.
//
symmetry_ghost_zone& patch::setup_mirror_symmetry_ghost_zone
	(const patch_edge& my_edge)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

symmetry_ghost_zone *temp = new symmetry_ghost_zone(my_edge);
set_ghost_zone(temp, my_edge);
return *temp;
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then sets it up as a periodic-symmetry
// ghost zone.  It returns a reference to the newly-set-up symmetry
// ghost zone object.
//
symmetry_ghost_zone& patch::setup_periodic_symmetry_ghost_zone
	(const patch_edge& my_edge, const patch_edge& symmetry_edge,
	 bool is_ipar_map_plus)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

int my_sample_ipar = my_edge.min_ipar_without_corners();
int symmetry_sample_ipar
	= is_ipar_map_plus ? symmetry_edge.min_ipar_without_corners()
			   : symmetry_edge.max_ipar_without_corners();

symmetry_ghost_zone *temp
	= new symmetry_ghost_zone(my_edge,        symmetry_edge,
				  my_sample_ipar, symmetry_sample_ipar,
				  is_ipar_map_plus);
ghost_zone_ptr(my_edge) = temp;
return *temp;
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then sets it up as an interpatch ghost
// zone.  It returns a reference to the newly-set-up interpatch ghost
// zone object.
//
interpatch_ghost_zone& patch::setup_interpatch_ghost_zone
	(const patch_edge& my_edge, const patch_edge& other_edge,
	 int N_overlap_points)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

interpatch_ghost_zone *temp
	= new interpatch_ghost_zone(my_edge, other_edge,
				    N_overlap_points);
ghost_zone_ptr(my_edge) = temp;
return *temp;
}

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

//
// This function finds which patch edge is adjacent to a neighboring
// patch q, or does an error_exit() if q isn't actually a neighboring patch.
// The computation is done using only (rho,sigma) coordinate sets and
// min/max dang bounds ==> it's ok to use this function in setting up
// interpatch ghost zones.
//
// Arguments:
// q = The (supposedly) neighboring patch.
// N_overlap_points = The number of grid points these patches overlap.
//		      If this is nonzero, then these patches must have the
//		      same grid spacing in the perpendicular direction.
//
const patch_edge& patch::edge_adjacent_to_patch(const patch& q,
						int N_overlap_points = 0)
	const
{
const patch& p = *this;

// which (rho,sigma) coordinate do the patches have in common?
// ... this is the perp coordinate for the border
const local_coords::coords_set common_coord_set
	= p.coords_set_rho_sigma() & q.coords_set_rho_sigma();

// is this coordinate rho or sigma in each patch?
const bool common_is_p_rho   = common_coord_set == p.coords_set_rho  ();
const bool common_is_p_sigma = common_coord_set == p.coords_set_sigma();
if ((common_is_p_rho ^ common_is_p_sigma) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of p.{rho,sigma}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_p_sigma=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_p_sigma));
								/*NOTREACHED*/
const bool common_is_q_rho   = common_coord_set == q.coords_set_rho  ();
const bool common_is_q_sigma = common_coord_set == q.coords_set_sigma();
if ((common_is_q_rho ^ common_is_q_sigma) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of q.{rho,sigma}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_q_rho=%d common_is_q_sigma=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_q_rho), int(common_is_q_sigma));
								/*NOTREACHED*/

// how much do the patches overlap?
if ((N_overlap_points != 0) && fuzzy<fp>::NE(p.delta_dang(common_is_p_rho),
					     q.delta_dang(common_is_q_rho)))
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        N_overlap_points != 0 must have same perp grid spacing in both patches!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_q_rho=%d\n"
"        p.delta_dang(common_is_p_rho)=%g\n"
"        q.delta_dang(common_is_q_rho)=%g\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_q_rho),
		   double(p.delta_dang(common_is_p_rho)),
		   double(q.delta_dang(common_is_q_rho)));	/*NOTREACHED*/
const fp doverlap = fp(N_overlap_points) * p.delta_dang(common_is_p_rho);

// where is the common boundary relative to the min/max sides of each patch?
const bool common_is_p_min_q_max
    = local_coords::fuzzy_EQ_dang(p.min_dang(common_is_p_rho),
				  q.max_dang(common_is_q_rho) - doverlap);
const bool common_is_p_max_q_min
    = local_coords::fuzzy_EQ_dang(p.max_dang(common_is_p_rho),
				  q.min_dang(common_is_q_rho) + doverlap);
if ((common_is_p_min_q_max ^ common_is_p_max_q_min) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of {pmax/qmin, pmin/qmax}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_q_rho=%d\n"
"        p.delta_dang(common_is_p_rho)=%g\n"
"        q.delta_dang(common_is_q_rho)=%g\n"
"        N_overlap_points=%d doverlap=%g\n"
"        common_is_p_min_q_max=%d common_is_p_max_q_min=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_q_rho),
		   double(p.delta_dang(common_is_p_rho)),
		   double(q.delta_dang(common_is_q_rho)),
		   N_overlap_points, double(doverlap),
		   int(common_is_p_min_q_max), int(common_is_p_max_q_min));
								/*NOTREACHED*/

return p.minmax_ang_patch_edge(common_is_p_min_q_max, common_is_p_rho);
}

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

//
// This function verifies that the ghost zone on a specified edge
// is indeed interpatch, and returns a reference to it as an
//  interpatch_ghost_zone .
//
interpatch_ghost_zone& patch::interpatch_ghost_zone_on_edge
	(const patch_edge &e)
	const
{
ghost_zone &g = ghost_zone_on_edge(e);
assert(g.is_interpatch());
return static_cast<interpatch_ghost_zone &>(g);
}

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

//
// This function verifies (via assert()) that all ghost zones of this
// patch, and all frontiers to which they should point, have been fully
// set up.
//
void patch::assert_all_ghost_zones_fully_setup() const
{
assert(min_rho_ghost_zone_ != NULL);
assert(max_rho_ghost_zone_ != NULL);
assert(min_sigma_ghost_zone_ != NULL);
assert(max_sigma_ghost_zone_ != NULL);

// these calls are no-ops for non-interpatch ghost zones
min_rho_ghost_zone().assert_fully_setup();
max_rho_ghost_zone().assert_fully_setup();
min_sigma_ghost_zone().assert_fully_setup();
max_sigma_ghost_zone().assert_fully_setup();
}

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

//
// This function prints a gridfn in ASCII format to an already-open
// stdio stream.  The output format is suitable for a gnuplot 'splot'
// command.  Individual patches may be selected with the  select.patch
// program (perl script).
//
void patch::print_gridfn(int gfn, bool want_ghost_zones,
			 FILE *output_fp = stdout)
	const
{
fprintf(output_fp, "### %s patch\n", name());
fprintf(output_fp, "# gfn=%d\n", gfn);
fprintf(output_fp, "# dpx = %s\n", name_of_dpx());
fprintf(output_fp, "# dpy = %s\n", name_of_dpy());
fprintf(output_fp, "#\n");
fprintf(output_fp, "# dpx\tdpy\tgridfn\tirho\tisigma\n");
	for (int irho = effective_min_irho(want_ghost_zones) ;
	     irho <= effective_max_irho(want_ghost_zones) ;
	     ++irho)
	{
		for (int isigma = effective_min_isigma(want_ghost_zones) ;
		     isigma <= effective_max_isigma(want_ghost_zones) ;
		     ++isigma)
		{
		const fp rho = rho_of_irho(irho);
		const fp sigma = sigma_of_isigma(isigma);
		const fp dpx = dpx_of_rho_sigma(rho, sigma);
		const fp dpy = dpy_of_rho_sigma(rho, sigma);
		fprintf(output_fp,
			"%g\t%g\t%.15g\t%d\t%d\n",
			dpx, dpy, gridfn(gfn,irho,isigma), irho, isigma);
		}
	printf("\n");
	}
}

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

//
// This function computes, and returns a reference to, a
//  struct grid_arrays::grid_array_pars  from the info in a
//  struct patch_info  and the additional information in the arguments.
//
// The result refers to an internal static buffer in this function; the
// usual caveats about lifetimes/overwriting apply.
//
// Arguments:
// N_ghost_points = Width in grid points of all ghost zones.
// N_extend_points = Number of grid points to extend each patch past
//		     "just touching" so as to overlap neighboring patches.
//		     Thus patches overlap by
//			N_overlap_points = 2*N_extend_points + 1
//		     grid points.  For example, with N_extend_points == 2,
//		     here are the grid points of two neighboring patches:
//			x   x   x   x   x   X   X
//                                      |
//			        O   O   o   o   o   o   o
//		     Here | marks the "just touching" boundary,
//		     x and o the grid points before this extension,
//		     and X and O the extra grid points added by this
//		     extension.
// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
//
const grid_arrays::grid_array_pars& patch_info::grid_array_pars
	(int N_ghost_points,
	 int N_extend_points,
	 fp delta_drho_dsigma)
	const
{
static struct grid_arrays::grid_array_pars grid_array_pars_buffer;

grid_array_pars_buffer.min_irho
	= jtutil::round<fp>::to_integer(min_drho  /delta_drho_dsigma);
grid_array_pars_buffer.min_isigma
	= jtutil::round<fp>::to_integer(min_dsigma/delta_drho_dsigma);
grid_array_pars_buffer.max_irho
	= array_pars.min_irho
	  + jtutil::round<fp>::to_integer(
		(max_drho  -min_drho  ) / delta_drho_dsigma
					 );
grid_array_pars_buffer.max_isigma
	= array_pars.min_isigma
	  + jtutil::round<fp>::to_integer(
		(max_dsigma-min_dsigma) / delta_drho_dsigma
					 );
grid_array_pars_buffer.min_irho   -= N_extend_points;
grid_array_pars_buffer.min_isigma -= N_extend_points;
grid_array_pars_buffer.max_irho   += N_extend_points;
grid_array_pars_buffer.max_isigma += N_extend_points;

grid_array_pars_buffer.min_rho_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.max_rho_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.min_sigma_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.max_sigma_N_ghost_points = N_ghost_points;

return grid_array_pars_buffer;
}
	  }

//******************************************************************************
//
//
// This function computes, and returns a reference to, a
//  struct grid_arrays::grid_pars  from the info in a  struct patch_info
// and the additional information in the arguments.
//
// The result refers to an internal static buffer in this function; the
// usual caveats about lifetimes/overwriting apply.
//
// Arguments:
// N_extend_points = Number of grid points to extend each patch past
//		     "just touching" so as to overlap neighboring patches.
//		     Thus patches overlap by  2*N_extend_points + 1  grid
//		     points.  For example, with N_extend_points == 2, here
//		     are the grid points of two neighboring patches:
//			x   x   x   x   x   X   X
//                                      |
//			        O   O   o   o   o   o   o
//		     Here | marks the "just touching" boundary,
//		     x and o the grid points before this extension,
//		     and X and O the extra grid points added by this
//		     extension.
// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
//
const grid::grid_pars& patch_info::grid_pars(int N_extend_points,
					     fp delta_drho_dsigma)
	const
{
static struct grid::grid_pars grid_pars_buffer;

const fp extend_drho_dsigma = fp(N_extend_points) * delta_drho_dsigma;

grid_pars_buffer.  min_drho   = min_drho   - extend_drho_dsigma;
grid_pars_buffer.delta_drho   = delta_drho_dsigma;
grid_pars_buffer.  max_drho   = max_drho   + extend_drho_dsigma;
grid_pars_buffer.  min_dsigma = min_dsigma - extend_drho_dsigma;
grid_pars_buffer.delta_dsigma = delta_drho_dsigma;
grid_pars_buffer.  max_dsigma = max_dsigma + extend_drho_dsigma;

return grid_pars_buffer;
}
	  }