E3MBKFT4GEFDAGZQQW4OROY5F6FWC46G6MRH54GDYTGO7O5YSRIAC
VIK5E6DBUCP5HGDVHDP6SWTHH7ODEUIGKOGSXR7VLVK46LDD4W6QC
BSMJ4V7GV3EOGY4KCSTOJQUOFE2OOCIKQETE4WC2WRNLWBQIBW3QC
GSGI6HIWST43XFEORCMKH6VPEGXHEZG5EK4JUUNWP3XFYUKGNA3AC
KCIWCVZOHG44WBOLKI2XK33WPHPRI5FWCETF4AOGTPZISKCW3CLQC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
BPRNUTY7MHK7LK4EY5MY5OFFG3ABOL7LWXD574L35M74YSQPULFAC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
TOBGHRPKEPSXDN56WGSZNWOMCBVJ4KUSLWYWI56MC2RR3MM3KLZAC
PQB3EKQ6MBCXPTW4HB7SGMSTOTYMB3EFZX2573OFCQI6PSOEKCSQC
2DD222JSYRPHTXKSRXLSOMSCQPZUORNFLLO2P3GMIDELAAMD5MEQC
JD6PQOJ6YYNQYEEWEXO2NM7NVYNBUI6V7ZU6Q3FNHGAT2VYOF5WAC
UUGQGVC4WEKN64WAP7F5QPS2UHGQB5ZLMFRIYNWKMIEBDO3QRX4AC
WE6MDRN5SPK3THM4COLQFE3IUWBCQ5ZYUIAUCBJAZVEMMOVTNBOAC
WASO7G5FJXRXWNH2U2FLUNEKU6VE63OI3HUYP64BVD4LMD6KE7OQC
UTHNLH3J3VG7BBJPOKGE6DY7Z2QHDLY72TR2VMEDQYP73SIWZGKAC
ULVYXPG2IOJ724MTFVKJLQMPMVL5VAJUOFH7BIXPXVJ2SUGOFRFAC
NGD6QWRDHWBMQXL2YXZVF3BFQALSII3LH6BMIOFL7VFMGPLHCZEQC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
for (auto &restrict leveldata : ghext->leveldata) {
const int min_level = current_level == -1 ? 0 : current_level;
const int max_level =
current_level == -1 ? ghext->leveldata.size() : current_level + 1;
for (int level = min_level; level < max_level; ++level) {
auto &restrict leveldata = ghext->leveldata.at(level);
# SCHEDULE HydroToyAMReX_EstimateError AT postinitial
# {
# LANG: C
# WRITES: AMReX::regrid_error(interior)
# } "Estimate local error for regridding initial conditions"
SCHEDULE HydroToyAMReX_Pressure AT initial AFTER HydroToyAMReX_Boundaries
{
LANG: C
READS: conserved(everywhere)
WRITES: primitive(everywhere)
} "Calculate pressure"
SCHEDULE HydroToyAMReX_Fluxes AT initial AFTER HydroToyAMReX_Pressure
{
LANG: C
READS: conserved(everywhere)
READS: primitive(everywhere)
WRITES: flux_x(interior) flux_y(interior) flux_z(interior)
SYNC: flux_x flux_y flux_z
} "Calculate the hydro fluxes"
SCHEDULE HydroToyAMReX_EstimateError AT postinitial AFTER HydroToyAMReX_Boundaries
{
LANG: C
READS: conserved(everywhere)
WRITES: AMReX::regrid_error(interior)
} "Estimate local error for regridding initial conditions"
# SCHEDULE HydroToyAMReX_EstimateError AT poststep
SCHEDULE HydroToyAMReX_Fluxes AT poststep AFTER HydroToyAMReX_Pressure
{
LANG: C
READS: conserved(everywhere)
READS: primitive(everywhere)
WRITES: flux_x(interior) flux_y(interior) flux_z(interior)
SYNC: flux_x flux_y flux_z
} "Calculate the hydro fluxes"
# SCHEDULE HydroToyAMReX_Output AT poststep
# WRITES: AMReX::regrid_error(interior)
# } "Estimate local error for regridding during evolution"
# READS: conserved(everywhere)
# } "Output grid data"
SCHEDULE HydroToyAMReX_EstimateError AT poststep
{
LANG: C
READS: conserved(everywhere)
WRITES: AMReX::regrid_error(interior)
} "Estimate local error for regridding during evolution"
# SCHEDULE HydroToyAMReX_Evolve1 AT evol
# {
# LANG: C
# READS: conserved_p(everywhere)
# READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
# WRITES: conserved(everywhere)
# } "Evolve the hydro equations, step 1"
#
# SCHEDULE HydroToyAMReX_Pressure AT evol AFTER HydroToyAMReX_Evolve1
# {
# LANG: C
# READS: conserved(everywhere)
# WRITES: primitive(everywhere)
# } "Calculate pressure"
#
# SCHEDULE HydroToyAMReX_Fluxes AT evol AFTER HydroToyAMReX_Pressure
# {
# LANG: C
# READS: conserved(everywhere)
# READS: primitive(everywhere)
# WRITES: flux_x(interior) flux_y(interior) flux_z(interior)
# SYNC: flux_x flux_y flux_z
# } "Calculate the hydro fluxes"
#
# SCHEDULE HydroToyAMReX_Evolve AT evol AFTER HydroToyAMReX_Fluxes
# {
# LANG: C
# READS: conserved_p(everywhere)
# READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
# WRITES: conserved(everywhere)
# } "Evolve the hydro equations, step 2"
# Restriction
SCHEDULE HydroToyAMReX_Boundaries AT postrestrict
{
LANG: C
SYNC: conserved
} "Apply boundary conditions after restricting"
fxmomx_(p.I) =
mkflux(momx_, [&](auto I) { return momx_(I) * velx_(I) + press_(I); });
fxmomy_(p.I) = mkflux(momy_, [&](auto I) { return momy_(I) * velx_(I); });
fxmomz_(p.I) = mkflux(momz_, [&](auto I) { return momz_(I) * velx_(I); });
fxmomx_(p.I) = calcflux(
momx_, [&](auto I) { return momx_(I) * velx_(I) + press_(I); });
fxmomy_(p.I) = calcflux(momy_, [&](auto I) { return momy_(I) * velx_(I); });
fxmomz_(p.I) = calcflux(momz_, [&](auto I) { return momz_(I) * velx_(I); });
fymomx_(p.I) = mkflux(momx_, [&](auto I) { return momx_(I) * vely_(I); });
fymomy_(p.I) =
mkflux(momy_, [&](auto I) { return momy_(I) * vely_(I) + press_(I); });
fymomz_(p.I) = mkflux(momz_, [&](auto I) { return momz_(I) * vely_(I); });
fymomx_(p.I) = calcflux(momx_, [&](auto I) { return momx_(I) * vely_(I); });
fymomy_(p.I) = calcflux(
momy_, [&](auto I) { return momy_(I) * vely_(I) + press_(I); });
fymomz_(p.I) = calcflux(momz_, [&](auto I) { return momz_(I) * vely_(I); });
fzmomx_(p.I) = mkflux(momx_, [&](auto I) { return momx_(I) * velz_(I); });
fzmomy_(p.I) = mkflux(momy_, [&](auto I) { return momy_(I) * velz_(I); });
fzmomz_(p.I) =
mkflux(momz_, [&](auto I) { return momz_(I) * velz_(I) + press_(I); });
fzmomx_(p.I) = calcflux(momx_, [&](auto I) { return momx_(I) * velz_(I); });
fzmomy_(p.I) = calcflux(momy_, [&](auto I) { return momy_(I) * velz_(I); });
fzmomz_(p.I) = calcflux(
momz_, [&](auto I) { return momz_(I) * velz_(I) + press_(I); });
auto mkupdate{[&](auto &fx, auto &fy, auto &fz) {
return 0.5 * dt *
((fx(p.I + p.DI(0)) - fx(p.I)) / dx +
(fy(p.I + p.DI(1)) - fy(p.I)) / dy +
(fz(p.I + p.DI(2)) - fz(p.I)) / dz);
auto calcupdate{[&](auto &fx, auto &fy, auto &fz) {
return dt * ((fx(p.I + p.DI(0)) - fx(p.I)) / dx +
(fy(p.I + p.DI(1)) - fy(p.I)) / dy +
(fz(p.I + p.DI(2)) - fz(p.I)) / dz);
momx_(p.I) = momx_p_(p.I) - mkupdate(fxmomx_, fymomx_, fzmomx_);
momy_(p.I) = momy_p_(p.I) - mkupdate(fxmomy_, fymomy_, fzmomy_);
momz_(p.I) = momz_p_(p.I) - mkupdate(fxmomz_, fymomz_, fzmomz_);
momx_(p.I) = momx_p_(p.I) - calcupdate(fxmomx_, fymomx_, fzmomx_);
momy_(p.I) = momy_p_(p.I) - calcupdate(fxmomy_, fymomy_, fzmomy_);
momz_(p.I) = momz_p_(p.I) - calcupdate(fxmomz_, fymomz_, fzmomz_);
const CCTK_REAL dt = CCTK_DELTA_TIME;
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, 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<const CCTK_REAL, 1, 1, 1> rho_p_(cctkGH, rho_p);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momx_p_(cctkGH, momx_p);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momy_p_(cctkGH, momy_p);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momz_p_(cctkGH, momz_p);
const Loop::GF3D<const CCTK_REAL, 1, 1, 1> etot_p_(cctkGH, etot_p);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxrho_(cctkGH, fxrho);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomx_(cctkGH, fxmomx);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomy_(cctkGH, fxmomy);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomz_(cctkGH, fxmomz);
const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxetot_(cctkGH, fxetot);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fyrho_(cctkGH, fyrho);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomx_(cctkGH, fymomx);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomy_(cctkGH, fymomy);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomz_(cctkGH, fymomz);
const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fyetot_(cctkGH, fyetot);
const Loop::GF3D<CCTK_REAL, 1, 1, 1> regrid_error_(cctkGH, regrid_error);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzrho_(cctkGH, fzrho);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomx_(cctkGH, fzmomx);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomy_(cctkGH, fzmomy);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomz_(cctkGH, fzmomz);
const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzetot_(cctkGH, fzetot);
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);
// Transport
// dt rho + d_i (rho vel^i) = 0
// dt mom_j + d_i (mom_j vel^i) = 0
// dt etot + d_i (etot vel^i) = 0
Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
auto mkupdate{[&](auto &fx, auto &fy, auto &fz) {
return dt * ((fx(p.I + p.DI(0)) - fx(p.I)) / dx +
(fy(p.I + p.DI(1)) - fy(p.I)) / dy +
(fz(p.I + p.DI(2)) - fz(p.I)) / dz);
Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
auto calcerr{[&](auto &var_) {
CCTK_REAL err{0};
for (int d = 0; d < dim; ++d) {
auto varm = var_(p.I - p.DI(d));
auto var0 = var_(p.I);
auto varp = var_(p.I + p.DI(d));
err = fmax(err, fabs(varm - 2 * var0 + varp));
}
return err;
rho_(p.I) = rho_p_(p.I) - mkupdate(fxrho_, fyrho_, fzrho_);
momx_(p.I) = momx_p_(p.I) - mkupdate(fxmomx_, fymomx_, fzmomx_);
momy_(p.I) = momy_p_(p.I) - mkupdate(fxmomy_, fymomy_, fzmomy_);
momz_(p.I) = momz_p_(p.I) - mkupdate(fxmomz_, fymomz_, fzmomz_);
etot_(p.I) = etot_p_(p.I) - mkupdate(fxetot_, fyetot_, fzetot_);
regrid_error_(p.I) =
fmax(fmax(fmax(fmax(calcerr(rho_), calcerr(momx_)), calcerr(momy_)),
calcerr(momz_)),
calcerr(etot_));