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.

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

thread lambda through pte #466

Merged
merged 11 commits into from
Feb 13, 2025
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
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
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
- [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS

### Fixed (Repair bugs, etc)
- [[PR466]](https://github.com/lanl/singularity-eos/pull/466) Fully thread lambda inputs through the PTE solver
- [[PR449]](https://github.com/lanl/singularity-eos/pull/449) Ensure that DensityEnergyFromPressureTemperature works for all equations of state and is properly tested
- [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API
- [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test
Expand Down
73 changes: 50 additions & 23 deletions singularity-eos/closure/mixed_cell_models.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -478,10 +478,11 @@ class PTESolverBase {

} // namespace mix_impl

template <typename EOSIndexer, typename RealIndexer>
template <typename EOSIndexer, typename RealIndexer, typename LambdaIndexer = NullIndexer>
PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU(
const std::size_t nmat, EOSIndexer &&eos, const Real u_tot, RealIndexer &&rho,
RealIndexer &&vfrac, const Real Tguess = 0.0) {
RealIndexer &&vfrac, const Real Tguess = 0.0,
LambdaIndexer &&lambda = NullIndexer()) {
// should these be passed in?
constexpr Real minimum_temperature = 1.e-9;
constexpr Real maximum_temperature = 1.e9;
Expand All @@ -493,7 +494,8 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU(
auto ufunc = [&](const Real T) {
Real usum = 0.0;
for (std::size_t m = 0; m < nmat; m++) {
usum += rho[m] * vfrac[m] * eos[m].InternalEnergyFromDensityTemperature(rho[m], T);
usum += rho[m] * vfrac[m] *
eos[m].InternalEnergyFromDensityTemperature(rho[m], T, lambda[m]);
}
return usum;
};
Expand Down Expand Up @@ -593,10 +595,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
dedv = AssignIncrement(scratch, nmat);
dpdT = AssignIncrement(scratch, nmat);
vtemp = AssignIncrement(scratch, nmat);
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (!variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
// JMM: Copy lambda data into cache. Since the cache uses
// pointers, it must be done this way.
for (std::size_t m = 0; m < nmat; ++m) {
if (!variadic_utils::is_nullptr(lambda[m])) {
for (std::size_t l = 0; l < eos[m].nlambda(); ++l) {
Cache[m][l] = lambda[m][l];
}
}
}
}

PORTABLE_INLINE_FUNCTION
Expand Down Expand Up @@ -819,10 +826,15 @@ class PTESolverPT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
: mix_impl::PTESolverBase<EOSIndexer, RealIndexer>(nmat, 2, eos, vfrac_tot, sie_tot,
rho, vfrac, sie, temp, press,
scratch, Tnorm, params) {
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (!variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
// JMM: Copy lambda data into cache. Since the cache uses
// pointers, it must be done this way.
for (std::size_t m = 0; m < nmat; ++m) {
if (!variadic_utils::is_nullptr(lambda[m])) {
for (std::size_t l = 0; l < eos[m].nlambda(); ++l) {
Cache[m][l] = lambda[m][l];
}
}
}
}

PORTABLE_INLINE_FUNCTION
Expand Down Expand Up @@ -1066,10 +1078,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
Tequil = T_true;
Ttemp = T_true;
Tnorm = 1.0;
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
// JMM: Copy lambda data into cache. Since the cache uses
// pointers, it must be done this way.
for (std::size_t m = 0; m < nmat; ++m) {
if (!variadic_utils::is_nullptr(lambda[m])) {
for (std::size_t l = 0; l < eos[m].nlambda(); ++l) {
Cache[m][l] = lambda[m][l];
}
}
}
}

PORTABLE_INLINE_FUNCTION
Expand Down Expand Up @@ -1275,10 +1292,15 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
dpdT = AssignIncrement(scratch, nmat);
vtemp = AssignIncrement(scratch, nmat);
Pequil = P;
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
// JMM: Copy lambda data into cache. Since the cache uses
// pointers, it must be done this way.
for (std::size_t m = 0; m < nmat; ++m) {
if (!variadic_utils::is_nullptr(lambda[m])) {
for (std::size_t l = 0; l < eos[m].nlambda(); ++l) {
Cache[m][l] = lambda[m][l];
}
}
}
}

PORTABLE_INLINE_FUNCTION
Expand Down Expand Up @@ -1514,10 +1536,15 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
dtde = AssignIncrement(scratch, nmat);
vtemp = AssignIncrement(scratch, nmat);
utemp = AssignIncrement(scratch, nmat);
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
// JMM: Copy lambda data into cache. Since the cache uses
// pointers, it must be done this way.
for (std::size_t m = 0; m < nmat; ++m) {
if (!variadic_utils::is_nullptr(lambda[m])) {
for (std::size_t l = 0; l < eos[m].nlambda(); ++l) {
Cache[m][l] = lambda[m][l];
}
}
}
}

PORTABLE_INLINE_FUNCTION
Expand Down
26 changes: 16 additions & 10 deletions test/test_pte_ideal.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,20 @@
#include <test/eos_unit_test_helpers.hpp>

using singularity::ApproxTemperatureFromRhoMatU;
using singularity::IdealGas;
using singularity::IdealElectrons;
using singularity::MeanAtomicProperties;
using singularity::PTESolverPT;
using singularity::PTESolverPTRequiredScratch;
using singularity::PTESolverRhoT;
using singularity::PTESolverRhoTRequiredScratch;
using singularity::Variant;
using EOS = Variant<IdealGas>;
using EOS = Variant<IdealElectrons>;

struct LambdaIndexer {
Real *operator[](const int i) { return &z; }
const Real *operator[](const int i) const { return &z; }
Real z = 0.9;
};

template <typename RealIndexer, typename EOSIndexer>
PORTABLE_INLINE_FUNCTION Real set_state(Real rho_nom, Real sie_nom, RealIndexer &&rho,
Expand Down Expand Up @@ -78,8 +85,9 @@ SCENARIO("PT space PTE solver for two ideal gases", "[PTESolverPT][IdealGas]") {
EOS_Indexer_t eoss;

for (std::size_t m = 0; m < NEOS; ++m) {
constexpr Real denom = NEOS + 1;
eoss[m] = IdealGas((m + 1) / denom, m + 1);
constexpr Real Abar = NEOS + 1;
MeanAtomicProperties AZBar(Abar, m + 1);
eoss[m] = IdealElectrons(AZBar);
}

WHEN("We create a thermodynamic state") {
Expand Down Expand Up @@ -110,10 +118,8 @@ SCENARIO("PT space PTE solver for two ideal gases", "[PTESolverPT][IdealGas]") {
Real *sies = (Real *)malloc(sizeof(Real) * NEOS);
Real *Ts = (Real *)malloc(sizeof(Real) * NEOS);
Real *Ps = (Real *)malloc(sizeof(Real) * NEOS);
Real *lambda[NEOS];
for (std::size_t i = 0; i < NEOS; i++) {
lambda[i] = nullptr;
}
LambdaIndexer lambda;

Real *scratch_rt = (Real *)malloc(sizeof(Real) * nscratch_rt);
Real *scratch_pt = (Real *)malloc(sizeof(Real) * nscratch_pt);

Expand All @@ -123,8 +129,8 @@ SCENARIO("PT space PTE solver for two ideal gases", "[PTESolverPT][IdealGas]") {
for (std::size_t m = 0; m < NEOS; ++m) {
rho_tot += rhos[m] * vfracs[m];
}
const Real Tguess =
ApproxTemperatureFromRhoMatU(NEOS, eoss, rho_tot * sie_tot, rhos, vfracs);
const Real Tguess = ApproxTemperatureFromRhoMatU(
NEOS, eoss, rho_tot * sie_tot, rhos, vfracs, 0.0, lambda);

auto method_rt = PTESolverRhoT<EOS_Indexer_t, Real *, decltype(lambda)>(
NEOS, eoss, 1.0, sie_tot, rhos, vfracs, sies, Ts, Ps, lambda, scratch_rt,
Expand Down
Loading