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
|
#include <assert.h>
#include <stdlib.h>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <cctk_Functions.h>
#include <loopcontrol.h>
#include "bits.h"
void
CoordBase_SetupMask (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_INT nboundaryzones[2*cctk_dim];
CCTK_INT is_internal[2*cctk_dim];
CCTK_INT is_staggered[2*cctk_dim];
CCTK_INT shiftout[2*cctk_dim];
CCTK_INT is_periodic[2*cctk_dim];
int bnd_points[2*cctk_dim]; /* points outside the domain */
int int_points[cctk_dim]; /* global interior points */
int gmin[cctk_dim], gmax[cctk_dim]; /* global domain extent */
int imin[cctk_dim], imax[cctk_dim]; /* domain extent */
int bmin[cctk_dim], bmax[cctk_dim]; /* boundary extent */
int ierr;
int const reflevel = GetRefinementLevel(cctkGH);
if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
int const m = MultiPatch_GetMap (cctkGH);
assert (m >= 0);
ierr = MultiPatch_GetBoundarySpecification
(m, 2*cctk_dim, nboundaryzones, is_internal, is_staggered, shiftout);
} else {
ierr = GetBoundarySpecification
(2*cctk_dim, nboundaryzones, is_internal, is_staggered, shiftout);
}
if (ierr != 0) {
CCTK_ERROR ("Could not get boundary specification");
}
for (int d=0; d<cctk_dim; ++d) {
for (int f=0; f<2; ++f) {
is_periodic[2*d+f] = 0;
}
}
CCTK_INT const *const periodic =
CCTK_ParameterGet("periodic", "PeriodicCarpet", NULL);
CCTK_INT const *const periodic_x =
CCTK_ParameterGet("periodic_x", "PeriodicCarpet", NULL);
CCTK_INT const *const periodic_y =
CCTK_ParameterGet("periodic_x", "PeriodicCarpet", NULL);
CCTK_INT const *const periodic_z =
CCTK_ParameterGet("periodic_x", "PeriodicCarpet", NULL);
is_periodic[0] = is_periodic[1] =
(periodic && *periodic) || (periodic_x && *periodic_x);
is_periodic[2] = is_periodic[3] =
(periodic && *periodic) || (periodic_y && *periodic_y);
is_periodic[4] = is_periodic[5] =
(periodic && *periodic) || (periodic_z && *periodic_z);
/* Calculate the number of boundary points. This excludes points
that are directly on the boundary. */
for (int d=0; d<cctk_dim; ++d) {
for (int f=0; f<2; ++f) {
if (is_internal[2*d+f]) {
/* The boundary extends inwards */
bnd_points[2*d+f] = shiftout[2*d+f];
} else {
/* The boundary extends outwards */
if (is_periodic[2*d+f]) {
bnd_points[2*d+f] = nboundaryzones[2*d+f];
} else{
bnd_points[2*d+f] =
nboundaryzones[2*d+f] + shiftout[2*d+f] + is_staggered[2*d+f] - 1;
}
}
}
}
/* Calculate the global extent of the domain */
for (int d=0; d<cctk_dim; ++d) {
gmin[d] = 0;
gmax[d] = cctk_gsh[d];
}
/* Ensure that the domain specification is consistent */
for (int d=0; d<cctk_dim; ++d) {
int_points[d] = (gmax[d] - bnd_points[2*d+1]) - (gmin[d] + bnd_points[2*d]);
if (int_points[d] < 0) {
CCTK_ERROR ("Number of internal grid points is negative");
}
}
/* Loop over all dimensions and faces */
for (int d=0; d<cctk_dim; ++d) {
for (int f=0; f<2; ++f) {
/* If this processor has the outer boundary */
if (cctk_bbox[2*d+f]) {
if (bnd_points[2*d+f] < 0 || bnd_points[2*d+f] > cctk_lsh[d]) {
CCTK_ERROR ("Illegal number of boundary points");
}
/* Calculate the extent of the local part of the domain */
for (int dd=0; dd<cctk_dim; ++dd) {
imin[dd] = 0;
imax[dd] = cctk_lsh[dd];
}
/* Calculate the extent of the boundary */
for (int dd=0; dd<cctk_dim; ++dd) {
bmin[dd] = imin[dd];
bmax[dd] = imax[dd];
}
if (f==0) {
/* lower face */
bmax[d] = imin[d] + bnd_points[2*d+f];
} else {
/* upper face */
bmin[d] = imax[d] - bnd_points[2*d+f];
}
/* Loop over the boundary */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting boundary points in direction %d face %d to weight 0 on level %d", d, f, reflevel);
}
#pragma omp parallel
CCTK_LOOP3(CoordBase_SetupMask_boundary,
i,j,k,
bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2],
cctk_ash[0],cctk_ash[1],cctk_ash[2])
{
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
iweight[ind] = 0;
} CCTK_ENDLOOP3(CoordBase_SetupMask_boundary);
/* When the boundary is neither staggered nor periodic, then
give the points directly on the boundary the weight 1/2 */
if (! is_staggered[2*d+f] && ! is_periodic[2*d+f]) {
/* Since we are going to cut off 1/2 at each boundary, the
domain size must be at least 1 to begin with */
assert (int_points[d] >= 1);
/* Check whether the domain is empty */
if (int_points[d] == 1) {
/* The domain is empty. The correct thing to do would be
to set the weights to 0. But this is boring, because
then there are no interior points left, and all
reductions become trivial. Instead, leave the
non-staggered boundary points alone, and assume that
the user wants to perform a reduction in one fewer
dimensions. */
/* do nothing */
} else {
/* Adapt the extent of the boundary */
if (f==0) {
/* lower face */
bmin[d] = bmax[d];
bmax[d] = bmin[d] + 1;
} else {
/* upper face */
bmax[d] = bmin[d];
bmin[d] = bmax[d] - 1;
}
/* Loop over the points next to boundary */
unsigned const bits = BMSK(cctk_dim);
unsigned bmask = 0;
for (unsigned b=0; b<bits; ++b) {
if (BGET(b,d) == !f) {
bmask = BSET(bmask, b);
}
}
assert (BCNT(bmask) == bits/2);
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting non-staggered boundary points on level %d in direction %d face %d to bmask 0x%x", reflevel, d, f, bmask);
}
#pragma omp parallel
CCTK_LOOP3(CoordBase_SetupMask_boundary2,
i,j,k,
bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2],
cctk_ash[0],cctk_ash[1],cctk_ash[2])
{
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
iweight[ind] &= bmask;
} CCTK_ENDLOOP3(CoordBase_SetupMask_boundary2);
} /* if the domain is not empty */
} /* if the boundary is not staggered */
} else { /* if this is a ghost boundary */
/* Calculate the extent of the local part of the domain */
for (int dd=0; dd<cctk_dim; ++dd) {
imin[dd] = 0;
imax[dd] = cctk_lsh[dd];
}
/* Calculate the extent of the boundary */
for (int dd=0; dd<cctk_dim; ++dd) {
bmin[dd] = imin[dd];
bmax[dd] = imax[dd];
}
if (f==0) {
/* lower face */
bmax[d] = imin[d] + cctk_nghostzones[d];
} else {
/* upper face */
bmin[d] = imax[d] - cctk_nghostzones[d];
}
/* Loop over the boundary */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting ghost points in direction %d face %d to weight 0 on level %d", d, f, reflevel);
}
#pragma omp parallel
CCTK_LOOP3(CoordBase_SetupMask_ghost,
i,j,k,
bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2],
cctk_ash[0],cctk_ash[1],cctk_ash[2])
{
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
iweight[ind] = 0;
} CCTK_ENDLOOP3(CoordBase_SetupMask_ghost);
} /* if this is a ghost boundary */
} /* loop over faces */
} /* loop over directions */
}
|
