aboutsummaryrefslogtreecommitdiff
path: root/src/patch/patch_interp.hh
blob: c0c6dcbd644afcac6c55abff2f56d77a9e806562 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
// patch_frontier.hh -- interpolate data for another patch's ghost zone
// $Id$

//
// prerequisites:
//	<stdio.h>
//	<assert.h>
//	<math.h>
//
//	"cctk.h"
//
//	"jt/util.hh"
//	"jt/array.hh"
//	"jt/linear_map.hh"
//	"jt/interpolate.hh"
//	fp.hh
//	coords.hh
//	grid.hh
//	fd_grid.hh
//	patch.hh
//	patch_edge.hh
//

//*****************************************************************************

//
// patch_frontier - interpolation from just inside patch edge
//

//
// A patch_frontier 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).
//

//
// 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.
//
// 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.
//
// 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
	{
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_; }

	// min/max (iperp,ipar) of the gridfn data we'll use
	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_; }

	//
	// ***** 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
	// ... result buffer indices are (gfn, iperp,parindex)
	void interpolate(int ghosted_min_gfn_to_interp,
			 int ghosted_max_gfn_to_interp,
			 jtutil::array3d<fp>& result_buffer);


	//
	// ***** internal functions *****
	//
private:
	// range of ipar in this patch from which we may use gridfn data
	// in interpolation -- depends on adjacent ghost zones' types as
	// described in comments to constructor (above), but needs only
	// my_patch_ and my_edge_ members of this object to be valid
	int min_ipar_for_gridfn_data() const;
	int max_ipar_for_gridfn_data() const;


	//
	// ***** constructor, destructor *****
	//
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).
	// [min,max]_parindex_used_in(iperp)
	//	= [min,max] range of parindex actually used at each iperp.
	// 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).
	//	  We keep a reference to the passed-in  interp_par_in
	//	  array, so this should have a lifetime at last as long
	//	  as that of this object.
	// 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.
	//	  This table is not modified by any actions of this
	//	  class.
	//
	// 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,
		       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 interp_handle_in, int interp_par_table_handle_in);
	~patch_frontier();

private:
        // we forbid copying and passing by value
        // by declaring the copy constructor and assignment operator
        // private, but never defining them
	patch_frontier(const patch_frontier& rhs);
	patch_frontier& operator=(const patch_frontier& rhs);

private:
	const patch& my_patch_;
	const patch_edge& my_edge_;

	// range of (iperp,ipar) in this patch from which
	// we will use gridfn data in interpolation
	int min_iperp_, max_iperp_;
	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_used_;
	const jtutil::array1d<int>& max_parindex_used_;

	// 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
	//     frontier; it's shared across all iperp
	// ... we own this table
	int interp_par_table_handle_;

	// (par) origin and delta values of the gridfn data
	fp gridfn_coord_origin_, gridfn_coord_delta_;

	// gridfn type codes for interpolator
	// ... must be CCTK_INT so we can pass by reference to interpolator
	// ... must be mutable so we can set values in ctor
	// ... index is (gfn)
	mutable jtutil::array1d<CCTK_INT> gridfn_type_codes_;

	// --> 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<fp *> gridfn_data_ptrs_;

	// --> start of result buffer for interpolation
	//     (reset for each iperp)
	// ... we do *not* own the pointed-to data!
	// ... index is (gfn)
	mutable jtutil::array1d<fp *> result_buffer_ptrs_;
	}