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 gridCCTK_REAL time = 0.0; // dummy timeghext->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 gridCCTK_REAL time = 0.0; // dummy timeghext->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 wordwhile (isspace(*p))++p;if (!*p)break;// Read grid function / group nameassert(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;elseassert(0);assert(*p == ')');++p;} else {assert(0); // missing regionvalid.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 twiceseen.insert(res);}return result;}
// Output domain informationif (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 timeassert(saved_cctkGH == nullptr);saved_cctkGH = cctkGH;ghext->amrcore->MakeNewGrids(time);saved_cctkGH = nullptr;}// Output domain informationif (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);elseassert(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 restrictedgrid.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);elseis_restricted = false; // be conservativeif (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);elseassert(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_phiSYNC: potential_ax potential_ay potential_azSYNC: field_ex field_ey field_ezSYNC: field_bx field_by field_bzSYNC: current_rhoSYNC: 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)