From 09a823a2175e081adbea6d2d265b73d2de61a9dc Mon Sep 17 00:00:00 2001 From: jthorn Date: Sat, 30 Mar 2002 17:35:32 +0000 Subject: further work on new design to use Cactus local interp git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@403 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/patch/patch_interp.cc | 96 ++++++++++++++------------------------ src/patch/patch_interp.hh | 115 ++++++++++++++++++++++++++++++++-------------- 2 files changed, 115 insertions(+), 96 deletions(-) diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 3335897..e21e09c 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -50,18 +50,33 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, int ghosted_min_gfn_in, int ghosted_max_gfn_in, int interp_handle_in, int interp_par_table_handle_in) - : my_edge_(my_edge_in), - my_patch_(my_edge_in.my_patch()), + : my_patch_(my_edge_in.my_patch()), + my_edge_(my_edge_in), min_iperp_(min_iperp_in), max_iperp_(max_iperp_in), + min_ipar_(my_patch() + .ghost_zone_on_edge(my_edge_in.min_par_adjacent_edge()) + .is_symmetry() + ? my_edge_in.min_ipar_with_corners() + : my_edge_in.min_ipar_without_corners()), + max_ipar_(my_patch() + .ghost_zone_on_edge(my_edge_in.max_par_adjacent_edge()) + .is_symmetry() + ? my_edge_in.max_ipar_with_corners() + : my_edge_in.max_ipar_without_corners()) + min_parindex_used_(min_parindex_used_in), + max_parindex_used_(max_parindex_used_in), + interp_par_(interp_par_in), ghosted_min_gfn_(ghosted_min_gfn_in), ghosted_max_gfn_(ghosted_max_gfn_in), interp_handle_(interp_handle_in), - interp_par_table_handle_(interp_par_table_handle_in), + interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)), interp_coord_origin_(0.0), // computed later interp_coord_delta_(my_patch().delta_ang(my_edge().par_is_rho())), - interp_par_coords_(ghosted_min_gfn_, ghosted_max_gfn_, - min_iperp_, max_iperp_, - min_parindex_in, max_parindex_in) + gridfn_data_ptrs_(ghosted_min_gfn_in, ghosted_max_gfn_in, + min_iperp_in, max_iperp_in), + gridfn_offsets_(ghosted_min_gfn_in, ghosted_max_gfn_in, + min_iperp_in, max_iperp_in), + gridfn_par_stride_(0) // computed later { } @@ -72,67 +87,24 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, // patch_frontier::~patch_frontier() { -// delete the interpolator objects - for (int gfn = my_patch().min_gfn() ; - gfn <= my_patch().max_gfn() ; - ++gfn) - { - for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) - { - delete (*interpolator_)(gfn,iperp); - } - } - -delete interpolator_; -delete gridfn_buffer_; -delete ripar_map_; +Util_TableDestroy(interp_par_table_handle_); } //***************************************************************************** // -// This function sets up a specified gridfn's interpolation state, i.e. -// it copies the gridfn data to the contiguous-in-ipar buffer and sets -// up any interpolator state +// This function interpolates the specified range of ghosted gridfns +// at all the coordinates specified when we were constructed. +// It stores the interpolated data in result_buffer(gfn, iperp,parindex) . +// It calls error_exit() if the point is out of range or any other +// interpolator error is detected. // -void patch_frontier::setup_interpolation_for_gridfn(int gfn) -{ - for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) - { - // - // copy this row of gridfn data from the (probably non-contiguous) - // gridfn array, into our (contiguous-in-ipar) scratch array - // - for (int ipar = min_ipar() ; ipar <= max_ipar() ; ++ipar) - { - int irho = my_edge(). irho_of_iperp_ipar(iperp, ipar); - int isigma = my_edge().isigma_of_iperp_ipar(iperp, ipar); - (*gridfn_buffer_)(gfn,iperp,ipar) - = my_patch().gridfn(gfn,irho,isigma); - } - - // - // Now we want to set up the interpolator for the just-copied row - // of data. The interpolator interface requires that we pass a - // pointer to the [0] data, where the [] subscripting refers to - // the interpolator's integer coordinate. Since ipar=0 may not - // even exist, we introduce ripar := ipar - min_ipar() and use - // this as the interpolator's integer coordinate, so that [0] - // refers to the start of this ipar row in the gridfn buffer. - // - const fp* row_ptr = & (*gridfn_buffer_)(gfn,iperp,min_ipar()); - interpolator_t_* I = (*interpolator_)(gfn,iperp); - I->setup_data_y(row_ptr); - } -} - -//****************************************************************************** - -// -// This function does the actual interpolation of data for a specified -// gridfn to a specified (iperp,par) point. It calls error_exit() -// if the point is out of range or any other interpolator error is found. -// -fp patch_frontier::interpolate(int gfn, int iperp, fp par) const +fp patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn, + jtutil::array3d& result_buffer) const { const int N_dims = 1; + +const int status + = CCTK_InterpLocalUniform(N_dims, + interp_handle_, + interp_par_table_handle_, diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index 19d6f54..39bdb13 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -39,9 +39,9 @@ // 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. +// However, we do assume that its range is not too large, and is roughly +// comparable at all values of iperp, so it's ok to treat (iperp,parindex) +// as a 2-D rectangular index space for the interpolation points. // // We use the Cactus local interpolator CCTK_InterpLocalUniform() // to do the interpolation. To minimize interpolator overheads, we @@ -81,36 +81,47 @@ private: public: // // 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) + // 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; + // = [min,max] range of parindex actually used at each iperp. + // See the note below about the required lifetime of this + // array. + // 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 (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 + // par coordinates to be independent of ghosted_gfn). + // See the note below about the required lifetime of this + // array. + // 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 + // 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 + // parameters (eg order) for the interpolation operator. + // This table is not modified by any actions of this + // class. // - // ... 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! + // This constructor keeps references to the passed-in + // [min,max]_parindex_used_in and interp_par_in arrays, so + // these should have lifetimes at least as long as that of + // this patch_frontier object. + // + // Note that this 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 + // Thus the 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, @@ -131,11 +142,29 @@ private: patch_frontier& operator=(const patch_frontier& rhs); private: - const patch_edge& my_edge_; const patch& my_patch_; + const patch_edge& my_edge_; int min_iperp_, max_iperp_; + // range of ipar in this patch from which we may use gridfn + // data in interpolation -- depends on adjacent ghost zones' + // type as described in comments to constructor (above) + int min_ipar_, max_ipar_; + + // [min,max]_parindex_used_(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). + // ... these are references to arrays passed in constructor + // ==> we do *not* own these + const jtutil::array1d& min_parindex_used_; + const jtutil::array1d& max_parindex_used_; + const jtutil::array2d& interp_par_; // indices (iperp,parindex) + // 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_; @@ -143,16 +172,34 @@ private: // 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 + // Cactus handle to our private Cactus key/value table + // giving parameters for the interpolation operator + // (we *do* own this -- it's a copy of the passed-in table) 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 interp_par_coords_; + // + // the interpolator accesses the gridfn data with the generic + // subscripting expression + // data[offset + parindex*stride] + // where we get to choose the pointer data and the offset + // value for each (gridfn,iperp), and the stride for each iperp. + // + // our strategy will be to choose offset so that the [] + // expression is 0 for parindex = min_parindex_used_(iperp) + // + + // --> start of gridfn data we can use for each iperp + // subscripts are (gfn,iperp) + jtutil::array2d gridfn_data_ptrs_; + + // offsets for gridfn subscripting by interpolator + // subscripts are (gfn,iperp) + jtutil::array2d gridfn_offsets_; + + // stride of par coordinate in gridfn data + // (same for all gfn and iperp) + int gridfn_par_stride_; }; -- cgit v1.2.3