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
|
#include <assert.h>
#include <vector>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "gh.hh"
#include "carpet.hh"
#include "regrid.hh"
extern "C" {
static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpoints.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetRegrid_manualgridpoints_cc);
}
namespace CarpetRegrid {
using namespace std;
using namespace Carpet;
int ManualGridpoints (cGH const * const cctkGH,
gh<dim> const & hh,
int const reflevel,
int const map,
int const size,
jjvect const & nboundaryzones,
jjvect const & is_internal,
jjvect const & is_staggered,
jjvect const & shiftout,
gh<dim>::rexts & bbsss,
gh<dim>::rbnds & obss,
gh<dim>::rprocs & pss)
{
DECLARE_CCTK_PARAMETERS;
assert (reflevel>=0 && reflevel<maxreflevels);
assert (map>=0 && map<maps);
if (refinement_levels > 4) {
CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
}
assert (refinement_levels >= 1 && refinement_levels <= 4);
// do nothing if the levels already exist
if (bbsss.size() == refinement_levels) return 0;
assert (bbsss.size() >= 1);
bbsss.resize (refinement_levels);
obss.resize (refinement_levels);
pss.resize (refinement_levels);
vector<ivect> ilower(3), iupper(3);
ilower.at(0) = ivect (l1ixmin, l1iymin, l1izmin);
iupper.at(0) = ivect (l1ixmax, l1iymax, l1izmax);
ilower.at(1) = ivect (l2ixmin, l2iymin, l2izmin);
iupper.at(1) = ivect (l2ixmax, l2iymax, l2izmax);
ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin);
iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax);
for (size_t rl=1; rl<bbsss.size(); ++rl) {
bbvect const ob (false);
vector<ibbox> bbs;
gh<dim>::cbnds obs;
ManualGridpoints_OneLevel
(cctkGH, hh, rl, refinement_levels,
ilower.at(rl-1), iupper.at(rl-1), ob, bbs, obs);
// make multiprocessor aware
gh<dim>::cprocs ps;
SplitRegions (cctkGH, bbs, obs, ps);
// make multigrid aware
vector<vector<ibbox> > bbss;
MakeMultigridBoxes
(cctkGH,
size, nboundaryzones, is_internal, is_staggered, shiftout,
bbs, obs, bbss);
bbsss.at(rl) = bbss;
obss.at(rl) = obs;
pss.at(rl) = ps;
} // for rl
return 1;
}
void ManualGridpoints_OneLevel (const cGH * const cctkGH,
const gh<dim> & hh,
const int reflevel,
const int reflevels,
const ivect ilower,
const ivect iupper,
const bbvect obound,
vector<ibbox> & bbs,
vector<bbvect> & obs)
{
if (reflevel >= reflevels) return;
const ivect rstr = hh.baseextent.stride();
const ivect rlb = hh.baseextent.lower();
const ivect rub = hh.baseextent.upper();
const int levfac = ipow(hh.reffact, reflevel);
assert (all (rstr % levfac == 0));
const ivect str (rstr / levfac);
const ivect lb (ilower);
const ivect ub (iupper);
if (! all(lb>=rlb && ub<=rub)) {
ostringstream buf;
buf << "The refinement region boundaries for refinement level #" << reflevel << " are not within the main grid. Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
CCTK_WARN (0, buf.str().c_str());
}
if (! all(lb<=ub)) {
ostringstream buf;
buf << "The refinement region boundaries for refinement level #" << reflevel << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends;
CCTK_WARN (0, buf.str().c_str());
}
if (! all(lb%str==0 && ub%str==0)) {
CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
"The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", reflevel);
}
assert (all(lb>=rlb && ub<=rub));
assert (all(lb<=ub));
assert (all(lb%str==0 && ub%str==0));
bbs.push_back (ibbox(lb, ub, str));
obs.push_back (obound);
}
} // namespace CarpetRegrid
|
