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
|
#include <cassert>
#include <istream>
#include <sstream>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <dh.hh>
#include <carpet.hh>
#include <loopcontrol.h>
#include "bits.h"
#define LOOP_OVER_NEIGHBOURS(dir) \
{ \
ivect dir_(-1); \
do { \
ivect const& dir = dir_; \
{
#define END_LOOP_OVER_NEIGHBOURS \
} \
for (int d_=0; d_<dim; ++d_) { \
if (dir_[d_] < +1) { \
++dir_[d_]; \
break; \
} \
dir_[d_] = -1; \
} \
} while (not all (dir_ == -1)); \
}
namespace CarpetMask {
using namespace std;
using namespace Carpet;
/**
* Reduce the weight on the current and the next coarser level to
* make things consistent. Set the weight to 0 inside the active
* region of the next coarser level, maybe to 1/2 near the boundary
* of that region, and also to 1/2 near the prolongation boundary of
* this level.
*/
extern "C" {
void
CarpetMaskSetup (CCTK_ARGUMENTS);
}
void
CarpetMaskSetup (CCTK_ARGUMENTS)
{
DECLARE_CCTK_PARAMETERS;
if (not is_singlemap_mode()) {
CCTK_WARN (CCTK_WARN_ABORT,
"This routine may only be called in singlemap mode");
}
gh const& hh = *vhh.AT(Carpet::map);
dh const& dd = *vdd.AT(Carpet::map);
ibbox const& base = hh.baseextent(mglevel, reflevel);
if (reflevel > 0) {
ivect const reffact =
spacereffacts.AT(reflevel) / spacereffacts.AT(reflevel-1);
assert (all (reffact == 2));
}
// Set prolongation boundaries and restricted region of this level
BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
DECLARE_CCTK_ARGUMENTS;
ibbox const& ext =
dd.light_boxes.AT(mglevel).AT(reflevel).AT(component).exterior;
ibbox const& owned =
dd.light_boxes.AT(mglevel).AT(reflevel).AT(component).owned;
ibbox const world = owned;
ibset const& active =
dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).active;
// There are no prolongation boundaries on the coarsest level;
// the outer boundary is treated elsewhere
ibset const not_active = reflevel==0 ? ibset() : world - active;
ibset const& fine_active =
dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).fine_active;
if (verbose) {
ostringstream buf;
buf << "Setting prolongation region " << not_active << " on level " << reflevel;
CCTK_INFO (buf.str().c_str());
buf << "Setting restriction region " << fine_active << " on level " << reflevel;
CCTK_INFO (buf.str().c_str());
}
// Set the weight in the interior of the not_active and the
// fine_active regions to zero, and set the weight on the
// boundary of the not_active and fine_active regions to 1/2.
//
// For the prolongation region, the "boundary" is the first
// layer outside of the region. For the restricted region, the
// "boundary" is the outermost layer of grid points if this
// layer is aligned with the next coarser (i.e. the current)
// grid; otherwise, the boundary is empty.
ibset test_boxes;
ibset test_cfboxes;
LOOP_OVER_NEIGHBOURS (shift) {
// In this loop, shift [1,1,1] denotes a convex corner of the
// region which should be masked out, i.e. a region where only
// a small bit (1/8) of the region should be masked out.
// Concave edges are treated implicitly (sequentially), i.e.
// larger chunks are cut out multiple times: a concave xy edge
// counts as both an x face and a y face.
ibset boxes = not_active;
ibset fboxes = fine_active;
for (int d=0; d<dim; ++d) {
ivect const dir = ivect::dir(d);
fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
}
for (int d=0; d<dim; ++d) {
// Calculate the boundary in direction d
ivect const dir = ivect::dir(d);
switch (shift[d]) {
case -1: {
// left boundary
boxes = boxes.shift (-dir) - boxes;
fboxes = fboxes.shift(-dir) - fboxes;
break;
}
case 0: {
// interior
// do nothing
break;
}
case +1: {
// right boundary
boxes = boxes.shift (+dir) - boxes;
fboxes = fboxes.shift(+dir) - fboxes;
break;
}
default:
assert (0);
}
}
boxes &= ext;
ibset const cfboxes = fboxes.contracted_for(base) & ext;
test_boxes |= boxes;
test_cfboxes |= cfboxes;
if (verbose) {
ostringstream buf;
buf << "Setting boundary " << shift << ": prolongation region " << boxes;
buf << "Setting boundary " << shift << ": restriction region " << fboxes;
CCTK_INFO (buf.str().c_str());
}
// Set up a bit mask that keeps the upper (when dir[d]=-1) or
// lower (when dir[d]=+1) half of the bits in each direction d
unsigned const bits = BMSK(dim);
unsigned bmask = 0;
for (int d=0; d<dim; ++d) {
for (unsigned b=0; b<bits; ++b) {
if ((shift[d]==-1 and not BGET(b,d)) or
(shift[d]==+1 and BGET(b,d)))
{
bmask = BSET(bmask, b);
}
}
}
// Handle prolongation region
LOOP_OVER_BSET (cctkGH, boxes, box, imin, imax) {
if (verbose) {
ostringstream buf;
buf << "Setting prolongation region "<< imin << ":" << imax-ivect(1) << " on level " << reflevel << " boundary " << shift << " to bmask " << bmask;
CCTK_INFO (buf.str().c_str());
}
#pragma omp parallel
CCTK_LOOP3(CarpetMaskSetup_prolongation,
i,j,k,
imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
{
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
iweight[ind] &= bmask;
} CCTK_ENDLOOP3(CarpetMaskSetup_prolongation);
} END_LOOP_OVER_BSET;
// Handle restricted region
LOOP_OVER_BSET (cctkGH, cfboxes, box, imin, imax) {
if (verbose) {
ostringstream buf;
buf << "Setting restricted region "<< imin << ":" << imax-ivect(1) << " on level " << reflevel << " boundary " << shift << " to bmask " << bmask;
CCTK_INFO (buf.str().c_str());
}
#pragma omp parallel
CCTK_LOOP3(CarpetMaskSetup_restriction,
i,j,k,
imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
{
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
iweight[ind] &= bmask;
} CCTK_ENDLOOP3(CarpetMaskSetup_restriction);
} END_LOOP_OVER_BSET;
} END_LOOP_OVER_NEIGHBOURS;
{
ibset const boxes = not_active.expand(ivect(1), ivect(1)) & ext;
ibset fboxes = fine_active;
for (int d=0; d<dim; ++d) {
ivect const dir = ivect::dir(d);
fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
}
fboxes = fboxes.expand(ivect(1), ivect(1));
ibset const cfboxes = fboxes.contracted_for(base) & ext;
if (not (test_boxes == boxes ) or
not (test_cfboxes == cfboxes))
{
cout << "boxes=" << boxes << "\n"
<< "test_boxes=" << test_boxes << "\n"
<< "b-t=" << (boxes-test_boxes) << "\n"
<< "t-b=" << (test_boxes-boxes) << "\n"
<< "cfboxes=" << cfboxes << "\n"
<< "test_cfboxes=" << test_cfboxes << "\n"
<< "b-t=" << (cfboxes-test_cfboxes) << "\n"
<< "t-b=" << (test_cfboxes-cfboxes) << "\n";
}
assert (test_boxes == boxes );
assert (test_cfboxes == cfboxes);
}
} END_LOCAL_COMPONENT_LOOP;
}
} // namespace CarpetMask
|
