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
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
|
// patch_system.hh -- describes the (an) entire system of interlinked patches
// $Id$
//
// patch_system - describes a system of interlinked patches
//
//
// prerequisites:
// <stdio.h>
// <assert.h>
// <math.h>
// "stdc.h"
// "config.hh"
// "../jtutil/util.hh"
// "../jtutil/array.hh"
// "../jtutil/linear_map.hh"
// "coords.hh"
// "grid.hh"
// "fd_grid.hh"
// "patch.hh"
// "patch_edge.hh"
// "ghost_zone.hh"
// "patch_interp.hh"
// "patch_info.hh"
//
//******************************************************************************
//
// A patch_system object describes a system of interlinked patches.
//
// Its const qualifiers refer (only) to the gridfn data. Notably, this
// means that synchronize() is a non-const function (it modifies gridfn
// data), while synchronize_Jacobian() et al are const functions (they
// don't modify gridfn data) even though they may update other internal
// state in the patch_system object and its subobjects.
//
class patch_system
{
//
// ***** static data & functions describing patch systems *****
//
public:
// what patch-system type are supported?
enum patch_system_type {
full_sphere_patch_system,
plus_z_hemisphere_patch_system,
plus_xy_quadrant_patch_system,
plus_xz_quadrant_patch_system,
plus_xyz_octant_patch_system
};
// maximum number of patches in any patch-system type
static const int max_N_patches = 6;
// decode patch system type into N_patches
static
int N_patches_of_type(enum patch_system_type type_in);
// patch system type <--> human-readable character-string name
static
const char* name_of_type(enum patch_system_type type_in);
static
enum patch_system_type type_of_name(const char* name_in);
//
// ***** coordinates *****
//
public:
// global (x,y,z) --> local (x,y,z)
fp local_x_of_global_x(fp global_x) const
{ return global_coords_.local_x_of_global_x(global_x); }
fp local_y_of_global_y(fp global_y) const
{ return global_coords_.local_y_of_global_y(global_y); }
fp local_z_of_global_z(fp global_z) const
{ return global_coords_.local_z_of_global_z(global_z); }
// local (x,y,z) --> global (x,y,z)
fp global_x_of_local_x(fp local_x) const
{ return global_coords_.global_x_of_local_x(local_x); }
fp global_y_of_local_y(fp local_y) const
{ return global_coords_.global_y_of_local_y(local_y); }
fp global_z_of_local_z(fp local_z) const
{ return global_coords_.global_z_of_local_z(local_z); }
// ... get global (x,y,z) coordinates of local origin point
fp origin_x() const { return global_coords_.origin_x(); }
fp origin_y() const { return global_coords_.origin_y(); }
fp origin_z() const { return global_coords_.origin_z(); }
//
// ***** meta-info about the entire patch system *****
//
public:
// patch-system type
enum patch_system_type type() const { return type_; }
// total number of patches
int N_patches() const { return N_patches_; }
// get patches by patch number
patch &ith_patch(int pn) const
{ return * all_patches_[pn]; }
// find a patch by name, return patch number; error_exit() if not found
int patch_number_of_name(const char* name) const;
// total number of grid points
int N_grid_points() const { return N_grid_points_; }
int ghosted_N_grid_points() const { return ghosted_N_grid_points_; }
//
// ***** meta-info about gridfns *****
//
public:
int min_gfn() const { return ith_patch(0).min_gfn(); }
int max_gfn() const { return ith_patch(0).max_gfn(); }
int N_gridfns() const { return ith_patch(0).N_gridfns(); }
int ghosted_min_gfn() const { return ith_patch(0).ghosted_min_gfn(); }
int ghosted_max_gfn() const { return ith_patch(0).ghosted_max_gfn(); }
int ghosted_N_gridfns() const
{ return ith_patch(0).ghosted_N_gridfns(); }
//
// ***** synchronize() and its Jacobian *****
//
public:
// "synchronize" all ghost zones of all patches,
// i.e. update the ghost-zone values of the specified gridfns
// via the appropriate sequence of symmetry operations
// and interpatch interpolations
void synchronize(int ghosted_min_gfn_to_sync,
int ghosted_max_gfn_to_sync);
// ... do this for all ghosted gridfns
void synchronize()
{
synchronize(ghosted_min_gfn(),
ghosted_max_gfn());
}
//
// compute Jacobian of synchronize() into internal buffers,
// taking into account synchronize()'s full 3-phase algorithm
//
void compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync,
int ghosted_max_gfn_to_sync)
const;
// ... do this for all ghosted gridfns
void compute_synchronize_Jacobian()
const
{
compute_synchronize_Jacobian(ghosted_min_gfn(),
ghosted_max_gfn());
}
//
// The following functions access the Jacobian computed by
// compute_synchronize_Jacobian() . Note this API is rather
// different than that of ghost_zone::comute_Jacobian() et al:
// here we must take into account synchronize()'s full 3-phase
// algorithm, and this may lead to a more general Jacobian
// structure.
//
// This API still implicitly assumes that the Jacobian is
// independent of ghosted_gfn , and that the set of y points
// (with nonzero Jacobian values) in a single row of the Jacobian
// matrix (i.e. the set of points on which a single ghost-zone
// point depends),
// - lies entirely within a single y patch
// - has a single yiperp value
// - have a contiguous interval of yipar; we parameterize this
// interval as yipar = posn+m
//
// what are the global min/max m over all ghost zone points?
// (this is useful for sizing the buffer for synchronize_Jacobian())
void synchronize_Jacobian_global_minmax_y_ipar_m
(int& min_m, int& max_m)
const;
// to which patch/edge do the y points in a Jacobian row belong?
patch& synchronize_Jacobian_y_patch(const ghost_zone& xgz,
int x_iperp, int x_ipar)
const;
const patch_edge& synchronize_Jacobian_y_edge(const ghost_zone& xgz,
int x_iperp, int x_ipar)
const;
// what is the iperp of the Jacobian y points in their (y) patch?
int synchronize_Jacobian_y_iperp(const ghost_zone& xgz,
int x_iperp, int x_ipar)
const;
// what are the posn and min/max m for a given x point?
int synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz,
int x_iperp, int x_ipar)
const;
void synchronize_Jacobian_minmax_y_ipar_m(const ghost_zone& xgz,
int x_iperp, int x_ipar,
int& min_m, int& max_m)
const;
public:
// store the Jacobian
// partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
// -------------------------------------------------------------
// partial gridfn(ghosted_gfn, py, y_iperp, y_posn+y_ipar_m)
// (taking into account synchronize()'s full 3-phase algorithm)
// in the caller-supplied buffer
// Jacobian_buffer(m)
// for each m , where
// y_iperp = Jacobian_y_iperp(x_iperp)
// y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
void synchronize_Jacobian(const ghost_zone& xgz,
int x_iperp, int x_ipar,
jtutil::array1d<fp>& Jacobian_buffer)
const;
//
// ***** gridfn operations *****
//
// dst = src
void gridfn_copy(int src_gfn, int dst_gfn);
// dst += delta
void add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn);
//
// ***** I/O *****
//
public:
// print to a named file (newly (re)created)
// output format is
// dpx dpy gridfn
void print_gridfn(int gfn, const char output_file_name[]) const
{
print_unknown_gridfn(false, gfn,
false, false, 0,
output_file_name, false);
}
void print_ghosted_gridfn(int ghosted_gfn,
const char output_file_name[],
bool want_ghost_zones = true)
const
{
print_unknown_gridfn(true, ghosted_gfn,
false, false, 0,
output_file_name, want_ghost_zones);
}
// print to a named file (newly (re)created)
// output format is
// dpx dpy gridfn global_x global_y global_z
// where global_[xyz} are derived from the angular position
// and a specified (unknown-grid) radius gridfn
void print_gridfn_with_xyz
(int gfn,
bool radius_is_ghosted_flag, int unknown_radius_gfn,
const char output_file_name[])
const
{
print_unknown_gridfn(false, gfn,
true, radius_is_ghosted_flag,
unknown_radius_gfn,
output_file_name, false);
}
void print_ghosted_gridfn_with_xyz
(int ghosted_gfn,
bool radius_is_ghosted_flag, int unknown_radius_gfn,
const char output_file_name[],
bool want_ghost_zones = true)
const
{
print_unknown_gridfn(true, ghosted_gfn,
true, radius_is_ghosted_flag,
unknown_radius_gfn,
output_file_name, want_ghost_zones);
}
public:
// read from a named file
void read_gridfn(int gfn, const char input_file_name[])
{ read_unknown_gridfn(false, gfn, input_file_name, false); }
void read_ghosted_gridfn(int ghosted_gfn,
const char input_file_name[],
bool want_ghost_zones = true)
{
read_unknown_gridfn(true, ghosted_gfn,
input_file_name, want_ghost_zones);
}
private:
// ... internal worker functions
void print_unknown_gridfn
(bool ghosted_flag, int unknown_gfn,
bool print_xyz_flag, bool radius_is_ghosted_flag,
int unknown_radius_gfn,
const char output_file_name[], bool want_ghost_zones)
const;
void read_unknown_gridfn(bool ghosted_flag, int unknown_gfn,
const char input_file_name[],
bool want_ghost_zones);
//
// ***** access to gridfns as 1-D arrays *****
//
// ... n.b. this interface implicitly assumes that gridfn data
// arrays are contiguous across patches; this is ensured by
// setup_gridfn_storage() (called by our constructor)
//
public:
// convert (patch,irho,isigma) <--> 1-D 0-origin grid point number (gpn)
int gpn_of_patch_irho_isigma(const patch& p, int irho, int isigma)
const
{
return starting_gpn_[p.patch_number()]
+ p.gpn_of_irho_isigma(irho,isigma);
}
int ghosted_gpn_of_patch_irho_isigma(const patch& p,
int irho, int isigma)
const
{
return ghosted_starting_gpn_[p.patch_number()]
+ p.ghosted_gpn_of_irho_isigma(irho,isigma);
}
// ... n.b. we return patch as a reference via the function result;
// an alternative would be to have a patch*& argument
patch& patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
const;
patch& ghosted_patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
const;
// access actual gridfn data arrays
// (low-level, dangerous, use with caution)
const fp* gridfn_data(int gfn) const
{ return ith_patch(0).gridfn_data_array(gfn); }
fp* gridfn_data(int gfn)
{ return ith_patch(0).gridfn_data_array(gfn); }
const fp* ghosted_gridfn_data(int ghosted_gfn) const
{ return ith_patch(0).ghosted_gridfn_data_array(ghosted_gfn); }
fp* ghosted_gridfn_data(int ghosted_gfn)
{ return ith_patch(0).ghosted_gridfn_data_array(ghosted_gfn); }
//
// ***** constructor, destructor *****
//
// This constructor doesn't support the full generality of the
// patch data structures (which would, eg, allow N_ghost_points
// and N_extend_points and the interpolator parameters to vary
// from ghost zone to ghost zone, and the grid spacings to vary
// from patch to patch. But in practice we'd probably never
// use that generality...
//
public:
patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
enum patch_system_type type_in,
int N_ghost_points, int N_overlap_points,
fp delta_drho_dsigma,
int min_gfn_in, int max_gfn_in,
int ghosted_min_gfn_in, int ghosted_max_gfn_in,
int interp_handle_in, int interp_par_table_handle_in);
~patch_system();
//
// ***** helper functions for constructor *****
//
private:
// construct patches as described by patch_info[] array,
// and link them into the patch system
// does *NOT* create ghost zones
// does *NOT* set up gridfns
void create_patches(const struct patch_info patch_info_in[],
int N_ghost_points, int N_extend_points,
fp delta_drho_dsigma);
// setup all gridfns with contiguous-across-patches storage
void setup_gridfn_storage
(int min_gfn_in, int max_gfn_in,
int ghosted_min_gfn_in, int ghosted_max_gfn_in);
// create/interlink ghost zones
void interlink_full_sphere_patch_system
(int N_overlap_points,
int interp_handle, int interp_par_table_handle);
void interlink_plus_z_hemisphere_patch_system
(int N_overlap_points,
int interp_handle, int interp_par_table_handle);
void interlink_plus_xy_quadrant_patch_system
(int N_overlap_points,
int interp_handle, int interp_par_table_handle);
void interlink_plus_xz_quadrant_patch_system
(int N_overlap_points,
int interp_handle, int interp_par_table_handle);
void interlink_plus_xyz_octant_patch_system
(int N_overlap_points,
int interp_handle, int interp_par_table_handle);
// create/interlink a pair of periodic-symmetry ghost zones
static
void create_periodic_symmetry_ghost_zones
(const patch_edge& ex, const patch_edge& ey,
bool ipar_map_is_plus);
// construct a pair of interpatch ghost zones
// ... automagically figures out which edges are adjacent
static
void create_interpatch_ghost_zones(patch &px, patch &py,
int N_overlap_points);
// finish setup of a pair of interpatch ghost zones
// ... automagically figures out which edges are adjacent
static
void finish_interpatch_setup
(patch &px, patch &py,
int N_overlap_points,
int interp_handle, int interp_par_table_handle);
// assert() that all ghost zones of all patches are fully setup
void assert_all_ghost_zones_fully_setup() const;
private:
// we forbid copying and passing by value
// by declaring the copy constructor and assignment operator
// private, but never defining them
patch_system(const patch_system &rhs);
patch_system& operator=(const patch_system &rhs);
private:
// local <--> global coordinate mapping
global_coords global_coords_;
// meta-info about patch system
enum patch_system_type type_;
int N_patches_;
int N_grid_points_, ghosted_N_grid_points_;
// [pn] = --> individual patches
std::vector<patch *> all_patches_;
// [pn] = starting grid point number of individual patches
// ... arrays are actually of size N_patches_+1, the [N_patches_]
// entries are == N_grid_points_ and ghosted_N_grid_points_
std::vector<int> starting_gpn_;
std::vector<int> ghosted_starting_gpn_;
// pointers to storage blocks for all gridfns
// ... patches point into these, but we own the storage blocks
fp* gridfn_storage_;
fp* ghosted_gridfn_storage_;
// min/max m over all ghost zone points
mutable int global_min_m_, global_max_m_;
};
|