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
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
|
// patch_interp.cc -- patch interpolation region
// $Id$
//
// patch_interp::patch_interp
// patch_interp::~patch_interp
// patch_interp::[min,max]_ipar_for_gridfn_data
// patch_interp::interpolate
//
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "util_Table.h"
#include "cctk.h"
#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;
#include "fp.hh"
#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "patch_interp.hh"
#include "ghost_zone.hh"
//*****************************************************************************
//
// This function constructs an patch_interp object.
//
// Note this 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.
//
// It also requires that the patch's gridfns exist, since we size various
// arrays based on the patch's min/max ghosted gfn.
//
patch_interp::patch_interp(const patch_edge& my_edge_in,
int min_iperp_in, int max_iperp_in,
const jtutil::array1d<int>& min_parindex_used_in,
const jtutil::array1d<int>& max_parindex_used_in,
const jtutil::array2d<fp>& interp_par_in,
int interp_handle_in,
int interp_par_table_handle_in)
: my_patch_(my_edge_in.my_patch()),
my_edge_(my_edge_in),
min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
min_ipar_(min_ipar_for_gridfn_data()),
max_ipar_(max_ipar_for_gridfn_data()),
min_parindex_used_(min_parindex_used_in),
max_parindex_used_(max_parindex_used_in),
interp_par_(interp_par_in),
interp_handle_(interp_handle_in),
interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
gridfn_coord_origin_(my_edge().par_of_ipar(min_ipar())),
gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
gridfn_type_codes_ (my_patch().ghosted_min_gfn(),
my_patch().ghosted_max_gfn()),
gridfn_data_ptrs_ (my_patch().ghosted_min_gfn(),
my_patch().ghosted_max_gfn()),
result_buffer_ptrs_(my_patch().ghosted_min_gfn(),
my_patch().ghosted_max_gfn())
{
int status;
// did interpolator-parameter-table cloning work ok?
status = interp_par_table_handle_;
if (status < 0)
then error_exit(ERROR_EXIT,
"***** patch_interp::patch_interp():\n"
" can't clone interpolator parmameter table (handle=%d)!\n"
" error status=%d\n"
,
interp_par_table_handle_in,
status); /*NOTREACHED*/
// set up gridfn storage addressing in interpolator parameter table
const int N_dims = 1;
const CCTK_INT stride = my_patch().ghosted_iang_stride(my_edge().is_rho());
status = Util_TableSetIntArray(interp_par_table_handle_,
N_dims, &stride,
"input_array_strides");
if (status < 0)
then error_exit(ERROR_EXIT,
"***** patch_interp::patch_interp():\n"
" can't set stride in interpolator parmameter table (handle=%d)!\n"
" error status=%d\n"
,
interp_par_table_handle_,
status); /*NOTREACHED*/
// set up gridfn type codes
for (int gfn = my_patch().ghosted_min_gfn() ;
gfn <= my_patch().ghosted_max_gfn() ;
++gfn)
{
gridfn_type_codes_(gfn) = CCTK_VARIABLE_REAL; // fp == CCTK_REAL
}
}
//*****************************************************************************
//
// This function destroys an patch_interp object
//
patch_interp::~patch_interp()
{
Util_TableDestroy(interp_par_table_handle_);
}
//*****************************************************************************
//
// These functions compute the [min,max] range of ipar in this patch
// from which we may use gridfn data in interpolation. This computation
// depends on our adjacent ghost zones' types as described in the comments
// to our constructor, but needs only my_patch_ and my_edge_ members of
// this object to be valid.
//
int patch_interp::min_ipar_for_gridfn_data() const
{
const ghost_zone& min_par_adjacent_ghost_zone
= my_patch()
.ghost_zone_on_edge(my_edge().min_par_adjacent_edge());
return min_par_adjacent_ghost_zone.is_symmetry()
? my_edge().min_ipar_with_corners()
: my_edge().min_ipar_without_corners();
}
int patch_interp::max_ipar_for_gridfn_data() const
{
const ghost_zone& max_par_adjacent_ghost_zone
= my_patch()
.ghost_zone_on_edge(my_edge().max_par_adjacent_edge());
return max_par_adjacent_ghost_zone.is_symmetry()
? my_edge().max_ipar_with_corners()
: my_edge().max_ipar_without_corners();
}
//*****************************************************************************
//
// 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_interp::interpolate(int ghosted_min_gfn_to_interp,
int ghosted_max_gfn_to_interp,
jtutil::array3d<fp>& result_buffer_array)
const
{
//
// We do a separate interpolation for each iperp.
//
// The interpolator accesses the gridfn data with the generic
// subscripting expression
// data[offset + i*stride]
// where i is a 0-origin point index for each interpolation, and where
// we get to choose the pointer data and the offset value for each
// (gridfn,iperp), and the stride for each iperp. Our strategy is to
// leave offset unspecified so it defaults to 0, use a common stride
// for all iperp (already set up in our constructor), and specify the
// appropriate data pointers (--> the start of the gridfn data we should
// use) for each (gfn,iperp). With this strategy we don't need to modify
// the parameter table at all after the constructor.
//
const int N_dims = 1;
const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
ghosted_max_gfn_to_interp);
const CCTK_INT N_gridfn_data_points
= jtutil::how_many_in_range(min_ipar(), max_ipar());
const int interp_coords_type_code = CCTK_VARIABLE_REAL; // fp == CCTK_REAL
const CCTK_INT *gridfn_type_codes_ptr
= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
//
// do the interpolations at each iperp
//
for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
{
const int min_parindex = min_parindex_used_(iperp);
const int max_parindex = max_parindex_used_(iperp);
// interpolation-point coordinates
const CCTK_INT N_interp_points
= jtutil::how_many_in_range(min_parindex, max_parindex);
const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex);
const void* const interp_coords[N_dims]
= { static_cast<const void *>(interp_coords_ptr) };
// pointers to gridfn data to interpolate, buffer for results
for (int gfn = ghosted_min_gfn_to_interp ;
gfn <= ghosted_max_gfn_to_interp ;
++gfn)
{
const int start_irho
= my_edge(). irho_of_iperp_ipar(iperp, min_ipar());
const int start_isigma
= my_edge().isigma_of_iperp_ipar(iperp, min_ipar());
gridfn_data_ptrs_(gfn)
= static_cast<const void *>(
& my_patch()
.ghosted_gridfn(gfn, start_irho,start_isigma)
);
result_buffer_ptrs_(gfn)
= static_cast<void *>(
& result_buffer_array(gfn, iperp,min_parindex)
);
}
const void* const* const gridfn_data
= & gridfn_data_ptrs_(ghosted_min_gfn_to_interp);
void* const* const result_buffer
= & result_buffer_ptrs_(ghosted_min_gfn_to_interp);
// do the actual interpolation
const int status
= CCTK_InterpLocalUniform(N_dims,
interp_handle_,
interp_par_table_handle_,
&gridfn_coord_origin_,
&gridfn_coord_delta_,
N_interp_points,
interp_coords_type_code,
interp_coords,
N_gridfns,
&N_gridfn_data_points,
gridfn_type_codes_ptr,
gridfn_data,
N_gridfns,
gridfn_type_codes_ptr,
result_buffer);
if (status < 0)
then error_exit(ERROR_EXIT,
"***** patch_interp::interpolate():\n"
" error return %d from interpolator at iperp=%d of [%d,%d]!\n"
,
status, iperp, min_iperp(), max_iperp());
/*NOTREACHED*/
}
}
|
