294 lines
11 KiB
C++
294 lines
11 KiB
C++
#ifndef TPATCH_INTERP_H
|
|
#define TPATCH_INTERP_H
|
|
namespace AHFinderDirect
|
|
{
|
|
|
|
//
|
|
// patch_interp - interpolation from a patch
|
|
//
|
|
|
|
//
|
|
// A patch_interp object is responsible for interpolating gridfn data
|
|
// from its owning patch for use by another patch's ghost_zone object
|
|
// (in setting up the gridfn in the other ghost zone). A patch_interp
|
|
// object deals only in its own patch's coordinates; other code elsewhere
|
|
// (in practice in interpatch_ghost_zone::) is responsible for translating
|
|
// other patch's coordinates into our coordinates.
|
|
//
|
|
|
|
//
|
|
// A patch_interp defines a "patch interpolation region", the region of
|
|
// its owning patch from which this interpolation will use gridfn data.
|
|
//
|
|
|
|
//
|
|
// The way the patch coordnates are constructed, any two adjacent patches
|
|
// share a common (perpendicular) coordinate. Thus we only have to do
|
|
// 1-dimensional interpolation here (in the parallel direction). In
|
|
// other words, for each iperp we interpolate in par.
|
|
//
|
|
// In general we interpolate each gridfn at a number of distinct par
|
|
// for each iperp; the integer "parindex" indexes these points. We
|
|
// attach no particular semantics to parindex, and it need not be
|
|
// 0-origin or have the same range for each iperp. [In practice,
|
|
// parindex will be the other patch's ipar coordinate.] However,
|
|
// we assume that the range of parindex is roughly similar for each
|
|
// iperp, so it's ok to use (iperp,parindex) as a 2-D rectangular
|
|
// index space.
|
|
//
|
|
// For example, we might interpolate at the points
|
|
// ipar ipar ipar ipar ipar ipar ipar ipar ipar
|
|
// 1 2 3 4 5 6 7 8 9
|
|
// iperp=10 (2a) (3b) (4c)
|
|
// iperp=11 (2d) (3e) (4f) (5g)
|
|
// where the (2a)-(5g) are the interpolation points, with 2-5 being the
|
|
// parindex values and a-g being unique identifiers used in our description
|
|
// below. In terms of our member data, this interpolation region would
|
|
// be described by
|
|
// [min,max]_iperp_=[10,11]
|
|
// [min,max]_ipar_=[1,9]
|
|
// [min,max]_parindex_array_(10)=[2,5]
|
|
// [min,max]_parindex_array_(11)=[2,6]
|
|
// interp_par_(10,2) = x[a]
|
|
// interp_par_(10,3) = x[b]
|
|
// interp_par_(10,4) = x[c]
|
|
// interp_par_(11,2) = x[d]
|
|
// interp_par_(11,3) = x[e]
|
|
// interp_par_(11,4) = x[f]
|
|
// interp_par_(11,5) = x[g]
|
|
//
|
|
|
|
//
|
|
// We use the Cactus local interpolator CCTK_InterpLocalUniform()
|
|
// to do the interpolation. To minimize interpolator overheads, we
|
|
// interpolate all the gridfns at each iperp in a single interpolator
|
|
// call. [Different iperp values involve different sets of (1-D)
|
|
// gridfn data, and so inherently require distinct interpolator calls.]
|
|
//
|
|
// Setting up the array subscripting for the interpolator to access
|
|
// the gridfn data is a bit tricky: The interpolator accesses the
|
|
// gridfn data using the generic (1-D) subscripting expression
|
|
// data[offset + i*stride]
|
|
// where i is the data array index. However, we'd rather not use
|
|
// offset , because it has to be supplied in the parameter table as
|
|
// an array subscripted by gfn , and so would require changing the
|
|
// parameter table for each call on interpolate() (with potentially
|
|
// different numbers of gridfns being interpolated). Instead, at each
|
|
// iperp we use i = ipar-min_ipar , so the default offset=0 makes
|
|
// the subscripting expression zero for ipar = min_ipar . This also
|
|
// makes the interpolator's min_i = 0 and max_i be dims-1 (both
|
|
// the defaults), so those also don't have to be set in the parameter
|
|
// table either. We set the interpolator's data coordinate origin to
|
|
// the par coordinate for min_ipar , so it correctly maps i --> par .
|
|
// With this strategy we can share the interpolator parameter table
|
|
// across all the iperp values, and we don't need to modify the
|
|
// parameter table at all after the initial setup in our constructor.
|
|
// However, we do have to adjust the molecule positions in
|
|
// patch_interp::molecule_posn() , since the interpolator will return
|
|
// i values, while molecule_posn() needs ipar values.
|
|
//
|
|
|
|
class patch_interp
|
|
{
|
|
public:
|
|
// to which patch/edge do we belong?
|
|
const patch& my_patch() const { return my_patch_; }
|
|
const patch_edge& my_edge() const { return my_edge_; }
|
|
|
|
|
|
public:
|
|
//
|
|
// ***** main client interface *****
|
|
//
|
|
// interpolate specified range of ghosted gridfns
|
|
// at all the coordinates specified when we were constructed,
|
|
// store interpolated data in
|
|
// data_buffer(ghosted_gfn, iperp, parindex)
|
|
void interpolate(int ghosted_min_gfn_to_interp,
|
|
int ghosted_max_gfn_to_interp,
|
|
jtutil::array3d<fp>& data_buffer)
|
|
const;
|
|
void interpolate(int ghosted_min_gfn_to_interp,
|
|
int ghosted_max_gfn_to_interp,
|
|
jtutil::array3d<fp>& data_buffer,
|
|
jtutil::array2d<CCTK_INT>& posn_buffer,
|
|
jtutil::array3d<fp>& Jacobian_buffe)
|
|
const;
|
|
|
|
public:
|
|
//
|
|
// ***** Jacobian of interpolate() *****
|
|
//
|
|
|
|
// verify (no-op if ok, error_exit() if not) that interpolator
|
|
// has a Jacobian sparsity pattern which we grok: at present this
|
|
// means molecules are fixed-sized hypercubes, with size/shape
|
|
// independent of interpolation coordinates and the floating-point
|
|
// values in the input arrays
|
|
void verify_Jacobian_sparsity_pattern_ok() const;
|
|
|
|
//
|
|
// The API for the remaining Jacobian functions implicitly
|
|
// assumes that the Jacobian sparsity pattern is "ok" as
|
|
// verified by verify_Jacobian_sparsity_pattern_ok() ,
|
|
// and in particular that [min,max]_ipar_m are independent
|
|
// of iperp and parindex.
|
|
//
|
|
|
|
// get [min,max] ipar m coordinates of interpolation molecules
|
|
void molecule_minmax_ipar_m(int& min_ipar_m, int& max_ipar_m) const;
|
|
|
|
// get interpolation molecule ipar positions in
|
|
// molecule_posn_buffer(iperp, parindex)
|
|
// ... array type is CCTK_INT so we can pass by reference
|
|
// to interpolator
|
|
void molecule_posn(jtutil::array2d<CCTK_INT>& posn_buffer) const;
|
|
|
|
// get Jacobian of interpolated data with respect to this patch's
|
|
// ghosted gridfns,
|
|
// partial interpolate() data_buffer(ghosted_gfn, iperp, parindex)
|
|
// ---------------------------------------------------------------
|
|
// partial ghosted_gridfn(ghosted_gfn, iperp, posn+ipar_m)
|
|
// store Jacobian in
|
|
// Jacobian_buffer(iperp, parindex, ipar_m)
|
|
// where we implicitly assume the Jacobian to be independent of
|
|
// ghosted_gfn, and where
|
|
// posn = posn_buffer(iperp, parindex)
|
|
// as returned by molecule_posn()
|
|
void Jacobian(jtutil::array3d<fp>& Jacobian_buffer) const;
|
|
|
|
//
|
|
// ***** internal functions *****
|
|
//
|
|
private:
|
|
// [min,max] iperp for interpolation and gridfn data
|
|
int min_iperp() const { return min_iperp_; }
|
|
int max_iperp() const { return max_iperp_; }
|
|
|
|
// min/max (iperp,ipar) of the gridfn data to use for interpolation
|
|
int min_ipar() const { return min_ipar_; }
|
|
int max_ipar() const { return max_ipar_; }
|
|
|
|
//
|
|
// ***** constructor, destructor, et al *****
|
|
//
|
|
public:
|
|
//
|
|
// Constructor arguments:
|
|
// my_edge_in = Identifies the patch/edge to which this
|
|
// interpolation region is to belong.
|
|
// [min,max]_iperp_in = The range of iperp for this interpolation
|
|
// region
|
|
// [min,max]_parindex_array_in(iperp)
|
|
// = [min,max] range of parindex actually used at each iperp.
|
|
// We keep references to these arrays, so they should have
|
|
// lifetimes at last as long as that of this object.
|
|
// interp_par_in(iperp,parindex)
|
|
// = Gives the par coordinates at which we will interpolate;
|
|
// array entries outside the range [min,max]_parindex_in
|
|
// are unused. We keep a reference to this array, so it
|
|
// should have a lifetime at last as long as that of this
|
|
// object.
|
|
// ok_to_use_[min,max]_par_ghost_zone
|
|
// = Boolean flags saying whether or not we should use gridfn
|
|
// data from the [min,max]_par ghost zones in the interpolation.
|
|
// interp_handle_in = Cactus handle to the interpatch interpolation
|
|
// operator.
|
|
// interp_par_table_handle_in
|
|
// = Cactus handle to a Cactus key/value table giving
|
|
// parameters (eg order) for the interpatch interpolation
|
|
// operator. This class internally clones this table and
|
|
// modifies the clone, so the original table is not modified
|
|
// by any actions of this class.
|
|
//
|
|
// This constructor requires that this patch's gridfns already
|
|
// exist, since we size various arrays based on the patch's min/max
|
|
// ghosted gfn.
|
|
//
|
|
patch_interp(const patch_edge& my_edge_in,
|
|
int min_iperp_in, int max_iperp_in,
|
|
const jtutil::array1d<int>& min_parindex_array_in,
|
|
const jtutil::array1d<int>& max_parindex_array_in,
|
|
const jtutil::array2d<fp>& interp_par_in,
|
|
bool ok_to_use_min_par_ghost_zone,
|
|
bool ok_to_use_max_par_ghost_zone,
|
|
int interp_handle_in, int interp_par_table_handle_in);
|
|
~patch_interp();
|
|
|
|
private:
|
|
// we forbid copying and passing by value
|
|
// by declaring the copy constructor and assignment operator
|
|
// private, but never defining them
|
|
patch_interp(const patch_interp& rhs);
|
|
patch_interp& operator=(const patch_interp& rhs);
|
|
|
|
|
|
//
|
|
// ***** data members *****
|
|
//
|
|
private:
|
|
const patch& my_patch_;
|
|
const patch_edge& my_edge_;
|
|
|
|
// range of gfn we can handle
|
|
// (any given interpolate() call may specify a subrange)
|
|
const int min_gfn_, max_gfn_;
|
|
|
|
// these are strictly speaking redundant
|
|
// but we keep them for use in debugging
|
|
bool ok_to_use_min_par_ghost_zone_, ok_to_use_max_par_ghost_zone_;
|
|
|
|
// patch interpolation region,
|
|
// i.e. range of (iperp,ipar) in this patch from which
|
|
// we will use gridfn data in interpolation
|
|
const int min_iperp_, max_iperp_;
|
|
const int min_ipar_, max_ipar_;
|
|
|
|
// [min,max] parindex at each iperp
|
|
// ... these are references to arrays passed in to our constructor
|
|
// ==> we do *not* own them!
|
|
// ... indices are (iperp)
|
|
const jtutil::array1d<int>& min_parindex_array_;
|
|
const jtutil::array1d<int>& max_parindex_array_;
|
|
|
|
// interp_par_(iperp,parindex)
|
|
// = Gives the par coordinates at which we will interpolate;
|
|
// array entries outside the range [min,max]_parindex_in
|
|
// are unused (n.b. this interface implicitly takes the
|
|
// par coordinates to be independent of ghosted_gfn).
|
|
// ... this is a reference to an array passed in to our constructor
|
|
// ==> we do *not* own this!
|
|
const jtutil::array2d<fp>& interp_par_; // indices (iperp,parindex)
|
|
|
|
// Cactus handle to the interpolation operator
|
|
int interp_handle_;
|
|
|
|
// Cactus handle to our private Cactus key/value table
|
|
// giving parameters for the interpolation operator
|
|
// ... this starts out as a copy of the passed-in table,
|
|
// then gets extra stuff added to it specific to this
|
|
// interpolation region; it's shared across all iperp
|
|
// ... we own this table
|
|
const int interp_par_table_handle_;
|
|
|
|
// (par) origin and delta values of the gridfn data
|
|
const fp gridfn_coord_origin_, gridfn_coord_delta_;
|
|
|
|
// --> start of gridfn data to use for interpolation
|
|
// (reset for each iperp)
|
|
// ... we do *not* own the pointed-to data!
|
|
// ... index is (gfn)
|
|
mutable jtutil::array1d<const void*> gridfn_data_ptrs_;
|
|
|
|
// --> start of interpolation data buffer for each gridfn
|
|
// (reset for each iperp)
|
|
// ... we do *not* own the pointed-to data!
|
|
// ... index is (gfn)
|
|
mutable jtutil::array1d<void*> interp_data_buffer_ptrs_;
|
|
};
|
|
|
|
//******************************************************************************
|
|
|
|
} // namespace AHFinderDirect
|
|
#endif /* TPATCH_INTERP_H */
|