YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC 722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC 2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC 33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC 5IAXY3XZJTRMMVT2OVIJ6OXQJI6OJPTPCHHA4IVLVMHANCCC5NKAC WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC QN2UTSQP4IMCMVXZNR24J22N24ASGXF6EWXQ2P3VWVQERUPFAFDQC RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC 3KQ5ZQE45TNGATMW4R4NHXVBULGOXHF7JYC53BZL2LNDQ2V7C2CAC OJZWEAJRLOA5FZAC6FMGBQ6XFQFDNO55IRWLFWDAYWM54MUIU5LAC UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC KG47IF4CPCUT3BHS34WDRHTH5HYMBTY4OSTB3X7APR2E5ZJ32IYQC 2DD222JSYRPHTXKSRXLSOMSCQPZUORNFLLO2P3GMIDELAAMD5MEQC M6DALF7V7J377IXGE5ABJRJ25PMBDOBGQX5ZXTQZHTGOHQSUSZKAC AEVGZIZEUIC52MCK3J4V547YEV2R4YQL3JUJW7FSP4R34PSZ43DAC # AMReX::out_tsv = yes
grid.loop_int([&](int i, int j, int k, int idx) {tag(grid.cactus_offset.x + i, grid.cactus_offset.y + j,grid.cactus_offset.z + k) =err[idx] >= regrid_error_threshold ? TagBox::SET : TagBox::CLEAR;
grid.loop_int(groupdata.indextype, [&](const Loop::PointDesc &p) {tag(grid.cactus_offset.x + p.i, grid.cactus_offset.y + p.j,grid.cactus_offset.z + p.k) =err[p.idx] >= regrid_error_threshold ? TagBox::SET : TagBox::CLEAR;
}}namespace {array<int, dim> get_group_indextype(const int gi) {assert(gi >= 0);const int tags = CCTK_GroupTagsTableI(gi);assert(tags >= 0);array<CCTK_INT, dim> index;int iret = Util_TableGetIntArray(tags, dim, index.data(), "index");if (iret == UTIL_ERROR_TABLE_NO_SUCH_KEY) {index = {1, 1, 1}; // default: cell-centred} else {assert(iret == dim);
if (leveldata.level == 0) {// Coarsest level: Copy from same level// This should not happenassert(0);PhysBCFunctNoOp fphysbc;for (int tl = 0; tl < ntls; ++tl) {auto mfab =make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);// 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(*mfab, 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);grid.loop_all([&](int i, int j, int k, int idx) { ptr[idx] = 0.0 / 0.0; });}}if (tl < prolongate_tl)FillPatchSingleLevel(*mfab, 0.0, {&*groupdata.mfab.at(tl)}, {0.0}, 0,0, groupdata.numvars,ghext->amrcore->Geom(level), fphysbc, 0);groupdata.mfab.at(tl) = move(mfab);}
auto &coarseleveldata = ghext->leveldata.at(level - 1);auto &restrict coarsegroupdata = coarseleveldata.groupdata.at(gi);assert(coarsegroupdata.numvars == groupdata.numvars);
// Copy from same level and/or prolongate from next coarser levelauto &coarseleveldata = ghext->leveldata.at(level - 1);auto &restrict coarsegroupdata = coarseleveldata.groupdata.at(gi);assert(coarsegroupdata.numvars == groupdata.numvars);
PhysBCFunctNoOp cphysbc;PhysBCFunctNoOp fphysbc;const IntVect reffact{2, 2, 2};// boundary conditionsconst BCRec bcrec(periodic_x ? BCType::int_dir : BCType::reflect_odd,periodic_y ? BCType::int_dir : BCType::reflect_odd,periodic_z ? BCType::int_dir : BCType::reflect_odd,periodic_x ? BCType::int_dir : BCType::reflect_odd,periodic_y ? BCType::int_dir : BCType::reflect_odd,periodic_z ? BCType::int_dir : BCType::reflect_odd);const Vector<BCRec> bcs(groupdata.numvars, bcrec);for (int tl = 0; tl < ntls; ++tl) {auto mfab =make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);// Set grid functions to nanauto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling({max_tile_size_x, max_tile_size_y, max_tile_size_z});
PhysBCFunctNoOp cphysbc;PhysBCFunctNoOp fphysbc;const IntVect reffact{2, 2, 2};// boundary conditionsconst BCRec bcrec(periodic_x ? BCType::int_dir : BCType::reflect_odd,periodic_y ? BCType::int_dir : BCType::reflect_odd,periodic_z ? BCType::int_dir : BCType::reflect_odd,periodic_x ? BCType::int_dir : BCType::reflect_odd,periodic_y ? BCType::int_dir : BCType::reflect_odd,periodic_z ? BCType::int_dir : BCType::reflect_odd);const Vector<BCRec> bcs(groupdata.numvars, bcrec);for (int tl = 0; tl < ntls; ++tl) {auto mfab = make_unique<MultiFab>(ba, dm, groupdata.numvars, ghost_size);// Set grid functions to nanauto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling({max_tile_size_x, max_tile_size_y, max_tile_size_z});
for (MFIter mfi(*mfab, 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);grid.loop_all([&](int i, int j, int k, int idx) { ptr[idx] = 0.0 / 0.0; });}
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);grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {ptr[p.idx] = 0.0 / 0.0;});
if (tl < prolongate_tl)FillPatchTwoLevels(*mfab, 0.0, {&*coarsegroupdata.mfab.at(tl)}, {0.0},{&*groupdata.mfab.at(tl)}, {0.0}, 0, 0,groupdata.numvars, ghext->amrcore->Geom(level - 1),ghext->amrcore->Geom(level), cphysbc, 0, fphysbc,0, reffact, &cell_bilinear_interp, bcs, 0);groupdata.mfab.at(tl) = move(mfab);
if (tl < prolongate_tl)FillPatchTwoLevels(*mfab, 0.0, {&*coarsegroupdata.mfab.at(tl)}, {0.0},{&*groupdata.mfab.at(tl)}, {0.0}, 0, 0,groupdata.numvars, ghext->amrcore->Geom(level - 1),ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0,reffact, &cell_bilinear_interp, bcs, 0);groupdata.mfab.at(tl) = move(mfab);
void OutputPlotfile(const cGH *restrict cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;const int numgroups = CCTK_NumGroups();for (int gi = 0; gi < numgroups; ++gi) {auto &restrict groupdata0 = ghext->leveldata.at(0).groupdata.at(gi);if (groupdata0.mfab.size() > 0) {const int tl = 0;string groupname = unique_ptr<char>(CCTK_GroupName(gi)).get();groupname = regex_replace(groupname, regex("::"), "-");for (auto &c : groupname)c = tolower(c);ostringstream buf;buf << out_dir << "/" << groupname;buf << ".it" << setw(6) << setfill('0') << cctk_iteration;string filename = buf.str();Vector<string> varnames(groupdata0.numvars);for (int vi = 0; vi < groupdata0.numvars; ++vi) {ostringstream buf;buf << CCTK_VarName(groupdata0.firstvarindex + vi);for (int i = 0; i < tl; ++i)buf << "_p";varnames.at(vi) = buf.str();}Vector<const MultiFab *> mfabs(ghext->leveldata.size());Vector<Geometry> geoms(ghext->leveldata.size());Vector<int> iters(ghext->leveldata.size());Vector<IntVect> reffacts(ghext->leveldata.size());for (const auto &restrict leveldata : ghext->leveldata) {mfabs.at(leveldata.level) = &*leveldata.groupdata.at(gi).mfab.at(tl);geoms.at(leveldata.level) = ghext->amrcore->Geom(leveldata.level);iters.at(leveldata.level) = cctk_iteration;reffacts.at(leveldata.level) = IntVect{2, 2, 2};}
#warning "TODO: Get these from schedule.cxx"struct TileBox {array<int, dim> tile_min;array<int, dim> tile_max;};extern vector<TileBox> thread_local_tilebox;void enter_level_mode(cGH *restrict cctkGH,const GHExt::LevelData &restrict leveldata);void leave_level_mode(cGH *restrict cctkGH,const GHExt::LevelData &restrict leveldata);void enter_local_mode(cGH *restrict cctkGH, TileBox &restrict tilebox,const GHExt::LevelData &restrict leveldata,const MFIter &mfi);void leave_local_mode(cGH *restrict cctkGH, TileBox &restrict tilebox,const GHExt::LevelData &restrict leveldata,const MFIter &mfi);
const MultiFab &mfab = *groupdata.mfab.at(tl);for (MFIter mfi(mfab); mfi.isValid(); ++mfi) {
for (MFIter mfi(*leveldata.mfab0); mfi.isValid(); ++mfi) {TileBox &tilebox = thread_local_tilebox.at(omp_get_thread_num());enter_local_mode(const_cast<cGH *>(cctkGH), tilebox, leveldata, mfi);
write_arrays(file, cctkGH, leveldata.level, mfi.index(), ptrs, grid);
// write_arrays(file, cctkGH, leveldata.level, mfi.index(), ptrs, grid);Loop::loop_all(cctkGH, groupdata.indextype, [&](const Loop::PointDesc &p) {file << cctkGH->cctk_iteration << "\t" << cctkGH->cctk_time << "\t"<< leveldata.level << "\t"<< mfi.index()// << "\t" << isghost<< "\t" << (grid.lbnd[0] + p.i) << "\t" << (grid.lbnd[1] + p.j)<< "\t" << (grid.lbnd[2] + p.k) << "\t" << p.x << "\t" << p.y<< "\t" << p.z;for (const auto &ptr : ptrs)file << "\t" << ptr[p.idx];file << "\n";});leave_local_mode(const_cast<cGH *>(cctkGH), tilebox, leveldata, mfi);
Vector<const MultiFab *> mfabs(ghext->leveldata.size());Vector<Geometry> geoms(ghext->leveldata.size());Vector<int> iters(ghext->leveldata.size());Vector<IntVect> reffacts(ghext->leveldata.size());for (const auto &restrict leveldata : ghext->leveldata) {mfabs.at(leveldata.level) = &*leveldata.groupdata.at(gi).mfab.at(tl);geoms.at(leveldata.level) = ghext->amrcore->Geom(leveldata.level);iters.at(leveldata.level) = cctk_iteration;reffacts.at(leveldata.level) = IntVect{2, 2, 2};}
WriteASCII(cctkGH, filename, gi, varnames);}}}
// TODO: Output all groups into a single fileWriteMultiLevelPlotfile(filename, mfabs.size(), mfabs, varnames, geoms,cctk_time, iters, reffacts);if (out_tsv) {ostringstream procbuf;procbuf << out_dir << "/" << groupname;procbuf << ".it" << setw(6) << setfill('0') << cctk_iteration;procbuf << ".p" << setw(4) << setfill('0') << CCTK_MyProc(nullptr);string procfilename = procbuf.str();
void OutputNorms(const cGH *restrict cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;
WriteASCII(cctkGH, procfilename, gi, varnames);}count_vars += ghext->leveldata.at(0).groupdata.at(gi).numvars;}}
const bool is_root = CCTK_MyProc(nullptr) == 0;
return count_vars;
int OutputGH(const cGH *restrict cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;const bool is_root = CCTK_MyProc(nullptr) == 0;if (is_root)cout << "OutputGH: iteration " << cctk_iteration << ", time " << cctk_time<< ", run time " << CCTK_RunTime() << " s\n";if (out_every <= 0 || cctk_iteration % out_every != 0)return 0;OutputPlotfile(cctkGH);OutputASCII(cctkGH);OutputNorms(cctkGH);// TODO: This should be the number of variables outputreturn 0;
template <typename T, int CI, int CJ, int CK> struct GF3D {static_assert(CI == 0 || CI == 1, "");static_assert(CJ == 0 || CJ == 1, "");static_assert(CK == 0 || CK == 1, "");typedef T value_type;T *restrict ptr;static constexpr int di = 1;int dj, dk;int ni, nj, nk;constexpr array<int, dim> indextype() const { return {CI, CJ, CK}; }inline GF3D(const cGH *restrict cctkGH, T *restrict ptr): ptr(ptr), dj(di * (cctkGH->cctk_ash[0] + 1 - CI)),dk(dj * (cctkGH->cctk_ash[1] + 1 - CJ)),ni(cctkGH->cctk_lsh[0] + 1 - CI), nj(cctkGH->cctk_lsh[1] + 1 - CJ),nk(cctkGH->cctk_lsh[2] + 1 - CK) {}inline int offset(int i, int j, int k) const {assert(i >= 0 && i < ni);assert(j >= 0 && j < nj);assert(k >= 0 && k < nk);return i * di + j * dj + k * dk;}inline T &restrict operator()(int i, int j, int k) const {return ptr[offset(i, j, k)];}};////////////////////////////////////////////////////////////////////////////////struct PointDesc {int idx;int i, j, k;CCTK_REAL x, y, z;static constexpr int di = 1;int dj, dk;CCTK_REAL dx, dy, dz;};
}}template <typename F>void loop_all(const array<int, dim> &indextype, const F &f) const {// typedef void (GridDescBase::*funptr)(const F &f) const;// constexpr array<funptr, 8> funptrs{// &GridDescBase::loop_all<0, 0, 0, F>,// &GridDescBase::loop_all<1, 0, 0, F>,// &GridDescBase::loop_all<0, 1, 0, F>,// &GridDescBase::loop_all<1, 1, 0, F>,// &GridDescBase::loop_all<0, 0, 1, F>,// &GridDescBase::loop_all<1, 0, 1, F>,// &GridDescBase::loop_all<0, 1, 1, F>,// &GridDescBase::loop_all<1, 1, 1, F>,// };// return (this->*funptrs[indextype[0] + 2 * indextype[1] + 4 *// indextype[2]])(// f);switch (indextype[0] + 2 * indextype[1] + 4 * indextype[2]) {case 0b000:return loop_all<0, 0, 0>(f);case 0b001:return loop_all<1, 0, 0>(f);case 0b010:return loop_all<0, 1, 0>(f);case 0b011:return loop_all<1, 1, 0>(f);case 0b100:return loop_all<0, 0, 1>(f);case 0b101:return loop_all<1, 0, 1>(f);case 0b110:return loop_all<0, 1, 1>(f);case 0b111:return loop_all<1, 1, 1>(f);default:assert(0);}}template <typename F>void loop_int(const array<int, dim> &indextype, const F &f) const {switch (indextype[0] + 2 * indextype[1] + 4 * indextype[2]) {case 0b000:return loop_int<0, 0, 0>(f);case 0b001:return loop_int<1, 0, 0>(f);case 0b010:return loop_int<0, 1, 0>(f);case 0b011:return loop_int<1, 1, 0>(f);case 0b100:return loop_int<0, 0, 1>(f);case 0b101:return loop_int<1, 0, 1>(f);case 0b110:return loop_int<0, 1, 1>(f);case 0b111:return loop_int<1, 1, 1>(f);default:assert(0);}}template <typename F>void loop_bnd(const array<int, dim> &indextype, const F &f) const {switch (indextype[0] + 2 * indextype[1] + 4 * indextype[2]) {case 0b000:return loop_bnd<0, 0, 0>(f);case 0b001:return loop_bnd<1, 0, 0>(f);case 0b010:return loop_bnd<0, 1, 0>(f);case 0b011:return loop_bnd<1, 1, 0>(f);case 0b100:return loop_bnd<0, 0, 1>(f);case 0b101:return loop_bnd<1, 0, 1>(f);case 0b110:return loop_bnd<0, 1, 1>(f);case 0b111:return loop_bnd<1, 1, 1>(f);default:assert(0);
template <typename F> void loop_all(const cGH *cctkGH, const F &f) {GridDescBase(cctkGH).loop_all(f);
template <int CI, int CJ, int CK, typename F>void loop_all(const cGH *cctkGH, const F &f) {GridDescBase(cctkGH).loop_all<CI, CJ, CK>(f);}template <int CI, int CJ, int CK, typename F>void loop_int(const cGH *cctkGH, const F &f) {GridDescBase(cctkGH).loop_int<CI, CJ, CK>(f);
template <typename F>void loop_int(const cGH *cctkGH, const array<int, dim> &indextype, const F &f) {GridDescBase(cctkGH).loop_int(indextype, f);}template <typename F>void loop_bnd(const cGH *cctkGH, const array<int, dim> &indextype, const F &f) {GridDescBase(cctkGH).loop_bnd(indextype, f);}
reduction<T> reduce_array(const Geometry &geom, const T *restrict ptr,const GridDesc &grid) {
reduction<T> reduce_array(const Geometry &geom,const GHExt::LevelData::GroupData &groupdata,const T *restrict ptr, const GridDesc &grid) {
const cGH *restrict threadGH = &AMReX::thread_local_cctkGH.at(thread_num);// Check whether this is the correct cGH structureassert(cctkGH == threadGH);
if (omp_in_parallel()) {const cGH *restrict threadGH = &AMReX::thread_local_cctkGH.at(thread_num);// Check whether this is the correct cGH structureassert(cctkGH == threadGH);}
const Geometry &geom = ghext->amrcore->Geom(0);const CCTK_REAL *restrict const global_x0 = geom.ProbLo();const CCTK_REAL *restrict const global_dx = geom.CellSize();for (int d = 0; d < dim; ++d) {const int levfac = 1 << leveldata.level;const int levoff = 1 - 2 * nghostzones[d];const int levoffdenom = 2;const CCTK_REAL origin_space = global_x0[d];const CCTK_REAL delta_space = global_dx[d];dx[d] = delta_space / levfac;x0[d] = origin_space + dx[d] * levoff / levoffdenom;}
const cGH *restrict threadGH = &thread_local_cctkGH.at(thread_num);// Check whether this is the correct cGH structureassert(cctkGH == threadGH);
if (omp_in_parallel()) {const cGH *restrict threadGH = &thread_local_cctkGH.at(thread_num);// Check whether this is the correct cGH structureassert(cctkGH == threadGH);}
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {rho[idx] = 1.0;momx[idx] = 0.0;momy[idx] = 0.0;momz[idx] = 0.0;etot[idx] = 1.0;
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {rho[p.idx] = 1.0;momx[p.idx] = 0.0;momy[p.idx] = 0.0;momz[p.idx] = 0.0;etot[p.idx] = 1.0;
rho[idx] = rho_p[idx] - dt * ((momx_p[idx + di] - momx_p[idx]) / dx +(momy_p[idx + dj] - momy_p[idx]) / dy +(momz_p[idx + dk] - momz_p[idx]) / dz);
rho[p.idx] =rho_p[p.idx] - dt * ((momx_p[p.idx + di] - momx_p[p.idx]) / p.dx +(momy_p[p.idx + dj] - momy_p[p.idx]) / p.dy +(momz_p[p.idx + dk] - momz_p[p.idx]) / p.dz);
CCTK_REAL potential TYPE=gf TIMELEVELS=2{# d*dA = mu J A^a,^b_b - A^b,^a_b = mu J^a# d*A = 0 A^a,a = 0
# d*dA = mu J A^a,^b_b - A^b,^a_b = mu J^a# d*A = 0 A^a,a = 0# -d*(dphi + A,t) = rho/epsilon# d*dA + dt*(dphi+A,t) = mu J# *phi,t + d*A = 0 phi,t + A^i,i = 0# F = dA F_ab = A_a,b - A_b,a# dF = 0 F_ab,c + F_bc,a + F_ca,b = 0# d*F = mu J F^ab,a = mu J^b
# -d*(dphi + A,t) = rho/epsilon# d*dA + dt*(dphi+A,t) = mu J# *phi,t + d*A = 0 phi,t + A^i,i = 0
# E = -dphi - A,t E_i = - phi,i - A_i,t# B = dA B_ij = A_j,i - A_i,j# dE + B,t = 0 E_j,i - E_i,j + B_ij,t = 0# dB = 0 B_jk,i + B_ki,j + B_ijk = 0# d*E = rho/epsilon E^i,i = rho/epsilon# d*B - dt*E = mu J - B^ij,i - E^j,t = mu J^j# CCTK_REAL potential TYPE=gf TIMELEVELS=2a# {# phi# ax ay az# } "Electromagnetic potential"## CCTK_REAL field TYPE=gf TIMELEVELS=2# {# ex ey ez# bx by bz # B_yz B_zx B_xy# } "Electromagnetic field"## CCTK_REAL current TYPE=gf TIMELEVELS=2# {# rho# jx jy jz# } "Electric charge current"## CCTK_REAL errors TYPE=gf# {# phierr# axerr ayerr azerr# exerr eyerr ezerr# bxerr byerr bzerr# rhoerr# jxerr jyerr jzerr# } "Errors"
CCTK_REAL field TYPE=gf TIMELEVELS=2
CCTK_REAL potential_ax TYPE=gf TIMELEVELS=2 TAGS='index={1 0 0}'{ax} "Electromagnetic potential"CCTK_REAL potential_ay TYPE=gf TIMELEVELS=2 TAGS='index={0 1 0}'{ay} "Electromagnetic potential"CCTK_REAL potential_az TYPE=gf TIMELEVELS=2 TAGS='index={0 0 1}'
# E = -dphi - A,t E_i = - phi,i - A_i,t# B = dA B_ij = A_j,i - A_i,j# dE + B,t = 0 E_j,i - E_i,j + B_ij,t = 0# dB = 0 B_jk,i + B_ki,j + B_ijk = 0# d*E = rho/epsilon E^i,i = rho/epsilon# d*B - dt*E = mu J - B^ij,i - E^j,t = mu J^j
CCTK_REAL field_ex TYPE=gf TIMELEVELS=2 TAGS='index={1 0 0}'{ex} "Electromagnetic field"CCTK_REAL field_ey TYPE=gf TIMELEVELS=2 TAGS='index={0 1 0}'{ey} "Electromagnetic field"CCTK_REAL field_ez TYPE=gf TIMELEVELS=2 TAGS='index={0 0 1}'{ez} "Electromagnetic field"
ex ey ezbx by bz # B_yz B_zx B_xy
CCTK_REAL field_bx TYPE=gf TIMELEVELS=2 TAGS='index={0 1 1}'{bx} "Electromagnetic field"CCTK_REAL field_by TYPE=gf TIMELEVELS=2 TAGS='index={1 0 1}'{by} "Electromagnetic field"CCTK_REAL field_bz TYPE=gf TIMELEVELS=2 TAGS='index={1 1 0}'{bz
CCTK_REAL errors TYPE=gf
CCTK_REAL current_jx TYPE=gf TIMELEVELS=2 TAGS='index={1 0 0}'{jx} "Electric charge current"CCTK_REAL current_jy TYPE=gf TIMELEVELS=2 TAGS='index={0 1 0}'{jy} "Electric charge current"CCTK_REAL current_jz TYPE=gf TIMELEVELS=2 TAGS='index={0 0 1}'{jz} "Electric charge current"CCTK_REAL errors_potential_phi TYPE=gf TAGS='index={0 0 0}'
axerr ayerr azerrexerr eyerr ezerrbxerr byerr bzerr
} "Electromagnetic potential"CCTK_REAL errors_potential_ax TYPE=gf TAGS='index={1 0 0}'{axerr} "Electromagnetic potential"CCTK_REAL errors_potential_ay TYPE=gf TAGS='index={0 1 0}'{ayerr} "Electromagnetic potential"CCTK_REAL errors_potential_az TYPE=gf TAGS='index={0 0 1}'{azerr} "Electromagnetic potential"CCTK_REAL errors_field_ex TYPE=gf TAGS='index={1 0 0}'{exerr} "Electromagnetic field"CCTK_REAL errors_field_ey TYPE=gf TAGS='index={0 1 0}'{eyerr} "Electromagnetic field"CCTK_REAL errors_field_ez TYPE=gf TAGS='index={0 0 1}'{ezerr} "Electromagnetic field"CCTK_REAL errors_field_bx TYPE=gf TAGS='index={0 1 1}'{bxerr} "Electromagnetic field"CCTK_REAL errors_field_by TYPE=gf TAGS='index={1 0 1}'{byerr} "Electromagnetic field"CCTK_REAL errors_field_bz TYPE=gf TAGS='index={1 1 0}'{bzerr} "Electromagnetic field"CCTK_REAL errors_current_rho TYPE=gf TAGS='index={0 0 0}'{
jxerr jyerr jzerr} "Errors"
} "Electric charge current"CCTK_REAL errors_current_jx TYPE=gf TAGS='index={1 0 0}'{jxerr} "Electric charge current"CCTK_REAL errors_current_jy TYPE=gf TAGS='index={0 1 0}'{jyerr} "Electric charge current"CCTK_REAL errors_current_jz TYPE=gf TAGS='index={0 0 1}'{jzerr} "Electric charge current"
return {.phi = x.phi / a, .ax = x.ax / a, .ay = x.ay / a, .az = x.az / a};
return {.phi = f(x.phi, y.phi),.ax = f(x.ax, y.ax),.ay = f(x.ay, y.ay),.az = f(x.az, y.az)};}potential operator+() const { return *this; }potential operator-() const { return map1(std::negate<T>()); }potential operator+(const potential &y) const {return map2(std::plus<T>(), y);}potential operator-(const potential &y) const {return map2(std::minus<T>(), y);}potential operator*(const potential &y) const {return map2(std::multiplies<T>(), y);}potential operator/(const potential &y) const {return map2(std::divides<T>(), y);
const CCTK_REAL dx = CCTK_DELTA_SPACE(0);const CCTK_REAL dy = CCTK_DELTA_SPACE(1);const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
const Loop::GF3D<CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);const Loop::GF3D<CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);const Loop::GF3D<CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);const Loop::GF3D<CCTK_REAL, 0, 0, 1> az_(cctkGH, az);const Loop::GF3D<CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);const Loop::GF3D<CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);const Loop::GF3D<CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);const Loop::GF3D<CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);const Loop::GF3D<CCTK_REAL, 1, 0, 1> by_(cctkGH, by);const Loop::GF3D<CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);const Loop::GF3D<CCTK_REAL, 0, 0, 0> rho_(cctkGH, rho);const Loop::GF3D<CCTK_REAL, 1, 0, 0> jx_(cctkGH, jx);const Loop::GF3D<CCTK_REAL, 0, 1, 0> jy_(cctkGH, jy);const Loop::GF3D<CCTK_REAL, 0, 0, 1> jz_(cctkGH, jz);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0 + (cctk_lbnd[0] + i) * dx;CCTK_REAL y = y0 + (cctk_lbnd[1] + j) * dy;CCTK_REAL z = z0 + (cctk_lbnd[2] + k) * dz;
Loop::loop_all<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {phi_(p.i, p.j, p.k) = plane_wave(t + dt / 2, p.x, p.y, p.z).phi;});
phi[idx] = plane_wave(t + dt / 2, x, y, z).phi;ax[idx] = plane_wave(t, x + dx / 2, y, z).ax;ay[idx] = plane_wave(t, x, y + dy / 2, z).ay;az[idx] = plane_wave(t, x, y, z + dz / 2).az;
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {ax_(p.i, p.j, p.k) = plane_wave(t, p.x, p.y, p.z).ax;});Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {ay_(p.i, p.j, p.k) = plane_wave(t, p.x, p.y, p.z).ay;});Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {az_(p.i, p.j, p.k) = plane_wave(t, p.x, p.y, p.z).az;});
ex[idx] =-xderiv(plane_wave<CCTK_REAL>, dx)(t + dt / 2, x + dx / 2, y, z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x + dx / 2, y, z).ax;ey[idx] =-yderiv(plane_wave<CCTK_REAL>, dy)(t + dt / 2, x, y + dy / 2, z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x, y + dy / 2, z).ay;ez[idx] =-zderiv(plane_wave<CCTK_REAL>, dz)(t + dt / 2, x, y, z + dz / 2).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x, y, z + dz / 2).az;bx[idx] =yderiv(plane_wave<CCTK_REAL>, dy)(t, x, y + dy / 2, z + dz / 2).az -zderiv(plane_wave<CCTK_REAL>, dz)(t, x, y + dy / 2, z + dz / 2).ay;by[idx] =zderiv(plane_wave<CCTK_REAL>, dz)(t, x + dx / 2, y, z + dz / 2).ax -xderiv(plane_wave<CCTK_REAL>, dx)(t, x + dx / 2, y, z + dz / 2).az;bz[idx] =xderiv(plane_wave<CCTK_REAL>, dx)(t, x, y + dy / 2, z + dz / 2).ay -yderiv(plane_wave<CCTK_REAL>, dy)(t, x, y + dy / 2, z + dz / 2).ax;
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {ex_(p.i, p.j, p.k) =-xderiv(plane_wave<CCTK_REAL>, p.dx)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ax;});Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {ey_(p.i, p.j, p.k) =-yderiv(plane_wave<CCTK_REAL>, p.dy)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ay;});Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {ez_(p.i, p.j, p.k) =-zderiv(plane_wave<CCTK_REAL>, p.dz)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).az;});
rho[idx] = 0.0;jx[idx] = 0.0;jy[idx] = 0.0;jz[idx] = 0.0;
Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {bx_(p.i, p.j, p.k) =yderiv(plane_wave<CCTK_REAL>, p.dy)(t, p.x, p.y, p.z).az -zderiv(plane_wave<CCTK_REAL>, p.dz)(t, p.x, p.y, p.z).ay;});Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {by_(p.i, p.j, p.k) =zderiv(plane_wave<CCTK_REAL>, p.dz)(t, p.x, p.y, p.z).ax -xderiv(plane_wave<CCTK_REAL>, p.dx)(t, p.x, p.y, p.z).az;});Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {bz_(p.i, p.j, p.k) =xderiv(plane_wave<CCTK_REAL>, p.dx)(t, p.x, p.y, p.z).ay -yderiv(plane_wave<CCTK_REAL>, p.dy)(t, p.x, p.y, p.z).ax;
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) { jx_(p.i, p.j, p.k) = 0.0; });Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) { jy_(p.i, p.j, p.k) = 0.0; });Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) { jz_(p.i, p.j, p.k) = 0.0; });
constexpr int di = 1;const int dj = di * cctk_ash[0];const int dk = dj * cctk_ash[1];
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_p_(cctkGH, ax_p);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_p_(cctkGH, ay_p);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_p_(cctkGH, az_p);
Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {jx[idx] = jx_p[idx];jy[idx] = jy_p[idx];jz[idx] = jz_p[idx];
const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_p_(cctkGH, ex_p);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_p_(cctkGH, ey_p);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_p_(cctkGH, ez_p);
bx[idx] = bx_p[idx] + dt * ((ey_p[idx + dk] - ey_p[idx]) / dz -(ez_p[idx + dj] - ez_p[idx]) / dy);by[idx] = by_p[idx] + dt * ((ez_p[idx + di] - ez_p[idx]) / dx -(ex_p[idx + dk] - ex_p[idx]) / dz);bz[idx] = bz_p[idx] + dt * ((ex_p[idx + dj] - ex_p[idx]) / dy -(ey_p[idx + di] - ey_p[idx]) / dx);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> bx_p_(cctkGH, bx_p);const Loop::GF3D<const CCTK_REAL, 1, 0, 1> by_p_(cctkGH, by_p);const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bz_p_(cctkGH, bz_p);
ax[idx] =ax_p[idx] - dt * ((phi_p[idx + di] - phi_p[idx]) / dx + ex_p[idx]);ay[idx] =ay_p[idx] - dt * ((phi_p[idx + dj] - phi_p[idx]) / dy + ey_p[idx]);az[idx] =az_p[idx] - dt * ((phi_p[idx + dk] - phi_p[idx]) / dz + ez_p[idx]);
const Loop::GF3D<const CCTK_REAL, 0, 0, 0> rho_p_(cctkGH, rho_p);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> jx_p_(cctkGH, jx_p);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> jy_p_(cctkGH, jy_p);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> jz_p_(cctkGH, jz_p);const Loop::GF3D<CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);const Loop::GF3D<CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);const Loop::GF3D<CCTK_REAL, 0, 0, 1> az_(cctkGH, az);const Loop::GF3D<CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);const Loop::GF3D<CCTK_REAL, 1, 0, 1> by_(cctkGH, by);const Loop::GF3D<CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);const Loop::GF3D<CCTK_REAL, 1, 0, 0> jx_(cctkGH, jx);const Loop::GF3D<CCTK_REAL, 0, 1, 0> jy_(cctkGH, jy);const Loop::GF3D<CCTK_REAL, 0, 0, 1> jz_(cctkGH, jz);Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {jx_(p.i, p.j, p.k) = jx_p_(p.i, p.j, p.k);});Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {jy_(p.i, p.j, p.k) = jy_p_(p.i, p.j, p.k);});Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {jz_(p.i, p.j, p.k) = jz_p_(p.i, p.j, p.k);});Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {bx_(p.i, p.j, p.k) =bx_p_(p.i, p.j, p.k) +dt * ((ey_p_(p.i, p.j, p.k + 1) - ey_p_(p.i, p.j, p.k)) / p.dz -(ez_p_(p.i, p.j + 1, p.k) - ez_p_(p.i, p.j, p.k)) / p.dy);});Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {by_(p.i, p.j, p.k) =by_p_(p.i, p.j, p.k) +dt * ((ez_p_(p.i + 1, p.j, p.k) - ez_p_(p.i, p.j, p.k)) / p.dx -(ex_p_(p.i, p.j, p.k + 1) - ex_p_(p.i, p.j, p.k)) / p.dz);});Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {bz_(p.i, p.j, p.k) =bz_p_(p.i, p.j, p.k) +dt * ((ex_p_(p.i, p.j + 1, p.k) - ex_p_(p.i, p.j, p.k)) / p.dy -(ey_p_(p.i + 1, p.j, p.k) - ey_p_(p.i, p.j, p.k)) / p.dx);});Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {ax_(p.i, p.j, p.k) =ax_p_(p.i, p.j, p.k) -dt * ((phi_p_(p.i + 1, p.j, p.k) - phi_p_(p.i, p.j, p.k)) / p.dx +ex_p_(p.i, p.j, p.k));});Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {ay_(p.i, p.j, p.k) =ay_p_(p.i, p.j, p.k) -dt * ((phi_p_(p.i, p.j + 1, p.k) - phi_p_(p.i, p.j, p.k)) / p.dy +ey_p_(p.i, p.j, p.k));});Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {az_(p.i, p.j, p.k) =az_p_(p.i, p.j, p.k) -dt * ((phi_p_(p.i, p.j, p.k + 1) - phi_p_(p.i, p.j, p.k)) / p.dz +ez_p_(p.i, p.j, p.k));
const CCTK_REAL dx = CCTK_DELTA_SPACE(0);const CCTK_REAL dy = CCTK_DELTA_SPACE(1);const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
const Loop::GF3D<const CCTK_REAL, 0, 0, 0> phi_p_(cctkGH, phi_p);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_p_(cctkGH, ex_p);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_p_(cctkGH, ey_p);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_p_(cctkGH, ez_p);const Loop::GF3D<const CCTK_REAL, 0, 0, 0> rho_p_(cctkGH, rho_p);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_(cctkGH, az);
constexpr int di = 1;const int dj = di * cctk_ash[0];const int dk = dj * cctk_ash[1];
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);const Loop::GF3D<const CCTK_REAL, 1, 0, 1> by_(cctkGH, by);const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> jx_(cctkGH, jx);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> jy_(cctkGH, jy);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> jz_(cctkGH, jz);const Loop::GF3D<CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);const Loop::GF3D<CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);const Loop::GF3D<CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);const Loop::GF3D<CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
ex[idx] = ex_p[idx] - dt * ((by[idx] - by[idx - dk]) / dz -(bz[idx] - bz[idx - dj]) / dy + jx[idx]);ey[idx] = ey_p[idx] - dt * ((bz[idx] - bz[idx - di]) / dx -(bx[idx] - bx[idx - dk]) / dz + jy[idx]);ez[idx] = ez_p[idx] - dt * ((bx[idx] - bx[idx - dj]) / dy -(by[idx] - by[idx - di]) / dx + jz[idx]);
Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {rho_(p.i, p.j, p.k) =rho_p_(p.i, p.j, p.k) -dt * ((jx_(p.i, p.j, p.k) - jx_(p.i - 1, p.j, p.k)) / p.dx +(jy_(p.i, p.j, p.k) - jy_(p.i, p.j - 1, p.k)) / p.dy +(jz_(p.i, p.j, p.k) - jz_(p.i, p.j, p.k - 1)) / p.dz);});
phi[idx] = phi_p[idx] - dt * ((ax[idx] - ax[idx - di]) / dx +(ay[idx] - ay[idx - dj]) / dy +(ay[idx] - az[idx - dk]) / dz);
Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {ex_(p.i, p.j, p.k) =ex_p_(p.i, p.j, p.k) -dt * ((by_(p.i, p.j, p.k) - by_(p.i, p.j, p.k - 1)) / p.dz -(bz_(p.i, p.j, p.k) - bz_(p.i, p.j - 1, p.k)) / p.dy +jx_(p.i, p.j, p.k));});Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {ey_(p.i, p.j, p.k) =ey_p_(p.i, p.j, p.k) -dt * ((bz_(p.i, p.j, p.k) - bz_(p.i - 1, p.j, p.k)) / p.dx -(bx_(p.i, p.j, p.k) - bx_(p.i, p.j, p.k - 1)) / p.dz +jy_(p.i, p.j, p.k));});Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {ez_(p.i, p.j, p.k) =ez_p_(p.i, p.j, p.k) -dt * ((bx_(p.i, p.j, p.k) - bx_(p.i, p.j - 1, p.k)) / p.dy -(by_(p.i, p.j, p.k) - by_(p.i - 1, p.j, p.k)) / p.dx +jz_(p.i, p.j, p.k));
Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {phi_(p.i, p.j, p.k) =phi_p_(p.i, p.j, p.k) -dt * ((ax_(p.i, p.j, p.k) - ax_(p.i - 1, p.j, p.k)) / p.dx +(ay_(p.i, p.j, p.k) - ay_(p.i, p.j - 1, p.k)) / p.dy +(ay_(p.i, p.j, p.k) - az_(p.i, p.j, p.k - 1)) / p.dz);});
const CCTK_REAL dx = CCTK_DELTA_SPACE(0);const CCTK_REAL dy = CCTK_DELTA_SPACE(1);const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
const Loop::GF3D<CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_(cctkGH, az);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);const Loop::GF3D<const CCTK_REAL, 0, 1, 1> bx_(cctkGH, bx);const Loop::GF3D<const CCTK_REAL, 1, 0, 1> by_(cctkGH, by);const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bz_(cctkGH, bz);const Loop::GF3D<const CCTK_REAL, 0, 0, 0> rho_(cctkGH, rho);const Loop::GF3D<const CCTK_REAL, 1, 0, 0> jx_(cctkGH, jx);const Loop::GF3D<const CCTK_REAL, 0, 1, 0> jy_(cctkGH, jy);const Loop::GF3D<const CCTK_REAL, 0, 0, 1> jz_(cctkGH, jz);const Loop::GF3D<CCTK_REAL, 0, 0, 0> phierr_(cctkGH, phierr);const Loop::GF3D<CCTK_REAL, 1, 0, 0> axerr_(cctkGH, axerr);const Loop::GF3D<CCTK_REAL, 0, 1, 0> ayerr_(cctkGH, ayerr);const Loop::GF3D<CCTK_REAL, 0, 0, 1> azerr_(cctkGH, azerr);const Loop::GF3D<CCTK_REAL, 1, 0, 0> exerr_(cctkGH, exerr);const Loop::GF3D<CCTK_REAL, 0, 1, 0> eyerr_(cctkGH, eyerr);const Loop::GF3D<CCTK_REAL, 0, 0, 1> ezerr_(cctkGH, ezerr);const Loop::GF3D<CCTK_REAL, 0, 1, 1> bxerr_(cctkGH, bxerr);const Loop::GF3D<CCTK_REAL, 1, 0, 1> byerr_(cctkGH, byerr);const Loop::GF3D<CCTK_REAL, 1, 1, 0> bzerr_(cctkGH, bzerr);const Loop::GF3D<CCTK_REAL, 0, 0, 0> rhoerr_(cctkGH, rhoerr);const Loop::GF3D<CCTK_REAL, 1, 0, 0> jxerr_(cctkGH, jxerr);const Loop::GF3D<CCTK_REAL, 0, 1, 0> jyerr_(cctkGH, jyerr);const Loop::GF3D<CCTK_REAL, 0, 0, 1> jzerr_(cctkGH, jzerr);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0 + (cctk_lbnd[0] + i) * dx;CCTK_REAL y = y0 + (cctk_lbnd[1] + j) * dy;CCTK_REAL z = z0 + (cctk_lbnd[2] + k) * dz;
Loop::loop_all<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {phierr_(p.i, p.j, p.k) =phi_(p.i, p.j, p.k) - plane_wave(t + dt / 2, p.x, p.y, p.z).phi;});
phierr[idx] = phi[idx] - plane_wave(t + dt / 2, x, y, z).phi;axerr[idx] = ax[idx] - plane_wave(t, x + dx / 2, y, z).ax;ayerr[idx] = ay[idx] - plane_wave(t, x, y + dy / 2, z).ay;azerr[idx] = az[idx] - plane_wave(t, x, y, z + dz / 2).az;
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {axerr_(p.i, p.j, p.k) =ax_(p.i, p.j, p.k) - plane_wave(t, p.x, p.y, p.z).ax;});Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {ayerr_(p.i, p.j, p.k) =ay_(p.i, p.j, p.k) - plane_wave(t, p.x, p.y, p.z).ay;});Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {azerr_(p.i, p.j, p.k) =az_(p.i, p.j, p.k) - plane_wave(t, p.x, p.y, p.z).az;});
exerr[idx] =ex[idx] -(-xderiv(plane_wave<CCTK_REAL>, dx)(t + dt / 2, x + dx / 2, y, z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x + dx / 2, y, z).ax);eyerr[idx] =ey[idx] -(-yderiv(plane_wave<CCTK_REAL>, dy)(t + dt / 2, x, y + dy / 2, z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x, y + dy / 2, z).ay);ezerr[idx] =ez[idx] -(-zderiv(plane_wave<CCTK_REAL>, dz)(t + dt / 2, x, y, z + dz / 2).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, x, y, z + dz / 2).az);bxerr[idx] =bx[idx] -(yderiv(plane_wave<CCTK_REAL>, dy)(t, x, y + dy / 2, z + dz / 2).az -zderiv(plane_wave<CCTK_REAL>, dz)(t, x, y + dy / 2, z + dz / 2).ay);byerr[idx] =by[idx] -(zderiv(plane_wave<CCTK_REAL>, dz)(t, x + dx / 2, y, z + dz / 2).ax -xderiv(plane_wave<CCTK_REAL>, dx)(t, x + dx / 2, y, z + dz / 2).az);bzerr[idx] =bz[idx] -(xderiv(plane_wave<CCTK_REAL>, dx)(t, x, y + dy / 2, z + dz / 2).ay -yderiv(plane_wave<CCTK_REAL>, dy)(t, x, y + dy / 2, z + dz / 2).ax);
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {exerr_(p.i, p.j, p.k) =ex_(p.i, p.j, p.k) -(-xderiv(plane_wave<CCTK_REAL>, p.dx)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ax);});Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {eyerr_(p.i, p.j, p.k) =ey_(p.i, p.j, p.k) -(-yderiv(plane_wave<CCTK_REAL>, p.dy)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).ay);});Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {ezerr_(p.i, p.j, p.k) =ez_(p.i, p.j, p.k) -(-zderiv(plane_wave<CCTK_REAL>, p.dz)(t + dt / 2, p.x, p.y, p.z).phi -tderiv(plane_wave<CCTK_REAL>, dt)(t + dt / 2, p.x, p.y, p.z).az);});Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {bxerr_(p.i, p.j, p.k) =bx_(p.i, p.j, p.k) -(yderiv(plane_wave<CCTK_REAL>, p.dy)(t, p.x, p.y, p.z).az -zderiv(plane_wave<CCTK_REAL>, p.dz)(t, p.x, p.y, p.z).ay);});Loop::loop_all<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {byerr_(p.i, p.j, p.k) =by_(p.i, p.j, p.k) -(zderiv(plane_wave<CCTK_REAL>, p.dz)(t, p.x, p.y, p.z).ax -xderiv(plane_wave<CCTK_REAL>, p.dx)(t, p.x, p.y, p.z).az);});Loop::loop_all<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {bzerr_(p.i, p.j, p.k) =bz_(p.i, p.j, p.k) -(xderiv(plane_wave<CCTK_REAL>, p.dx)(t, p.x, p.y, p.z).ay -yderiv(plane_wave<CCTK_REAL>, p.dy)(t, p.x, p.y, p.z).ax);});
Loop::loop_all<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) { jxerr_(p.i, p.j, p.k) = 0.0; });Loop::loop_all<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) { jyerr_(p.i, p.j, p.k) = 0.0; });Loop::loop_all<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) { jzerr_(p.i, p.j, p.k) = 0.0; });
Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phi[idx] = standing(t, x, y, z);// phi_p[idx] = standing(t - dt, x, y, z);psi[idx] = timederiv(standing, dt)(t, x, y, z);
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phi[p.idx] = standing(t, p.x, p.y, p.z);// phi_p[p.idx] = standing(t - dt, p.x, p.y, p.z);psi[p.idx] = timederiv(standing, dt)(t, p.x, p.y, p.z);
Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phi[idx] = periodic_gaussian(t, x, y, z);// phi_p[idx] = periodic_gaussian(t - dt, x, y, z);psi[idx] = timederiv(periodic_gaussian, dt)(t, x, y, z);
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phi[p.idx] = periodic_gaussian(t, p.x, p.y, p.z);// phi_p[p.idx] = periodic_gaussian(t - dt, p.x, p.y, p.z);psi[p.idx] = timederiv(periodic_gaussian, dt)(t, p.x, p.y, p.z);
Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phi[idx] = gaussian(t, x, y, z);// phi_p[idx] = gaussian(t - dt, x, y, z);psi[idx] = timederiv(gaussian, dt)(t, x, y, z);
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phi[p.idx] = gaussian(t, p.x, p.y, p.z);// phi_p[p.idx] = gaussian(t - dt, p.x, p.y, p.z);psi[p.idx] = timederiv(gaussian, dt)(t, p.x, p.y, p.z);
Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phi[idx] = central_potential(t, x, y, z);// phi_p[idx] = central_potential(t - dt, x, y, z);psi[idx] = timederiv(central_potential, dt)(t, x, y, z);
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phi[p.idx] = central_potential(t, p.x, p.y, p.z);// phi_p[p.idx] = central_potential(t - dt, p.x, p.y, p.z);psi[p.idx] = timederiv(central_potential, dt)(t, p.x, p.y, p.z);
constexpr int di = 1;const int dj = di * cctk_ash[0];const int dk = dj * cctk_ash[1];Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
psi[idx] = psi_p[idx] +dt * (ddx_phi + ddy_phi + ddz_phi - pow(mass, 2) * phi_p[idx] +4 * M_PI * central_potential(t, x, y, z));phi[idx] = phi_p[idx] + dt * psi[idx];
psi[p.idx] =psi_p[p.idx] +dt * (ddx_phi + ddy_phi + ddz_phi - pow(mass, 2) * phi_p[p.idx] +4 * M_PI * central_potential(t, p.x, p.y, p.z));phi[p.idx] = phi_p[p.idx] + dt * psi[p.idx];
// (phi_p[idx] - phi_p[idx - ndir * ddir]) / dx[dir];// phi[idx] =// phi_p[idx] -// dt * phi_vel * (r * gradphi + (phi_p[idx] - phi_inf) / r);
// (phi_p[p.idx] - phi_p[p.idx - ndir * ddir]) / dx[dir];// phi[p.idx] =// phi_p[p.idx] -// dt * phi_vel * (r * gradphi + (phi_p[p.idx] - phi_inf) / r);
// (psi_p[idx] - psi_p[idx - ndir * ddir]) / dx[dir];// psi[idx] =// psi_p[idx] -// dt * psi_vel * (r * gradpsi + (psi_p[idx] - psi_inf) / r);
// (psi_p[p.idx] - psi_p[p.idx - ndir * ddir]) / dx[dir];// psi[p.idx] =// psi_p[p.idx] -// dt * psi_vel * (r * gradpsi + (psi_p[p.idx] - psi_inf) / r);
(phi_p[idx] - phi_p[idx - ndir * ddir]) / dx[dir];psi[idx] = -phi_vel * (r * gradphi + (phi_p[idx] - phi_inf) / r);phi[idx] = phi_p[idx] + dt * psi[idx];
(phi_p[p.idx] - phi_p[p.idx - ndir * ddir]) / dx[dir];psi[p.idx] = -phi_vel * (r * gradphi + (phi_p[p.idx] - phi_inf) / r);phi[p.idx] = phi_p[p.idx] + dt * psi[p.idx];
Loop::loop_bnd(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phi[idx] = central_potential(t, x, y, z);psi[idx] = timederiv(central_potential, dt)(t, x, y, z);
Loop::loop_bnd<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phi[p.idx] = central_potential(t, p.x, p.y, p.z);psi[p.idx] = timederiv(central_potential, dt)(t, p.x, p.y, p.z);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {assert(CCTK_isfinite(phi_p[idx]));assert(CCTK_isfinite(psi_p[idx]));
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {assert(CCTK_isfinite(phi_p[p.idx]));assert(CCTK_isfinite(psi_p[p.idx]));
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {if (!CCTK_isfinite(phi[idx]) || !CCTK_isfinite(psi[idx])) {CCTK_VINFO("level %d idx=[%d,%d,%d] gidx=[%d,%d,%d] (%g,%g)", lev, i, j,k, cctk_lbnd[0] + i, cctk_lbnd[1] + j, cctk_lbnd[2] + k,phi[idx], psi[idx]);
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {if (!CCTK_isfinite(phi[p.idx]) || !CCTK_isfinite(psi[p.idx])) {CCTK_VINFO("level %d idx=[%d,%d,%d] gidx=[%d,%d,%d] (%g,%g)", lev, p.i,p.j, p.k, cctk_lbnd[0] + p.i, cctk_lbnd[1] + p.j,cctk_lbnd[2] + p.k, phi[p.idx], psi[p.idx]);
(phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / pow(dx[2], 2);CCTK_REAL psi_n = psi[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
(phi[p.idx - p.dk] - 2 * phi[p.idx] + phi[p.idx + p.dk]) / pow(p.dz, 2);CCTK_REAL psi_n = psi[p.idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
CCTK_REAL dt_phi_m = psi_p[idx];CCTK_REAL dx_phi_p = (phi[idx + di] - phi[idx]) / dx[0];CCTK_REAL dx_phi_m = (phi[idx] - phi[idx - di]) / dx[0];CCTK_REAL dy_phi_p = (phi[idx + dj] - phi[idx]) / dx[1];CCTK_REAL dy_phi_m = (phi[idx] - phi[idx - dj]) / dx[1];CCTK_REAL dz_phi_p = (phi[idx + dk] - phi[idx]) / dx[2];CCTK_REAL dz_phi_m = (phi[idx] - phi[idx - dk]) / dx[2];eps[idx] = (pow(dt_phi_p, 2) + pow(dt_phi_m, 2) + pow(dx_phi_p, 2) +pow(dx_phi_m, 2) + pow(dy_phi_p, 2) + pow(dy_phi_m, 2) +pow(dz_phi_p, 2) + pow(dz_phi_m, 2)) /4;
CCTK_REAL dt_phi_m = psi_p[p.idx];CCTK_REAL dx_phi_p = (phi[p.idx + p.di] - phi[p.idx]) / p.dx;CCTK_REAL dx_phi_m = (phi[p.idx] - phi[p.idx - p.di]) / p.dx;CCTK_REAL dy_phi_p = (phi[p.idx + p.dj] - phi[p.idx]) / p.dy;CCTK_REAL dy_phi_m = (phi[p.idx] - phi[p.idx - p.dj]) / p.dy;CCTK_REAL dz_phi_p = (phi[p.idx + p.dk] - phi[p.idx]) / p.dz;CCTK_REAL dz_phi_m = (phi[p.idx] - phi[p.idx - p.dk]) / p.dz;eps[p.idx] = (pow(dt_phi_p, 2) + pow(dt_phi_m, 2) + pow(dx_phi_p, 2) +pow(dx_phi_m, 2) + pow(dy_phi_p, 2) + pow(dy_phi_m, 2) +pow(dz_phi_p, 2) + pow(dz_phi_m, 2)) /4;
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phierr[idx] = phi[idx] - standing(t, x, y, z);psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phierr[p.idx] = phi[p.idx] - standing(t, p.x, p.y, p.z);psierr[p.idx] = psi[p.idx] - timederiv(standing, dt)(t, p.x, p.y, p.z);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phierr[idx] = phi[idx] - periodic_gaussian(t, x, y, z);psierr[idx] = psi[idx] - timederiv(periodic_gaussian, dt)(t, x, y, z);
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phierr[p.idx] = phi[p.idx] - periodic_gaussian(t, p.x, p.y, p.z);psierr[p.idx] =psi[p.idx] - timederiv(periodic_gaussian, dt)(t, p.x, p.y, p.z);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phierr[idx] = phi[idx] - gaussian(t, x, y, z);psierr[idx] = psi[idx] - timederiv(gaussian, dt)(t, x, y, z);
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phierr[p.idx] = phi[p.idx] - gaussian(t, p.x, p.y, p.z);psierr[p.idx] = psi[p.idx] - timederiv(gaussian, dt)(t, p.x, p.y, p.z);
Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];phierr[idx] = phi[idx] - central_potential(t, x, y, z);psierr[idx] = psi[idx] - timederiv(central_potential, dt)(t, x, y, z);
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {phierr[p.idx] = phi[p.idx] - central_potential(t, p.x, p.y, p.z);psierr[p.idx] =psi[p.idx] - timederiv(central_potential, dt)(t, p.x, p.y, p.z);