path: root/src/patch/ghost_zone.cc

// ghost_zone.cc -- fill in gridfn data in patch ghost zones
// $Id$

//
// symmetry_ghost_zone::symmetry_ghost_zone (mirror symmetry)
// symmetry_ghost_zone::symmetry_ghost_zone (periodic BC)
// symmetry_ghost_zone::~symmetry_ghost_zone
// symmetry_ghost_zone::setup_scalar_gridfn_in_ghost_zone
//
// interpatch_ghost_zone::interpatch_ghost_zone
// interpatch_ghost_zone::~interpatch_ghost_zone
// interpatch_ghost_zone::assert_fully_setup
// interpatch_ghost_zone::setup_scalar_gridfn_in_ghost_zone
//

#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 "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "ghost_zone.hh"
#include "patch_frontier.hh"

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

//
// This function constructs a mirror-symmetry ghost zone object
//
symmetry_ghost_zone::symmetry_ghost_zone(const patch_edge& my_edge_in)
	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
	  symmetry_patch_(my_edge_in.my_patch()),
	  symmetry_edge_(my_edge_in)
{
// iperp_map: i --> (i of ghost zone) - i
iperp_map_ = new cpm_map(min_iperp(), max_iperp(),
			 my_edge_in.fp_grid_outer_iperp());

// ipar_map_: identity map
ipar_map_ = new cpm_map(my_edge_in.min_ipar_with_corners(),
			my_edge_in.max_ipar_with_corners());
}

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

//
// This function constructs a periodic-symmetry ghost zone object.
//
symmetry_ghost_zone::symmetry_ghost_zone
	(const patch_edge& my_edge_in, const patch_edge& symmetry_edge_in,
	 int my_edge_sample_ipar,      int symmetry_edge_sample_ipar,
	 bool ipar_map_sign)
	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
	  symmetry_patch_(symmetry_edge_in.my_patch()),
	  symmetry_edge_(symmetry_edge_in)
{
//
// perpendicular map
//
fp fp_my_period_plane_iperp = my_edge().fp_grid_outer_iperp();
fp fp_symmetry_period_plane_iperp = symmetry_edge().fp_grid_outer_iperp();

fp fp_iperp_shift_amount
	= fp_symmetry_period_plane_iperp - fp_my_period_plane_iperp;
if (! fuzzy<fp>::is_integer(fp_iperp_shift_amount))
   then error_exit(ERROR_EXIT,
"***** symmetry_ghost_zone::symmetry_ghost_zone (periodic symmetry):\n"
"        fp_symmetry_period_plane_iperp=%g - fp_my_period_plane_iperp=%g\n"
"        isn't fuzzily integral!\n"
"        ==> periodic-symmetry periodicity is incommensurate with grid!\n"
,
		   double(fp_symmetry_period_plane_iperp),
		   double(fp_my_period_plane_iperp));		/*NOTREACHED*/
int iperp_shift_amount = round<fp>::to_integer(fp_iperp_shift_amount);

// iperp mapping must be outside --> inside
// i.e. if both edges have iperp as the same min/max "direction",
//	   then the mapping is  iperp increasing --> iperp decreasing
bool iperp_map_sign = (my_edge().is_min() == symmetry_edge().is_min())
		      ? cpm_map::minus_map
		      : cpm_map::plus_map;
iperp_map_ = new
	     cpm_map(min_iperp(), max_iperp(),
		     fp_my_period_plane_iperp, fp_symmetry_period_plane_iperp,
		     iperp_map_sign);

//
// parallel map
//
ipar_map_ = new cpm_map(my_edge().min_ipar_with_corners(),
			my_edge().max_ipar_with_corners(),
			my_edge_sample_ipar, symmetry_edge_sample_ipar,
			ipar_map_sign);
}

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

//
// This function destroys a  symmetry_ghost_zone  object.
//
symmetry_ghost_zone::~symmetry_ghost_zone()
{
delete ipar_map_;
delete iperp_map_;
}

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

//
// This function symmetry-maps a scalar gridfn in a specified part of
// the ghost zone.  at a specified grid point.  ("Scalar" here refers
// solely to the symmetry properties of the gridfn under coordinate
// reflections, it need not be an actual scalar under more general
// coordinate transformations.)
//
void symmetry_ghost_zone::setup_scalar_gridfn_in_ghost_zone
	(int gfn,
	 bool want_min_par_corner,
	 bool want_non_corner,
	 bool want_max_par_corner)
	const
{
	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
	{
	// do we want to do this point?
	if (! my_edge().ipar_is_in_selected_corner(want_min_par_corner,
						   want_non_corner,
						   want_max_par_corner,
						   ipar) )
	   then continue;				// *** LOOP CONTROL ***

	const int sym_iperp = iperp_map_of_iperp(iperp);
	const int sym_ipar  = ipar_map_of_ipar  (ipar );
	const int sym_irho
		= symmetry_edge().irho_of_iperp_ipar  (sym_iperp,sym_ipar);
	const int sym_isigma
		= symmetry_edge().isigma_of_iperp_ipar(sym_iperp,sym_ipar);
	const fp sym_gridfn = symmetry_patch().gridfn(gfn, sym_irho,sym_isigma);

	const int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
	const int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
	my_patch().gridfn(gfn, irho,isigma) = sym_gridfn;
	}
	}
}

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

//
// This function constructs an  interpatch_ghost_zone  object.
//
interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
					     const patch_edge& other_edge_in)
	: ghost_zone(my_edge_in, ghost_zone_is_interpatch),
	  other_patch_(other_edge_in.my_patch()),
	  other_edge_(other_edge_in),
	  other_frontier_(NULL)		// other frontier object
					// may not exist yet
{
//
// verify that we have the expected relationships between
// this and the other patch's (mu,nu,phi) coordinates:
//

// perp coordinate is common to us and the other patch, so
// ghost zone/frontier must be min in one patch, max in the other
if (my_edge().is_min() == other_edge().is_min())
   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        ghost zone/frontier must be min in one patch, max in the other!\n"
,
		   my_patch().name(), my_edge().name(),
		   other_patch().name(), other_edge().name());	/*NOTREACHED*/

// coord in common between the two patches must be perp coord in both patches
// and this patch's tau coordinate must be other edge's parallel coordinate
local_coords::coords_set common_coords_set
	= local_coords::coords_set_not(my_patch().coords_set_rho_sigma()
				       ^
				       other_patch().coords_set_rho_sigma());
if (! (    (common_coords_set == my_edge().coords_set_perp())
	&& (common_coords_set == other_edge().coords_set_perp())
	&& (my_patch().coords_set_tau() == other_edge().coords_set_par())    ) )
   then error_exit(PANIC_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        (rho,sigma,tau) coordinates don't match up properly\n"
"        between this patch/edge and the other patch/edge!\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        my_patch().coords_set_{rho,sigma,tau}={%s,%s,%s}\n"
"        my_edge().coords_set_{perp,par}={%s,%s}\n"
"        other_patch().coords_set_{rho,sigma,tau}={%s,%s,%s}\n"
"        other_edge().coords_set_{perp,par}={%s,%s}\n"
,
	my_patch().name(), my_edge().name(),
	other_patch().name(), other_edge().name(),
	local_coords::name_of_coords_set(my_patch().coords_set_rho()),
	local_coords::name_of_coords_set(my_patch().coords_set_sigma()),
	local_coords::name_of_coords_set(my_patch().coords_set_tau()),
	local_coords::name_of_coords_set(my_edge().coords_set_perp()),
	local_coords::name_of_coords_set(my_edge().coords_set_par()),
	local_coords::name_of_coords_set(other_patch().coords_set_rho()),
	local_coords::name_of_coords_set(other_patch().coords_set_sigma()),
	local_coords::name_of_coords_set(other_patch().coords_set_tau()),
	local_coords::name_of_coords_set(other_edge().coords_set_perp()),
	local_coords::name_of_coords_set(other_edge().coords_set_par()));
								/*NOTREACHED*/

// perp coordinate must match across the two patches
if (fuzzy<fp>::NE(my_edge().grid_outer_perp(), other_edge().grid_outer_perp()))
   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        perp coordinate doesn't match across the two patches!\n"
"        my_edge().grid_outer_perp()=%g\n"
"        other_edge.grid_outer_perp()=%g\n"
,
		   my_patch().name(), my_edge().name(),
		   other_patch().name(), other_edge().name(),
		   double(my_edge().grid_outer_perp()),
		   double(other_edge().grid_outer_perp()));	/*NOTREACHED*/


//
// how long will this ghost zone be?
//
// ... note that we can't yet count on any other ghost zone existing,
//     so we don't yet know whether or not we'll want our corners
//     (since that depends on the type of our patch's adjacent ghost zones)
//     ==> we include the corners on the chance we may want them later
//
int ghost_zone_min_ipar = my_edge().min_ipar_with_corners();
int ghost_zone_max_ipar = my_edge().max_ipar_with_corners();


//
// set up the interpatch coordinate mappings
//
other_iperp_ = new array1d<int>(min_iperp(), max_iperp());
other_par_ = new array2d<fp>(min_iperp(), max_iperp(),
			     ghost_zone_min_ipar, ghost_zone_max_ipar);

	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	// compute the  other_iperp  corresponding to  iperp
	// using the fact that  perp  is the same coordinate in both cases
	const fp perp = my_edge().perp_of_iperp(iperp);
	const fp fp_other_iperp = other_edge().fp_iperp_of_perp(perp);

	// verify that  other_iperp  is fuzzily a grid point
	if (! fuzzy<fp>::is_integer(fp_other_iperp))
	   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        iperp=%d perp=%g\n"
"        ==> fp_other_iperp=%g isn't fuzzily an integer\n"
"        ==> patches aren't commensurate in the perpendicular coordinate!\n"
,
			   my_patch().name(), my_edge().name(),
			   other_patch().name(), other_edge().name(),
			   iperp, double(perp),
			   double(fp_other_iperp));		/*NOTREACHED*/

	(*other_iperp_)(iperp) = round<fp>::to_integer(fp_other_iperp);

		for (int ipar = ghost_zone_min_ipar ;
		     ipar <= ghost_zone_max_ipar ;
		     ++ipar)
		{
		// compute the  other_par corresponding to  (iperp,ipar)
		// ... here we use the fact (which we verified above) that
		//     other edge's parallel coordinate == our tau coordinate
		const fp par   = my_edge().par_of_ipar(ipar);

		const fp rho   = my_edge().  rho_of_perp_par(perp, par);
		const fp sigma = my_edge().sigma_of_perp_par(perp, par);

		const fp tau   = my_patch().tau_of_rho_sigma(rho, sigma);

		(*other_par_)(iperp,ipar) = tau;
		}
	}
}

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

//
// this function destroys an  interpatch_ghost_zone  object.
//
interpatch_ghost_zone::~interpatch_ghost_zone()
{
delete other_par_;
delete other_iperp_;
}

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

//
// This function assert()s that
// * our frontier pointer is non-NULL,
// * the other patch has an interpatch ghost zone on this edge
// * that interpatch ghost zone points back to us
//
void interpatch_ghost_zone::assert_fully_setup() const
{
assert(other_frontier_ != NULL);

const ghost_zone& ogz = other_patch().ghost_zone_on_edge(other_edge());
assert(ogz.is_interpatch());

const interpatch_ghost_zone& oigz
	= static_cast<const interpatch_ghost_zone &>(ogz);
assert(& oigz.other_patch() == & my_patch());
}

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

//
// This function sets up a scalar gridfn in a specified part of the
// ghost zone, by interpatch-interpolating it from the neighboring
// patch's frontier.  ("Scalar" here refers solely to the symmetry
// properties of the gridfn under coordinate reflections, it need not
// be an actual scalar under more general coordinate transformations.)
//
void interpatch_ghost_zone::setup_scalar_gridfn_in_ghost_zone
	(int gfn,
	 bool want_min_par_corner,
	 bool want_non_corner,
	 bool want_max_par_corner)
	const
{
other_frontier_->setup_interpolation_for_gridfn(gfn);

	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
	{
	// do we want to do this point?
	if (! my_edge().ipar_is_in_selected_corner(want_min_par_corner,
						   want_non_corner,
						   want_max_par_corner,
						   ipar) )
	   then continue;				// *** LOOP CONTROL ***

	int oiperp = (*other_iperp_)(iperp);
	fp opar = (*other_par_)(iperp,ipar);
	fp ogridfn = other_frontier_->interpolate(gfn, oiperp, opar);

	int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
	int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
	my_patch().gridfn(gfn, irho,isigma) = ogridfn;
	}
	}
}