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 nanauto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling({max_tile_size_x, max_tile_size_y, max_tile_size_z});#pragma omp parallelfor (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 validbool 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 validbool 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 nanvoid 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 parallelfor (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 levelsconst 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 nanconst int tl = 0;auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling({max_tile_size_x, max_tile_size_y, max_tile_size_z});#pragma omp parallelfor (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 variablesif (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);elseassert(0);
// Only restrict valid grid functionsbool 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);elseassert(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 parallelfor (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 boundariesgrid.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 0Loop::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 0Loop::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 0Loop::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