KCIWCVZOHG44WBOLKI2XK33WPHPRI5FWCETF4AOGTPZISKCW3CLQC
TFYA5EBLWIMWQ4ULT6MK4Y2XWAFJH5JPLS2MGM7IR4ZWLNJHTNXQC
IIDW4VXRBD7WM6WBZ7JNZ4GMGCIWF2NMYM2W75SYV5XFLRTV63RAC
TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC
RTNZAS3UPI6GG3KY4Z5WVXJ4R2YF5427BB6WAV3GHRS5W7XPOSUQC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC
33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC
BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC
R6B5YHARRV34JO6W4GAWBOQ4BW5H6VEKXWMVKJKPWC7MUO4ZYZVQC
QNV4LD7UGYSSNYDXGJC6SRP6Y3USUVDQLEERBMBWRC7LILWREB7AC
2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
UUGQGVC4WEKN64WAP7F5QPS2UHGQB5ZLMFRIYNWKMIEBDO3QRX4AC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC
RF3SDQ3ELEBCGJSFQL2VOSHNN34RAH6KM3LEITEYOHKJMMHHKMFAC
2XYZZL42IEZHGDJA6NDKGSQKGJP24LOTLFJ6RNHOKWHHSUYIHGKQC
UTHNLH3J3VG7BBJPOKGE6DY7Z2QHDLY72TR2VMEDQYP73SIWZGKAC
NE2O3IMQBR3OAPLJC22SPJ7AX2ABULA7XNSIZYFHCQLWC45T6GJAC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC
5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC
UD2O4RHXGNMJII7BYM75WDELMJFIAEJMMZM72VJMY34FGAB2DT4AC
ORAAQXS3UNKXYDJBCTHNEIDJIZDHWMM5EVCZKVMFRKQK2NIQNJGAC
AEVGZIZEUIC52MCK3J4V547YEV2R4YQL3JUJW7FSP4R34PSZ43DAC
if (poison_undefined_values) {
// Set grid functions to nan
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
for (int vi = 0; vi < groupdata.numvars; ++vi) {
CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.0;
});
}
}
}
for (int vi = 0; vi < groupdata.numvars; ++vi)
poison_invalid(leveldata, groupdata, vi, tl);
#warning "TODO: only interpolate if coarse grid data are valid"
InterpFromCoarseLevel(*groupdata.mfab.at(tl), 0.0,
*coarsegroupdata.mfab.at(tl), 0, 0,
// Only interpolate if coarse grid data are valid
bool all_invalid = true;
for (int vi = 0; vi < groupdata.numvars; ++vi)
all_invalid &= !coarsegroupdata.valid.at(tl).at(vi).valid_int &&
!coarsegroupdata.valid.at(tl).at(vi).valid_bnd;
if (all_invalid) {
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi).valid_int = false;
groupdata.valid.at(tl).at(vi).valid_bnd = false;
}
} else {
for (int vi = 0; vi < groupdata.numvars; ++vi)
assert(coarsegroupdata.valid.at(tl).at(vi).valid_int &&
coarsegroupdata.valid.at(tl).at(vi).valid_bnd);
InterpFromCoarseLevel(
*groupdata.mfab.at(tl), 0.0, *coarsegroupdata.mfab.at(tl), 0, 0,
ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0,
reffact, &cell_bilinear_interp, bcs, 0);
ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0, reffact,
&cell_bilinear_interp, bcs, 0);
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi).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_bnd = false;
}
}
// Only interpolate if coarse grid data are valid
bool all_invalid = true;
for (int vi = 0; vi < groupdata.numvars; ++vi)
all_invalid &= !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;
if (all_invalid) {
valid.at(tl).valid_int = false;
valid.at(tl).valid_bnd = false;
} else {
for (int vi = 0; vi < groupdata.numvars; ++vi) {
const bool cond = 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;
if (!cond)
CCTK_VERROR("Found invalid input data: RemakeLevel level %d, "
"variable %s%s: need everything defined, have coarse "
"%s, have current %s",
leveldata.level,
CCTK_FullVarName(groupdata.firstvarindex + vi),
string("_p", tl).c_str(),
string(coarsegroupdata.valid.at(tl).at(vi)).c_str(),
string(groupdata.valid.at(tl).at(vi)).c_str());
}
// Set grid functions to nan
void poison_invalid(const GHExt::LevelData &leveldata,
const GHExt::LevelData::GroupData &groupdata, int vi,
int tl) {
DECLARE_CCTK_PARAMETERS;
if (!poison_undefined_values)
return;
const valid_t &valid = groupdata.valid.at(tl).at(vi);
if (valid.valid_int && valid.valid_bnd)
return;
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
if (!valid.valid_int) {
if (!valid.valid_bnd) {
grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.0;
});
} else {
grid.loop_int(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.0;
});
}
} else {
if (!valid.valid_bnd) {
grid.loop_bnd(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.0;
});
} else {
assert(0);
}
}
}
}
}
bool get_group_checkpoint_flag(const int gi) {
int tags = CCTK_GroupTagsTableI(gi);
assert(tags >= 0);
char buf[100];
int iret = Util_TableGetString(tags, sizeof buf, buf, "checkpoint");
if (iret == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
return true;
} else if (iret >= 0) {
string str(buf);
for (auto &c : str)
c = tolower(c);
if (str == "yes")
return true;
if (str == "no")
return false;
assert(0);
} else {
assert(0);
}
}
void InvalidateTimelevels(cGH *restrict const cctkGH) {
DECLARE_CCTK_PARAMETERS;
for (auto &restrict leveldata : ghext->leveldata) {
for (auto &restrict groupdata : leveldata.groupdata) {
const bool checkpoint = get_group_checkpoint_flag(groupdata.groupindex);
if (!checkpoint) {
// Invalidate all time levels
const int ntls = groupdata.mfab.size();
for (int tl = 0; tl < ntls; ++tl) {
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi) = valid_t();
poison_invalid(leveldata, groupdata, vi, tl);
}
}
}
}
}
groupdata.valid.at(0) = vector<valid_t>(groupdata.numvars);
if (poison_undefined_values) {
// Set grid functions to nan
const int tl = 0;
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
for (int vi = 0; vi < groupdata.numvars; ++vi) {
CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.0;
});
}
}
}
groupdata.valid.at(0) = move(vtmp);
for (int vi = 0; vi < groupdata.numvars; ++vi)
groupdata.valid.at(0).at(vi) = valid_t();
for (int vi = 0; vi < groupdata.numvars; ++vi)
poison_invalid(leveldata, groupdata, vi, 0);
#warning "TODO: poison those output variables that are not input variables"
// Poison those output variables that are not input variables
if (poison_undefined_values) {
map<clause_t, valid_t> isread;
const vector<clause_t> &reads = decode_clauses(
attribute, attribute->n_ReadsClauses, attribute->ReadsClauses);
for (const auto &rd : reads) {
clause_t cl = rd;
cl.valid = valid_t();
assert(isread.count(cl) == 0);
isread[cl] = rd.valid;
}
const vector<clause_t> &writes = decode_clauses(
attribute, attribute->n_WritesClauses, attribute->WritesClauses);
for (const auto &wr : writes) {
clause_t cl = wr;
cl.valid = valid_t();
valid_t need;
if (isread.count(cl) > 0)
need = isread[cl];
const valid_t &provided = wr.valid;
for (int level = min_level; level < max_level; ++level) {
auto &restrict leveldata = ghext->leveldata.at(level);
auto &restrict groupdata = leveldata.groupdata.at(wr.gi);
valid_t &have = groupdata.valid.at(wr.tl).at(wr.vi);
have.valid_int &= need.valid_int || !provided.valid_int;
have.valid_bnd &= need.valid_bnd || !provided.valid_bnd;
}
}
}
for (int vi = 0; vi < groupdata.numvars; ++vi)
assert(groupdata.valid.at(tl).at(vi).valid_int);
for (int vi = 0; vi < groupdata.numvars; ++vi)
groupdata.valid.at(tl).at(vi).valid_bnd = false;
for (int vi = 0; vi < groupdata.numvars; ++vi)
poison_invalid(leveldata, groupdata, vi, tl);
for (int vi = 0; vi < groupdata.numvars; ++vi)
assert(coarsegroupdata.valid.at(tl).at(vi).valid_int &&
coarsegroupdata.valid.at(tl).at(vi).valid_bnd &&
groupdata.valid.at(tl).at(vi).valid_int);
for (int vi = 0; vi < groupdata.numvars; ++vi)
groupdata.valid.at(tl).at(vi).valid_bnd = false;
for (int vi = 0; vi < groupdata.numvars; ++vi)
poison_invalid(leveldata, groupdata, vi, tl);
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);
// Only restrict valid grid functions
bool all_invalid = true;
for (int vi = 0; vi < groupdata.numvars; ++vi)
all_invalid &= !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_int;
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 (all_invalid) {
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi).valid_int = false;
groupdata.valid.at(tl).at(vi).valid_bnd = false;
}
} else {
for (int vi = 0; vi < groupdata.numvars; ++vi)
assert(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_int);
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);
if (poison_undefined_values) {
MultiFab &mfab = *groupdata.mfab.at(tl);
const MultiFab &fine_mfab = *fine_groupdata.mfab.at(tl);
const IntVect reffact{2, 2, 2};
const iMultiFab finemask =
makeFineMask(mfab, fine_mfab.boxArray(), reffact);
auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
{max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
for (MFIter mfi(*leveldata.mfab0, mfitinfo); mfi.isValid(); ++mfi) {
GridPtrDesc grid(leveldata, mfi);
const Array4<CCTK_REAL> &vars = mfab.array(mfi);
for (int vi = 0; vi < groupdata.numvars; ++vi) {
CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
// Poison boundaries
grid.loop_bnd(groupdata.indextype, [&](const Loop::PointDesc &p) {
ptr[p.idx] = 0.0 / 0.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;
// Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
// regrid_error[p.idx] = fabs(p.x) < 0.5 && fabs(p.y) < 0.5 && fabs(p.z) <
// 0.5;
// });
#if 0
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
regrid_error[p.idx] = fabs(p.x) < 0.5 && fabs(p.y) < 0.5 && fabs(p.z) <
0.5;
});
#endif
#if 0
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
CCTK_REAL err = 0;
for (int b = 0; b < 2; ++b)
for (int a = 0; a < 2; ++a)
err += fabs(ax_(p.i - 1, p.j + a, p.k + b) -
2 * ax_(p.i, p.j + a, p.k + b) +
ax_(p.i + 1, p.j + a, p.k + b)) +
fabs(ay_(p.i + a, p.j - 1, p.k + b) -
2 * ay_(p.i + a, p.j, p.k + b) +
ay_(p.i + a, p.j + 1, p.k + b)) +
fabs(az_(p.i + a, p.j + b, p.k - 1) -
2 * az_(p.i + a, p.j + b, p.k) +
az_(p.i + a, p.j + b, p.k + 1));
regrid_error[p.idx] = err;
});
#endif
// Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
// CCTK_REAL err = 0;
// for (int b = 0; b < 2; ++b)
// for (int a = 0; a < 2; ++a)
// err += fabs(ax_(p.i - 1, p.j + a, p.k + b) -
// 2 * ax_(p.i, p.j + a, p.k + b) +
// ax_(p.i + 1, p.j + a, p.k + b)) +
// fabs(ay_(p.i + a, p.j - 1, p.k + b) -
// 2 * ay_(p.i + a, p.j, p.k + b) +
// ay_(p.i + a, p.j + 1, p.k + b)) +
// fabs(az_(p.i + a, p.j + b, p.k - 1) -
// 2 * az_(p.i + a, p.j + b, p.k) +
// az_(p.i + a, p.j + b, p.k + 1));
// regrid_error[p.idx] = err;
// });
#if 0
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
auto closeto = [&](CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL r) {
return fabs(p.x - x) <= r && fabs(p.y - y) <= r && fabs(p.z - z) <= r;
};
constexpr CCTK_REAL x = 0.25;
constexpr CCTK_REAL y = 0.25;
constexpr CCTK_REAL z = 0.25;
constexpr CCTK_REAL r = 0.2;
regrid_error[p.idx] = // closeto(-x, -y, -z, r) ||
closeto(+x, -y, -z, r) || closeto(-x, +y, -z, r) ||
closeto(+x, +y, -z, r) || closeto(-x, -y, +z, r) ||
closeto(+x, -y, +z, r) || closeto(-x, +y, +z, r) ||
closeto(+x, +y, +z, r);
});
#endif