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
|
#include <assert.h>
#include <sstream>
#include <vector>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "gh.hh"
#include "carpet.hh"
#include "regrid.hh"
namespace CarpetRegrid {
using namespace std;
using namespace Carpet;
int ManualGridpoints (cGH const * const cctkGH,
gh const & hh,
gh::rexts & bbsss,
gh::rbnds & obss,
gh::rprocs & pss)
{
DECLARE_CCTK_PARAMETERS;
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 (reflevel == 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);
assert (! smart_outer_boundaries);
for (size_t rl=1; rl<bbsss.size(); ++rl) {
bbvect const ob (false);
vector<ibbox> bbs;
gh::cbnds obs;
ManualGridpoints_OneLevel
(cctkGH, hh, rl,refinement_levels,
ilower.at(rl-1), iupper.at(rl-1), ob, bbs, obs);
// make multiprocessor aware
gh::cprocs ps;
SplitRegions (cctkGH, bbs, obs, ps);
// make multigrid aware
vector<vector<ibbox> > bbss;
MakeMultigridBoxes (cctkGH, 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 & hh,
const int rl,
const int numrl,
const ivect ilower,
const ivect iupper,
const bbvect obound,
vector<ibbox> & bbs,
vector<bbvect> & obs)
{
const ivect rstr = hh.baseextent.stride();
const ivect rlb = hh.baseextent.lower();
const ivect rub = hh.baseextent.upper();
const int levfac = ipow(hh.reffact, rl);
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 #" << rl << " 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 #" << rl << " 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", rl);
}
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
|
