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

thread lambda through pte #466

Merged
merged 11 commits into from
Feb 13, 2025
Merged

thread lambda through pte #466

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
5a5829b
thread lambda through pte
jonahm-LANL Feb 13, 2025
9dbd27e
changelog
jonahm-LANL Feb 13, 2025
c6d8aec
test_pte_ideal now tests ideal electrons and ideal gas
jonahm-LANL Feb 13, 2025
c405d38
fix get_sg_eos_rho_p_bug
jonahm-LANL Feb 13, 2025
2845164
nlambda returns an int
jonahm-LANL Feb 13, 2025
7e0599d
eos_v should be cache[0]
jonahm-LANL Feb 13, 2025
67c5b5a
size_t
jonahm-LANL Feb 13, 2025
7e57e14
CC
jonahm-LANL Feb 13, 2025
26092e3
missed a spot
jonahm-LANL Feb 13, 2025
88f8c29
no structured binding on HIP
jonahm-LANL Feb 13, 2025
59b0034
portable inline
jonahm-LANL Feb 13, 2025
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
View file
Open in desktop
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
View file
Open in desktop
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
12 changes: 6 additions & 6 deletions singularity-eos/eos/get_sg_eos_rho_p.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
Real *ptemp_pte = &temp_pte(tid, 0);
Real *ppress_pte = &press_pte(tid, 0);
Real *pscratch = &solver_scratch(tid, 0);
PTESolverFixedP<singularity::EOSAccessor_, Real *, Real *> method(
PTESolverFixedP<singularity::EOSAccessor_, Real *, Real **> method(
npte, eos_inx, vfrac_sum, press_pte(tid, 0), prho_pte, pvfrac_pte, psie_pte,
ptemp_pte, ppress_pte, cache[0], pscratch);
ptemp_pte, ppress_pte, cache, pscratch);
auto status = PTESolver(method);
pte_converged = status.converged;
// calculate total sie
Expand All @@ -76,7 +76,7 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
// calculate sie from single eos
auto p_from_t = [&](const Real &t_i) {
return eos_v(pte_idxs(tid, 0))
.PressureFromDensityTemperature(rho_pte(tid, 0), t_i, cache[0]);
.PressureFromDensityTemperature(rho_pte(tid, 0), t_i, cache);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't this be cache[0] since this is the single material case and it is expecting the indexer of a material (the only material in this case)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oops yes that one should be. Good catch.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually I have been thinking about this a bit more. While this is correct in the vast majority of cases where npte == nmat == 1, there could be cases where 1 == npte < nmat. In this case we need to use the correct offsetting index for the density and the cache. Maybe this should be its own issue and in a separate PR, but thought I should point it out.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll create an issue. This code will also need to change in the same way I changed the C++ code to support non-trivial lambdas for, e.g., partial ionization at some point.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

all of the pte solver files likely suffer from this bug for the single material case. Wonder how many times its caused issues?

These are my fault, sorry about that this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's harmless right now, since the cache isn't filled with anything anyway for the EOS's used. But when threading through real lambdas, it will need to change. #467

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

all of the pte solver files likely suffer from this bug for the single material case. Wonder how many times its caused issues?

These are my fault, sorry about that this.

I will double check---I think it's not so bad right now because the cache is only used as an initial guess.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the rho index is also potentially incorrect though as well.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah actually no you're right it could be more serious. The eos index too?

};
// calculate sie root bounds
Real r_rho{}, r_temp{}, r_sie{}, r_press{}, r_cv{}, r_bmod{};
Expand All @@ -87,9 +87,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
Real temp_i;
RootFinding1D::findRoot(p_from_t, press_v(i), r_temp, 1.e-10 * r_temp,
1.e10 * r_temp, 1.e-12, 1.e-12, temp_i);
sie_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
.InternalEnergyFromDensityTemperature(rho_pte(tid, 0),
temp_i, cache[0]);
sie_pte(tid, 0) =
eos_v(pte_idxs(tid, 0))
.InternalEnergyFromDensityTemperature(rho_pte(tid, 0), temp_i, cache);
sie_tot_true = sie_pte(tid, 0);
// set temperature
temp_pte(tid, 0) = temp_i;
Expand Down
Loading
Loading