Skip to content

Commit e210832

Browse files
zingalezingale
authored
wrap some of the entropy calcs in if constexpr (#1912)
this avoids a log and sqrt and a few divides note: there is one place where I swapped a divide by multiply by inverse, so some small roundoff might happen
1 parent 64cb627 commit e210832

File tree

1 file changed

+74
-63
lines changed

1 file changed

+74
-63
lines changed

‎EOS/helmholtz/actual_eos.H‎

Lines changed: 74 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -117,32 +117,19 @@ void fwt (const amrex::Real* AMREX_RESTRICT fi,
117117

118118
// cubic hermite polynomial functions
119119
// psi0 & derivatives
120-
AMREX_GPU_HOST_DEVICE AMREX_INLINE
120+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
121121
amrex::Real xpsi0(amrex::Real z) noexcept
122122
{
123123
return z * z * (2.0e0_rt * z - 3.0e0_rt) + 1.0_rt;
124124
}
125125

126-
AMREX_GPU_HOST_DEVICE AMREX_INLINE
127-
amrex::Real xdpsi0(amrex::Real z) noexcept
128-
{
129-
return z * (6.0e0_rt * z - 6.0e0_rt);
130-
}
131-
132126
// psi1 & derivatives
133-
AMREX_GPU_HOST_DEVICE AMREX_INLINE
127+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
134128
amrex::Real xpsi1(amrex::Real z) noexcept
135129
{
136130
return z * (z * (z - 2.0e0_rt) + 1.0e0_rt);
137131
}
138132

139-
AMREX_GPU_HOST_DEVICE AMREX_INLINE
140-
amrex::Real xdpsi1(amrex::Real z) noexcept
141-
{
142-
return z * (3.0e0_rt * z - 4.0e0_rt) + 1.0e0_rt;
143-
}
144-
145-
146133

147134
namespace electron_table_indexing {
148135

@@ -334,15 +321,16 @@ void apply_electrons (T& state)
334321
// giving the j offset, i offset, and component
335322

336323
// pressure derivative with density
337-
amrex::Real dpepdd = 0.0e0_rt;
338-
amrex::constexpr_for<0, 16>([&] (auto K) {
339-
constexpr int k = K;
340-
constexpr auto m = electron_table_indexing::map[k];
341-
342-
dpepdd += dpdf[jat+m.dj][iat+m.di][m.n] * wdt[k];
343-
});
344-
dpepdd = amrex::max(state.y_e * dpepdd, 0.0e0_rt);
324+
[[maybe_unused]] amrex::Real dpepdd = 0.0e0_rt;
325+
if constexpr (has_pressure<T>::value) {
326+
amrex::constexpr_for<0, 16>([&] (auto K) {
327+
constexpr int k = K;
328+
constexpr auto m = electron_table_indexing::map[k];
345329

330+
dpepdd += dpdf[jat+m.dj][iat+m.di][m.n] * wdt[k];
331+
});
332+
dpepdd = amrex::max(state.y_e * dpepdd, 0.0e0_rt);
333+
}
346334

347335
// electron chemical potential etaele
348336
[[maybe_unused]] amrex::Real etaele{};
@@ -377,13 +365,15 @@ void apply_electrons (T& state)
377365
amrex::Real dpepdt = x * df_dt;
378366
[[maybe_unused]] amrex::Real dpepda{};
379367
[[maybe_unused]] amrex::Real dpepdz{};
380-
if constexpr (has_dpdA<T>::value || has_dpdZ<T>::value) {
381-
amrex::Real s = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
382-
if constexpr (has_dpdA<T>::value) {
383-
dpepda = -ytot1 * (2.0e0_rt * pele + s * din);
384-
}
385-
if constexpr (has_dpdZ<T>::value) {
386-
dpepdz = state.rho*ytot1*(2.0e0_rt * din * df_d + s);
368+
if constexpr (has_pressure<T>::value) {
369+
if constexpr (has_dpdA<T>::value || has_dpdZ<T>::value) {
370+
amrex::Real s = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
371+
if constexpr (has_dpdA<T>::value) {
372+
dpepda = -ytot1 * (2.0e0_rt * pele + s * din);
373+
}
374+
if constexpr (has_dpdZ<T>::value) {
375+
dpepdz = state.rho*ytot1*(2.0e0_rt * din * df_d + s);
376+
}
387377
}
388378
}
389379

@@ -487,16 +477,21 @@ void apply_ions (T& state)
487477
[[maybe_unused]] amrex::Real deionda = 1.5e0_rt * dpionda * deni;
488478
[[maybe_unused]] amrex::Real deiondz = 0.0e0_rt;
489479

490-
amrex::Real x = state.abar * state.abar * std::sqrt(state.abar) * deni / avo_eos;
491-
amrex::Real s = sioncon * state.T;
492-
amrex::Real z = x * s * std::sqrt(s);
493-
amrex::Real y = std::log(z);
494-
amrex::Real sion = (pion * deni + eion) * tempi + kergavo * ytot1 * y;
495-
amrex::Real dsiondd = (dpiondd * deni - pion * deni * deni + deiondd) * tempi -
496-
kergavo * deni * ytot1;
497-
amrex::Real dsiondt = (dpiondt * deni + deiondt) * tempi -
498-
(pion * deni + eion) * tempi * tempi +
499-
1.5e0_rt * kergavo * tempi * ytot1;
480+
[[maybe_unused]] amrex::Real sion{};
481+
[[maybe_unused]] amrex::Real dsiondd{};
482+
[[maybe_unused]] amrex::Real dsiondt{};
483+
if constexpr (has_entropy<T>::value) {
484+
amrex::Real x = state.abar * state.abar * std::sqrt(state.abar) * deni / avo_eos;
485+
amrex::Real s = sioncon * state.T;
486+
amrex::Real z = x * s * std::sqrt(s);
487+
amrex::Real y = std::log(z);
488+
sion = (pion * deni + eion) * tempi + kergavo * ytot1 * y;
489+
dsiondd = (dpiondd * deni - pion * deni * deni + deiondd) * tempi -
490+
kergavo * deni * ytot1;
491+
dsiondt = (dpiondt * deni + deiondt) * tempi -
492+
(pion * deni + eion) * tempi * tempi +
493+
1.5e0_rt * kergavo * tempi * ytot1;
494+
}
500495

501496
if constexpr (has_pressure<T>::value) {
502497
state.p = state.p + pion;
@@ -570,9 +565,14 @@ void apply_radiation (T& state)
570565
amrex::Real deraddd = -erad * deni;
571566
amrex::Real deraddt = 3.0e0_rt * dpraddt * deni;
572567

573-
amrex::Real srad = (prad * deni + erad) * tempi;
574-
amrex::Real dsraddd = (dpraddd * deni - prad * deni * deni + deraddd) * tempi;
575-
amrex::Real dsraddt = (dpraddt * deni + deraddt - srad) * tempi;
568+
[[maybe_unused]] amrex::Real srad{};
569+
[[maybe_unused]] amrex::Real dsraddd{};
570+
[[maybe_unused]] amrex::Real dsraddt{};
571+
if constexpr (has_entropy<T>::value) {
572+
srad = (prad * deni + erad) * tempi;
573+
dsraddd = (dpraddd * deni - prad * deni * deni + deraddd) * tempi;
574+
dsraddt = (dpraddt * deni + deraddt - srad) * tempi;
575+
}
576576

577577
// Note that unlike the other terms, radiation
578578
// sets these terms instead of adding to them,
@@ -647,7 +647,7 @@ void apply_coulomb_corrections (T& state)
647647
[[maybe_unused]] amrex::Real decoulda = 0.e0_rt;
648648
[[maybe_unused]] amrex::Real decouldz = 0.e0_rt;
649649

650-
amrex::Real s, x, y, z; // temporary variables
650+
amrex::Real s, z; // temporary variables
651651

652652
// uniform background corrections only
653653
// from yakovlev & shalybkov 1989
@@ -683,18 +683,21 @@ void apply_coulomb_corrections (T& state)
683683
// .yakovlev & shalybkov 1989 equations 82, 85, 86, 87
684684
if (plasg >= 1.0e0_rt)
685685
{
686-
x = std::pow(plasg, 0.25e0_rt);
687-
y = avo_eos * ytot1 * kerg;
686+
amrex::Real x = std::pow(plasg, 0.25e0_rt);
687+
amrex::Real y = avo_eos * ytot1 * kerg;
688688
ecoul = y * state.T * (a1 * plasg + b1 * x + c1 / x + d1);
689689
pcoul = onethird * state.rho * ecoul;
690-
scoul = -y * (3.0e0_rt * b1 * x - 5.0e0_rt*c1 / x +
691-
d1 * (std::log(plasg) - 1.0e0_rt) - e1);
690+
691+
if constexpr (has_entropy<T>::value) {
692+
scoul = -y * (3.0e0_rt * b1 * x - 5.0e0_rt*c1 / x +
693+
d1 * (std::log(plasg) - 1.0e0_rt) - e1);
694+
}
692695

693696
y = avo_eos*ytot1*kt*(a1 + 0.25e0_rt/plasg*(b1*x - c1/x));
694697
decouldd = y * plasgdd;
695698
decouldt = y * plasgdt + ecoul/state.T;
696699

697-
decoulda = y * plasgda - ecoul/state.abar;
700+
decoulda = y * plasgda - ecoul * ytot1;
698701
decouldz = y * plasgdz;
699702

700703
y = onethird * state.rho;
@@ -704,10 +707,12 @@ void apply_coulomb_corrections (T& state)
704707
dpcoulda = y * decoulda;
705708
dpcouldz = y * decouldz;
706709

707-
y = -avo_eos * kerg / (state.abar * plasg) *
708-
(0.75e0_rt * b1 * x + 1.25e0_rt * c1 / x + d1);
709-
dscouldd = y * plasgdd;
710-
dscouldt = y * plasgdt;
710+
if constexpr (has_entropy<T>::value) {
711+
y = -avo_eos * kerg / (state.abar * plasg) *
712+
(0.75e0_rt * b1 * x + 1.25e0_rt * c1 / x + d1);
713+
dscouldd = y * plasgdd;
714+
dscouldt = y * plasgdt;
715+
}
711716

712717
// yakovlev & shalybkov 1989 equations 102, 103, 104
713718
}
@@ -720,12 +725,14 @@ void apply_coulomb_corrections (T& state)
720725
amrex::Real dpionda = dxnida * kt;
721726
amrex::Real dpiondz = 0.0e0_rt;
722727

723-
x = plasg * std::sqrt(plasg);
724-
y = std::pow(plasg, b2);
728+
amrex::Real x = plasg * std::sqrt(plasg);
729+
amrex::Real y = std::pow(plasg, b2);
725730
z = c2 * x - onethird * a2 * y;
726731
pcoul = -pion * z;
727732
ecoul = 3.0e0_rt * pcoul / state.rho;
728-
scoul = -avo_eos / state.abar * kerg * (c2 * x -a2 * (b2 - 1.0e0_rt) / b2 * y);
733+
if constexpr (has_entropy<T>::value) {
734+
scoul = -avo_eos * ytot1 * kerg * (c2 * x -a2 * (b2 - 1.0e0_rt) / b2 * y);
735+
}
729736

730737
s = 1.5e0_rt * c2 * x / plasg - onethird * a2 * b2 * y / plasg;
731738
dpcouldd = -dpiondd * z - pion * s * plasgdd;
@@ -739,10 +746,12 @@ void apply_coulomb_corrections (T& state)
739746
decoulda = s * dpcoulda;
740747
decouldz = s * dpcouldz;
741748

742-
s = -avo_eos * kerg / (state.abar * plasg) *
743-
(1.5e0_rt * c2 * x - a2 * (b2 - 1.0e0_rt) * y);
744-
dscouldd = s * plasgdd;
745-
dscouldt = s * plasgdt;
749+
if constexpr (has_entropy<T>::value) {
750+
s = -avo_eos * kerg / (state.abar * plasg) *
751+
(1.5e0_rt * c2 * x - a2 * (b2 - 1.0e0_rt) * y);
752+
dscouldd = s * plasgdd;
753+
dscouldt = s * plasgdt;
754+
}
746755
}
747756

748757
// Disable Coulomb corrections if they cause
@@ -766,9 +775,11 @@ void apply_coulomb_corrections (T& state)
766775
ecoul = 0.0e0_rt;
767776
decouldd = 0.0e0_rt;
768777
decouldt = 0.0e0_rt;
769-
scoul = 0.0e0_rt;
770-
dscouldd = 0.0e0_rt;
771-
dscouldt = 0.0e0_rt;
778+
if constexpr (has_entropy<T>::value) {
779+
scoul = 0.0e0_rt;
780+
dscouldd = 0.0e0_rt;
781+
dscouldt = 0.0e0_rt;
782+
}
772783

773784
dpcoulda = 0.0e0_rt;
774785
dpcouldz = 0.0e0_rt;

0 commit comments

Comments
 (0)