Y7FSCHC3E3KIGVSB2KQRIPKBD2HL6JPUNBBSVZQAJKPXJKKRVNQQC Cactus Code Thorn HydroAuthor(s) : Erik Schnetter <schnetter@gmail.com>Maintainer(s): Erik Schnetter <schnetter@gmail.com>Licence : LGPL--------------------------------------------------------------------------1. PurposeSolve the non-relativistic hydrodynamics equations
# Configuration definitions for thorn HydroREQUIRES Vectors
# Interface definition for thorn HydroIMPLEMENTS: HydroUSES INCLUDE HEADER: loop.hxxUSES INCLUDE HEADER: vectors.hvoid FUNCTION GetTileExtent(CCTK_POINTER_TO_CONST IN cctkGH,CCTK_INT ARRAY OUT tile_min,CCTK_INT ARRAY OUT tile_max)REQUIRES FUNCTION GetTileExtentCCTK_REAL conserved TYPE=gf TIMELEVELS=2 TAGS='index={1 1 1}'{rhomomx 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"'{pressvelx vely velzeint} "Primitive variables"CCTK_REAL flux_x TYPE=gf TAGS='index={0 1 1} checkpoint="no" restrict="no" nghostzones={0 0 0}'{fxrhofxmomx fxmomy fxmomzfxetot} "Fluxes in x direction"CCTK_REAL flux_y TYPE=gf TAGS='index={1 0 1} checkpoint="no" restrict="no" nghostzones={0 0 0}'{fyrhofymomx fymomy fymomzfyetot} "Fluxes in y direction"CCTK_REAL flux_z TYPE=gf TAGS='index={1 1 0} checkpoint="no" restrict="no" nghostzones={0 0 0}'{fzrhofzmomx fzmomy fzmomzfzetot} "Fluxes in z direction"
# Parameter definitions for thorn HydroKEYWORD setup "Initial setup" STEERABLE=never{"equilibrium" :: """sound wave" :: """shock tube" :: """spherical shock" :: ""} "equilibrium"CCTK_REAL amplitude "Wave amplitude" STEERABLE=never{0.0:* :: ""} 1.0e-3CCTK_REAL shock_radius "Shock radius" STEERABLE=never{0.0:* :: ""} 0.1CCTK_REAL gamma "EOS parameter" STEERABLE=always{0.0:* :: ""} 1.6666666666666667CCTK_INT output_every "Calculate output quantities every that many iterations" STEERABLE=recover{0:* :: ""} 1
# Schedule definitions for thorn HydroSTORAGE: conserved[2]STORAGE: conserved_rhsSTORAGE: primitiveSTORAGE: flux_x flux_y flux_z# Initial conditionsSCHEDULE Hydro_Initialize AT initial{LANG: CWRITES: conserved(interior)INVALIDATES: primitive(everywhere)INVALIDATES: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)SYNC: conserved} "Set up hydro initial conditions"SCHEDULE Hydro_EstimateError AT postinitial{LANG: CREADS: conserved(everywhere)WRITES: CarpetX::regrid_error(interior)} "Estimate local error for regridding initial conditions"# RegriddingSCHEDULE Hydro_Boundaries AT postregrid{LANG: CSYNC: conserved} "Apply boundary conditions after regridding"# Dependent quantities# TODO: Use a more accurate conditionif (output_every > 0) {SCHEDULE GROUP Hydro_OutputGroup AT poststep{} "Output HydroToy quantities"SCHEDULE Hydro_Con2prim IN Hydro_OutputGroup{LANG: CREADS: conserved(everywhere)WRITES: primitive(everywhere)INVALIDATES: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)} "Calculate pressure"SCHEDULE Hydro_Fluxes IN Hydro_OutputGroup AFTER Hydro_Con2prim{LANG: CREADS: conserved(everywhere)READS: primitive(everywhere)WRITES: flux_x(interior) flux_y(interior) flux_z(interior)} "Calculate the hydro fluxes"}SCHEDULE Hydro_EstimateError AT poststep{LANG: CREADS: conserved(everywhere)WRITES: CarpetX::regrid_error(interior)} "Estimate local error for regridding during evolution"# Time steppingSCHEDULE Hydro_Con2prim AT evol{LANG: CREADS: conserved(everywhere)WRITES: primitive(everywhere)INVALIDATES: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)} "Calculate pressure"SCHEDULE Hydro_Fluxes AT evol AFTER Hydro_Con2prim{LANG: CREADS: conserved(everywhere)READS: primitive(everywhere)WRITES: flux_x(interior) flux_y(interior) flux_z(interior)SYNC: flux_x flux_y flux_z # restrict} "Calculate the hydro fluxes"SCHEDULE Hydro_RHS AT evol AFTER Hydro_Fluxes{LANG: CREADS: flux_x(interior) flux_y(interior) flux_z(interior)WRITES: conserved_rhs(interior)SYNC: conserved_rhs # sync and restrict} "Calculate RHS of the hydro equations"
#include "defs.hxx"#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>namespace Hydro {using namespace std;extern "C" void Hydro_Boundaries(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_Boundaries;DECLARE_CCTK_PARAMETERS;// do nothing}} // namespace Hydro
#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_Con2prim(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_Con2prim;DECLARE_CCTK_PARAMETERS;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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC rho1 = vloadu(rho[ind]);CCTK_REALVEC momx1 = vloadu(momx[ind]);CCTK_REALVEC momy1 = vloadu(momy[ind]);CCTK_REALVEC momz1 = vloadu(momz[ind]);CCTK_REALVEC etot1 = vloadu(etot[ind]);CCTK_REALVEC rho_inv = CCTK_REAL(1.0) / rho1;CCTK_REALVEC ekin = CCTK_REAL(0.5) * rho_inv *sqrt(pow2(momx1) + pow2(momy1) + pow2(momz1));CCTK_REALVEC eint1 = etot1 - ekin;// Equation of state: p = (gamma - 1) eCCTK_REALVEC press1 = CCTK_REAL(gamma - 1) * eint1;// vel^j = delta^j_i mom_i / rhoCCTK_REALVEC velx1 = rho_inv * momx1;CCTK_REALVEC vely1 = rho_inv * momy1;CCTK_REALVEC velz1 = rho_inv * momz1;press1.storeu_partial(press[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]);}}}}} // namespace Hydro
#ifndef DEFS_HH#define DEFS_HH#include <vectors.h>#include <array>#include <cmath>namespace Hydro {using namespace std;constexpr int dim = 3;template <typename T> T pow2(T x) { return x * x; }template <typename T> inline T fmax3(T x0, T x1, T x2) {T x01 = fmax(x0, x1);T x012 = fmax(x2, x01);return x012;}template <typename T> inline T fmax5(T x0, T x1, T x2, T x3, T x4) {T x01 = fmax(x0, x1);T x23 = fmax(x2, x3);T x014 = fmax(x4, x01);T x01234 = fmax(x23, x014);return x01234;}typedef vectype<CCTK_REAL> CCTK_REALVEC;constexpr int VS = vecprops<CCTK_REAL>::size();inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REALVECvloadu(const CCTK_REAL &restrict x) {return CCTK_REALVEC::loadu(x);}inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REALVEC viota() {array<CCTK_REAL, VS> iota_;for (int i = 0; i < VS; ++i)iota_[i] = i;return CCTK_REALVEC::loadu(iota_[0]);}inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REALVEC vsin(CCTK_REALVEC x) {return ksin(x);}} // namespace Hydro#endif // #ifndef DEFS_HH
#include "defs.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>#include <array>namespace Hydro {using namespace std;////////////////////////////////////////////////////////////////////////////////extern "C" void Hydro_EstimateError(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_EstimateError;DECLARE_CCTK_PARAMETERS;const Loop::vect<int, dim> lsh{{cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]}};const Loop::vect<int, dim> nghostzones{{cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[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 = max(tmin, nghostzones);const Loop::vect<int, dim> imax = min(tmax, lsh - nghostzones);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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;auto calcerr{[&](const CCTK_REAL *restrict var) {CCTK_REALVEC varxx = vloadu(var[ind - 1]) -CCTK_REAL(2.0) * vloadu(var[ind]) +vloadu(var[ind + 1]);CCTK_REALVEC varyy = vloadu(var[ind - dj]) -CCTK_REAL(2.0) * vloadu(var[ind]) +vloadu(var[ind + dj]);CCTK_REALVEC varzz = vloadu(var[ind - dk]) -CCTK_REAL(2.0) * vloadu(var[ind]) +vloadu(var[ind + dk]);return fmax3(varxx, varyy, varzz);}};CCTK_REALVEC regrid_error1 =fmax5(calcerr(rho), calcerr(momx), calcerr(momy), calcerr(momz),calcerr(etot));regrid_error1.storeu_partial(regrid_error[ind], i, imin[0], imax[0]);}}}}} // namespace Hydro
#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;////////////////////////////////////////////////////////////////////////////////namespace {// LLF (local Lax-Friedrichs) Riemann solvertemplate <typename T>inline CCTK_ATTRIBUTE_ALWAYS_INLINE T llf(const array<T, 2> &var,const array<T, 2> &flux) {return T(0.5) * ((flux[0] + flux[1]) - (var[1] - var[0]));}} // namespace////////////////////////////////////////////////////////////////////////////////extern "C" void Hydro_Fluxes(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_Fluxes;DECLARE_CCTK_PARAMETERS;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 CCTK_REAL dAx = dy * dz;const CCTK_REAL dAy = dx * dz;const CCTK_REAL dAz = dx * dy;const Loop::vect<int, dim> zero{{0, 0, 0}};const Loop::vect<int, dim> nx{{1, 0, 0}};const Loop::vect<int, dim> ny{{0, 1, 0}};const Loop::vect<int, dim> nz{{0, 0, 1}};const Loop::vect<int, dim> lsh{{cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]}};const Loop::vect<int, dim> nghostzones{{cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[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> ash{{cctk_ash[0], cctk_ash[1], cctk_ash[2]}};constexpr ptrdiff_t di = 1;const ptrdiff_t dj = di * ash[0];const ptrdiff_t dk = dj * ash[1];const auto calcflux{[&](const int dir, const int ddir, const CCTK_REAL dAdir,const Loop::vect<int, dim> &ndir,const CCTK_REAL *restrict const veldir,CCTK_REAL *restrict const fdirrho,CCTK_REAL *restrict const fdirmomx,CCTK_REAL *restrict const fdirmomy,CCTK_REAL *restrict const fdirmomz,CCTK_REAL *restrict constfdiretot) CCTK_ATTRIBUTE_ALWAYS_INLINE {const Loop::vect<int, dim> imindir = max(tmin, nghostzones);const Loop::vect<int, dim> imaxdir =min(tmax + (tmax >= lsh).ifelse(ndir, zero), lsh + ndir - nghostzones);// fluxes: face-centred without ghostsconst Loop::vect<int, dim> ashdir = ash + ndir - 2 * nghostzones;constexpr ptrdiff_t didir = 1;const ptrdiff_t djdir = didir * ashdir[0];const ptrdiff_t dkdir = djdir * ashdir[1];for (int k = imindir[2]; k < imaxdir[2]; ++k) {for (int j = imindir[1]; j < imaxdir[1]; ++j) {for (int i = imindir[0]; i < imaxdir[0]; i += VS) {ptrdiff_t ind = i + dj * j + dk * k;ptrdiff_t inddir = i + djdir * j + dkdir * k;// Read conserved and primitive variablesarray<CCTK_REALVEC, 2> rho2, momx2, momy2, momz2, etot2, press2,veldir2;for (int n = 0; n < 2; ++n) {rho2[n] = vloadu(rho[ind + (n - 1) * ddir]);momx2[n] = vloadu(momx[ind + (n - 1) * ddir]);momy2[n] = vloadu(momy[ind + (n - 1) * ddir]);momz2[n] = vloadu(momz[ind + (n - 1) * ddir]);etot2[n] = vloadu(etot[ind + (n - 1) * ddir]);press2[n] = vloadu(press[ind + (n - 1) * ddir]);veldir2[n] = vloadu(veldir[ind + (n - 1) * ddir]);}// Calculate centred fluxesarray<CCTK_REALVEC, 2> fdirrho2, fdirmomx2, fdirmomy2, fdirmomz2,fdiretot2;for (int n = 0; n < 2; ++n) {fdirrho2[n] = rho2[n] * veldir2[n];fdirmomx2[n] = momx2[n] * veldir2[n] + (dir == 0 ? press2[n] : 0);fdirmomy2[n] = momy2[n] * veldir2[n] + (dir == 1 ? press2[n] : 0);fdirmomz2[n] = momz2[n] * veldir2[n] + (dir == 2 ? press2[n] : 0);fdiretot2[n] = (etot2[n] + press2[n]) * veldir2[n];}// Calculate face fluxesCCTK_REALVEC fdirrho1 = dAdir * llf(rho2, fdirrho2);CCTK_REALVEC fdirmomx1 = dAdir * llf(momx2, fdirmomx2);CCTK_REALVEC fdirmomy1 = dAdir * llf(momy2, fdirmomy2);CCTK_REALVEC fdirmomz1 = dAdir * llf(momz2, fdirmomz2);CCTK_REALVEC fdiretot1 = dAdir * llf(etot2, fdiretot2);// Store fluxesfdirrho1.storeu_partial(fdirrho[inddir], i, imindir[0], imaxdir[0]);fdirmomx1.storeu_partial(fdirmomx[inddir], i, imindir[0], imaxdir[0]);fdirmomy1.storeu_partial(fdirmomy[inddir], i, imindir[0], imaxdir[0]);fdirmomz1.storeu_partial(fdirmomz[inddir], i, imindir[0], imaxdir[0]);fdiretot1.storeu_partial(fdiretot[inddir], i, imindir[0], imaxdir[0]);}}}}};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);}} // namespace Hydro
#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_Initialize(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_Initialize;DECLARE_CCTK_PARAMETERS;const Loop::vect<CCTK_REAL, dim> x0{{CCTK_ORIGIN_SPACE(0), CCTK_ORIGIN_SPACE(1), CCTK_ORIGIN_SPACE(2)}};const Loop::vect<CCTK_REAL, dim> dx{{CCTK_DELTA_SPACE(0), CCTK_DELTA_SPACE(1), CCTK_DELTA_SPACE(2)}};const Loop::vect<int, dim> lsh{{cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]}};const Loop::vect<int, dim> nghostzones{{cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[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 = max(tmin, nghostzones);const Loop::vect<int, dim> imax = min(tmax, lsh - nghostzones);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];if (CCTK_EQUALS(setup, "equilibrium")) {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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC rho1 = CCTK_REAL(1.0);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);rho1.storeu_partial(rho[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]);}}}} else if (CCTK_EQUALS(setup, "sound wave")) {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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC x1 = x0[0] + (i + CCTK_REAL(0.5) + viota()) * dx[0];CCTK_REALVEC rho1 = CCTK_REAL(1.0);CCTK_REALVEC momx1 = CCTK_REAL(0.0) + CCTK_REAL(amplitude) *vsin(CCTK_REAL(M_PI) * x1);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 hererho1.storeu_partial(rho[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]);}}}} else if (CCTK_EQUALS(setup, "shock tube")) {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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;CCTK_REALVEC x1 = x0[0] + (i + CCTK_REAL(0.5) + viota()) * dx[0];// left stateCCTK_REALVEC rho1l = CCTK_REAL(2.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);// right stateCCTK_REALVEC rho1r = CCTK_REAL(1.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);// chooseCCTK_REALVEC rho1 = ifthen(x1 <= CCTK_REAL(0.0), rho1l, rho1r);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);rho1.storeu_partial(rho[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]);}}}} else if (CCTK_EQUALS(setup, "spherical shock")) {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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;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_REAL sr2 = pow2(shock_radius);// inner stateCCTK_REALVEC rho1i = CCTK_REAL(2.0);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);// outer stateCCTK_REALVEC rho1o = CCTK_REAL(1.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);// chooseCCTK_REALVEC rho1 = ifthen(r2 <= sr2, rho1i, rho1o);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);rho1.storeu_partial(rho[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]);}}}} else {assert(0);}}} // namespace Hydro
# Main make.code.defn file for thorn Hydro# Source files in this directorySRCS = boundaries.cxx \con2prim.cxx \estimateerror.cxx \fluxes.cxx \initialize.cxx \rhs.cxx# Subdirectories containing source filesSUBDIRS =
#include "defs.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>namespace Hydro {using namespace std;////////////////////////////////////////////////////////////////////////////////extern "C" void Hydro_RHS(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Hydro_RHS;DECLARE_CCTK_PARAMETERS;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 CCTK_REAL dV1 = 1 / (dx * dy * dz);const Loop::vect<int, dim> zero{{0, 0, 0}};const Loop::vect<int, dim> nx{{1, 0, 0}};const Loop::vect<int, dim> ny{{0, 1, 0}};const Loop::vect<int, dim> nz{{0, 0, 1}};const Loop::vect<int, dim> lsh{{cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]}};const Loop::vect<int, dim> nghostzones{{cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[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 = max(tmin, nghostzones);const Loop::vect<int, dim> imax = min(tmax, lsh - nghostzones);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 * ash[0];const ptrdiff_t dk = dj * ash[1];// x// fluxes: face-centred without ghostsconst Loop::vect<int, dim> ashx = ash + nx - 2 * nghostzones;constexpr ptrdiff_t dix = 1;const ptrdiff_t djx = dix * ashx[0];const ptrdiff_t dkx = djx * ashx[1];// y// fluxes: face-centred without ghostsconst Loop::vect<int, dim> ashy = ash + ny - 2 * nghostzones;constexpr ptrdiff_t diy = 1;const ptrdiff_t djy = diy * ashy[0];const ptrdiff_t dky = djy * ashy[1];// z// fluxes: face-centred without ghostsconst Loop::vect<int, dim> ashz = ash + nz - 2 * nghostzones;constexpr ptrdiff_t diz = 1;const ptrdiff_t djz = diz * ashz[0];const ptrdiff_t dkz = djz * ashz[1];const auto calcrhs{[&](const CCTK_REAL *restrict const fxvar,const CCTK_REAL *restrict const fyvar,const CCTK_REAL *restrict const fzvar,CCTK_REAL *restrict const dtvar) CCTK_ATTRIBUTE_ALWAYS_INLINE {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 += VS) {ptrdiff_t ind = i + dj * j + dk * k;ptrdiff_t indx = i + djx * j + dkx * k;ptrdiff_t indy = i + djy * j + dky * k;ptrdiff_t indz = i + djz * j + dkz * k;CCTK_REALVEC dtvar1 =dV1 * ((vloadu(fxvar[indx + dix]) - vloadu(fxvar[indx])) +(vloadu(fyvar[indy + djy]) - vloadu(fyvar[indy])) +(vloadu(fzvar[indz + dkz]) - vloadu(fzvar[indz])));dtvar1.store_partial(dtvar[ind], i, imin[0], imax[0]);}}}}};calcrhs(fxrho, fyrho, fzrho, dtrho);calcrhs(fxmomx, fymomx, fzmomx, dtmomx);calcrhs(fxmomy, fymomy, fzmomy, dtmomy);calcrhs(fxmomz, fymomz, fzmomz, dtmomz);calcrhs(fxetot, fyetot, fzetot, dtetot);}} // namespace Hydro