605 lines
24 KiB
C
605 lines
24 KiB
C
#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
|