ULVYXPG2IOJ724MTFVKJLQMPMVL5VAJUOFH7BIXPXVJ2SUGOFRFAC } "Density"
momx momy momzetot} "Conserved variables"CCTK_REAL primitive TYPE=gf TIMELEVELS=1 TAGS='index={1 1 1}'{pressvelx vely velz} "Primitive variables"CCTK_REAL flux_x TYPE=gf TIMELEVELS=1 TAGS='index={0 1 1}'{fxrhofxmomx fxmomy fxmomzfxetot} "Fluxes in x direction"
ActiveThorns = "AMReXFormalineIOUtilHydroToyAMReXSystemTopologyTimerReport"$nlevels = 1$ncells = 32Cactus::cctk_show_schedule = no# Cactus::terminate = "iteration"# Cactus::cctk_itlast = 0Cactus::terminate = "time"Cactus::cctk_final_time = 2.0AMReX::verbose = noAMReX::xmin = -1.0AMReX::ymin = -1.0AMReX::zmin = -1.0AMReX::xmax = +1.0AMReX::ymax = +1.0AMReX::zmax = +1.0AMReX::ncells_x = $ncellsAMReX::ncells_y = $ncellsAMReX::ncells_z = $ncellsAMReX::max_num_levels = $nlevelsAMReX::regrid_every = 16AMReX::regrid_error_threshold = 0.1AMReX::dtfac = 0.5HydroToyAMReX::setup = "sound wave"IO::out_dir = $parfileIO::out_every = $ncells * 2 ** ($nlevels - 1) / 32# AMReX::out_tsv = yesTimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32TimerReport::out_filename = "TimerReport"TimerReport::output_all_timers_together = yesTimerReport::output_all_timers_readable = yesTimerReport::n_top_timers = 50
ActiveThorns = "AMReXFormalineIOUtilHydroToyAMReXSystemTopologyTimerReport"$nlevels = 1$ncells = 32Cactus::cctk_show_schedule = no# Cactus::terminate = "iteration"# Cactus::cctk_itlast = 0Cactus::terminate = "time"Cactus::cctk_final_time = 2.0AMReX::verbose = noAMReX::xmin = -1.0AMReX::ymin = -1.0AMReX::zmin = -1.0AMReX::xmax = +1.0AMReX::ymax = +1.0AMReX::zmax = +1.0AMReX::ncells_x = $ncellsAMReX::ncells_y = $ncellsAMReX::ncells_z = $ncellsAMReX::max_num_levels = $nlevelsAMReX::regrid_every = 16AMReX::regrid_error_threshold = 0.1AMReX::dtfac = 0.5IO::out_dir = $parfileIO::out_every = $ncells * 2 ** ($nlevels - 1) / 32# AMReX::out_tsv = yesTimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32TimerReport::out_filename = "TimerReport"TimerReport::output_all_timers_together = yesTimerReport::output_all_timers_readable = yesTimerReport::n_top_timers = 50
KEYWORD setup "Initial setup" STEERABLE=never{"equilibrium" :: """sound wave" :: ""} "equilibrium"CCTK_REAL amplitude "Wave amplitude" STEERABLE=never{0.0:* :: ""} 1.0e-3CCTK_REAL gamma "EOS parameter" STEERABLE=always{0.0:* :: ""} 1.6666666666666667
# SCHEDULE HydroToyAMReX_EstimateError AT postinitial# {# LANG: C# WRITES: AMReX::regrid_error(interior)# } "Estimate local error for regridding initial conditions"SCHEDULE HydroToyAMReX_Pressure AT poststep{LANG: CREADS: conserved(everywhere)WRITES: primitive(everywhere)} "Calculate pressure"
SCHEDULE HydroToyAMReX_Fluxes AT poststep AFTER HydroToyAMReX_Pressure{LANG: CREADS: 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 poststep# {# LANG: C# WRITES: AMReX::regrid_error(interior)# } "Estimate local error for regridding during evolution"
DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;// const CCTK_REAL t = cctk_time;// const CCTK_REAL dt = CCTK_DELTA_TIME;if (CCTK_EQUALS(setup, "equilibrium")) {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;});} else if (CCTK_EQUALS(setup, "sound wave")) {Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {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;});} else {assert(0);}}extern "C" void HydroToyAMReX_Pressure(CCTK_ARGUMENTS) {
// Equation of state: p = (gamma - 1) e// vel^j = delta^j_i mom_i / rhoLoop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {CCTK_REAL ekin =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;
const CCTK_REAL t = cctk_time;
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];});}// Minmod slope limitertemplate <typename T> T minmod(T a, T b) {// Same sign: return value with smaller magnitude// Different signs: return zeroreturn (copysign(T(0.5), a) + copysign(T(0.5), b)) * fmin(fabs(a), fabs(b));}// Van Leer slope limitertemplate <typename T> T vanleer(T a, T b) {return ifthen(a * b > 0, 2 * a * b / (a + b), T(0));}// MC2 slope limitertemplate <typename T> T mc2(T a, T b) {return (copysign(T(0.25), a) + copysign(T(0.25), b)) *fmin(4 * fmin(fabs(a), fabs(b)), fabs(a) + fabs(b));}// Superbee slope limitertemplate <typename T> T superbee(T a, T b) {return (copysign(T(0.5), a) + copysign(T(0.5), b)) *fmax(fmin(2 * fabs(a), fabs(b)), fmin(fabs(a), 2 * fabs(b)));}// LLF (local Lax-Friedrichs) Riemann solvertemplate <typename T>T llf(T lambda_m, T lambda_p, T var_m, T var_p, T flux_m, T flux_p) {return T(0.5) * ((flux_m + flux_p) -fmax(fabs(lambda_m), fabs(lambda_p)) * (var_p - var_m));}// HLLE (Harten, Lax, van Leer, Einfeldt) Riemann solvertemplate <typename T>T hlle(T lambda_m, T lambda_p, T var_m, T var_p, T flux_m, T flux_p) {// var_m and var_p are probably swappedassert(0);return (lambda_p * flux_p - lambda_m * flux_m +lambda_p * lambda_m * (var_m - var_p)) /(lambda_p - lambda_m);}extern "C" void HydroToyAMReX_Fluxes(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;
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;
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> press_(cctkGH, press);const Loop::GF3D<const CCTK_REAL, 1, 1, 1> velx_(cctkGH, velx);const Loop::GF3D<const CCTK_REAL, 1, 1, 1> vely_(cctkGH, vely);const Loop::GF3D<const CCTK_REAL, 1, 1, 1> velz_(cctkGH, velz);const Loop::GF3D<CCTK_REAL, 0, 1, 1> fxrho_(cctkGH, fxrho);const Loop::GF3D<CCTK_REAL, 0, 1, 1> fxmomx_(cctkGH, fxmomx);const Loop::GF3D<CCTK_REAL, 0, 1, 1> fxmomy_(cctkGH, fxmomy);const Loop::GF3D<CCTK_REAL, 0, 1, 1> fxmomz_(cctkGH, fxmomz);const Loop::GF3D<CCTK_REAL, 0, 1, 1> fxetot_(cctkGH, fxetot);const Loop::GF3D<CCTK_REAL, 1, 0, 1> fyrho_(cctkGH, fyrho);const Loop::GF3D<CCTK_REAL, 1, 0, 1> fymomx_(cctkGH, fymomx);const Loop::GF3D<CCTK_REAL, 1, 0, 1> fymomy_(cctkGH, fymomy);const Loop::GF3D<CCTK_REAL, 1, 0, 1> fymomz_(cctkGH, fymomz);const Loop::GF3D<CCTK_REAL, 1, 0, 1> fyetot_(cctkGH, fyetot);const Loop::GF3D<CCTK_REAL, 1, 1, 0> fzrho_(cctkGH, fzrho);const Loop::GF3D<CCTK_REAL, 1, 1, 0> fzmomx_(cctkGH, fzmomx);const Loop::GF3D<CCTK_REAL, 1, 1, 0> fzmomy_(cctkGH, fzmomy);const Loop::GF3D<CCTK_REAL, 1, 1, 0> fzmomz_(cctkGH, fzmomz);const Loop::GF3D<CCTK_REAL, 1, 1, 0> fzetot_(cctkGH, fzetot);// frho^i = rho vel^i// fmom^i_j = mom_j vel^i + delta^i_j press// fetot^i = (etot + press) vel^iLoop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {auto mkflux{[&](auto &u, auto f) {auto lambda_m = -dx / dt;auto lambda_p = dx / dt;auto var_m = u(p.I - p.DI(0));auto var_p = u(p.I);auto flux_m = f(p.I - p.DI(0));auto flux_p = f(p.I);return llf(lambda_m, lambda_p, var_m, var_p, flux_m, flux_p);}};fxrho_(p.I) = mkflux(rho_, [&](auto I) { return rho_(I) * velx_(I); });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); });fxetot_(p.I) = mkflux(etot_, [&](auto I) { return (etot_(I) + press_(I)) * velx_(I); });});Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {auto mkflux{[&](auto &u, auto f) {auto lambda_m = -dy / dt;auto lambda_p = dy / dt;auto var_m = u(p.I - p.DI(1));auto var_p = u(p.I);auto flux_m = f(p.I - p.DI(1));auto flux_p = f(p.I);return llf(lambda_m, lambda_p, var_m, var_p, flux_m, flux_p);}};fyrho_(p.I) = mkflux(rho_, [&](auto I) { return rho_(I) * vely_(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); });fyetot_(p.I) = mkflux(etot_, [&](auto I) { return (etot_(I) + press_(I)) * vely_(I); });});Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {auto mkflux{[&](auto &u, auto f) {auto lambda_m = -dz / dt;auto lambda_p = dz / dt;auto var_m = u(p.I - p.DI(2));auto var_p = u(p.I);auto flux_m = f(p.I - p.DI(2));auto flux_p = f(p.I);return llf(lambda_m, lambda_p, var_m, var_p, flux_m, flux_p);}};fzrho_(p.I) = mkflux(rho_, [&](auto I) { return rho_(I) * velz_(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); });fzetot_(p.I) = mkflux(etot_, [&](auto I) { return (etot_(I) + press_(I)) * velz_(I); });
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, 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<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) = 0Loop::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) {// dt rho + d_i (rho vel^i) = 0rho[p.idx] =rho_p[p.idx] - dt * ((momx_p[p.idx + di] - momx_p[p.idx]) / dx +(momy_p[p.idx + dj] - momy_p[p.idx]) / dy +(momz_p[p.idx + dk] - momz_p[p.idx]) / dz);
rho_(p.I) = rho_p_(p.I) - mkupdate(fxrho_, fyrho_, fzrho_);
// // dt mom^i + d_j (mom^i vel^j) = 0// momx[p.idx] = momx_p[p.idx] - dt * (// momx_p[p.idx] * momx_p[p.idx]// );
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_);