BJDGFYBMECTJG7BHLNHLSCUCBVYHAY6OGY37FIJP6JDGNDXQNQVAC
2D74GU7KGIN5GTSMCY7752YYQN62T5UY2FIKMWKYVMGKS2CDAJKAC
BBUXGUGVTYMZXMDU2KTRUCACCBWTRFVVVMJUJFCQJPDVFU3DH53QC
OHAARXLBTFOIV5FOTY7IDMOIAE275CKXCE6XQMWV3EC63YCSMIWAC
GQVQJCNQNO2KD7ZMC7RESCUAMUAP7OED6CTA6SYLZKQGXKXZ6T3QC
BZ7T742UYYUIK2X3Y6UWJ2G33ER4PJJBA7YLCVAQKAJQUUECHVRQC
A7ETPFXEHA2RM4LINSBVMJJ3G62NF7Q5ZQOKJPNJK3YOQ5WS5HKAC
B2MBQPT3B3N2EVCTHX34PKGIQSLUU622HEEVZNHQFL2PAJFFX3SQC
GHWN22UDWSZ32PM6RFYH6IS5VHJELGXKI2S2HTUINKYSU25NFFPQC
XFO76KCYC4GCJ353GCS3SSFNMMHUPUOIEHBBS26JKGZZHJSQS2TQC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC
2DKSL6DKZAIYQUJGDULORCKU5K4Z5Z3W4RIKQYDSLKMCNQNDZFBAC
5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC
KCIWCVZOHG44WBOLKI2XK33WPHPRI5FWCETF4AOGTPZISKCW3CLQC
33IC3UHCEPZLGS5ACS2JXHGT6CRU5LXU6PM6RDHCERPOIELVRVXQC
5IAXY3XZJTRMMVT2OVIJ6OXQJI6OJPTPCHHA4IVLVMHANCCC5NKAC
BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
QNV4LD7UGYSSNYDXGJC6SRP6Y3USUVDQLEERBMBWRC7LILWREB7AC
U77PE56ICORZNQW33NXGSEMW7GDHCSSZ4EXB6OHBJSHEG6WHYSSQC
RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC
TVBD244E7Q7WV44CRBTFST535NUP3JAZH6OLL4IKDR3OWEXSU7HAC
VAF66DTVLDWNG7N2PEQYEH4OH5SPSMFBXKPR2PU67IIM6CVPCJ7AC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
BSMJ4V7GV3EOGY4KCSTOJQUOFE2OOCIKQETE4WC2WRNLWBQIBW3QC
GECUITHDXKCWB7HBCM7EA5Q56JDDWUVUWHMW2K6OM7UW36DFAZ3QC
2DD222JSYRPHTXKSRXLSOMSCQPZUORNFLLO2P3GMIDELAAMD5MEQC
UZAKARMGORRQG733ZUPJOEGL5FG243I32NCC2SRSFDCZKUQ5A52QC
PQB3EKQ6MBCXPTW4HB7SGMSTOTYMB3EFZX2573OFCQI6PSOEKCSQC
IVHURSHY4636OGIF3PNDO5CWOVRLJ75M4LP65J6I2E6KAM4QKF4AC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
UTHNLH3J3VG7BBJPOKGE6DY7Z2QHDLY72TR2VMEDQYP73SIWZGKAC
ULVYXPG2IOJ724MTFVKJLQMPMVL5VAJUOFH7BIXPXVJ2SUGOFRFAC
NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC
WJFK2JGGVKQN4LVHZNB4SOKH2GAE5BODBU2D6Y4CO76EIHSW6EVAC
GSGI6HIWST43XFEORCMKH6VPEGXHEZG5EK4JUUNWP3XFYUKGNA3AC
E3MBKFT4GEFDAGZQQW4OROY5F6FWC46G6MRH54GDYTGO7O5YSRIAC
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
CCTK_REAL *restrict const err = grid.ptr(vars, vi);
const Array4<const CCTK_REAL> &err_array4 =
groupdata.mfab.at(tl)->array(mfi);
const GF3D1<const CCTK_REAL> err_ = grid.gf3d(err_array4, vi);
grid.loop_int(groupdata.indextype, [&](const Loop::PointDesc &p) {
tags_array4(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;
});
grid.loop_idx(
where_t::interior, groupdata.indextype, [&](const Loop::PointDesc &p) {
tags_array4(grid.cactus_offset.x + p.i, grid.cactus_offset.y + p.j,
grid.cactus_offset.z + p.k) =
err_(p.I) >= regrid_error_threshold ? TagBox::SET : TagBox::CLEAR;
});
// grid.loop_bnd(groupdata.indextype, [&](const Loop::PointDesc &p) {
// tags_array4(grid.cactus_offset.x + p.i, grid.cactus_offset.y + p.j,
// grid.cactus_offset.z + p.k) = TagBox::CLEAR;
// });
if (false) {
grid.loop_idx(where_t::boundary, groupdata.indextype,
[&](const Loop::PointDesc &p) {
tags_array4(grid.cactus_offset.x + p.i,
grid.cactus_offset.y + p.j,
grid.cactus_offset.z + p.k) = TagBox::CLEAR;
});
}
array<int, dim> get_group_nghostzones(const int gi) {
DECLARE_CCTK_PARAMETERS;
assert(gi >= 0);
const int tags = CCTK_GroupTagsTableI(gi);
assert(tags >= 0);
array<CCTK_INT, dim> nghostzones;
int iret =
Util_TableGetIntArray(tags, dim, nghostzones.data(), "nghostzones");
if (iret == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
// default: use driver parameter
nghostzones = {ghost_size, ghost_size, ghost_size};
} else if (iret >= 0) {
assert(iret == dim);
} else {
assert(0);
}
return nghostzones;
}
groupdata.mfab.at(tl) =
make_unique<MultiFab>(gba, dm, groupdata.numvars, ghost_size);
groupdata.mfab.at(tl) = make_unique<MultiFab>(
gba, dm, groupdata.numvars, IntVect(groupdata.nghostzones));
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;
});
const GF3D1<CCTK_REAL> ptr_ = grid.gf3d(vars, vi);
grid.loop_idx(
where_t::everywhere, groupdata.indextype, groupdata.nghostzones,
[&](const Loop::PointDesc &p) { ptr_(p.I) = 0.0 / 0.0; });
string banner =
"AMR driver provided by CarpetX, using AMReX " + amrex::Version();
string banner = "AMR driver provided by CarpetX, using AMReX " +
amrex::Version() +
" ("
#ifdef AMREX_USE_MPI
"MPI, "
#else
"no MPI, "
#endif
#ifdef AMREX_USE_OMP
"OpenMP, "
#else
"no OpenMP, "
#endif
#ifdef AMREX_USE_GPU
"Accelerators, "
#else
"no Accelerators, "
#endif
#ifdef AMREX_USE_ASSERTION
"DEBUG, "
#else
"OPTIMIZED, "
#endif
")";
GridPtrDesc grid(leveldata, mfi);
const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
vector<const CCTK_REAL *> ptrs(groupdata.numvars);
GridPtrDesc1 grid(leveldata, groupdata, mfi);
const Array4<const CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
vector<GF3D1<const CCTK_REAL> > ptrs_;
ptrs_.reserve(groupdata.numvars);
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";
});
Loop::loop_idx(cctkGH, where_t::everywhere, 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.I);
file << "\n";
});
Loop::loop_all<1, 1, 1>(
cctkGH, [&](const Loop::PointDesc &p) { regrid_error[p.idx] = 0; });
const Loop::GF3D<CCTK_REAL, 1, 1, 1> regrid_error_(cctkGH, regrid_error);
Loop::loop<1, 1, 1>(
cctkGH, where_t::everywhere,
[&](const Loop::PointDesc &p) { regrid_error_(p.I) = 0; });
template <typename T> struct GF3D1 {
typedef T value_type;
T *restrict ptr;
#ifdef CCTK_DEBUG
array<int, dim> imin, imax;
array<int, dim> ash;
#endif
static constexpr int di = 1;
int dj, dk;
int off;
inline GF3D1(T *restrict ptr, const array<int, dim> &imin,
const array<int, dim> &imax, const array<int, dim> &ash)
: ptr(ptr),
#ifdef CCTK_DEBUG
imin(imin), imax(imax), ash(ash),
#endif
dj(di * ash[0]), dk(dj * ash[1]),
off(imin[0] * di + imin[1] * dj + imin[2] * dk) {
}
inline GF3D1(const cGH *restrict cctkGH, const array<int, dim> &indextype,
const array<int, dim> &nghostzones, T *restrict ptr) {
for (int d = 0; d < dim; ++d)
assert(indextype[d] == 0 || indextype[d] == 1);
for (int d = 0; d < dim; ++d) {
assert(nghostzones[d] >= 0);
assert(nghostzones[d] <= cctkGH->cctk_nghostzones[d]);
}
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
imin[d] = cctkGH->cctk_nghostzones[d] - nghostzones[d];
imax[d] = cctkGH->cctk_lsh[d] + (1 - indextype[d]) -
(cctkGH->cctk_nghostzones[d] - nghostzones[d]);
}
array<int, dim> ash;
for (int d = 0; d < dim; ++d)
ash[d] = cctkGH->cctk_ash[d] + (1 - indextype[d]) -
2 * (cctkGH->cctk_nghostzones[d] - nghostzones[d]);
*this = GF3D1(ptr, imin, imax, ash);
}
inline int offset(int i, int j, int k) const {
// These index checks prevent vectorization. We thus only enable
// them in debug mode.
#ifdef CCTK_DEBUG
assert(i >= imin[0] && i < imax[0]);
assert(j >= imin[1] && j < imax[1]);
assert(k >= imin[2] && k < imax[2]);
#endif
return i * di + j * dj + k * dk - off;
}
inline T &restrict operator()(int i, int j, int k) const {
return ptr[offset(i, j, k)];
}
inline T &restrict operator()(const vect<int, dim> &I) const {
return ptr[offset(I[0], I[1], I[2])];
}
};
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);
template <int CI, int CJ, int CK, typename F>
void loop(where_t where, const array<int, dim> &group_nghostzones,
const F &f) const {
switch (where) {
case where_t::everywhere:
return loop_all<CI, CJ, CK>(group_nghostzones, f);
case where_t::interior:
return loop_int<CI, CJ, CK>(group_nghostzones, f);
case where_t::boundary:
return loop_bnd<CI, CJ, CK>(group_nghostzones, 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);
}
void loop_idx(where_t where, const array<int, dim> &indextype,
const F &f) const {
loop_idx(where, indextype, nghostzones, 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 <typename F>
void loop_idx(const cGH *cctkGH, where_t where,
const array<int, dim> &indextype,
const array<int, dim> &nghostzones, const F &f) {
GridDescBase(cctkGH).loop_idx(where, indextype, nghostzones, 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_idx(const cGH *cctkGH, where_t where,
const array<int, dim> &indextype, const F &f) {
GridDescBase(cctkGH).loop_idx(where, indextype, f);
template <typename F>
void loop_all(const cGH *cctkGH, const array<int, dim> &indextype, const F &f) {
GridDescBase(cctkGH).loop_all(indextype, f);
// Keep these for convenience
template <int CI, int CJ, int CK, typename F>
void loop_all(const cGH *cctkGH, const F &f) {
loop<CI, CJ, CK>(cctkGH, where_t::everywhere, 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 <int CI, int CJ, int CK, typename F>
void loop_int(const cGH *cctkGH, const F &f) {
loop<CI, CJ, CK>(cctkGH, where_t::interior, f);
template <typename F>
void loop_bnd(const cGH *cctkGH, const array<int, dim> &indextype, const F &f) {
GridDescBase(cctkGH).loop_bnd(indextype, f);
template <int CI, int CJ, int CK, typename F>
void loop_bnd(const cGH *cctkGH, const F &f) {
loop<CI, CJ, CK>(cctkGH, where_t::boundary, f);
grid.loop_int(groupdata.indextype, [&](const Loop::PointDesc &p) {
if (!finemask_array4 || !(*finemask_array4)(grid.cactus_offset.x + p.i,
grid.cactus_offset.y + p.j,
grid.cactus_offset.z + p.k))
red = reduction<T>(red, reduction<T>(dV, ptr[p.idx]));
});
grid.loop_idx(
where_t::interior, groupdata.indextype, [&](const Loop::PointDesc &p) {
if (!finemask_array4 || !(*finemask_array4)(grid.cactus_offset.x + p.i,
grid.cactus_offset.y + p.j,
grid.cactus_offset.z + p.k))
red = reduction<T>(red, reduction<T>(dV, ptr_(p.I)));
});
const CCTK_REAL *restrict ptr = grid.ptr(vars, vi);
red += reduce_array(geom, groupdata, grid, ptr, finemask_array4.get());
const GF3D1<const CCTK_REAL> ptr_ = grid.gf3d(vars, vi);
red += reduce_array(geom, groupdata, grid, ptr_, finemask_array4.get());
vector<TileBox> thread_local_tilebox;
struct thread_local_info_t {
// TODO: store only MFIter here; recalculate other things from it
cGH cctkGH;
TileBox tilebox;
MFIter *mfiter;
unsigned char padding[128]; // Prevent false sharing
};
vector<unique_ptr<thread_local_info_t> > thread_local_info;
vector<unique_ptr<thread_local_info_t> > saved_thread_local_info;
for (int d = 0; d < dim; ++d) {
assert(groupdata.nghostzones[d] >= 0);
assert(groupdata.nghostzones[d] <= nghostzones[d]);
}
for (int d = 0; d < dim; ++d)
gimin[d] = nghostzones[d] - groupdata.nghostzones[d];
for (int d = 0; d < dim; ++d)
gimax[d] = lsh[d] + (1 - groupdata.indextype[d]) -
(nghostzones[d] - groupdata.nghostzones[d]);
for (int d = 0; d < dim; ++d)
gash[d] = ash[d] + (1 - groupdata.indextype[d]) -
2 * (nghostzones[d] - groupdata.nghostzones[d]);
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;
});
}
if (!valid.valid_bnd)
where = where_t::everywhere;
else
where = where_t::interior;
if (valid.valid_bnd) {
grid.loop_all(groupdata.indextype, [&](const Loop::PointDesc &p) {
if (CCTK_BUILTIN_EXPECT(isnan(ptr[p.idx]), false))
found_nan = true;
});
} else {
grid.loop_int(groupdata.indextype, [&](const Loop::PointDesc &p) {
if (CCTK_BUILTIN_EXPECT(isnan(ptr[p.idx]), false))
found_nan = true;
});
}
if (valid.valid_bnd)
where = where_t::everywhere;
else
where = where_t::interior;
if (valid.valid_bnd) {
grid.loop_bnd(groupdata.indextype, [&](const Loop::PointDesc &p) {
if (CCTK_BUILTIN_EXPECT(isnan(ptr[p.idx]), false))
found_nan = true;
});
} else {
if (valid.valid_bnd)
where = where_t::boundary;
else
for (int vi = 0; vi < groupdata.numvars; ++vi) {
cctkGH->data[groupdata.firstvarindex + vi][tl] = grid.ptr(vars, vi);
}
for (int vi = 0; vi < groupdata.numvars; ++vi)
cctkGH->data[groupdata.firstvarindex + vi][tl] = grid1.ptr(vars, vi);
if (groupdata.indextype == array<int, dim>{0, 0, 0})
#warning \
"TODO: Allow different restriction operators, and ensure this is conservative"
// rank: 0: vertex, 1: edge, 2: face, 3: volume
int rank = 0;
for (int d = 0; d < dim; ++d)
rank += groupdata.indextype[d];
switch (rank) {
case 0:
}
};
struct GridPtrDesc1 : GridDesc {
Dim3 cactus_offset;
array<int, dim> gimin, gimax;
array<int, dim> gash;
GridPtrDesc1() = delete;
GridPtrDesc1(const GHExt::LevelData &leveldata,
const GHExt::LevelData::GroupData &groupdata, const MFIter &mfi);
template <typename T> T *ptr(const Array4<T> &vars, int vi) const {
return vars.ptr(cactus_offset.x + gimin[0], cactus_offset.y + gimin[1],
cactus_offset.z + gimin[2], vi);
}
template <typename T>
T &idx(const Array4<T> &vars, int i, int j, int k, int vi) const {
return vars(cactus_offset.x + gimin[0] + i, cactus_offset.y + gimin[1] + i,
cactus_offset.z + gimin[2] + j, vi);
}
template <typename T> GF3D1<T> gf3d(const Array4<T> &vars, int vi) const {
return GF3D1<T>(ptr(vars, vi), gimin, gimax, gash);
void check_valid(
const GHExt::LevelData &leveldata,
const GHExt::LevelData::GroupData &groupdata, int vi, int tl,
const function<string()> &msg);
void check_valid(const GHExt::LevelData &leveldata,
const GHExt::LevelData::GroupData &groupdata, int vi, int tl,
const function<string()> &msg);
SCHEDULE HydroToyAMReX_Output IN HydroToyAMReX_OutputGroup AFTER HydroToyAMReX_Fluxes
{
LANG: C
READS: conserved(interior)
READS: flux_x(interior) flux_y(interior) flux_z(interior)
READS: AMReX::regrid_error(interior)
} "Output grid data"
} # if output
if (output) {
SCHEDULE HydroToyAMReX_Output IN HydroToyAMReX_OutputGroup AFTER HydroToyAMReX_Fluxes
{
LANG: C
READS: conserved(interior)
READS: flux_x(interior) flux_y(interior) flux_z(interior)
READS: AMReX::regrid_error(interior)
} "Output grid data"
} # if output
}
const Loop::GF3D<CCTK_REAL, 1, 1, 1> rho_(cctkGH, rho);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> momx_(cctkGH, momx);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> momy_(cctkGH, momy);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> momz_(cctkGH, momz);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> etot_(cctkGH, etot);
rho[p.idx] = 1.0;
momx[p.idx] = 0.0 + amplitude * sin(M_PI * p.x);
momy[p.idx] = 0.0;
momz[p.idx] = 0.0;
etot[p.idx] = 1.0; // should add kinetic energy here
rho_(p.I) = 1.0;
momx_(p.I) = 0.0 + amplitude * sin(M_PI * p.x);
momy_(p.I) = 0.0;
momz_(p.I) = 0.0;
etot_(p.I) = 1.0; // should add kinetic energy here
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> rho_(cctkGH, rho);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momx_(cctkGH, momx);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momy_(cctkGH, momy);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momz_(cctkGH, momz);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> etot_(cctkGH, etot);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> press_(cctkGH, press);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> velx_(cctkGH, velx);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> vely_(cctkGH, vely);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> velz_(cctkGH, velz);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> eint_(cctkGH, eint);
sqrt(pow(momx[p.idx], 2) + pow(momy[p.idx], 2) + pow(momz[p.idx], 2)) /
(2 * rho[p.idx]);
CCTK_REAL eint = etot[p.idx] - ekin;
press[p.idx] = (gamma - 1) * eint;
sqrt(pow(momx_(p.I), 2) + pow(momy_(p.I), 2) + pow(momz_(p.I), 2)) /
(2 * rho_(p.I));
eint_(p.I) = etot_(p.I) - ekin;
velx[p.idx] = momx[p.idx] / rho[p.idx];
vely[p.idx] = momy[p.idx] / rho[p.idx];
velz[p.idx] = momz[p.idx] / rho[p.idx];
press_(p.I) = (gamma - 1) * eint_(p.I);
velx_(p.I) = momx_(p.I) / rho_(p.I);
vely_(p.I) = momy_(p.I) / rho_(p.I);
velz_(p.I) = momz_(p.I) / rho_(p.I);
#else
const array<int, dim> nghosts{0, 0, 0};
const Loop::GF3D1<CCTK_REAL> fxrho_(cctkGH, {0, 1, 1}, nghosts, fxrho);
const Loop::GF3D1<CCTK_REAL> fxmomx_(cctkGH, {0, 1, 1}, nghosts, fxmomx);
const Loop::GF3D1<CCTK_REAL> fxmomy_(cctkGH, {0, 1, 1}, nghosts, fxmomy);
const Loop::GF3D1<CCTK_REAL> fxmomz_(cctkGH, {0, 1, 1}, nghosts, fxmomz);
const Loop::GF3D1<CCTK_REAL> fxetot_(cctkGH, {0, 1, 1}, nghosts, fxetot);
const Loop::GF3D1<CCTK_REAL> fyrho_(cctkGH, {1, 0, 1}, nghosts, fyrho);
const Loop::GF3D1<CCTK_REAL> fymomx_(cctkGH, {1, 0, 1}, nghosts, fymomx);
const Loop::GF3D1<CCTK_REAL> fymomy_(cctkGH, {1, 0, 1}, nghosts, fymomy);
const Loop::GF3D1<CCTK_REAL> fymomz_(cctkGH, {1, 0, 1}, nghosts, fymomz);
const Loop::GF3D1<CCTK_REAL> fyetot_(cctkGH, {1, 0, 1}, nghosts, fyetot);
const Loop::GF3D1<CCTK_REAL> fzrho_(cctkGH, {1, 1, 0}, nghosts, fzrho);
const Loop::GF3D1<CCTK_REAL> fzmomx_(cctkGH, {1, 1, 0}, nghosts, fzmomx);
const Loop::GF3D1<CCTK_REAL> fzmomy_(cctkGH, {1, 1, 0}, nghosts, fzmomy);
const Loop::GF3D1<CCTK_REAL> fzmomz_(cctkGH, {1, 1, 0}, nghosts, fzmomz);
const Loop::GF3D1<CCTK_REAL> fzetot_(cctkGH, {1, 1, 0}, nghosts, fzetot);
#endif
#else
const array<int, dim> nghosts{0, 0, 0};
const Loop::GF3D1<const CCTK_REAL> fxrho_(cctkGH, {0, 1, 1}, nghosts, fxrho);
const Loop::GF3D1<const CCTK_REAL> fxmomx_(cctkGH, {0, 1, 1}, nghosts,
fxmomx);
const Loop::GF3D1<const CCTK_REAL> fxmomy_(cctkGH, {0, 1, 1}, nghosts,
fxmomy);
const Loop::GF3D1<const CCTK_REAL> fxmomz_(cctkGH, {0, 1, 1}, nghosts,
fxmomz);
const Loop::GF3D1<const CCTK_REAL> fxetot_(cctkGH, {0, 1, 1}, nghosts,
fxetot);
const Loop::GF3D1<const CCTK_REAL> fyrho_(cctkGH, {1, 0, 1}, nghosts, fyrho);
const Loop::GF3D1<const CCTK_REAL> fymomx_(cctkGH, {1, 0, 1}, nghosts,
fymomx);
const Loop::GF3D1<const CCTK_REAL> fymomy_(cctkGH, {1, 0, 1}, nghosts,
fymomy);
const Loop::GF3D1<const CCTK_REAL> fymomz_(cctkGH, {1, 0, 1}, nghosts,
fymomz);
const Loop::GF3D1<const CCTK_REAL> fyetot_(cctkGH, {1, 0, 1}, nghosts,
fyetot);
const Loop::GF3D1<const CCTK_REAL> fzrho_(cctkGH, {1, 1, 0}, nghosts, fzrho);
const Loop::GF3D1<const CCTK_REAL> fzmomx_(cctkGH, {1, 1, 0}, nghosts,
fzmomx);
const Loop::GF3D1<const CCTK_REAL> fzmomy_(cctkGH, {1, 1, 0}, nghosts,
fzmomy);
const Loop::GF3D1<const CCTK_REAL> fzmomz_(cctkGH, {1, 1, 0}, nghosts,
fzmomz);
const Loop::GF3D1<const CCTK_REAL> fzetot_(cctkGH, {1, 1, 0}, nghosts,
fzetot);
#endif
#else
const array<int, dim> nghosts{0, 0, 0};
const Loop::GF3D1<const CCTK_REAL> fxrho_(cctkGH, {0, 1, 1}, nghosts, fxrho);
const Loop::GF3D1<const CCTK_REAL> fxmomx_(cctkGH, {0, 1, 1}, nghosts,
fxmomx);
const Loop::GF3D1<const CCTK_REAL> fxmomy_(cctkGH, {0, 1, 1}, nghosts,
fxmomy);
const Loop::GF3D1<const CCTK_REAL> fxmomz_(cctkGH, {0, 1, 1}, nghosts,
fxmomz);
const Loop::GF3D1<const CCTK_REAL> fxetot_(cctkGH, {0, 1, 1}, nghosts,
fxetot);
#endif