diff options
Diffstat (limited to 'src/patch/patch_interp.hh')
| -rw-r--r-- | src/patch/patch_interp.hh | 146 |
1 files changed, 87 insertions, 59 deletions
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index b1aea5f..19d6f54 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -34,15 +34,20 @@ // The way our coordinates 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 do interpolation in par. -// We implement this in the obvious way, keeping (a pointer to) a -// separate interpolation object for each (gridfn,iperp). // -// FIXME: -// Right now the choice of interpolator is hard-wired, and the set of -// parameters it takes (the order) is hard-wired into our interface -// and thus propagates up to our caller's code. +// 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. +// However, we do assume that its range is roughly comparable at all +// values of iperp, so there's no major waste in treating (iperp,parindex) +// as a 2-D index space for the interpolation points. +// +// 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.] // class patch_frontier @@ -52,46 +57,70 @@ public: const patch& my_patch() const { return my_patch_; } const patch_edge& my_edge() const { return my_edge_; } - // set up a specified gridfn's interpolation state - // ... this copies gridfn data if needed, - // and sets up any interpolator state - // ... it looks at what type of ghost zones are on the - // adjacent edges to decide how to handle the corners: - // * if an adjacent ghost zone is a symmetry ghost zone - // we assume it's already been mirrored by the time - // we get set up ==> we can (and do) use the data from it - // * if an adjacent ghost zone is an interpatch ghost zone, - // we can't assume it's already been interpatch-interpolated - // by the time we're called ==> we don't use data from it - void setup_interpolation_for_gridfn(int gfn); - - // interpolate a specified gridfn - // ... setup_interpolation_for_gridfn(gfn) must have already been called - fp interpolate(int gfn, int iperp, fp par) const; - -private: // min/max (iperp,ipar) of the frontier int min_iperp() const { return min_iperp_; } int max_iperp() const { return max_iperp_; } int min_ipar() const { return min_ipar_; } int max_ipar() const { return max_ipar_; } - // for reasons described in setup_interpolation_for_gridfn(), - // we use 0-origin "relative ipar" coordinates - int ripar_of_ipar(int ipar) const { return ipar - min_ipar(); } - int ipar_of_ripar(int ripar) const { return ripar + min_ipar(); } - int min_ripar() const { return 0; } - int max_ripar() const { return ripar_of_ipar(max_ipar()); } - + // + // ***** main client interface ***** + // +public: + // interpolate specified range of ghosted gridfns + // at all the coordinates specified when we were constructed, + // store interpolated data in result_buffer(gfn, iperp,parindex) + void interpolate(int ghosted_min_gfn, ghosted_max_gfn, + jtutil::array3d<fp>& result_buffer); + +private: + // + // ***** constructor, destructor ***** + // public: - // constructor, destructor - // ... constructor requires that the adjacent-side ghost_zone - // objects already exist, since they're used in determining - // the ipar range over which the interpolation will be done + // + // Constructor arguments: + // my_edge_in = identifies the patch/edge to which this frontier + // is to belong + // [min,max]_iperp_in = the range of iperp for this frontier + // [min,max]_parindex_in = extreme [min,max] range of parindex + // (for sizing (iperp,parindex) arrays etc) + // [min,max]_parindex_used_in(iperp) + // = [min,max] range of parindex actually used at each iperp + // 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) + // ghosted_[min,max]_gfn = the largest ghosted_gfn range for this + // frontier; any given interpolate() call + // may specify a subrange + // interp_handle_in = Cactus handle to the interpolation operator + // interp_par_table_handle_in + // = Cactus handle to a Cactus key/value table giving + // parameters (eg order) for the interpolation operator + // + // ... constructor looks at what type of ghost zones are on the + // adjacent edges to decide how to handle the corners: + // * if an adjacent ghost zone is a symmetry ghost zone + // we assume it's already been mirrored by the time + // we get set up ==> we can (and do) use the data from it + // * if an adjacent ghost zone is an interpatch ghost zone, + // we can't assume it's already been interpatch-interpolated + // by the time we're called ==> we don't use data from it + // ==> constructor requires that the adjacent-side ghost_zone + // objects already exist! + // patch_frontier(const patch_edge& my_edge_in, int min_iperp_in, int max_iperp_in, - int interpolator_order); + int min_parindex_in, int max_parindex_in, + const jtutil::array1d<fp>& min_parindex_used_in, + const jtutil::array1d<fp>& max_parindex_used_in, + const jtutil::array2d<fp>& interp_par_in, + int ghosted_min_gfn_in, int ghosted_max_gfn_in, + int interp_handle_in, + int interp_par_table_handle_in); ~patch_frontier(); private: @@ -102,29 +131,28 @@ private: patch_frontier& operator=(const patch_frontier& rhs); private: - const patch& my_patch_; const patch_edge& my_edge_; + const patch& my_patch_; - // range of iperp,ipar from which we can use data int min_iperp_, max_iperp_; - int min_ipar_, max_ipar_; - - // maps ripar <--> par - linear_map<fp>* ripar_map_; - - // the interpolator interface we're using requires that - // the gridfn data to be interpolated (i.e. for any given - // gridfn and iperp) must be in a contiguous 1-D data array, - // so we have to copy from the patch data into this buffer - // ... subscripts are (gfn,iperp,ipar) - // (this is contiguous in ipar as required) - array3d<fp>* gridfn_buffer_; - - // --> array of --> to interpolator objects - // ... subscripts are (gfn,iperp) - // ... we own the pointed-to objects - // ... interpolator objects use ripar as integer coordinate; - // see setup_interpolation_for_gridfn() for discussion of this - typedef interpolation::Lagrange_uniform_1d<fp> interpolator_t_; - array2d<interpolator_t_*>* interpolator_; + + // range of ghosted gfns for which we may interpolate + // (an actual interpolation may cover a subrange of this) + int ghosted_min_gfn_, ghosted_max_gfn_; + + // Cactus handle to the interpolation operator + int interp_handle_; + + // Cactus handle to a Cactus key/value table + // giving parameters (eg order) for the interpolation operator + int interp_par_table_handle_; + + // (par) coordinate origin and delta values for the interpolator + fp interp_coord_origin_, interp_coord_delta_; + + // par coordinates at which we are to interpolate + // ... interpolation treats this as a 1-D array; + // we set up values in 3-D + // ... subscripts are (gfn, iperp,parindex) + jtutil::array3d<fp> interp_par_coords_; }; |