W6W3EC56H4PPREAK7DMX3DAHYKDAHSASGDP63DSU55TRCBBGJXBAC
momx momy momz
etot
} "Conserved variables"
CCTK_REAL conserved_rhs TYPE=gf TAGS='index={1 1 1} fluxes="HydroToyCarpetX::flux_x HydroToyCarpetX::flux_y HydroToyCarpetX::flux_z"'
{
dtrho
dtmomx dtmomy dtmomz
dtetot
} "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"'
{
dtdens
dtmomx dtmomy dtmomz
dtetot
} "Time derivatives of conserved variables"
CCTK_REAL conserved TYPE=gf TAGS='index={1 1 1} rhs="Hydro::conserved_rhs"'
{
dens
momx momy momz
etot
} "Conserved variables"
ActiveThorns = "
CarpetX
Formaline
IOUtil
Hydro
ODESolvers
SystemTopology
TimerReport
"
$nlevels = 3
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = no
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::max_num_levels = $nlevels
CarpetX::regrid_every = 16
CarpetX::regrid_error_threshold = 0.01
CarpetX::prolongation_type = "conservative"
CarpetX::do_reflux = no
CarpetX::do_restrict = no
CarpetX::restrict_during_sync = yes
CarpetX::dtfac = 0.25
Hydro::setup = "spherical shock"
Hydro::output_every = 1
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32
CarpetX::out_tsv = no
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50
SCHEDULE Hydro_EstimateError AT postinitial
SCHEDULE Hydro_Prim2Con AT postinitial
{
LANG: C
READS: 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 pressure
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 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;