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
|
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <cmath>
#include <carpet.hh>
#include <mpi_string.hh>
#include "boundary.hh"
namespace CarpetRegrid2 {
using namespace std;
using namespace Carpet;
void evaluate_level_mask (cGH const* restrict const cctkGH,
vector<ibset>& regions, int const rl)
{
DECLARE_CCTK_PARAMETERS;
assert (is_singlemap_mode());
int const m = Carpet::map;
assert (rl > 0);
BEGIN_GLOBAL_MODE (cctkGH) {
ENTER_LEVEL_MODE (cctkGH, rl-1) {
ENTER_SINGLEMAP_MODE (cctkGH, m, CCTK_GF) {
if (verbose or veryverbose) {
cout << "Refinement level " << rl << ":\n";
}
gh const& hh = *vhh.AT(Carpet::map);
ivect const gsh = ivect::ref(cctkGH->cctk_gsh);
ivect const nghostzones = ivect::ref(cctkGH->cctk_nghostzones);
if (veryverbose) {
cout << " gsh: " << gsh << "\n"
<< " nghostzones: " << nghostzones << "\n";
}
// Determine the block size
ivect block_size = adaptive_block_size;
if (adaptive_block_size_x > 0) block_size[0] = adaptive_block_size_x;
if (adaptive_block_size_y > 0) block_size[1] = adaptive_block_size_y;
if (adaptive_block_size_z > 0) block_size[2] = adaptive_block_size_z;
assert (all (block_size > 0));
// Determine the block offset, i.e. the grid point index of
// the coordinate origin. We don't start blocks at index
// zero, so that one can change the location of the outer
// boundary without influencing the refinement mechanism.
// Instead, we align the block such that one of the blocks
// begins at the coordinate origin.
ivect origin;
rvect rorigin;
{
DECLARE_CCTK_ARGUMENTS;
for (int d=0; d<dim; ++d) {
// We round up, because the origin is located between
// grid points for cell-centred grids, and we want the
// next upper grid point in this case. We subtract 0.25
// just for proper rounding.
origin[d] =
ceil(- CCTK_ORIGIN_SPACE(d) / CCTK_DELTA_SPACE(d) - 0.25);
rorigin[d] =
CCTK_ORIGIN_SPACE(d) + origin[d] * CCTK_DELTA_SPACE(d);
}
}
// block_index := ([0...gsh] + block_offset) / block_size
// [0...gsh] = block_index * block_size - block_offset
ivect const block_offset =
(block_size - origin % block_size) % block_size;
ivect const num_blocks =
(gsh + block_offset + block_size - 1) / block_size;
if (veryverbose) {
cout << " origin: " << origin << "\n"
<< " coordinates of origin: " << rorigin << "\n"
<< " offset: " << block_offset << "\n"
<< " num blocks: " << num_blocks << "\n";
}
assert (all (num_blocks > 0));
// Mask containing all refined block indices (in arbitrary
// order, and potentially even duplicated). This is really a
// set, but we need to communicate this data structure, and
// therefore a vector is a better choice.
vector<ivect> mask;
// For vertex centred grids, the blocks need to overlap by
// one grid point, so that e.g. an 8^3 block is turned into
// a 9^3 block. This will ensure symmetry of refined regions
// in the simulation domain.
int const overlap = hh.refcent == vertex_centered ? 1 : 0;
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
DECLARE_CCTK_ARGUMENTS;
if (veryverbose) {
cout << " component " << component << ":\n";
}
ivect const lbnd = ivect::ref(cctk_lbnd);
ivect const lsh = ivect::ref(cctk_lbnd);
ivect const bboxlo (cctk_bbox[0], cctk_bbox[2], cctk_bbox[4]);
ivect const bboxhi (cctk_bbox[1], cctk_bbox[3], cctk_bbox[5]);
ivect const imin = 0 + either(bboxlo, 0, nghostzones);
ivect const imax = lsh - either(bboxhi, 0, nghostzones);
ivect const bmin =
max (0,
((lbnd + imin + block_offset + block_size - overlap) /
block_size) -
1);
ivect const bmax = (lbnd + imax + block_offset) / block_size;
// Loop over all blocks
int nblocks = 0;
int nrefined = 0;
for (int bk=bmin[2]; bk<bmax[2]; ++bk) {
for (int bj=bmin[1]; bj<bmax[1]; ++bj) {
for (int bi=bmin[0]; bi<bmax[0]; ++bi) {
ivect const bind (bi,bj,bk);
ivect const bimin =
max(imin,
(bind ) * block_size - block_offset - lbnd );
ivect const bimax =
min(imax,
(bind+1) * block_size - block_offset - lbnd + overlap);
bool refine = false;
bool have_nan = false;
// Loop over all points in this block
for (int k=bimin[2]; k<bimax[2]; ++k) {
for (int j=bimin[1]; j<bimax[1]; ++j) {
for (int i=bimin[0]; i<bimax[0]; ++i) {
int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
if (not isnan(level_mask[ind3d])) {
refine = refine or level_mask[ind3d] >= rl;
} else {
have_nan = true;
}
}
}
}
// Refine this block if any point in this block requires
// refinement
if (have_nan) {
cout << " *** found nan in block " << bind << " ***\n";
}
if (refine) {
if (veryverbose) {
cout << " refining block " << bind << "\n";
}
mask.push_back(bind);
++ nrefined;
}
++ nblocks;
}
}
}
if (veryverbose) {
cout << " refining " << nrefined << " "
<< "of " << nblocks << " blocks on this component\n";
}
} END_LOCAL_COMPONENT_LOOP;
// Combine this mask from all processes, and to all processes
vector<ivect> const fullmask = allgatherv1 (dist::comm(), mask);
// Convert vector into ibset
const ibbox& cbaseext = hh.baseextent(mglevel, rl-1);
const ibbox& baseext = hh.baseextent(mglevel, rl);
ibset& region = regions.at(rl);
for (vector<ivect>::const_iterator
it = fullmask.begin(); it != fullmask.end(); ++it)
{
ivect const& bind = *it;
ivect const cstr = cbaseext.stride();
ivect const clo =
cbaseext.lower() + (bind * block_size - block_offset) * cstr;
ivect const cup = clo + (block_size + overlap - 1) * cstr;
ibbox const cblock (clo, cup, cstr);
ibbox const block = cblock.expanded_for(baseext);
region |= block;
}
if (verbose or veryverbose) {
if (veryverbose) {
cout << " refined region is " << region << "\n"
<< " domain is " << baseext << "\n"
<< " next coarser level is " << regions.at(rl-1) << "\n";
}
ibbox::size_type const nrefined = region.size();
ibbox::size_type const npoints = baseext.size();
ibbox::size_type const ncoarse =
regions.at(rl-1).expanded_for(baseext).size();
cout << " refining "
<< (100.0*nrefined/ncoarse) << "% of the next coarser level\n"
<< " refining "
<< (100.0*nrefined/npoints) << "% of the domain\n";
}
} LEAVE_SINGLEMAP_MODE;
} LEAVE_LEVEL_MODE;
} END_GLOBAL_MODE;
}
} // namespace CarpetRegrid2
|
