diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index f84e9e31..137f019c 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -304,7 +304,7 @@ class ConToPrim { KOKKOS_INLINE_FUNCTION ConToPrimStatus solve(const VarAccessor &v, const CellGeom &g, const singularity::EOS &eos, const Real x1, const Real x2, const Real x3) const { - int num_nans = std::isnan(v(crho)) + std::isnan(v(cmom_lo)) + std::isnan(ceng); + int num_nans = std::isnan(v(crho)) + std::isnan(v(cmom_lo)) + std::isnan(v(ceng)); if (num_nans > 0) return ConToPrimStatus::failure; const Real igdet = 1.0/g.gdet; diff --git a/src/fluid/riemann.hpp b/src/fluid/riemann.hpp index 05bb2f1d..12be3da6 100644 --- a/src/fluid/riemann.hpp +++ b/src/fluid/riemann.hpp @@ -159,7 +159,6 @@ class FluxState { bcovm += g.gcov[m+1][n] * (Bcon[n-1]/W + b0*(vcon[n-1]-g.beta[n-1]/g.alpha)); vcovm += g.gcov[m+1][n] * vcon[n-1]; } - U[cmom_lo+m] = (rhohWsq+bsqWsq)*vcovm - b0*bcovm; F[cmom_lo+m] = U[cmom_lo+m]*vtil + (P + 0.5*bsq)*Delta(dir,m) - bcovm*Bcon[dir]/W; } diff --git a/src/geometry/inchworm.cpp b/src/geometry/inchworm.cpp index d40363e0..49ec3517 100644 --- a/src/geometry/inchworm.cpp +++ b/src/geometry/inchworm.cpp @@ -36,8 +36,10 @@ void Initialize(ParameterInput *pin, StateDescriptor *geometry) { Params ¶ms = geometry->AllParams(); Real a = pin->GetOrAddReal("geometry", "a", 0.3); - Real k = pin->GetOrAddReal("geometry", "k", M_PI/2.); - params.Add("a", a); + Real k = pin->GetOrAddReal("geometry", "k", 2.*M_PI); + Real kmult = pin->GetOrAddReal("geometry", "kmult", 1); + k *= kmult; + params.Add("a", a); params.Add("k", k); } diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index 72ef749e..63dbc00f 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -43,25 +43,20 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const int ndim = pmb->pmy_mesh->ndim; PackIndexMap imap; - auto v = rc->PackVariables({"p.density", - "p.velocity", - "p.energy", - fluid_prim::bfield, - fluid_prim::ye, - "pressure", - "temperature", - "gamma1", - "cs"}, - imap); - - const int irho = imap["p.density"].first; - const int ivlo = imap["p.velocity"].first; - const int ivhi = imap["p.velocity"].second; - const int ieng = imap["p.energy"].first; + std::vector vars({fluid_prim::density, fluid_prim::velocity, + fluid_prim::energy, fluid_prim::bfield, + fluid_prim::pressure, fluid_prim::temperature, + fluid_prim::ye}); + auto v = rc->PackVariables(vars, imap); + + const int irho = imap[fluid_prim::density].first; + const int ivlo = imap[fluid_prim::velocity].first; + const int ivhi = imap[fluid_prim::velocity].second; + const int ieng = imap[fluid_prim::energy].first; const int ib_lo = imap[fluid_prim::bfield].first; const int ib_hi = imap[fluid_prim::bfield].second; - const int iprs = imap["pressure"].first; - const int itmp = imap["temperature"].first; + const int iprs = imap[fluid_prim::pressure].first; + const int itmp = imap[fluid_prim::temperature].first; const int iye = imap[fluid_prim::ye].second; const int nv = ivhi - ivlo + 1; @@ -226,8 +221,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real eos_lambda[2]; if (iye > 0) { - v(iye, k, j, i) = 0.5; - eos_lambda[0] = v(iye, k, j, i); + v(iye, k, j, i) = 0.5; + eos_lambda[0] = v(iye, k, j, i); } // This line causes NaNs and I don't know why //v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy(rho, v(ieng, k, j, i)/rho); @@ -257,33 +252,47 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } Real Gamma = sqrt(1. + vsq); - if (is_snake ||is_inchworm || is_boosted_minkowski) { + if (is_snake || is_inchworm || is_boosted_minkowski) { PARTHENON_REQUIRE(ivhi == 3, "Only works for 3D velocity!"); // Transform velocity Real gcov[NDFULL][NDFULL] = {0}; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); Real ucon[NDFULL] = {Gamma, // alpha = 1 in Minkowski - v(ivlo, k, j, i), // beta^i = 0 in Minkowski + v(ivlo, k, j, i), // beta^i = 0 in Minkowski v(ivlo+1, k, j, i), v(ivlo+2, k, j, i)}; + Real Bdotu = 0.0; + for(int d = ib_lo; d <= ib_hi; d++){ + Bdotu += v(d, k, j, i) * ucon[d-ib_lo+1]; + } + Real bcon[NDFULL] = {Bdotu, 0.0,0.0,0.0}; + for(int d =ib_lo; d <= ib_hi; d++){ + bcon[d-ib_lo+1] = (v(d, k, j, i) + alpha*bcon[0]*ucon[d-ib_lo+1]) / Gamma; + } Real J[NDFULL][NDFULL] = {0}; if (is_snake) { J[0][0] = 1/alpha; - J[2][0] = -betay/alpha; + J[2][0] = -betay/alpha; J[2][1] = a_snake*k_snake*cos(k_snake*x); - J[1][1] = J[2][2] = J[3][3] = 1; - } else if (is_boosted_minkowski) { - J[0][0] = J[1][1] = J[2][2] = J[3][3] = 1; - J[1][0] = -betax; - J[2][0] = -betay; - J[3][0] = -betaz; + J[1][1] = J[2][2] = J[3][3] = 1; + } else if (is_boosted_minkowski) { + J[0][0] = J[1][1] = J[2][2] = J[3][3] = 1; + J[1][0] = -betax; + J[2][0] = -betay; + J[3][0] = -betaz; } else if (is_inchworm) { - PARTHENON_FAIL("This geometry isn't supported with a J!"); + J[0][0] = J[2][2] = J[3][3] = 1; + J[1][1] = 1 + a_snake*k_snake*cos(k_snake*x); } Real ucon_transformed[NDFULL] = {0, 0, 0, 0}; SPACETIMELOOP(mu) SPACETIMELOOP(nu){ ucon_transformed[mu] += J[mu][nu]*ucon[nu]; } + Real bcon_transformed[NDFULL] = {0, 0, 0, 0}; + SPACETIMELOOP(mu) SPACETIMELOOP(nu){ + bcon_transformed[mu] += J[mu][nu]*bcon[nu]; + } + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); Gamma = lapse * ucon_transformed[0]; Real shift[NDSPACE]; @@ -291,12 +300,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo, k, j, i) = ucon_transformed[1] + Gamma*shift[0]/lapse; v(ivlo+1, k, j, i) = ucon_transformed[2] + Gamma*shift[1]/lapse; v(ivlo+2, k, j, i) = ucon_transformed[3] + Gamma*shift[2]/lapse; - - // Enforce zero B fields for now - if (ib_hi >= 3) { - v(ib_lo, k, j, i) = 0.; - v(ib_lo + 1, k, j, i) = 0.; - v(ib_lo + 2, k, j, i) = 0.; + for(int d = ib_lo; d <= ib_hi; d++){ + v(d, k, j, i) = bcon_transformed[d-ib_lo+1]*Gamma - lapse*bcon_transformed[0]*ucon_transformed[d-ib_lo+1]; } } });