From 0c157798cca9d49faa0b81d9359c118387e7f76c Mon Sep 17 00:00:00 2001 From: Janiris Altagracia Rodriguez Bueno - 358952 Date: Wed, 8 Dec 2021 16:01:50 -0700 Subject: [PATCH 1/9] fixing alfven and inchworm --- src/geometry/inchworm.cpp | 6 ++++-- src/pgen/linear_modes.cpp | 25 ++++++++++++++++++------- 2 files changed, 22 insertions(+), 9 deletions(-) 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 8a33c785..fa4852b7 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -263,6 +263,14 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Gamma*v(ivlo, k, j, i), Gamma*v(ivlo+1, k, j, i), Gamma*v(ivlo+2, k, j, i)}; + Real Bdotu = 0.0; + for(int d = 0; d < 3; d++){ + Bdotu += v(ib_lo+d, k, j, i) * ucon[d+1]; + } + Real bcon[NDFULL] = {Bdotu, 0.0,0.0,0.0}; + for(int d =1; d < 4; d++){ + bcon[d] = (v(ib_lo+d-1, k, j, i) + alpha*bcon[0]*ucon[d]) / Gamma; + } Real J[NDFULL][NDFULL] = {0}; if (is_snake) { J[0][0] = 1/alpha; @@ -275,12 +283,18 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { 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]; @@ -289,12 +303,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { 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 = 0; d < 3; d++){ + v(ib_lo+d, k, j, i) = bcon_transformed[d+1]*Gamma - lapse*bcon_transformed[0]*ucon_transformed[d+1]; + } } }); From ab1981060983a6452f20c644ece4cb200518f0b1 Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Thu, 9 Dec 2021 16:11:08 -0700 Subject: [PATCH 2/9] deal with segfault when there is no B field --- src/pgen/linear_modes.cpp | 40 +++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index fa4852b7..e4461a8f 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -263,28 +263,28 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Gamma*v(ivlo, k, j, i), Gamma*v(ivlo+1, k, j, i), Gamma*v(ivlo+2, k, j, i)}; - Real Bdotu = 0.0; - for(int d = 0; d < 3; d++){ - Bdotu += v(ib_lo+d, k, j, i) * ucon[d+1]; - } - Real bcon[NDFULL] = {Bdotu, 0.0,0.0,0.0}; - for(int d =1; d < 4; d++){ - bcon[d] = (v(ib_lo+d-1, k, j, i) + alpha*bcon[0]*ucon[d]) / Gamma; - } + 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) { J[0][0] = J[2][2] = J[3][3] = 1; - J[1][1] = 1 + a_snake*k_snake*cos(k_snake*x); + J[1][1] = 1 + a_snake*k_snake*cos(k_snake*x); } Real ucon_transformed[NDFULL] = {0, 0, 0, 0}; SPACETIMELOOP(mu) SPACETIMELOOP(nu){ @@ -294,7 +294,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { 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]; @@ -303,9 +303,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { 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; - for(int d = 0; d < 3; d++){ - v(ib_lo+d, k, j, i) = bcon_transformed[d+1]*Gamma - lapse*bcon_transformed[0]*ucon_transformed[d+1]; - } + 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]; + } } }); From c95ad09465ea5c0b4333112ef20204a84d811168 Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Tue, 8 Feb 2022 12:17:17 -0700 Subject: [PATCH 3/9] debugging --- external/singularity-eos | 2 +- src/fluid/con2prim_robust.hpp | 6 +++- src/fluid/fluid.cpp | 5 +++ src/geometry/snake.hpp | 6 +++- src/pgen/linear_modes.cpp | 42 ++++++++++++++------------ src/phoebus_utils/relativity_utils.hpp | 4 +-- 6 files changed, 41 insertions(+), 24 deletions(-) diff --git a/external/singularity-eos b/external/singularity-eos index 6a3fe45e..89480875 160000 --- a/external/singularity-eos +++ b/external/singularity-eos @@ -1 +1 @@ -Subproject commit 6a3fe45ec44c5cf7456a4e1fe2073b3d0a96058b +Subproject commit 89480875001f8d1c3fa0998a4f9a22887c0330b1 diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index e09394f1..2091d19c 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -300,7 +300,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; @@ -415,6 +415,10 @@ class ConToPrim { Real S[3]; Real bcons[3]; Real sig[3]; + /*if (fabs(vel[0]) > 0.0 || fabs(vel[1]) > 0.0) { + printf("Uninitialized? %g %g %g\n", vel[0], vel[1], vel[2]); + printf("%g %g %g %g %g %g %g %g %g %g\n", W, mu, x, rcon[0], rcon[1], rcon[2], bdotr, bu[0], bu[1], bu[2]); + }*/ prim2con::p2c(v(prho), vel, bu, v(peng), ye_prim, v(prs), v(gm1), g.gcov4, g.gcon, g.beta, g.lapse, g.gdet, v(crho), S, bcons, v(ceng), ye_cons, sig); diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 36e988be..70423aaf 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -330,6 +330,10 @@ TaskStatus PrimitiveToConservedRegion(MeshBlockData *rc, const IndexRange const Real vel[] = {v(b, pvel_lo, k, j, i), v(b, pvel_lo+1, k, j, i), v(b, pvel_hi, k, j, i)}; + if (fabs(vel[0]) > 0.0 || fabs(vel[1]) > 0.0) { + printf("Uninitialized? %d %d %d %d %g %g %g\n", b, k, j, i, vel[0], vel[1], vel[2]); + } + //printf("vel: %g %g %g\n", vel[0], vel[1], vel[2]); Real bcons[3]; Real bp[3] = {0.0, 0.0, 0.0}; if (pb_hi > 0) { @@ -407,6 +411,7 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, invert.NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + //printf("cell: %d %d %d %d\n", b, k, j, i); auto status = invert(geom, eos, coords, k, j, i); fail(k, j, i) = (status == con2prim_robust::ConToPrimStatus::success ? con2prim_robust::FailFlags::success diff --git a/src/geometry/snake.hpp b/src/geometry/snake.hpp index 4762d64c..1f34c1a8 100644 --- a/src/geometry/snake.hpp +++ b/src/geometry/snake.hpp @@ -63,6 +63,7 @@ class Snake { g[1][3] = g[3][1] = 0; g[2][2] = 1; g[2][3] = g[3][2] = 0; + g[3][3] = 1; } KOKKOS_INLINE_FUNCTION void SpacetimeMetricInverse(Real X0, Real X1, Real X2, Real X3, @@ -98,7 +99,7 @@ class Snake { gamma[2][0] = gamma[0][2] = 0; gamma[1][1] = d * d + 1; gamma[1][2] = gamma[2][1] = 0; - gamma[2][2] = 0; + gamma[2][2] = 1; } KOKKOS_INLINE_FUNCTION Real DetGamma(Real X0, Real X1, Real X2, Real X3) const { return 1.; } @@ -108,6 +109,8 @@ class Snake { KOKKOS_INLINE_FUNCTION void ConnectionCoefficient(Real X0, Real X1, Real X2, Real X3, Real Gamma[NDFULL][NDFULL][NDFULL]) const { + //Utils::SetConnectionCoeffByFD(*this, Gamma, X0, X1, X2, X3); + const Real a2 = a_ * a_; const Real k2 = k_ * k_; const Real k3 = k_ * k_ * k_; @@ -119,6 +122,7 @@ class Snake { KOKKOS_INLINE_FUNCTION void MetricDerivative(Real X0, Real X1, Real X2, Real X3, Real dg[NDFULL][NDFULL][NDFULL]) const { + //Utils::SetMetricGradientByFD(*this, 1.e-6, X0, X1, X2, X3, dg); const Real a2 = a_ * a_; const Real k2 = k_ * k_; const Real k3 = k_ * k_ * k_; diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index 135561ed..d22e5477 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -39,30 +39,32 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(is_minkowski || is_boosted_minkowski || is_snake || is_inchworm, "Problem \"linear_modes\" requires \"Minkowski\" geometry!"); + printf("minkowski/boosted_minkowski/snake/inchworm: %d/%d/%d/%d\n", + is_minkowski, is_boosted_minkowski, is_snake, is_inchworm); + auto &rc = pmb->meshblock_data.Get(); const int ndim = pmb->pmy_mesh->ndim; PackIndexMap imap; - auto v = rc->PackVariables({"p.density", - "p.velocity", - "p.energy", - fluid_prim::bfield, - "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}); + 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 nv = ivhi - ivlo + 1; + for (auto &field : vars) { + std::cout << field << ": " << imap[field].first << " " << imap[field].second << std::endl; + } + const Real gam = pin->GetReal("eos", "Gamma"); PARTHENON_REQUIRE_THROWS(nv == 3, "3 have 3 velocities"); @@ -124,6 +126,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_FAIL(msg); } } else if (physics == "mhd") { + if (ib_lo < 0 || ib_hi < 0) { + PARTHENON_THROW("physics = mhd but bfields are not present"); + } + printf("ib_lo = %d ib_hi = %d\n", ib_lo, ib_hi); B10 = 1.0; if (mode == "slow") { omega = complex(0., 2.41024185339); @@ -252,11 +258,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + ii, k, j, i) *= Gamma; } - 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}; - Real vcon[NDSPACE] = {v(ivlo, k, j, i)/Gamma, v(ivlo+1, k, j, i)/Gamma, v(ivlo+2, k, j, i)/Gamma}; 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 @@ -301,7 +306,6 @@ 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; - 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]; } diff --git a/src/phoebus_utils/relativity_utils.hpp b/src/phoebus_utils/relativity_utils.hpp index ba02006f..884cced0 100644 --- a/src/phoebus_utils/relativity_utils.hpp +++ b/src/phoebus_utils/relativity_utils.hpp @@ -32,7 +32,7 @@ GetLorentzFactor(const Real vcon[Geometry::NDSPACE], const Real gammacov[Geometr SPACELOOP2(ii, jj) { vsq += gammacov[ii][jj]*vcon[ii]*vcon[jj]; } - return sqrt(1. + vsq); + return std::sqrt(1. + vsq); } /* @@ -49,7 +49,7 @@ GetLorentzFactor(const Real vcon[Geometry::NDSPACE], const Real gcov[Geometry::N SPACELOOP2(ii, jj) { vsq += gcov[ii+1][jj+1]*vcon[ii]*vcon[jj]; } - return sqrt(1. + vsq); + return std::sqrt(1. + vsq); } /* From 05f9f33ea6b1aa4019b88d2dd0f506b10d3c0d2d Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Wed, 9 Feb 2022 11:15:48 -0700 Subject: [PATCH 4/9] fix gamma --- src/pgen/linear_modes.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index d22e5477..5d0b90d4 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -253,10 +253,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { SPACELOOP2(ii, jj) { vsq += v(ivlo + ii, k, j, i)*v(ivlo + jj, k, j, i); } - Real Gamma = 1./sqrt(1. - vsq); - SPACELOOP(ii) { - v(ivlo + ii, k, j, i) *= Gamma; - } + Real Gamma = sqrt(1. + vsq); if (is_snake || is_inchworm || is_boosted_minkowski) { PARTHENON_REQUIRE(ivhi == 3, "Only works for 3D velocity!"); From fca0f0c7c9d89c48a8a10b381b5fd44204aa1f73 Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Wed, 9 Feb 2022 11:16:30 -0700 Subject: [PATCH 5/9] clean up --- src/fluid/fluid.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 70423aaf..462e2bd4 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -330,10 +330,6 @@ TaskStatus PrimitiveToConservedRegion(MeshBlockData *rc, const IndexRange const Real vel[] = {v(b, pvel_lo, k, j, i), v(b, pvel_lo+1, k, j, i), v(b, pvel_hi, k, j, i)}; - if (fabs(vel[0]) > 0.0 || fabs(vel[1]) > 0.0) { - printf("Uninitialized? %d %d %d %d %g %g %g\n", b, k, j, i, vel[0], vel[1], vel[2]); - } - //printf("vel: %g %g %g\n", vel[0], vel[1], vel[2]); Real bcons[3]; Real bp[3] = {0.0, 0.0, 0.0}; if (pb_hi > 0) { From 513678e34ca85f6a6242bb81a125ff1fb4ca0ddc Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Wed, 9 Feb 2022 14:00:23 -0700 Subject: [PATCH 6/9] fixes --- src/fluid/con2prim_robust.hpp | 4 ---- src/fluid/riemann.hpp | 7 ++----- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/fluid/con2prim_robust.hpp b/src/fluid/con2prim_robust.hpp index 2091d19c..0783d3b5 100644 --- a/src/fluid/con2prim_robust.hpp +++ b/src/fluid/con2prim_robust.hpp @@ -415,10 +415,6 @@ class ConToPrim { Real S[3]; Real bcons[3]; Real sig[3]; - /*if (fabs(vel[0]) > 0.0 || fabs(vel[1]) > 0.0) { - printf("Uninitialized? %g %g %g\n", vel[0], vel[1], vel[2]); - printf("%g %g %g %g %g %g %g %g %g %g\n", W, mu, x, rcon[0], rcon[1], rcon[2], bdotr, bu[0], bu[1], bu[2]); - }*/ prim2con::p2c(v(prho), vel, bu, v(peng), ye_prim, v(prs), v(gm1), g.gcov4, g.gcon, g.beta, g.lapse, g.gdet, v(crho), S, bcons, v(ceng), ye_cons, sig); diff --git a/src/fluid/riemann.hpp b/src/fluid/riemann.hpp index 30ad9966..05bb2f1d 100644 --- a/src/fluid/riemann.hpp +++ b/src/fluid/riemann.hpp @@ -178,13 +178,10 @@ class FluxState { SPACETIMELOOP2(mu, nu) { ucov[mu] += g.gcov[mu][nu]*ucon[nu]; } - Real bcon[4] = {0., 0., 0., 0.}; + Real bcon[4] = {b0/g.alpha, 0., 0., 0.}; Real bcov0 = 0.; - SPACELOOP(ii) SPACETIMELOOP(mu) { - bcon[0] += Bcon[ii]*ucon[mu]*g.gcov[ii][mu]; - } SPACELOOP(ii) { - bcon[ii] = (Bcon[ii] + bcon[0]*ucon[ii])/ucon[0]; + bcon[ii+1] = (Bcon[ii] + b0*ucon[ii+1])/W; } SPACETIMELOOP(mu) { bcov0 += g.gcov[0][mu]*bcon[mu]; From 9db9e66947647597ed905239cf0d3afc88488c41 Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Fri, 11 Feb 2022 11:46:35 -0700 Subject: [PATCH 7/9] I think this works --- src/fluid/fluid.cpp | 3 +-- src/fluid/riemann.hpp | 1 - src/fluid/tmunu.hpp | 1 + src/geometry/snake.hpp | 3 --- src/phoebus_driver.cpp | 3 +++ 5 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 462e2bd4..c68cbb4c 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -203,10 +203,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { Metadata::Derived, Metadata::OneCopy}, five_vec); physics->AddField(diag::divf, mdiv); - std::vector seven_vec(1,5); Metadata mdiag = Metadata({Metadata::Cell, Metadata::Intensive, Metadata::Vector, Metadata::Derived, Metadata::OneCopy}, - seven_vec); + five_vec); physics->AddField(diag::src_terms, mdiag); #endif 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/fluid/tmunu.hpp b/src/fluid/tmunu.hpp index b37f3eb3..b8f37657 100644 --- a/src/fluid/tmunu.hpp +++ b/src/fluid/tmunu.hpp @@ -126,6 +126,7 @@ class StressEnergyTensorCon { u[l] = v_(l, std::forward(args)...) - u[0] * beta[l - 1]; b[l] = iW * (b_(l, std::forward(args)...) + alpha * b[0] * u[l]); } + //printf("b[] = %g %g %g %g\n", b[0], b[1], b[2], b[3]); bsq = (Bsq + alpha * alpha * b[0] * b[0]) * iW * iW; } diff --git a/src/geometry/snake.hpp b/src/geometry/snake.hpp index 1f34c1a8..9efdbf5a 100644 --- a/src/geometry/snake.hpp +++ b/src/geometry/snake.hpp @@ -109,8 +109,6 @@ class Snake { KOKKOS_INLINE_FUNCTION void ConnectionCoefficient(Real X0, Real X1, Real X2, Real X3, Real Gamma[NDFULL][NDFULL][NDFULL]) const { - //Utils::SetConnectionCoeffByFD(*this, Gamma, X0, X1, X2, X3); - const Real a2 = a_ * a_; const Real k2 = k_ * k_; const Real k3 = k_ * k_ * k_; @@ -122,7 +120,6 @@ class Snake { KOKKOS_INLINE_FUNCTION void MetricDerivative(Real X0, Real X1, Real X2, Real X3, Real dg[NDFULL][NDFULL][NDFULL]) const { - //Utils::SetMetricGradientByFD(*this, 1.e-6, X0, X1, X2, X3, dg); const Real a2 = a_ * a_; const Real k2 = k_ * k_; const Real k3 = k_ * k_ * k_; diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 8b80b1db..78364fef 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -101,6 +101,9 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { if (fluid_active) { src_names.push_back(fluid_cons::momentum); src_names.push_back(fluid_cons::energy); +#if SET_FLUX_SRC_DIAGS + src_names.push_back(diagnostic_variables::src_terms); +#endif } if (rad_active) { src_names.push_back(radmoment_cons::E); From 1daf20b72aed2825f6bdde8b51e85c1c8ca3ccc7 Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Mon, 14 Feb 2022 10:14:22 -0700 Subject: [PATCH 8/9] cleanup --- src/pgen/linear_modes.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index 2d6dae86..b6b244ec 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -127,7 +127,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (ib_lo < 0 || ib_hi < 0) { PARTHENON_THROW("physics = mhd but bfields are not present"); } - printf("ib_lo = %d ib_hi = %d\n", ib_lo, ib_hi); B10 = 1.0; if (mode == "slow") { omega = complex(0., 2.41024185339); From f20e264607cce43531c569c4df83d2b82b8b0e7b Mon Sep 17 00:00:00 2001 From: Josh Dolence Date: Mon, 14 Feb 2022 10:20:37 -0700 Subject: [PATCH 9/9] more cleanup --- src/fluid/fluid.cpp | 1 - src/fluid/tmunu.hpp | 1 - src/pgen/linear_modes.cpp | 8 +------- 3 files changed, 1 insertion(+), 9 deletions(-) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 2a264290..7f20133b 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -406,7 +406,6 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, invert.NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - //printf("cell: %d %d %d %d\n", b, k, j, i); auto status = invert(geom, eos, coords, k, j, i); fail(k, j, i) = (status == con2prim_robust::ConToPrimStatus::success ? con2prim_robust::FailFlags::success diff --git a/src/fluid/tmunu.hpp b/src/fluid/tmunu.hpp index dd0626a4..93ead1dc 100644 --- a/src/fluid/tmunu.hpp +++ b/src/fluid/tmunu.hpp @@ -126,7 +126,6 @@ class StressEnergyTensorCon { u[l] = v_(l, std::forward(args)...) - u[0] * beta[l - 1]; b[l] = iW * (b_(l, std::forward(args)...) + alpha * b[0] * u[l]); } - //printf("b[] = %g %g %g %g\n", b[0], b[1], b[2], b[3]); bsq = (Bsq + alpha * alpha * b[0] * b[0]) * iW * iW; } diff --git a/src/pgen/linear_modes.cpp b/src/pgen/linear_modes.cpp index b6b244ec..63dbc00f 100644 --- a/src/pgen/linear_modes.cpp +++ b/src/pgen/linear_modes.cpp @@ -39,9 +39,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(is_minkowski || is_boosted_minkowski || is_snake || is_inchworm, "Problem \"linear_modes\" requires \"Minkowski\" geometry!"); - printf("minkowski/boosted_minkowski/snake/inchworm: %d/%d/%d/%d\n", - is_minkowski, is_boosted_minkowski, is_snake, is_inchworm); - auto &rc = pmb->meshblock_data.Get(); const int ndim = pmb->pmy_mesh->ndim; @@ -60,7 +57,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const int ib_hi = imap[fluid_prim::bfield].second; const int iprs = imap[fluid_prim::pressure].first; const int itmp = imap[fluid_prim::temperature].first; - const int iye = imap[fluid_prim::ye].first; + const int iye = imap[fluid_prim::ye].second; const int nv = ivhi - ivlo + 1; const Real gam = pin->GetReal("eos", "Gamma"); @@ -124,9 +121,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_FAIL(msg); } } else if (physics == "mhd") { - if (ib_lo < 0 || ib_hi < 0) { - PARTHENON_THROW("physics = mhd but bfields are not present"); - } B10 = 1.0; if (mode == "slow") { omega = complex(0., 2.41024185339);