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::ghost_zone_on_edge
// patch::corner_ghost_zone_containing_point
// patch::symmetry_ghost_zone_on_edge
// patch::interpatch_ghost_zone_on_edge
// patch::create_mirror_symmetry_ghost_zone
// patch::create_periodic_symmetry_ghost_zone
// patch::create_interpatch_ghost_zone
// patch::set_ghost_zone
// patch::edge_adjacent_to_patch
// patch::assert_all_ghost_zones_fully_setup
//

#include <cstdio>
#include <cmath>
#include <assert.h>
using std::fprintf;
using std::printf;

#include "cctk.h"

#include "stdc.h"
#include "config.hh"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

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

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

//
// 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) // 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_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 ghost zone on a specified
// edge, after first assert()ing that the edge belongs to this patch.
//
// N.b. This function can't be inline in "patch.hh" because it needs
//	member functions of class patch_edge, which comes after class patch
//	in our #include order.
//
ghost_zone& patch::ghost_zone_on_edge(const patch_edge& e)
	const
{
assert(& e.my_patch() == this);
return minmax_ang_ghost_zone(e.is_min(), e.is_rho());
}

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

//
// This function determines which of the two adjacent ghost zones meeting
// at a specified corner, contains a specified point.  If the point isn't
// in either ghost zone, an abort(0) is done.
//
// Arguments:
// {rho,sigma}_is_min = Specify the corner (and implicitly the ghost zones).
// irho,isigma = Specify the point.
//
// Results:
// This function returns (a reference to) the desired ghost zone.
ghost_zone& patch::corner_ghost_zone_containing_point
	(bool rho_is_min, bool sigma_is_min,
	 int irho, int isigma)
	const
{
ghost_zone&   rho_gz =   minmax_rho_ghost_zone(  rho_is_min);
ghost_zone& sigma_gz = minmax_sigma_ghost_zone(sigma_is_min);

const patch_edge&   rho_edge =   rho_gz.my_edge();
const patch_edge& sigma_edge = sigma_gz.my_edge();

const int   rho_iperp =   rho_edge.iperp_of_irho_isigma(irho, isigma);
const int   rho_ipar  =   rho_edge. ipar_of_irho_isigma(irho, isigma);
const int sigma_iperp = sigma_edge.iperp_of_irho_isigma(irho, isigma);
const int sigma_ipar  = sigma_edge. ipar_of_irho_isigma(irho, isigma);

const bool is_in_rho_ghost_zone
	=   rho_gz.is_in_ghost_zone(  rho_iperp,   rho_ipar);
const bool is_in_sigma_ghost_zone
	= sigma_gz.is_in_ghost_zone(sigma_iperp, sigma_ipar);

// check that point is in exactly one ghost zone
assert(is_in_rho_ghost_zone ^ is_in_sigma_ghost_zone);

return is_in_rho_ghost_zone ? rho_gz : sigma_gz;
}

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

//
// These functions verify that the ghost zone on a specified edge is
// indeed of the specified type, and return a reference to it as the
// appropriate derived class.
//
symmetry_ghost_zone& patch::symmetry_ghost_zone_on_edge
	(const patch_edge &e)
	const
{
ghost_zone &g = ghost_zone_on_edge(e);
assert(g.is_symmetry());
return static_cast<symmetry_ghost_zone &>(g);
}

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

//
// 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 assert()s that a specified ghost zone of this patch
// hasn't already been set up, then constructs it as a mirror-symmetry
// ghost zone and properly links this to/from the patch.
//
void patch::create_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(my_edge, temp);
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then creates it as a periodic-symmetry
// ghost zone and properly links this to/from the patch.
//
void patch::create_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);
set_ghost_zone(my_edge, temp);
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then creates it as an interpatch ghost
// zone (with lots of NULL pointers for info we can't compute yet)
// and properly links this to/from the patch.
//
void patch::create_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);
set_ghost_zone(my_edge, temp);
}

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

//
// This is a helper function for  setup_*_ghost_zone().  This function
// assert()s that one of the ghost zone pointers (which one is selected
// by  edge ) is NULL, then stores a value in it.
//
void patch::set_ghost_zone(const patch_edge& edge, ghost_zone* gzp)
{
ghost_zone*& ghost_zone_ptr_to_set
	= 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(ghost_zone_ptr_to_set == NULL);
ghost_zone_ptr_to_set = gzp;
}

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

//
// 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?
// ... eg N_overlap_points = 3 would be
//	p   p   p   p   p
//		q   q   q   q   q
//     so the overlap would be (N_overlap_points-1) * delta = 2 * delta
if ( (N_overlap_points-1 != 0)
     && jtutil::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-1) * 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 (via assert()) that all ghost zones of this
// patch 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();
}