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
|
#include <assert.h>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
void
CoordBase_SetupMask (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_INT nboundaryzones[6];
CCTK_INT is_internal[6];
CCTK_INT is_staggered[6];
CCTK_INT shiftout[6];
int imin[3], imax[3]; /* domain extent */
int bmin[3], bmax[3]; /* boundary extent */
int i, j, k;
int d, f;
int dd;
int ierr;
ierr = GetBoundarySpecification
(6, nboundaryzones, is_internal, is_staggered, shiftout);
if (ierr != 0) {
CCTK_WARN (0, "Could not get boundary specification");
}
/* Loop over all dimensions and faces */
for (d=0; d<3; ++d) {
for (f=0; f<2; ++f) {
/* If this processor has the outer boundary */
if (cctk_bbox[2*d+f]) {
int npoints;
if (is_internal[2*d+f]) {
/* The boundary extends inwards */
npoints = - nboundaryzones[2*d+f] + shiftout[2*d+f];
} else {
/* The boundary extends outwards */
npoints = nboundaryzones[2*d+f] + shiftout[2*d+f] - 1;
}
/* If there are boundary that should be ignored */
if (npoints >= 0) {
if (npoints < 0 || npoints > cctk_lsh[d]) {
CCTK_WARN (0, "Illegal number of boundary points");
}
/* Calculate the extent of the domain */
for (dd=0; dd<3; ++dd) {
imin[dd] = 0;
imax[dd] = cctk_lsh[dd];
}
/* Calculate the extent of the boundary */
for (dd=0; dd<3; ++dd) {
bmin[dd] = imin[dd];
bmax[dd] = imax[dd];
}
if (f==0) {
/* lower face */
bmax[d] = imin[d] + npoints;
} else {
/* upper face */
bmin[d] = imax[d] - npoints;
}
/* Loop over the boundary */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting boundary points in direction %d face %d to weight 0", d, f);
}
for (k=bmin[2]; k<bmax[2]; ++k) {
for (j=bmin[1]; j<bmax[1]; ++j) {
for (i=bmin[0]; i<bmax[0]; ++i) {
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
weight[ind] = 0.0;
}
}
}
/* When the boundary is not staggered, then give the points
on the boundary the weight 1/2 */
if (! is_staggered[2*d+f]) {
/* Check whether the domain is empty */
/* TODO: This check is flawed, and therefore disabled */
if (0 && imin[d] == imax[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
dimension less. */
/* 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 */
if (verbose) {
CCTK_VInfo (CCTK_THORNSTRING,
"Setting non-staggered boundary points in direction %d face %d to weight 1/2", d, f);
}
for (k=bmin[2]; k<bmax[2]; ++k) {
for (j=bmin[1]; j<bmax[1]; ++j) {
for (i=bmin[0]; i<bmax[0]; ++i) {
int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
weight[ind] *= 0.5;
}
}
}
} /* if the domain is not empty */
} /* if the boundary is not staggered */
} /* if there are boundary points */
} /* if is outer boundary */
} /* loop over faces */
} /* loop over directions */
}
|
