BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC
JD6PQOJ6YYNQYEEWEXO2NM7NVYNBUI6V7ZU6Q3FNHGAT2VYOF5WAC
R6B5YHARRV34JO6W4GAWBOQ4BW5H6VEKXWMVKJKPWC7MUO4ZYZVQC
NE2O3IMQBR3OAPLJC22SPJ7AX2ABULA7XNSIZYFHCQLWC45T6GJAC
UD2O4RHXGNMJII7BYM75WDELMJFIAEJMMZM72VJMY34FGAB2DT4AC
MFSNKENJILMYLGGTLYSL4SH3K7KHM6WA5STADZGD3TKDTFWCWYAQC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC
MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC
2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC
33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
KG47IF4CPCUT3BHS34WDRHTH5HYMBTY4OSTB3X7APR2E5ZJ32IYQC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
QN2UTSQP4IMCMVXZNR24J22N24ASGXF6EWXQ2P3VWVQERUPFAFDQC
2DD222JSYRPHTXKSRXLSOMSCQPZUORNFLLO2P3GMIDELAAMD5MEQC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
TOBGHRPKEPSXDN56WGSZNWOMCBVJ4KUSLWYWI56MC2RR3MM3KLZAC
UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC
ORAAQXS3UNKXYDJBCTHNEIDJIZDHWMM5EVCZKVMFRKQK2NIQNJGAC
VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC
5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
UUGQGVC4WEKN64WAP7F5QPS2UHGQB5ZLMFRIYNWKMIEBDO3QRX4AC
IVHURSHY4636OGIF3PNDO5CWOVRLJ75M4LP65J6I2E6KAM4QKF4AC
PQB3EKQ6MBCXPTW4HB7SGMSTOTYMB3EFZX2573OFCQI6PSOEKCSQC
QNV4LD7UGYSSNYDXGJC6SRP6Y3USUVDQLEERBMBWRC7LILWREB7AC
2XYZZL42IEZHGDJA6NDKGSQKGJP24LOTLFJ6RNHOKWHHSUYIHGKQC
UTHNLH3J3VG7BBJPOKGE6DY7Z2QHDLY72TR2VMEDQYP73SIWZGKAC
NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC
for (int vi = 0; vi < groupdata.numvars; ++vi) {
valid.at(tl).valid_int =
coarsegroupdata.valid.at(tl).at(vi).valid_int &&
coarsegroupdata.valid.at(tl).at(vi).valid_bnd &&
groupdata.valid.at(tl).at(vi).valid_int &&
groupdata.valid.at(tl).at(vi).valid_bnd;
valid.at(tl).valid_bnd = false;
}
}
// Create coarse grid
CCTK_REAL time = 0.0; // dummy time
ghext->amrcore->MakeNewGrids(time);
void CreateRefinedGrid(int level) {
DECLARE_CCTK_PARAMETERS;
if (level > ghext->amrcore->maxLevel())
return;
if (verbose)
CCTK_VINFO("CreateRefinedGrid level %d", level);
// Create refined grid
CCTK_REAL time = 0.0; // dummy time
ghext->amrcore->regrid(0, time);
int numlevels = ghext->amrcore->finestLevel() + 1;
int maxnumlevels = ghext->amrcore->maxLevel() + 1;
assert(numlevels >= 0 && numlevels <= maxnumlevels);
assert(numlevels <= level + 1);
}
struct valid_t {
bool valid_int, valid_bnd;
valid_t() : valid_int(false), valid_bnd(false) {}
friend bool operator<(const valid_t &x, const valid_t &y) {
if (x.valid_int < y.valid_int)
return true;
return x.valid_bnd < y.valid_bnd;
}
friend ostream &operator<<(ostream &os, const valid_t v) {
auto str = [](bool v) { return v ? "VAL" : "INV"; };
return os << "[int:" << str(v.valid_int) << ",bnd:" << str(v.valid_bnd)
<< "]";
}
operator string() const {
ostringstream buf;
buf << *this;
return buf.str();
}
};
struct clause_t {
int gi, vi, tl;
valid_t valid;
friend bool operator<(const clause_t &x, const clause_t &y) {
if (x.gi < y.gi)
return true;
if (x.vi < y.vi)
return true;
if (x.tl < y.tl)
return true;
return x.valid < y.valid;
}
};
vector<clause_t> decode_clauses(const cFunctionData *restrict attribute,
int n_clauses, const char **clauses) {
vector<clause_t> result;
for (int n = 0; n < n_clauses; ++n) {
for (const char *restrict p = clauses[n]; *p;) {
// Find beginning of word
while (isspace(*p))
++p;
if (!*p)
break;
// Read grid function / group name
assert(isalnum(*p) || *p == '_' || *p == ':');
const char *const name0 = p;
while (isalnum(*p) || *p == '_' || *p == ':')
++p;
string name(name0, p);
valid_t valid;
if (*p == '(') {
++p;
const char *const reg0 = p;
while (isalpha(*p))
++p;
string regstr(reg0, p);
if (regstr == "interior")
valid.valid_int = true;
else if (regstr == "boundary")
valid.valid_bnd = true;
else if (regstr == "everywhere")
valid.valid_int = valid.valid_bnd = true;
else
assert(0);
assert(*p == ')');
++p;
} else {
assert(0); // missing region
valid.valid_int = valid.valid_bnd = true;
}
assert(!*p || isspace(*p));
int tl = 0;
while (name.size() >= 2 && name.substr(name.size() - 2) == "_p") {
name = name.substr(0, name.size() - 2);
++tl;
}
if (name.find(':') == string::npos)
name = string(attribute->thorn) + "::" + name;
const int gi0 = CCTK_GroupIndex(name.c_str());
if (gi0 >= 0) {
const int gi = gi0;
const int nv = CCTK_NumVarsInGroupI(gi);
for (int vi = 0; vi < nv; ++vi)
result.push_back({gi, vi, tl, valid});
} else {
const int var = CCTK_VarIndex(name.c_str());
if (var >= 0) {
const int gi = CCTK_GroupIndexFromVarI(var);
const int v0 = CCTK_FirstVarIndexI(gi);
const int vi = var - v0;
result.push_back({gi, vi, tl, valid});
} else {
CCTK_VERROR("Cannot decode group/variable name \"%s\" in %s: %s::%s",
name.c_str(), attribute->where, attribute->thorn,
attribute->routine);
}
}
}
}
set<clause_t> seen;
for (const auto &res : result) {
assert(seen.count(res) == 0); // variable listed twice
seen.insert(res);
}
return result;
}
// Output domain information
if (CCTK_MyProc(nullptr) == 0) {
enter_level_mode(cctkGH, ghext->leveldata.at(0));
const int *restrict gsh = cctkGH->cctk_gsh;
const int *restrict nghostzones = cctkGH->cctk_nghostzones;
CCTK_REAL x0[dim], x1[dim], dx[dim];
for (int d = 0; d < dim; ++d) {
dx[d] = cctkGH->cctk_delta_space[d];
x0[d] = cctkGH->cctk_origin_space[d];
x1[d] = x0[d] + (gsh[d] - 2 * nghostzones[d]) * dx[d];
}
CCTK_VINFO("Grid extent:");
CCTK_VINFO(" gsh=[%d,%d,%d]", gsh[0], gsh[1], gsh[2]);
CCTK_VINFO("Domain extent:");
CCTK_VINFO(" xmin=[%g,%g,%g]", x0[0], x0[1], x0[2]);
CCTK_VINFO(" xmax=[%g,%g,%g]", x1[0], x1[1], x1[2]);
CCTK_VINFO(" base dx=[%g,%g,%g]", dx[0], dx[1], dx[2]);
CCTK_VINFO("Time stepping:");
CCTK_VINFO(" t0=%g", cctkGH->cctk_time);
CCTK_VINFO(" dt=%g", cctkGH->cctk_delta_time);
leave_level_mode(cctkGH, ghext->leveldata.at(0));
}
// Create coarse grid
{
static Timer timer("InitialiseRegrid");
Interval interval(timer);
CCTK_REAL time = 0.0; // dummy time
assert(saved_cctkGH == nullptr);
saved_cctkGH = cctkGH;
ghext->amrcore->MakeNewGrids(time);
saved_cctkGH = nullptr;
}
// Output domain information
if (CCTK_MyProc(nullptr) == 0) {
enter_level_mode(cctkGH, ghext->leveldata.at(0));
const int *restrict gsh = cctkGH->cctk_gsh;
const int *restrict nghostzones = cctkGH->cctk_nghostzones;
CCTK_REAL x0[dim], x1[dim], dx[dim];
for (int d = 0; d < dim; ++d) {
dx[d] = cctkGH->cctk_delta_space[d];
x0[d] = cctkGH->cctk_origin_space[d];
x1[d] = x0[d] + (gsh[d] - 2 * nghostzones[d]) * dx[d];
}
CCTK_VINFO("Grid extent:");
CCTK_VINFO(" gsh=[%d,%d,%d]", gsh[0], gsh[1], gsh[2]);
CCTK_VINFO("Domain extent:");
CCTK_VINFO(" xmin=[%g,%g,%g]", x0[0], x0[1], x0[2]);
CCTK_VINFO(" xmax=[%g,%g,%g]", x1[0], x1[1], x1[2]);
CCTK_VINFO(" base dx=[%g,%g,%g]", dx[0], dx[1], dx[2]);
CCTK_VINFO("Time stepping:");
CCTK_VINFO(" t0=%g", cctkGH->cctk_time);
CCTK_VINFO(" dt=%g", cctkGH->cctk_delta_time);
leave_level_mode(cctkGH, ghext->leveldata.at(0));
}
const int min_level = current_level == -1 ? 0 : current_level;
const int max_level =
current_level == -1 ? ghext->leveldata.size() : current_level + 1;
// Check whether input variables have valid data
{
const vector<clause_t> reads = decode_clauses(
attribute, attribute->n_ReadsClauses, attribute->ReadsClauses);
for (const auto &rd : reads) {
for (int level = min_level; level < max_level; ++level) {
const auto &restrict leveldata = ghext->leveldata.at(level);
const auto &restrict groupdata = leveldata.groupdata.at(rd.gi);
const valid_t &need = rd.valid;
const valid_t &have = groupdata.valid.at(rd.tl).at(rd.vi);
// "x <= y" for booleans means "x implies y"
const bool cond = need.valid_int <= have.valid_int &&
need.valid_bnd <= have.valid_bnd;
if (!cond)
CCTK_VERROR(
"Found invalid input data: iteration %d %s: %s::%s, level %d, "
"variable %s%s: need %s, have %s",
cctkGH->cctk_iteration, attribute->where, attribute->thorn,
attribute->routine, leveldata.level,
CCTK_FullVarName(groupdata.firstvarindex + rd.vi),
string("_p", rd.tl).c_str(), string(need).c_str(),
string(have).c_str());
}
}
}
#warning "TODO: poison those output variables that are not input variables"
}
// Mark output variables as having valid data
{
const vector<clause_t> writes = decode_clauses(
attribute, attribute->n_WritesClauses, attribute->WritesClauses);
for (const auto &wr : writes) {
for (int level = min_level; level < max_level; ++level) {
auto &restrict leveldata = ghext->leveldata.at(level);
auto &restrict groupdata = leveldata.groupdata.at(wr.gi);
const valid_t &provided = wr.valid;
valid_t &have = groupdata.valid.at(wr.tl).at(wr.vi);
// Code cannot invalidate...
have.valid_int |= provided.valid_int;
have.valid_bnd |= provided.valid_bnd;
}
}
if (groupdata.indextype == array<int, dim>{0, 0, 0})
average_down_nodal(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 0, 0} ||
groupdata.indextype == array<int, dim>{0, 1, 0} ||
groupdata.indextype == array<int, dim>{0, 0, 1})
average_down_edges(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 1, 0} ||
groupdata.indextype == array<int, dim>{1, 0, 1} ||
groupdata.indextype == array<int, dim>{0, 1, 1})
average_down_faces(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 1, 1})
average_down(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl), 0,
groupdata.numvars, reffact);
else
assert(0);
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi).valid_int &=
fine_groupdata.valid.at(tl).at(vi).valid_int &&
fine_groupdata.valid.at(tl).at(vi).valid_bnd;
groupdata.valid.at(tl).at(vi).valid_bnd = false;
}
#if 0
#warning \
"TODO: Remove this loop; it only checks AMReX's internal correctness, and this is not a good way for doing so"
// Poison points that will be restricted
grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {
const int i = grid.cactus_offset.x + p.i;
const int j = grid.cactus_offset.y + p.j;
const int k = grid.cactus_offset.z + p.k;
bool is_restricted = true;
// For cell centred indices (indextype=1), check the
// mask directly. For vertex centred indices
// (indextype=0), check the two neighbouring cells. We
// assume restriction is possible if either of the cells
// is unmasked. At the boundary of the mask grid
// function, be conservative and don't poison.
for (int c = groupdata.indextype[2] - 1; c < 1; ++c)
for (int b = groupdata.indextype[1] - 1; b < 1; ++b)
for (int a = groupdata.indextype[0] - 1; a < 1; ++a)
if (i + a >= mask.begin.x && i + a < mask.end.x &&
j + b >= mask.begin.y && j + b < mask.end.y &&
k + c >= mask.begin.z && k + c < mask.end.z)
is_restricted &= mask(i + a, j + b, k + c);
else
is_restricted = false; // be conservative
if (is_restricted)
ptr[p.idx] = 0.0 / 0.0;
});
#endif
if (groupdata.indextype == array<int, dim>{0, 0, 0})
average_down_nodal(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 0, 0} ||
groupdata.indextype == array<int, dim>{0, 1, 0} ||
groupdata.indextype == array<int, dim>{0, 0, 1})
average_down_edges(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 1, 0} ||
groupdata.indextype == array<int, dim>{1, 0, 1} ||
groupdata.indextype == array<int, dim>{0, 1, 1})
average_down_faces(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl),
reffact);
else if (groupdata.indextype == array<int, dim>{1, 1, 1})
average_down(*fine_groupdata.mfab.at(tl), *groupdata.mfab.at(tl), 0,
groupdata.numvars, reffact);
else
assert(0);
WRITES: potential_phi(everywhere)
WRITES: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
WRITES: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
WRITES: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
WRITES: current_rho(everywhere)
WRITES: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
WRITES: potential_phi(everywhere)
WRITES: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
WRITES: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
WRITES: current_rho(everywhere)
WRITES: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
SYNC: potential_phi
SYNC: potential_ax potential_ay potential_az
SYNC: field_ex field_ey field_ez
SYNC: field_bx field_by field_bz
SYNC: current_rho
SYNC: current_jx current_jy current_jz
} "Apply boundary conditions for the Maxwell equations"
READS: potential_phi(everywhere)
READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
READS: potential_phi_p(everywhere)
READS: potential_ax_p(everywhere) potential_ay_p(everywhere) potential_az_p(everywhere)
READS: field_ex_p(everywhere) field_ey_p(everywhere) field_ez_p(everywhere)
READS: field_bx_p(everywhere) field_by_p(everywhere) field_bz_p(everywhere)
READS: current_rho_p(everywhere)
READS: current_jx_p(everywhere) current_jy_p(everywhere) current_jz_p(everywhere)
WRITES: potential_ax(interior) potential_ay(interior) potential_az(interior)
WRITES: field_bx(interior) field_by(interior) field_bz(interior)
WRITES: current_jx(interior) current_jy(interior) current_jz(interior)
READS: potential_phi_p(everywhere)
READS: potential_ax_p(everywhere) potential_ay_p(everywhere) potential_az_p(everywhere)
READS: field_ex_p(everywhere) field_ey_p(everywhere) field_ez_p(everywhere)
READS: current_rho_p(everywhere)
READS: current_jx_p(everywhere) current_jy_p(everywhere) current_jz_p(everywhere)
WRITES: potential_ax(interior) potential_ay(interior) potential_az(interior)
WRITES: current_jx(interior) current_jy(interior) current_jz(interior)
READS: potential_phi_p(everywhere)
READS: field_ex_p(everywhere) field_ey_p(everywhere) field_ez_p(everywhere)
READS: current_rho_p(everywhere)
READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
READS: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
WRITES: potential_phi(interior)
WRITES: field_ex(interior) field_ey(interior) field_ez(interior)
WRITES: current_rho(interior)
READS: potential_phi(everywhere)
READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
WRITES: AMReX::regrid_error(interior)
READS: potential_phi(everywhere)
READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
READS: current_rho(everywhere)
READS: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)
WRITES: errors_potential_phi(everywhere)
WRITES: errors_potential_ax(everywhere) errors_potential_ay(everywhere) errors_potential_az(everywhere)
WRITES: errors_field_ex(everywhere) errors_field_ey(everywhere) errors_field_ez(everywhere)
WRITES: errors_field_bx(everywhere) errors_field_by(everywhere) errors_field_bz(everywhere)
WRITES: errors_current_rho(everywhere)
WRITES: errors_current_jx(everywhere) errors_current_jy(everywhere) errors_current_jz(everywhere)
READS: potential_phi(everywhere)
READS: potential_ax(everywhere) potential_ay(everywhere) potential_az(everywhere)
READS: field_ex(everywhere) field_ey(everywhere) field_ez(everywhere)
READS: field_bx(everywhere) field_by(everywhere) field_bz(everywhere)
READS: current_rho(everywhere)
READS: current_jx(everywhere) current_jy(everywhere) current_jz(everywhere)