W6W3EC56H4PPREAK7DMX3DAHYKDAHSASGDP63DSU55TRCBBGJXBAC momx momy momzetot} "Conserved variables"CCTK_REAL conserved_rhs TYPE=gf TAGS='index={1 1 1} fluxes="HydroToyCarpetX::flux_x HydroToyCarpetX::flux_y HydroToyCarpetX::flux_z"'{dtrhodtmomx dtmomy dtmomzdtetot} "Time derivatives of conserved variables"CCTK_REAL primitive TYPE=gf TAGS='index={1 1 1} checkpoint="no" restrict="no"'{press
CCTK_REAL conserved_rhs TYPE=gf TAGS='index={1 1 1} checkpoint="no" fluxes="Hydro::flux_x Hydro::flux_y Hydro::flux_z"'{dtdensdtmomx dtmomy dtmomzdtetot} "Time derivatives of conserved variables"CCTK_REAL conserved TYPE=gf TAGS='index={1 1 1} rhs="Hydro::conserved_rhs"'{densmomx momy momzetot} "Conserved variables"
ActiveThorns = "CarpetXFormalineIOUtilHydroODESolversSystemTopologyTimerReport"$nlevels = 3$ncells = 32Cactus::cctk_show_schedule = no# Cactus::terminate = "iteration"# Cactus::cctk_itlast = 0Cactus::terminate = "time"Cactus::cctk_final_time = 1.0CarpetX::verbose = noCarpetX::xmin = -1.0CarpetX::ymin = -1.0CarpetX::zmin = -1.0CarpetX::xmax = +1.0CarpetX::ymax = +1.0CarpetX::zmax = +1.0CarpetX::ncells_x = $ncellsCarpetX::ncells_y = $ncellsCarpetX::ncells_z = $ncellsCarpetX::max_num_levels = $nlevelsCarpetX::regrid_every = 16CarpetX::regrid_error_threshold = 0.01CarpetX::prolongation_type = "conservative"CarpetX::do_reflux = noCarpetX::do_restrict = noCarpetX::restrict_during_sync = yesCarpetX::dtfac = 0.25Hydro::setup = "spherical shock"Hydro::output_every = 1IO::out_dir = $parfileIO::out_every = $ncells * 2 ** ($nlevels - 1) / 32CarpetX::out_tsv = noTimerReport::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
SCHEDULE Hydro_EstimateError AT postinitial
SCHEDULE Hydro_Prim2Con AT postinitial{LANG: CREADS: primitive(everywhere)WRITES: conserved(everywhere)INVALIDATES: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)INVALIDATES: conserved_rhs(everywhere)} "Calculate conserved variables"SCHEDULE Hydro_EstimateError AT postinitial AFTER Hydro_Prim2Con
CCTK_REALVEC ekin = CCTK_REAL(0.5) * rho_inv *sqrt(pow2(momx1) + pow2(momy1) + pow2(momz1));CCTK_REALVEC eint1 = etot1 - ekin;
CCTK_REALVEC ekin1 = CCTK_REAL(0.5) * rho1_inv *sqrt(pow2(momx1) + pow2(momy1) + pow2(momz1));CCTK_REALVEC eint1 = etot1 - ekin1;
calcflux(0, di, dAx, nx, velx, fxrho, fxmomx, fxmomy, fxmomz, fxetot);calcflux(1, dj, dAy, ny, vely, fyrho, fymomx, fymomy, fymomz, fyetot);calcflux(2, dk, dAz, nz, velz, fzrho, fzmomx, fzmomy, fzmomz, fzetot);
calcflux(0, di, dAx, nx, velx, fxdens, fxmomx, fxmomy, fxmomz, fxetot);calcflux(1, dj, dAy, ny, vely, fydens, fymomx, fymomy, fymomz, fyetot);calcflux(2, dk, dAz, nz, velz, fzdens, fzmomx, fzmomy, fzmomz, fzetot);
CCTK_REALVEC momx1 = CCTK_REAL(0.0);CCTK_REALVEC momy1 = CCTK_REAL(0.0);CCTK_REALVEC momz1 = CCTK_REAL(0.0);CCTK_REALVEC etot1 = CCTK_REAL(1.0);
CCTK_REALVEC velx1 = CCTK_REAL(0.0);CCTK_REALVEC vely1 = CCTK_REAL(0.0);CCTK_REALVEC velz1 = CCTK_REAL(0.0);CCTK_REALVEC eint1 = CCTK_REAL(1.0);
momx1.storeu_partial(momx[ind], i, imin[0], imax[0]);momy1.storeu_partial(momy[ind], i, imin[0], imax[0]);momz1.storeu_partial(momz[ind], i, imin[0], imax[0]);etot1.storeu_partial(etot[ind], i, imin[0], imax[0]);
velx1.storeu_partial(velx[ind], i, imin[0], imax[0]);vely1.storeu_partial(vely[ind], i, imin[0], imax[0]);velz1.storeu_partial(velz[ind], i, imin[0], imax[0]);eint1.storeu_partial(eint[ind], i, imin[0], imax[0]);
CCTK_REALVEC momy1 = CCTK_REAL(0.0);CCTK_REALVEC momz1 = CCTK_REAL(0.0);CCTK_REALVEC etot1 = CCTK_REAL(1.0); // should add kinetic energy here
CCTK_REALVEC vely1 = CCTK_REAL(0.0);CCTK_REALVEC velz1 = CCTK_REAL(0.0);CCTK_REALVEC eint1 = CCTK_REAL(1.0);
momx1.storeu_partial(momx[ind], i, imin[0], imax[0]);momy1.storeu_partial(momy[ind], i, imin[0], imax[0]);momz1.storeu_partial(momz[ind], i, imin[0], imax[0]);etot1.storeu_partial(etot[ind], i, imin[0], imax[0]);
velx1.storeu_partial(velx[ind], i, imin[0], imax[0]);vely1.storeu_partial(vely[ind], i, imin[0], imax[0]);velz1.storeu_partial(velz[ind], i, imin[0], imax[0]);eint1.storeu_partial(eint[ind], i, imin[0], imax[0]);
CCTK_REALVEC momx1l = CCTK_REAL(0.0);CCTK_REALVEC momy1l = CCTK_REAL(0.0);CCTK_REALVEC momz1l = CCTK_REAL(0.0);CCTK_REALVEC etot1l = CCTK_REAL(2.0);
CCTK_REALVEC velx1l = CCTK_REAL(0.0);CCTK_REALVEC vely1l = CCTK_REAL(0.0);CCTK_REALVEC velz1l = CCTK_REAL(0.0);CCTK_REALVEC eint1l = CCTK_REAL(2.0);
CCTK_REALVEC momx1r = CCTK_REAL(0.0);CCTK_REALVEC momy1r = CCTK_REAL(0.0);CCTK_REALVEC momz1r = CCTK_REAL(0.0);CCTK_REALVEC etot1r = CCTK_REAL(1.0);
CCTK_REALVEC velx1r = CCTK_REAL(0.0);CCTK_REALVEC vely1r = CCTK_REAL(0.0);CCTK_REALVEC velz1r = CCTK_REAL(0.0);CCTK_REALVEC eint1r = CCTK_REAL(1.0);
CCTK_REALVEC momx1 = ifthen(x1 <= CCTK_REAL(0.0), momx1l, momx1r);CCTK_REALVEC momy1 = ifthen(x1 <= CCTK_REAL(0.0), momy1l, momy1r);CCTK_REALVEC momz1 = ifthen(x1 <= CCTK_REAL(0.0), momz1l, momz1r);CCTK_REALVEC etot1 = ifthen(x1 <= CCTK_REAL(0.0), etot1l, etot1r);
CCTK_REALVEC velx1 = ifthen(x1 <= CCTK_REAL(0.0), velx1l, velx1r);CCTK_REALVEC vely1 = ifthen(x1 <= CCTK_REAL(0.0), vely1l, vely1r);CCTK_REALVEC velz1 = ifthen(x1 <= CCTK_REAL(0.0), velz1l, velz1r);CCTK_REALVEC eint1 = ifthen(x1 <= CCTK_REAL(0.0), eint1l, eint1r);
momx1.storeu_partial(momx[ind], i, imin[0], imax[0]);momy1.storeu_partial(momy[ind], i, imin[0], imax[0]);momz1.storeu_partial(momz[ind], i, imin[0], imax[0]);etot1.storeu_partial(etot[ind], i, imin[0], imax[0]);
velx1.storeu_partial(velx[ind], i, imin[0], imax[0]);vely1.storeu_partial(vely[ind], i, imin[0], imax[0]);velz1.storeu_partial(velz[ind], i, imin[0], imax[0]);eint1.storeu_partial(eint[ind], i, imin[0], imax[0]);
CCTK_REALVEC x1 = x0[0] + (i + CCTK_REAL(0.5) + viota()) * dx[0];CCTK_REALVEC y1 = x0[1] + (j + CCTK_REAL(0.5)) * dx[1];CCTK_REALVEC z1 = x0[2] + (k + CCTK_REAL(0.5)) * dx[2];CCTK_REALVEC r2 = pow2(x1) + pow2(y1) + pow2(z1);
CCTK_REALVEC x1 = x0[0] + (lbnd[0] + i) * dx[0] + viota() * dx[0];CCTK_REALVEC y1 = x0[1] + (lbnd[1] + j) * dx[1];CCTK_REALVEC z1 = x0[2] + (lbnd[2] + k) * dx[2];CCTK_REALVEC r2 = pow2(x1) + (pow2(y1) + pow2(z1));
CCTK_REALVEC momx1i = CCTK_REAL(0.0);CCTK_REALVEC momy1i = CCTK_REAL(0.0);CCTK_REALVEC momz1i = CCTK_REAL(0.0);CCTK_REALVEC etot1i = CCTK_REAL(2.0);
CCTK_REALVEC velx1i = CCTK_REAL(0.0);CCTK_REALVEC vely1i = CCTK_REAL(0.0);CCTK_REALVEC velz1i = CCTK_REAL(0.0);CCTK_REALVEC eint1i = CCTK_REAL(2.0);
CCTK_REALVEC momx1o = CCTK_REAL(0.0);CCTK_REALVEC momy1o = CCTK_REAL(0.0);CCTK_REALVEC momz1o = CCTK_REAL(0.0);CCTK_REALVEC etot1o = CCTK_REAL(1.0);
CCTK_REALVEC velx1o = CCTK_REAL(0.0);CCTK_REALVEC vely1o = CCTK_REAL(0.0);CCTK_REALVEC velz1o = CCTK_REAL(0.0);CCTK_REALVEC eint1o = CCTK_REAL(1.0);
CCTK_REALVEC momx1 = ifthen(r2 <= sr2, momx1i, momx1o);CCTK_REALVEC momy1 = ifthen(r2 <= sr2, momy1i, momy1o);CCTK_REALVEC momz1 = ifthen(r2 <= sr2, momz1i, momz1o);CCTK_REALVEC etot1 = ifthen(r2 <= sr2, etot1i, etot1o);
CCTK_REALVEC velx1 = ifthen(r2 <= sr2, velx1i, velx1o);CCTK_REALVEC vely1 = ifthen(r2 <= sr2, vely1i, vely1o);CCTK_REALVEC velz1 = ifthen(r2 <= sr2, velz1i, velz1o);CCTK_REALVEC eint1 = ifthen(r2 <= sr2, eint1i, eint1o);
momx1.storeu_partial(momx[ind], i, imin[0], imax[0]);momy1.storeu_partial(momy[ind], i, imin[0], imax[0]);momz1.storeu_partial(momz[ind], i, imin[0], imax[0]);etot1.storeu_partial(etot[ind], i, imin[0], imax[0]);
velx1.storeu_partial(velx[ind], i, imin[0], imax[0]);vely1.storeu_partial(vely[ind], i, imin[0], imax[0]);velz1.storeu_partial(velz[ind], i, imin[0], imax[0]);eint1.storeu_partial(eint[ind], i, imin[0], imax[0]);
}// Set an ignorable pressurefor (int k = imin[2]; k < imax[2]; ++k) {for (int j = imin[1]; j < imax[1]; ++j) {for (int i = imin[0]; i < imax[0]; i += vsize) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC press1 = CCTK_REAL(-1.0);press1.storeu_partial(press[ind], i, imin[0], imax[0]);}}
#include "defs.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>#include <array>#include <cmath>namespace Hydro {using namespace std;////////////////////////////////////////////////////////////////////////////////extern "C" void Hydro_Prim2Con(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_Prim2Con;DECLARE_CCTK_PARAMETERS;const Loop::vect<int, dim> lsh{{cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]}};array<CCTK_INT, dim> tmin_, tmax_;GetTileExtent(cctkGH, tmin_.data(), tmax_.data());const Loop::vect<int, dim> tmin(tmin_);const Loop::vect<int, dim> tmax(tmax_);const Loop::vect<int, dim> imin = tmin;const Loop::vect<int, dim> imax = tmax;const Loop::vect<int, dim> ash{{cctk_ash[0], cctk_ash[1], cctk_ash[2]}};constexpr ptrdiff_t di = 1;const ptrdiff_t dj = di * cctk_ash[0];const ptrdiff_t dk = dj * cctk_ash[1];for (int k = imin[2]; k < imax[2]; ++k) {for (int j = imin[1]; j < imax[1]; ++j) {for (int i = imin[0]; i < imax[0]; i += vsize) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC rho1 = vloadu(rho[ind]);CCTK_REALVEC velx1 = vloadu(velx[ind]);CCTK_REALVEC vely1 = vloadu(vely[ind]);CCTK_REALVEC velz1 = vloadu(velz[ind]);CCTK_REALVEC eint1 = vloadu(eint[ind]);CCTK_REALVEC ekin1 =CCTK_REAL(0.5) * rho1 * (pow2(velx1) + pow2(vely1) + pow2(velz1));CCTK_REALVEC dens1 = rho1;CCTK_REALVEC momx1 = rho1 * velx1;CCTK_REALVEC momy1 = rho1 * vely1;CCTK_REALVEC momz1 = rho1 * velz1;CCTK_REALVEC etot1 = ekin1 + eint1;dens1.storeu_partial(dens[ind], i, imin[0], imax[0]);momx1.storeu_partial(momx[ind], i, imin[0], imax[0]);momy1.storeu_partial(momy[ind], i, imin[0], imax[0]);momz1.storeu_partial(momz[ind], i, imin[0], imax[0]);etot1.storeu_partial(etot[ind], i, imin[0], imax[0]);}}}}} // namespace Hydro
ptrdiff_t indx = i + djx * j + dkx * k;ptrdiff_t indy = i + djy * j + dky * k;ptrdiff_t indz = i + djz * j + dkz * k;
ptrdiff_t indx = i + djx * j + dkx * k - offx;ptrdiff_t indy = i + djy * j + dky * k - offy;ptrdiff_t indz = i + djz * j + dkz * k - offz;