From 346737dcad51639c4499caedd249feeed40dd5ca Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 4 Mar 2024 13:56:43 -0500 Subject: [PATCH 1/9] entropy added in outputs --- src/fluid/fluid.cpp | 1 + src/phoebus_driver.cpp | 34 +++++++++++++++++++++++++++++++++ src/phoebus_utils/variables.hpp | 1 + 3 files changed, 36 insertions(+) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 82aa5826..cb4fe21f 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -204,6 +204,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { } physics->AddField(p::pressure::name(), mprim_scalar); physics->AddField(p::temperature::name(), mprim_scalar); + physics->AddField(p::entropy::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); if (ye) { physics->AddField(p::ye::name(), mprim_scalar); diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 432f9655..8e88f818 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -36,6 +36,7 @@ #include "phoebus_driver.hpp" #include "phoebus_package.hpp" #include "phoebus_utils/robust.hpp" +#include "phoebus_utils/variables.hpp" #include "progenitor/progenitordata.hpp" #include "radiation/radiation.hpp" #include "tov/tov.hpp" @@ -1019,6 +1020,39 @@ TaskListStatus PhoebusDriver::MonteCarloStep() { void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto tracer_pkg = pmb->packages.Get("tracers"); bool do_tracers = tracer_pkg->Param("active"); + bool is_stellarcollapse = (pin->GetString("eos", "type") == "StellarCollapse"); + + if (is_stellarcollapse) { + namespace p = fluid_prim; + auto rc = pmb->meshblock_data.Get().get(); + PackIndexMap imap; + auto v = rc->PackVariables( + {p::density::name(), p::ye::name(), p::temperature::name(), p::entropy::name()}, + imap); + + const int irho = imap[fluid_prim::density::name()].first; + const int iye = imap[fluid_prim::ye::name()].first; + const int itemp = imap[fluid_prim::temperature::name()].first; + const int is = imap[fluid_prim::entropy::name()].first; + + IndexRange ib = rc->GetBoundsI(IndexDomain::interior); + IndexRange jb = rc->GetBoundsJ(IndexDomain::interior); + IndexRange kb = rc->GetBoundsK(IndexDomain::interior); + + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + singularity::StellarCollapse eos_sc = + eos.GetUnmodifiedObject().get(); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "CoolingFunctionCalculateFourForce", DevExecSpace(), kb.s, + kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + Real lambda[2]; + lambda[0] = v(iye, k, j, i); + const Real s = eos_sc.EntropyFromDensityTemperature(v(irho, k, j, i), + v(itemp, k, j, i), lambda); + v(is, k, j, i) = s; + }); + } if (do_tracers) { auto &mbd = pmb->meshblock_data.Get(); tracers::FillTracers(mbd.get()); diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index d510af9e..a6dde283 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -49,6 +49,7 @@ VARIABLE(p, bfield); VARIABLE(p, ye); VARIABLE_NONS(pressure); VARIABLE_NONS(temperature); +VARIABLE_NONS(entropy); VARIABLE_NONS(gamma1); } // namespace fluid_prim From 3617d14a6f79958919e1b5030eb91d1d56818ca1 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 4 Mar 2024 15:58:26 -0500 Subject: [PATCH 2/9] cleaned and changed so that it's for all eos --- src/phoebus_driver.cpp | 61 +++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 33 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 8e88f818..ca788d48 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1020,39 +1020,34 @@ TaskListStatus PhoebusDriver::MonteCarloStep() { void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto tracer_pkg = pmb->packages.Get("tracers"); bool do_tracers = tracer_pkg->Param("active"); - bool is_stellarcollapse = (pin->GetString("eos", "type") == "StellarCollapse"); - - if (is_stellarcollapse) { - namespace p = fluid_prim; - auto rc = pmb->meshblock_data.Get().get(); - PackIndexMap imap; - auto v = rc->PackVariables( - {p::density::name(), p::ye::name(), p::temperature::name(), p::entropy::name()}, - imap); - - const int irho = imap[fluid_prim::density::name()].first; - const int iye = imap[fluid_prim::ye::name()].first; - const int itemp = imap[fluid_prim::temperature::name()].first; - const int is = imap[fluid_prim::entropy::name()].first; - - IndexRange ib = rc->GetBoundsI(IndexDomain::interior); - IndexRange jb = rc->GetBoundsJ(IndexDomain::interior); - IndexRange kb = rc->GetBoundsK(IndexDomain::interior); - - auto eos = pmb->packages.Get("eos")->Param("d.EOS"); - singularity::StellarCollapse eos_sc = - eos.GetUnmodifiedObject().get(); - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "CoolingFunctionCalculateFourForce", DevExecSpace(), kb.s, - kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int k, const int j, const int i) { - Real lambda[2]; - lambda[0] = v(iye, k, j, i); - const Real s = eos_sc.EntropyFromDensityTemperature(v(irho, k, j, i), - v(itemp, k, j, i), lambda); - v(is, k, j, i) = s; - }); - } + + namespace p = fluid_prim; + auto rc = pmb->meshblock_data.Get().get(); + PackIndexMap imap; + auto v = rc->PackVariables( + {p::density::name(), p::ye::name(), p::temperature::name(), p::entropy::name()}, + imap); + + const int irho = imap[p::density::name()].first; + const int iye = imap[p::ye::name()].first; + const int itemp = imap[p::temperature::name()].first; + const int is = imap[p::entropy::name()].first; + + IndexRange ib = rc->GetBoundsI(IndexDomain::interior); + IndexRange jb = rc->GetBoundsJ(IndexDomain::interior); + IndexRange kb = rc->GetBoundsK(IndexDomain::interior); + + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput", DevExecSpace(), kb.s, kb.e, jb.s, + jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { + Real lambda[2]; + lambda[0] = v(iye, k, j, i); + const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), + v(itemp, k, j, i), lambda); + v(is, k, j, i) = s; + }); + if (do_tracers) { auto &mbd = pmb->meshblock_data.Get(); tracers::FillTracers(mbd.get()); From f1ec0af67800968081aa62ed182d4cbde3afda39 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Mon, 4 Mar 2024 21:02:54 +0000 Subject: [PATCH 3/9] doc: some comments and loop label --- src/phoebus_driver.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index ca788d48..546b20fb 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1015,7 +1015,9 @@ TaskListStatus PhoebusDriver::MonteCarloStep() { /** * Gets called before output. - * Currently: Fills tracers. + * Currently: + * Fills Tracers + * Computes entropy for output **/ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto tracer_pkg = pmb->packages.Get("tracers"); @@ -1039,8 +1041,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput", DevExecSpace(), kb.s, kb.e, jb.s, - jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { + DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::Entropy", kb.s, kb.e, jb.s, jb.e, ib.s, + ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real lambda[2]; lambda[0] = v(iye, k, j, i); const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), From a7853a4c6f9d9894f37204d340ad0ffe3cf622f3 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 4 Mar 2024 16:13:53 -0500 Subject: [PATCH 4/9] turns out we do have to specify DevExecSpace() --- src/phoebus_driver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 546b20fb..0f5c717c 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1041,8 +1041,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::Entropy", kb.s, kb.e, jb.s, jb.e, ib.s, - ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { + DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::Entropy", DevExecSpace(), kb.s, kb.e, + jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real lambda[2]; lambda[0] = v(iye, k, j, i); const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), From da631dfcf5cab3af484ee4f1af1d9cac7468766a Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 4 Mar 2024 16:15:58 -0500 Subject: [PATCH 5/9] label changed --- src/phoebus_driver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 0f5c717c..84f0b47e 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1041,8 +1041,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::Entropy", DevExecSpace(), kb.s, kb.e, - jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { + DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::H5", DevExecSpace(), kb.s, kb.e, jb.s, + jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real lambda[2]; lambda[0] = v(iye, k, j, i); const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), From 3f1fb18fedd6b9cb952d90b169ac93b1c725896b Mon Sep 17 00:00:00 2001 From: mari2895 Date: Mon, 4 Mar 2024 18:08:14 -0500 Subject: [PATCH 6/9] added check iye>0 --- src/phoebus_driver.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 84f0b47e..0cb4f0b7 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1044,7 +1044,9 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::H5", DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real lambda[2]; - lambda[0] = v(iye, k, j, i); + if (iye > 0) { + lambda[0] = v(iye, k, j, i); + } const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), v(itemp, k, j, i), lambda); v(is, k, j, i) = s; From 4d9a971a078c11319941eebd670fa06ef8635b87 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Tue, 5 Mar 2024 23:08:36 -0500 Subject: [PATCH 7/9] added sound speed as an output. Also added divv/cs to find shock --- src/fluid/fluid.cpp | 2 + src/phoebus_driver.cpp | 201 +++++++++++++++++++++++++++++--- src/phoebus_utils/variables.hpp | 2 + 3 files changed, 192 insertions(+), 13 deletions(-) diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index cb4fe21f..2aade70d 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -205,6 +205,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { physics->AddField(p::pressure::name(), mprim_scalar); physics->AddField(p::temperature::name(), mprim_scalar); physics->AddField(p::entropy::name(), mprim_scalar); + physics->AddField(p::cs::name(), mprim_scalar); + physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar); physics->AddField(p::gamma1::name(), mprim_scalar); if (ye) { physics->AddField(p::ye::name(), mprim_scalar); diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 0cb4f0b7..5879b98f 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -41,9 +41,11 @@ #include "radiation/radiation.hpp" #include "tov/tov.hpp" #include "tracers/tracers.hpp" +#include using namespace parthenon::driver::prelude; using parthenon::AllReduce; +using namespace Geometry; namespace phoebus { @@ -1024,16 +1026,20 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { bool do_tracers = tracer_pkg->Param("active"); namespace p = fluid_prim; + namespace diag = diagnostic_variables; auto rc = pmb->meshblock_data.Get().get(); - PackIndexMap imap; - auto v = rc->PackVariables( - {p::density::name(), p::ye::name(), p::temperature::name(), p::entropy::name()}, - imap); + using parthenon::MakePackDescriptor; + Mesh *pmesh = rc->GetMeshPointer(); + const int ndim = pmesh->ndim; + auto &resolved_pkgs = pmesh->resolved_packages; - const int irho = imap[p::density::name()].first; - const int iye = imap[p::ye::name()].first; - const int itemp = imap[p::temperature::name()].first; - const int is = imap[p::entropy::name()].first; + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto v = desc.GetPack(rc); + auto coords = pmb->coords; + + auto geom = Geometry::GetCoordinateSystem(rc); IndexRange ib = rc->GetBoundsI(IndexDomain::interior); IndexRange jb = rc->GetBoundsJ(IndexDomain::interior); @@ -1044,12 +1050,181 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { DEFAULT_LOOP_PATTERN, "UserWorkBeforeOutput::H5", DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real lambda[2]; - if (iye > 0) { - lambda[0] = v(iye, k, j, i); + if (v.Contains(0, p::ye())) { + lambda[0] = v(0, p::ye(), k, j, i); + } + const Real s = eos.EntropyFromDensityTemperature( + v(0, p::density(), k, j, i), v(0, p::temperature(), k, j, i), lambda); + const Real p = eos.PressureFromDensityTemperature( + v(0, p::density(), k, j, i), v(0, p::temperature(), k, j, i), lambda); + const Real bmod = eos.BulkModulusFromDensityTemperature( + v(0, p::density(), k, j, i), v(0, p::temperature(), k, j, i), lambda); + const Real sie = eos.InternalEnergyFromDensityTemperature( + v(0, p::density(), k, j, i), v(0, p::temperature(), k, j, i), lambda); + const Real h = sie + p / v(0, p::density(), k, j, i) + 1; + const Real cs = std::sqrt(bmod / v(0, p::density(), k, j, i) / h); + Real divv; + Real gam[3][3]; + Real gam1[3][3]; + Real gam2[3][3]; + Real gam3[3][3]; + Real gam12[3][3]; + Real gam13[3][3]; + Real gam23[3][3]; + Real gam123[3][3]; + const Real vp[3] = {v(0, p::velocity(0), k, j, i), v(0, p::velocity(1), k, j, i), + v(0, p::velocity(2), k, j, i)}; + const Real vp1[3] = {v(0, p::velocity(0), k, j, i - 1), + v(0, p::velocity(1), k, j, i - 1), + v(0, p::velocity(2), k, j, i - 1)}; + geom.Metric(CellLocation::Cent, 0, k, j, i, gam); + geom.Metric(CellLocation::Cent, 0, k, j, i - 1, gam1); + Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i); + const Real W = phoebus::GetLorentzFactor(vp, gam); + const Real W1 = phoebus::GetLorentzFactor(vp1, gam1); + + if (ndim == 1) { + divv = 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(0), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(0), k, j, i - 1) / W1) / + coords.CellWidthFA(X1DIR, k, j, i); + } + + if (ndim == 2) { + Real vp2[3] = {v(0, p::velocity(0), k, j - 1, i), + v(0, p::velocity(1), k, j - 1, i), + v(0, p::velocity(2), k, j - 1, i)}; + Real vp3[3] = {v(0, p::velocity(0), k - 1, j, i), + v(0, p::velocity(1), k - 1, j, i), + v(0, p::velocity(2), k - 1, j, i)}; + Real vp12[3] = {v(0, p::velocity(0), k, j - 1, i - 1), + v(0, p::velocity(1), k, j - 1, i - 1), + v(0, p::velocity(2), k, j - 1, i - 1)}; + geom.Metric(CellLocation::Cent, 0, k, j - 1, i, gam2); + geom.Metric(CellLocation::Cent, 0, k - 1, j, i, gam3); + geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam12); + Real W2 = phoebus::GetLorentzFactor(vp2, gam2); + Real W3 = phoebus::GetLorentzFactor(vp3, gam3); + Real W12 = phoebus::GetLorentzFactor(vp12, gam12); + + divv = 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(0), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(0), k, j, i - 1) / W1 + + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(0), k, j - 1, i) / W2 - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * + v(0, p::velocity(0), k, j - 1, i - 1) / W12) / + coords.CellWidthFA(X1DIR, k, j, i) / 2 + + + 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(1), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(1), k, j - 1, i) / W2 + + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(1), k, j, i - 1) / W1 - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * + v(0, p::velocity(1), k, j - 1, i - 1) / W12) / + coords.CellWidthFA(X2DIR, k, j, i) / 2; } - const Real s = eos.EntropyFromDensityTemperature(v(irho, k, j, i), - v(itemp, k, j, i), lambda); - v(is, k, j, i) = s; + + if (ndim == 3) { + Real vp2[3] = {v(0, p::velocity(0), k, j - 1, i), + v(0, p::velocity(1), k, j - 1, i), + v(0, p::velocity(2), k, j - 1, i)}; + Real vp3[3] = {v(0, p::velocity(0), k - 1, j, i), + v(0, p::velocity(1), k - 1, j, i), + v(0, p::velocity(2), k - 1, j, i)}; + Real vp12[3] = {v(0, p::velocity(0), k, j - 1, i - 1), + v(0, p::velocity(1), k, j - 1, i - 1), + v(0, p::velocity(2), k, j - 1, i - 1)}; + geom.Metric(CellLocation::Cent, 0, k, j - 1, i, gam2); + geom.Metric(CellLocation::Cent, 0, k - 1, j, i, gam3); + geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam12); + Real W2 = phoebus::GetLorentzFactor(vp2, gam2); + Real W3 = phoebus::GetLorentzFactor(vp3, gam3); + Real W12 = phoebus::GetLorentzFactor(vp12, gam12); + Real vp13[3] = {v(0, p::velocity(0), k - 1, j, i - 1), + v(0, p::velocity(1), k - 1, j, i - 1), + v(0, p::velocity(2), k - 1, j, i - 1)}; + Real vp23[3] = {v(0, p::velocity(0), k - 1, j - 1, i), + v(0, p::velocity(1), k - 1, j - 1, i), + v(0, p::velocity(2), k - 1, j - 1, i)}; + Real vp123[3] = {v(0, p::velocity(0), k - 1, j - 1, i - 1), + v(0, p::velocity(1), k - 1, j - 1, i - 1), + v(0, p::velocity(2), k - 1, j - 1, i - 1)}; + geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam13); + geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam23); + geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam123); + Real W13 = phoebus::GetLorentzFactor(vp13, gam13); + Real W23 = phoebus::GetLorentzFactor(vp23, gam23); + Real W123 = phoebus::GetLorentzFactor(vp123, gam123); + + divv = 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(0), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(0), k, j, i - 1) / W1 + + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(0), k, j - 1, i) / W2 - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * + v(0, p::velocity(0), k, j - 1, i - 1) / W12 + + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * + v(0, p::velocity(0), k - 1, j, i) / W3 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * + v(0, p::velocity(0), k - 1, j, i - 1) / W13 + + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * + v(0, p::velocity(0), k - 1, j - 1, i) / W23 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * + v(0, p::velocity(0), k - 1, j - 1, i - 1) / W123) / + coords.CellWidthFA(X1DIR, k, j, i) / 4 + + + 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(1), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(1), k, j - 1, i) / W2 + + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(1), k, j, i - 1) / W1 - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * + v(0, p::velocity(1), k, j - 1, i - 1) / W12 + + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * + v(0, p::velocity(1), k - 1, j, i) / W3 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * + v(0, p::velocity(1), k - 1, j - 1, i) / W23 + + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * + v(0, p::velocity(1), k - 1, j, i - 1) / W13 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * + v(0, p::velocity(1), k - 1, j - 1, i - 1) / W123) / + coords.CellWidthFA(X2DIR, k, j, i) / 4 + + + 1 / gdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(2), k, j, i) / W - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * + v(0, p::velocity(2), k - 1, j, i) / W3 + + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(2), k, j - 1, i) / W2 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * + v(0, p::velocity(2), k - 1, j - 1, i) / W23 + + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(2), k, j, i - 1) / W1 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * + v(0, p::velocity(2), k - 1, j, i - 1) / W13 + + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * + v(0, p::velocity(2), k, j - 1, i - 1) / W12 - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * + v(0, p::velocity(2), k - 1, j - 1, i - 1) / W123) / + coords.CellWidthFA(X3DIR, k, j, i) / 4; + } + + v(0, p::entropy(), k, j, i) = s; + v(0, p::cs(), k, j, i) = cs; + v(0, diag::ratio_divv_cs(), k, j, i) = divv / cs; }); if (do_tracers) { diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index a6dde283..4746c5ba 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -50,6 +50,7 @@ VARIABLE(p, ye); VARIABLE_NONS(pressure); VARIABLE_NONS(temperature); VARIABLE_NONS(entropy); +VARIABLE_NONS(cs); VARIABLE_NONS(gamma1); } // namespace fluid_prim @@ -119,6 +120,7 @@ VARIABLE_CUSTOM(node_coords, g.n.coord); namespace diagnostic_variables { VARIABLE_NONS(divb); +VARIABLE_NONS(ratio_divv_cs); VARIABLE_CUSTOM(divf, flux_divergence); VARIABLE_NONS(src_terms); VARIABLE_CUSTOM(r_divf, r.flux_divergence); From 425141483683e8256d343ba9b6a6a18dd8c97503 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 6 Mar 2024 00:39:26 -0500 Subject: [PATCH 8/9] copy-paste error fixed that caused linear_modes test to fail --- src/phoebus_driver.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 5879b98f..2d17ead2 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1096,17 +1096,12 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { Real vp2[3] = {v(0, p::velocity(0), k, j - 1, i), v(0, p::velocity(1), k, j - 1, i), v(0, p::velocity(2), k, j - 1, i)}; - Real vp3[3] = {v(0, p::velocity(0), k - 1, j, i), - v(0, p::velocity(1), k - 1, j, i), - v(0, p::velocity(2), k - 1, j, i)}; Real vp12[3] = {v(0, p::velocity(0), k, j - 1, i - 1), v(0, p::velocity(1), k, j - 1, i - 1), v(0, p::velocity(2), k, j - 1, i - 1)}; geom.Metric(CellLocation::Cent, 0, k, j - 1, i, gam2); - geom.Metric(CellLocation::Cent, 0, k - 1, j, i, gam3); geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam12); Real W2 = phoebus::GetLorentzFactor(vp2, gam2); - Real W3 = phoebus::GetLorentzFactor(vp3, gam3); Real W12 = phoebus::GetLorentzFactor(vp12, gam12); divv = 1 / gdet * From 3e8f085a156ce9e4490cdee7d20ef4f91b3d1664 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 6 Mar 2024 13:30:05 -0500 Subject: [PATCH 9/9] no longer averaging --- src/phoebus_driver.cpp | 151 ++++++++--------------------------------- 1 file changed, 27 insertions(+), 124 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 2d17ead2..5c56bd01 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -1068,10 +1068,6 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { Real gam1[3][3]; Real gam2[3][3]; Real gam3[3][3]; - Real gam12[3][3]; - Real gam13[3][3]; - Real gam23[3][3]; - Real gam123[3][3]; const Real vp[3] = {v(0, p::velocity(0), k, j, i), v(0, p::velocity(1), k, j, i), v(0, p::velocity(2), k, j, i)}; const Real vp1[3] = {v(0, p::velocity(0), k, j, i - 1), @@ -1080,141 +1076,48 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { geom.Metric(CellLocation::Cent, 0, k, j, i, gam); geom.Metric(CellLocation::Cent, 0, k, j, i - 1, gam1); Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i); + Real igdet = robust::ratio(1., gdet); + Real lapse = geom.Lapse(CellLocation::Cent, 0, k, j, i); + Real lapse1 = geom.Lapse(CellLocation::Cent, 0, k, j, i - 1); const Real W = phoebus::GetLorentzFactor(vp, gam); const Real W1 = phoebus::GetLorentzFactor(vp1, gam1); - if (ndim == 1) { - divv = 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(0), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(0), k, j, i - 1) / W1) / - coords.CellWidthFA(X1DIR, k, j, i); - } + divv = igdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(0), k, j, i) * lapse / W - + geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * + v(0, p::velocity(0), k, j, i - 1) * lapse1 / W1) / + coords.CellWidthFA(X1DIR, k, j, i); - if (ndim == 2) { + if (ndim >= 2) { Real vp2[3] = {v(0, p::velocity(0), k, j - 1, i), v(0, p::velocity(1), k, j - 1, i), v(0, p::velocity(2), k, j - 1, i)}; - Real vp12[3] = {v(0, p::velocity(0), k, j - 1, i - 1), - v(0, p::velocity(1), k, j - 1, i - 1), - v(0, p::velocity(2), k, j - 1, i - 1)}; geom.Metric(CellLocation::Cent, 0, k, j - 1, i, gam2); - geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam12); Real W2 = phoebus::GetLorentzFactor(vp2, gam2); - Real W12 = phoebus::GetLorentzFactor(vp12, gam12); - - divv = 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(0), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(0), k, j, i - 1) / W1 + - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * - v(0, p::velocity(0), k, j - 1, i) / W2 - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * - v(0, p::velocity(0), k, j - 1, i - 1) / W12) / - coords.CellWidthFA(X1DIR, k, j, i) / 2 - - + 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(1), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * - v(0, p::velocity(1), k, j - 1, i) / W2 + - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(1), k, j, i - 1) / W1 - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * - v(0, p::velocity(1), k, j - 1, i - 1) / W12) / - coords.CellWidthFA(X2DIR, k, j, i) / 2; + Real lapse2 = geom.Lapse(CellLocation::Cent, 0, k, j - 1, i); + + divv += igdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(1), k, j, i) * lapse / W - + geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * + v(0, p::velocity(1), k, j - 1, i) * lapse2 / W2) / + coords.CellWidthFA(X2DIR, k, j, i); } - if (ndim == 3) { - Real vp2[3] = {v(0, p::velocity(0), k, j - 1, i), - v(0, p::velocity(1), k, j - 1, i), - v(0, p::velocity(2), k, j - 1, i)}; + if (ndim >= 3) { Real vp3[3] = {v(0, p::velocity(0), k - 1, j, i), v(0, p::velocity(1), k - 1, j, i), v(0, p::velocity(2), k - 1, j, i)}; - Real vp12[3] = {v(0, p::velocity(0), k, j - 1, i - 1), - v(0, p::velocity(1), k, j - 1, i - 1), - v(0, p::velocity(2), k, j - 1, i - 1)}; - geom.Metric(CellLocation::Cent, 0, k, j - 1, i, gam2); - geom.Metric(CellLocation::Cent, 0, k - 1, j, i, gam3); - geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam12); - Real W2 = phoebus::GetLorentzFactor(vp2, gam2); Real W3 = phoebus::GetLorentzFactor(vp3, gam3); - Real W12 = phoebus::GetLorentzFactor(vp12, gam12); - Real vp13[3] = {v(0, p::velocity(0), k - 1, j, i - 1), - v(0, p::velocity(1), k - 1, j, i - 1), - v(0, p::velocity(2), k - 1, j, i - 1)}; - Real vp23[3] = {v(0, p::velocity(0), k - 1, j - 1, i), - v(0, p::velocity(1), k - 1, j - 1, i), - v(0, p::velocity(2), k - 1, j - 1, i)}; - Real vp123[3] = {v(0, p::velocity(0), k - 1, j - 1, i - 1), - v(0, p::velocity(1), k - 1, j - 1, i - 1), - v(0, p::velocity(2), k - 1, j - 1, i - 1)}; - geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam13); - geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam23); - geom.Metric(CellLocation::Cent, 0, k, j - 1, i - 1, gam123); - Real W13 = phoebus::GetLorentzFactor(vp13, gam13); - Real W23 = phoebus::GetLorentzFactor(vp23, gam23); - Real W123 = phoebus::GetLorentzFactor(vp123, gam123); - - divv = 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(0), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(0), k, j, i - 1) / W1 + - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * - v(0, p::velocity(0), k, j - 1, i) / W2 - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * - v(0, p::velocity(0), k, j - 1, i - 1) / W12 + - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * - v(0, p::velocity(0), k - 1, j, i) / W3 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * - v(0, p::velocity(0), k - 1, j, i - 1) / W13 + - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * - v(0, p::velocity(0), k - 1, j - 1, i) / W23 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * - v(0, p::velocity(0), k - 1, j - 1, i - 1) / W123) / - coords.CellWidthFA(X1DIR, k, j, i) / 4 - - + 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(1), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * - v(0, p::velocity(1), k, j - 1, i) / W2 + - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(1), k, j, i - 1) / W1 - - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * - v(0, p::velocity(1), k, j - 1, i - 1) / W12 + - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * - v(0, p::velocity(1), k - 1, j, i) / W3 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * - v(0, p::velocity(1), k - 1, j - 1, i) / W23 + - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * - v(0, p::velocity(1), k - 1, j, i - 1) / W13 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * - v(0, p::velocity(1), k - 1, j - 1, i - 1) / W123) / - coords.CellWidthFA(X2DIR, k, j, i) / 4 - - + 1 / gdet * - (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * - v(0, p::velocity(2), k, j, i) / W - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * - v(0, p::velocity(2), k - 1, j, i) / W3 + - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i) * - v(0, p::velocity(2), k, j - 1, i) / W2 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i) * - v(0, p::velocity(2), k - 1, j - 1, i) / W23 + - geom.DetGamma(CellLocation::Cent, 0, k, j, i - 1) * - v(0, p::velocity(2), k, j, i - 1) / W1 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i - 1) * - v(0, p::velocity(2), k - 1, j, i - 1) / W13 + - geom.DetGamma(CellLocation::Cent, 0, k, j - 1, i - 1) * - v(0, p::velocity(2), k, j - 1, i - 1) / W12 - - geom.DetGamma(CellLocation::Cent, 0, k - 1, j - 1, i - 1) * - v(0, p::velocity(2), k - 1, j - 1, i - 1) / W123) / - coords.CellWidthFA(X3DIR, k, j, i) / 4; + Real lapse3 = geom.Lapse(CellLocation::Cent, 0, k - 1, j, i); + + divv += igdet * + (geom.DetGamma(CellLocation::Cent, 0, k, j, i) * + v(0, p::velocity(2), k, j, i) * lapse / W - + geom.DetGamma(CellLocation::Cent, 0, k - 1, j, i) * + v(0, p::velocity(2), k - 1, j, i) * lapse3 / W3) / + coords.CellWidthFA(X3DIR, k, j, i); } v(0, p::entropy(), k, j, i) = s;