#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include "cctk.h"
#include "config.h"
#include "stdc.h"
#include "util.h"
#include "array.h"
#include "cpm_map.h"
#include "linear_map.h"
#include "coords.h"
#include "tgrid.h"
#include "fd_grid.h"
#include "patch.h"
#include "patch_edge.h"
#include "patch_interp.h"
#include "ghost_zone.h"
namespace AHFinderDirect
{
using jtutil::error_exit;
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// These functions verify (assert()) that a ghost zone is indeed of
// the specified type, then static_cast to the appropriate derived class.
//
const symmetry_ghost_zone &ghost_zone::cast_to_symmetry_ghost_zone()
const
{
assert(is_symmetry());
return static_cast<const symmetry_ghost_zone &>(*this);
}
symmetry_ghost_zone &ghost_zone::cast_to_symmetry_ghost_zone()
{
assert(is_symmetry());
return static_cast<symmetry_ghost_zone &>(*this);
}
//**************************************
const interpatch_ghost_zone &ghost_zone::cast_to_interpatch_ghost_zone()
const
{
assert(is_interpatch());
return static_cast<const interpatch_ghost_zone &>(*this);
}
interpatch_ghost_zone &ghost_zone::cast_to_interpatch_ghost_zone()
{
assert(is_interpatch());
return static_cast<interpatch_ghost_zone &>(*this);
}
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// 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,
my_edge_in, // other edge == my edge
ghost_zone_is_symmetry)
{
// iperp_map: i --> (i of ghost zone) - i
iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
my_edge_in.fp_grid_outer_iperp());
// ipar_map_: identity map
ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar());
}
//******************************************************************************
//
// This function constructs a periodic-symmetry ghost zone object.
//
symmetry_ghost_zone::symmetry_ghost_zone(const patch_edge &my_edge_in, const patch_edge &other_edge_in,
int my_edge_sample_ipar, int other_edge_sample_ipar,
bool ipar_map_is_plus)
: ghost_zone(my_edge_in,
other_edge_in,
ghost_zone_is_symmetry)
{
//
// perpendicular map
//
const fp fp_my_period_plane_iperp = my_edge().fp_grid_outer_iperp();
const fp fp_other_period_plane_iperp = other_edge().fp_grid_outer_iperp();
// 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
// (i.e. the map's sign is -1)
const bool is_iperp_map_plus = !(my_edge().is_min() == other_edge().is_min());
iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
fp_my_period_plane_iperp,
fp_other_period_plane_iperp,
is_iperp_map_plus);
//
// parallel map
//
ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar(),
my_edge_sample_ipar, other_edge_sample_ipar,
ipar_map_is_plus);
}
//******************************************************************************
//
// This function destroys a symmetry_ghost_zone object.
//
symmetry_ghost_zone::~symmetry_ghost_zone()
{
delete ipar_map_;
delete iperp_map_;
}
//******************************************************************************
//
// This function "synchronizes" a ghost zone, i.e. it updates the
// ghost-zone values of the specified gridfns via the appropriate
// symmetry operations.The flags specify which part(s) of the ghost zone
// we want.
//
void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
bool want_corners /* = true */,
bool want_noncorner /* = true */)
{
// printf("*Sync sym ghost zone in %s patch\n", my_patch().name());
for (int gfn = ghosted_min_gfn; gfn <= ghosted_max_gfn; ++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_part(want_corners, want_noncorner,
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 = other_edge()
.irho_of_iperp_ipar(sym_iperp, sym_ipar);
const int sym_isigma = other_edge()
.isigma_of_iperp_ipar(sym_iperp, sym_ipar);
const fp sym_gridfn = other_patch()
.ghosted_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().ghosted_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,
int patch_overlap_width)
: ghost_zone(my_edge_in,
other_edge_in,
ghost_zone_is_interpatch),
// remaining pointers are all set up properly by finish_setup()
other_patch_interp_(NULL),
other_iperp_(NULL),
min_ipar_used_(NULL), max_ipar_used_(NULL),
other_par_(NULL),
interp_result_buffer_(NULL),
Jacobian_y_ipar_posn_(NULL), Jacobian_buffer_(NULL) // no comma
{
//
// 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 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 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
const 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 (mod 2*pi) across the two patches
// after taking into account any overlap
// ... eg patch_overlap_width = 3 would be
// p p p p p
// q q q q q
// so the overlap would be (patch_overlap_width-1) * delta
const fp other_overlap = (patch_overlap_width - 1) * other_edge().perp_map().delta_fp();
const fp other_outer_perp_minus_overlap // move back inwards into other patch
// by overlap distance, to get a value
// that should match our own
// grid_outer_perp() value
= other_edge().grid_outer_perp() + (other_edge().is_min() ? +other_overlap : -other_overlap);
if (!local_coords::fuzzy_EQ_ang(my_edge().grid_outer_perp(),
other_outer_perp_minus_overlap))
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 (mod 2*pi) across the two patches!\n"
" my_edge().grid_outer_perp()=%g <--(compare this)\n"
" patch_overlap_width=%d other_overlap=%g\n"
" other_edge.grid_outer_perp()=%g\n"
" other_outer_perp_minus_overlap=%g <--(against this)\n",
my_patch().name(), my_edge().name(),
other_patch().name(), other_edge().name(),
double(my_edge().grid_outer_perp()),
patch_overlap_width, double(other_overlap),
double(other_edge().grid_outer_perp()),
double(other_outer_perp_minus_overlap)); /*NOTREACHED*/
//
// set up the iperp interpatch coordinate mapping
// (gives other patch's iperp coordinate for interpolation)
//
// compute the iperp --> other_iperp mapping for a sample point;
// ... if the ghost zone is empty, then the sample point will necessarily
// be out-of-range in the ghost zone, so we use the *unchecked*
// conversions to avoid errors in this case
// ... we do the computation using the fact that perp is the same
// coordinate in both patches (modulo 2*pi radians = 360 degrees)
const int sample_iperp = outer_iperp();
const fp sample_perp = my_edge().perp_map().fp_of_int_unchecked(sample_iperp);
// unchecked conversion here!
const fp other_sample_perp = other_patch()
.modulo_reduce_ang(other_edge().perp_is_rho(),
sample_perp);
const fp fp_other_sample_iperp = other_edge()
.fp_iperp_of_perp(other_sample_perp);
// verify that this is fuzzily a grid point
if (!jtutil::fuzzy<fp>::is_integer(fp_other_sample_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"
" sample_iperp=%d sample_perp=%g\n"
" other_sample_perp=%g fp_other_sample_iperp=%g\n"
" ==> fp_other_sample_iperp 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(),
sample_iperp, double(sample_perp),
double(other_sample_perp),
double(fp_other_sample_iperp)); /*NOTREACHED*/
const int other_sample_iperp = jtutil::round<fp>::to_integer(fp_other_sample_iperp);
// compute the +/- sign (direction) of the iperp --> other_iperp mapping
//
// Since perp is the same in both patches (mod 2*pi radians = 360 degrees),
// the overall +/- sign is just the product of the signs of the two individual
// iperp <--> perp mappings.
//
// ... signs encoded as (floating-point) +/- 1.0
const double iperp_map_sign_pm1 = jtutil::signum(my_edge().perp_map().delta_fp()) * jtutil::signum(other_edge().perp_map().delta_fp());
// ... signs encoded as is_plus bool flag
const bool is_iperp_map_plus = (iperp_map_sign_pm1 > 0.0);
// now we finally know enough to set up the other_iperp(iperp)
// coordinate mapping
other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
sample_iperp, other_sample_iperp,
is_iperp_map_plus);
}
//******************************************************************************
//
// this function destroys an interpatch_ghost_zone object.
//
interpatch_ghost_zone::~interpatch_ghost_zone()
{
delete Jacobian_buffer_;
delete Jacobian_y_ipar_posn_;
delete interp_result_buffer_;
delete other_par_;
delete max_ipar_used_;
delete min_ipar_used_;
delete other_iperp_;
delete other_patch_interp_;
}
//******************************************************************************
//
// These functions compute the [min,max] ipar of the ghost zone for
// a given iperp, taking into account how we treat the corners
// (cf. the example in the header comments in "ghost_zone.hh"):
//
// If an adjacent ghost zone is symmetry,
// we do not include that corner;
// If an adjacent ghost zone is interpatch,
// we include up to the diagonal line, and if we are a rho ghost zone,
// then also the diagonal line itself. E.g. For the example in the
// header comments "ghost_zone.hh", the +x ghost zone includes (6,6),
// (7,6), and (7,7), while the +y ghost zone includes (6,7)
//
// ... in the following 2 functions,
// the iabs() term includes the diagonal,
// so we must remove the diagonal for !is_rho,
// i.e. add 1 to min_ipar and subtract 1 from max_ipar
//
int interpatch_ghost_zone::min_ipar(int iperp) const
{
return min_par_adjacent_ghost_zone().is_symmetry()
? my_edge().min_ipar_without_corners()
: my_edge().min_ipar_without_corners() - iabs(iperp - my_edge().nominal_grid_outer_iperp()) + (is_rho() ? 0 : 1);
}
int interpatch_ghost_zone::max_ipar(int iperp) const
{
return max_par_adjacent_ghost_zone().is_symmetry()
? my_edge().max_ipar_without_corners()
: my_edge().max_ipar_without_corners() + iabs(iperp - my_edge().nominal_grid_outer_iperp()) - (is_rho() ? 0 : 1);
}
//******************************************************************************
//
// This function finishes the construction/setup of an interpatch_ghost_zone
// object. It
// - sets up the par coordinate mapping information
// - sets up the interpatch interpolator data pointer and result arrays
// - constructs the patch_interp object to interpolate from the *other* patch
//
// We use our ipar as the patch_interp's parindex.
//
void interpatch_ghost_zone::finish_setup(int interp_handle,
int interp_par_table_handle)
{
min_other_iperp_ = min(other_iperp(min_iperp()),
other_iperp(max_iperp()));
max_other_iperp_ = max(other_iperp(min_iperp()),
other_iperp(max_iperp()));
//
// set up arrays giving actual [min,max] ipar that we'll use
// at each other_iperp (later on we will pass these arrays to the
// other patch's patch_interp object, with ipar being parindex there
//
min_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
max_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
{
for (int iperp = min_iperp(); iperp <= max_iperp(); ++iperp)
{
(*min_ipar_used_)(other_iperp(iperp)) = min_ipar(iperp);
(*max_ipar_used_)(other_iperp(iperp)) = max_ipar(iperp);
}
}
//
// set up array giving other patch's par coordinate for interpolation
//
other_par_ = new jtutil::array2d<fp>(min_other_iperp_, max_other_iperp_,
extreme_min_ipar(), extreme_max_ipar());
{
for (int iperp = min_iperp(); iperp <= max_iperp(); ++iperp)
{
for (int ipar = min_ipar(iperp); ipar <= max_ipar(iperp); ++ipar)
{
// compute the other_par corresponding to (iperp,ipar)
// ... here we use the fact (which we verified in our constructor)
// that other edge's parallel coordinate == our tau coordinate
// (at least modulo 2*pi radians = 360 degrees)
const fp perp = my_edge().perp_of_iperp(iperp);
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);
const fp other_par = other_patch()
.modulo_reduce_ang(other_edge().par_is_rho(), tau);
(*other_par_)(other_iperp(iperp), ipar) = other_par;
}
}
}
//
// set up interpolation result buffer
//
interp_result_buffer_ = new jtutil::array3d<fp>(my_patch().ghosted_min_gfn(),
my_patch().ghosted_max_gfn(),
min_other_iperp_, max_other_iperp_,
extreme_min_ipar(), extreme_max_ipar());
//
// construct the patch_interp object to interpolate from the *other* patch
// ... the patch_interp should use gridfn data from it's (the other patch's)
// min/max par ghost zones if those (adjacent) adjacent ghost zones
// are symmetry, but not if they're interpatch,
// cf the header comments in "ghost_zone.hh"
//
const ghost_zone &other_ghost_zone = other_patch()
.ghost_zone_on_edge(other_edge());
const bool ok_to_use_min_par_ghost_zone = other_ghost_zone.min_par_adjacent_ghost_zone()
.is_symmetry()
? true
: false;
const bool ok_to_use_max_par_ghost_zone = other_ghost_zone.max_par_adjacent_ghost_zone()
.is_symmetry()
? true
: false;
other_patch_interp_ = new patch_interp(other_edge(),
min_other_iperp_, max_other_iperp_,
*min_ipar_used_, *max_ipar_used_,
*other_par_,
ok_to_use_min_par_ghost_zone,
ok_to_use_max_par_ghost_zone,
interp_handle, interp_par_table_handle);
}
//******************************************************************************
//
// This function asserts() that
// - we have a patch_interp object
// - our and the patch_interp object's notions of the "other patch" agree
// - the other patch has an interpatch ghost zone on this edge
// - the other patch's interpatch ghost zone on this edge,
// points back to our patch
//
void interpatch_ghost_zone::assert_fully_setup() const
{
assert(other_patch_interp_ != NULL);
assert(other_patch() == other_patch_interp_->my_patch());
assert(other_patch()
.ghost_zone_on_edge(other_edge())
.is_interpatch());
assert(my_patch() == other_patch()
.ghost_zone_on_edge(other_edge())
.other_patch());
}
//******************************************************************************
//
// This function "synchronizes" a ghost zone, i.e. it updates the
// ghost-zone values of the specified gridfns via the appropriate
// interpatch interpolations.
//
// The flags specify which part(s) of the ghost zone we want, but
// the present implementation only supports the case where all the
// flags are true , i.e. we want the entire ghost zone.
//
void interpatch_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
bool want_corners /* = true */,
bool want_noncorner /* = true */)
{
#ifdef DEBUG_AHFD
printf("*Sync interpatch ghost zone in %s\n", my_patch().name());
#endif
// make sure the caller wants the entire ghost zone
if (!(want_corners && want_noncorner))
then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::synchronize():\n"
" we only support operating on the *entire* ghost zone,\n"
" but we were passed flags specifying a proper subset!\n"
" want_corners=(int)%d want_noncorner=(int)%d\n",
want_corners, want_noncorner); /*NOTREACHED*/
//
// move from 'Compute_Jacobian' below
//
assert(other_patch_interp_ != NULL);
other_patch_interp_->molecule_minmax_ipar_m(Jacobian_min_y_ipar_m_,
Jacobian_max_y_ipar_m_);
#ifdef DEBUG_AHFD
printf("%d %d %d %d %d %d \n", Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_,
min_other_iperp_, max_other_iperp_, extreme_min_ipar(), extreme_max_ipar());
getchar();
#endif
// /*
if (Jacobian_y_ipar_posn_ == NULL)
Jacobian_y_ipar_posn_ = new jtutil::array2d<CCTK_INT>(min_other_iperp_, max_other_iperp_,
extreme_min_ipar(), extreme_max_ipar());
if (Jacobian_buffer_ == NULL)
Jacobian_buffer_ = new jtutil::array3d<fp>(min_other_iperp_, max_other_iperp_,
extreme_min_ipar(), extreme_max_ipar(),
Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_);
// do the interpolation into our result buffer
other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
*interp_result_buffer_, //);
*Jacobian_y_ipar_posn_,
*Jacobian_buffer_);
// other_patch_interp_->molecule_posn(*Jacobian_y_ipar_posn_);
// store the results back into our gridfns
for (int gfn = ghosted_min_gfn; gfn <= ghosted_max_gfn; ++gfn)
{
for (int iperp = min_iperp(); iperp <= max_iperp(); ++iperp)
{
const int oiperp = other_iperp(iperp);
for (int ipar = min_ipar(iperp); ipar <= max_ipar(iperp); ++ipar)
{
int irho = my_edge().irho_of_iperp_ipar(iperp, ipar);
int isigma = my_edge().isigma_of_iperp_ipar(iperp, ipar);
my_patch().ghosted_gridfn(gfn, irho, isigma) = (*interp_result_buffer_)(gfn, oiperp, ipar);
}
}
}
}
//******************************************************************************
//
// This function allocates the internal buffers for the Jacobian, and
// computes that Jacobian
// partial synchronize gridfn(ghosted_gfn, iperp, ipar)
// ------------------------------------------------------------
// partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m)
// where
// oiperp = Jacobian_oiperp(iperp)
// posn = Jacobian_oipar_posn(iperp, ipar)
// into the internal buffers.
//
void interpatch_ghost_zone::compute_Jacobian(int ghosted_min_gfn, int ghosted_max_gfn,
bool want_corners /* = true */,
bool want_noncorner /* = true */)
const
{
// make sure the caller wants the entire ghost zone
if (!(want_corners && want_noncorner))
then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::compute_Jacobian():\n"
" we only support operating on the *entire* ghost zone,\n"
" but we were passed flags specifying a proper subset!\n"
" want_corners=(int)%d want_noncorner=(int)%d\n",
want_corners, want_noncorner); /*NOTREACHED*/
}
//******************************************************************************
//******************************************************************************
//******************************************************************************
} // namespace AHFinderDirect