Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

entropy added in outputs #201

Merged
merged 9 commits into from
Mar 6, 2024
Merged

entropy added in outputs #201

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
346737d
entropy added in outputs
mari2895 Mar 4, 2024
3617d14
cleaned and changed so that it's for all eos
mari2895 Mar 4, 2024
f1ec0af
doc: some comments and loop label
AstroBarker Mar 4, 2024
a7853a4
turns out we do have to specify DevExecSpace()
mari2895 Mar 4, 2024
da631df
label changed
mari2895 Mar 4, 2024
3f1fb18
added check iye>0
mari2895 Mar 4, 2024
4d9a971
added sound speed as an output. Also added divv/cs to find shock
mari2895 Mar 6, 2024
4251414
copy-paste error fixed that caused linear_modes test to fail
mari2895 Mar 6, 2024
3e8f085
no longer averaging
mari2895 Mar 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/fluid/fluid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,6 +204,9 @@ std::shared_ptr<StateDescriptor> 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);
Expand Down
205 changes: 204 additions & 1 deletion src/phoebus_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,16 @@
#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"
#include "tracers/tracers.hpp"
#include <interface/sparse_pack.hpp>

using namespace parthenon::driver::prelude;
using parthenon::AllReduce;
using namespace Geometry;

namespace phoebus {

Expand Down Expand Up @@ -1014,11 +1017,211 @@ 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");
bool do_tracers = tracer_pkg->Param<bool>("active");

namespace p = fluid_prim;
namespace diag = diagnostic_variables;
auto rc = pmb->meshblock_data.Get().get();
using parthenon::MakePackDescriptor;
Mesh *pmesh = rc->GetMeshPointer();
const int ndim = pmesh->ndim;
auto &resolved_pkgs = pmesh->resolved_packages;

static auto desc =
MakePackDescriptor<p::velocity, p::density, p::ye, p::temperature, p::entropy,
p::cs, diag::ratio_divv_cs>(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);
IndexRange kb = rc->GetBoundsK(IndexDomain::interior);

auto eos = pmb->packages.Get("eos")->Param<Microphysics::EOS::EOS>("d.EOS");
parthenon::par_for(
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 (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)};
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
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) /
mari2895 marked this conversation as resolved.
Show resolved Hide resolved
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 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)};
mari2895 marked this conversation as resolved.
Show resolved Hide resolved
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
mari2895 marked this conversation as resolved.
Show resolved Hide resolved

+ 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;
}

if (ndim == 3) {
mari2895 marked this conversation as resolved.
Show resolved Hide resolved
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) {
auto &mbd = pmb->meshblock_data.Get();
tracers::FillTracers(mbd.get());
Expand Down
3 changes: 3 additions & 0 deletions src/phoebus_utils/variables.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ VARIABLE(p, bfield);
VARIABLE(p, ye);
VARIABLE_NONS(pressure);
VARIABLE_NONS(temperature);
VARIABLE_NONS(entropy);
VARIABLE_NONS(cs);
VARIABLE_NONS(gamma1);
} // namespace fluid_prim

Expand Down Expand Up @@ -118,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);
Expand Down
Loading