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
|
/*@@
@file GetAllActive.cc
@date Sun Sep 19 20:34:27 EDT 2010
@author Erik Schnetter, Roland Haas
@desc
a piece of dh::regrid to compute the active region
@enddesc
@@*/
#include <cassert>
#include <vector>
#include "cctk.h"
#include "CarpetIOHDF5.hh"
namespace CarpetIOHDF5 {
/********************************************************************
********************* External Routines **********************
********************************************************************/
void
GetAllActive(const dh *dd, const gh *hh, int ml, int rl, ibset &allactive);
/********************************************************************
********************* Internal Routines ********************
********************************************************************/
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl,
char const * restrict const message);
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl, int const c,
char const * restrict const message);
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl, int const c, int const cc,
char const * restrict const message);
/********************************************************************
********************* Local Data *****************************
********************************************************************/
// accumulates any consistency errors
static bool there_was_an_error = false;
/*@@
@routine
@date Sun Sep 19 20:46:35 EDT 2010
@author Erik Schnetter, Roland Haas
@desc
Output warning for consistency checks for a refinement level.
@enddesc
@history
@endhistory
@var checkstring
@vdesc the check that failed
@vtype char const * restrict const
@vio in
@endvar
@var file
@vdesc name of source file where error occured
@vtype char const * restrict const
@vio in
@endvar
@var line
@vdesc source line in file where error occured
@vtype int const
@vio in
@endvar
@var ml
@vdesc meta level where error occured
@vtype int const
@vio in
@endvar
@var rl
@vdesc refinement level where error occured
@vtype int const
@vio in
@endvar
@var message
@vdesc error message to display
@vtype char const * restrict const
@vio in
@endvar
@returntype void
@@*/
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl,
char const * restrict const message)
{
CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n%s:%d:\n [ml=%d rl=%d] The following grid structure consistency check failed:\n %s\n %s",
file, line, ml, rl, message, checkstring);
there_was_an_error = true;
}
/*@@
@routine
@date Sun Sep 19 20:46:35 EDT 2010
@author Erik Schnetter, Roland Haas
@desc
Output warning for consistency checks for a single component.
@enddesc
@history
@endhistory
@var checkstring
@vdesc the check that failed
@vtype char const * restrict const
@vio in
@endvar
@var file
@vdesc name of source file where error occured
@vtype char const * restrict const
@vio in
@endvar
@var line
@vdesc source line in file where error occured
@vtype int const
@vio in
@endvar
@var ml
@vdesc meta level where error occured
@vtype int const
@vio in
@endvar
@var rl
@vdesc refinement level where error occured
@vtype int const
@vio in
@endvar
@var c
@vdesc component where error occured
@vtype int const
@vio in
@endvar
@var message
@vdesc error message to display
@vtype char const * restrict const
@vio in
@endvar
@returntype void
@@*/
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl, int const c,
char const * restrict const message)
{
CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n%s:%d:\n [ml=%d rl=%d c=%d] The following grid structure consistency check failed:\n %s\n %s",
file, line, ml, rl, c, message, checkstring);
there_was_an_error = true;
}
/*@@
@routine
@date Sun Sep 19 20:46:35 EDT 2010
@author Erik Schnetter, Roland Haas
@desc
Output warning for consistency checks between two components.
@enddesc
@history
@endhistory
@var checkstring
@vdesc the check that failed
@vtype char const * restrict const
@vio in
@endvar
@var file
@vdesc name of source file where error occured
@vtype char const * restrict const
@vio in
@endvar
@var line
@vdesc source line in file where error occured
@vtype int const
@vio in
@endvar
@var ml
@vdesc meta level where error occured
@vtype int const
@vio in
@endvar
@var rl
@vdesc refinement level where error occured
@vtype int const
@vio in
@endvar
@var c
@vdesc first component where error occured
@vtype int const
@vio in
@endvar
@var cc
@vdesc second component where error occured
@vtype int const
@vio in
@endvar
@var message
@vdesc error message to display
@vtype char const * restrict const
@vio in
@endvar
@returntype void
@@*/
static
void
assert_error (char const * restrict const checkstring,
char const * restrict const file, int const line,
int const ml, int const rl, int const c, int const cc,
char const * restrict const message)
{
CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n%s:%d:\n [ml=%d rl=%d c=%d cc=%d] The following grid structure consistency check failed:\n %s\n %s",
file, line, ml, rl, c, cc, message, checkstring);
there_was_an_error = true;
}
#ifdef CARPET_OPTIMISE
// For highest efficiency, omit all self-checks
#define ASSERT_rl(check, message)
#define ASSERT_c(check, message)
#define ASSERT_cc(check, message)
#else
#define ASSERT_rl(check, message) \
do { \
if (not (check)) { \
assert_error (#check, __FILE__, __LINE__, ml, rl, message); \
} \
} while (false)
#define ASSERT_c(check, message) \
do { \
if (not (check)) { \
assert_error (#check, __FILE__, __LINE__, ml, rl, c, message); \
} \
} while (false)
#define ASSERT_cc(check, message) \
do { \
if (not (check)) { \
assert_error (#check, __FILE__, __LINE__, ml, rl, c, cc, message); \
} \
} while (false)
#endif
/*@@
@routine
@date Sun Sep 19 20:46:35 EDT 2010
@author Erik Schnetter, Roland Haas
@desc
Recompute the union of all active regions. This routine is a
stripped down copy of dh::regrid.
@enddesc
@history
@endhistory
@var dh
@vdesc data hirarchy for this group
@vtype const dh *
@vio in
@endvar
@var hh
@vdesc grid hirarchy
@vtype const gh *
@vio in
@endvar
@var ml
@vdesc meta level to compute active region on
@vtype int
@vio in
@endvar
@var rl
@vdesc refinement level to compute active region on
@vtype int
@vio in
@endvar
@var allactive
@vdesc receives the list of active boxes
@vtype ibset &
@vio out
@endvar
@returntype void
@@*/
void
GetAllActive(const dh *dd, const gh *hh, int ml, int rl, ibset &allactive)
{
DECLARE_CCTK_PARAMETERS;
// All owned regions
ibset allowned;
// Domain:
ibbox const & domain_exterior = hh->baseextent(ml,rl);
// Variables may have size zero
// ASSERT_rl (not domain_exterior.empty(),
// "The exterior of the domain must not be empty");
i2vect const & buffer_width = dd->buffer_widths.AT(rl);
i2vect const & boundary_width = hh->boundary_width;
ASSERT_rl (all (all (boundary_width >= 0)),
"The gh boundary widths must not be negative");
ibbox const domain_active = domain_exterior.expand (- boundary_width);
// Variables may have size zero
// ASSERT_rl (not domain_active.empty(),
// "The active part of the domain must not be empty");
ASSERT_rl (domain_active <= domain_exterior,
"The active part of the domain must be contained in the exterior part of the domain");
for (int c = 0; c < hh->components(rl); ++ c) {
// Interior:
ibbox intr;
intr = ibbox::poison();
// The interior of the grid has the extent as specified by the
// regridding thorn
intr = hh->extent (ml,rl,c);
// (The interior must not be empty)
// Variables may have size zero
// ASSERT_c (not intr.empty(),
// "The interior must not be empty");
// The interior must be contained in the domain
ASSERT_c (intr <= hh->baseextent(ml,rl),
"The interior must be contained in the domain");
// All interiors must be disjunct
#ifdef CARPET_DEBUG
for (int cc = 0; cc < c; ++ cc) {
ASSERT_cc (not intr.intersects (full_level.AT(cc).interior),
"All interiors must be disjunct");
}
#endif
// Outer boundary faces:
b2vect is_outer_boundary;
// The outer boundary faces are where the interior extends up
// to the outer boundary of the domain. It is not possible to
// check whether it extends past the active part of the
// domain, since this would be wrong when the outer boundary
// width is zero.
is_outer_boundary[0] = intr.lower() == domain_exterior.lower();
is_outer_boundary[1] = intr.upper() == domain_exterior.upper();
// Owned region:
ibbox owned;
owned = ibbox::poison();
owned = intr.expand (i2vect (is_outer_boundary) * (- boundary_width));
// (The owned region must not be empty)
// Variables may have size zero
// ASSERT_c (not owned.empty(),
// "The owned region must not be empty");
// The owned region must be contained in the active part of
// the domain
ASSERT_c (owned <= domain_active,
"The owned region must be contained in the active part of the domain");
// All owned regions must be disjunct
#ifdef CARPET_DEBUG
for (int cc = 0; cc < c; ++ cc) {
ASSERT_cc (not owned.intersects (full_level.AT(cc).owned),
"All owned regions must be disjunct");
}
#endif
allowned |= owned;
ASSERT_rl (allowned <= domain_active,
"The owned regions must be contained in the active part of the domain");
} // for c
// Enlarge active part of domain
i2vect const safedist = i2vect (0);
ibbox const domain_enlarged = domain_active.expand (safedist);
// All not-owned regions
ibset const notowned = domain_enlarged - allowned;
// All not-active points
ibset const notactive = notowned.expand (buffer_width);
// All active points
allactive = allowned - notactive;
// Fail if any of the consistency checks failed
assert(there_was_an_error == false);
}
} // namespace CarpetIOHDF5
|