Skip to content

Commit 019405e

Browse files
bgclayton-lanljonahm-LANLjhp-lanl
authored
Fix for Carnahan-Starling EOS with DensityFromPressureTemperature (#563)
* Updated the DensityFromPressureTemperature function for Carn-Star EOS * Cleaned up DensityFromPressureTemperature, added comments, and updated the CHANGELOG * Didn't mean to touch the CMakeLists * formatting * Tighten tolerance on root finder * loosen tolerances on test. return to 1e-12 on regula falsi * Enforce bb > 0 strictly * Do not set bb to too small of a number * Apply suggestions from code review Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> * Updated the error tolerences with a reference value. * Cleaned up the testing some more. * Reduced a tolerance for one test. * Reduced the tolerances some. * Fixed a ill-defined recursive definition of * Fixed another bug with the use of . * Reduced the tolerance since the CI was failing. * format tests --------- Co-authored-by: Jonah Miller <jonahm@lanl.gov> Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com>
1 parent 2b7ef1b commit 019405e

File tree

3 files changed

+82
-38
lines changed

3 files changed

+82
-38
lines changed

‎CHANGELOG.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77

88
### Fixed (Repair bugs, etc)
99
- [[PR561]](https://github.com/lanl/singularity-eos/pull/561) Fix logic for kokkos-kernels in spackage so that it is only required for closure models on GPU
10+
- [[PR563]](https://github.com/lanl/singularity-eos/pull/563) Fixed DensityFromPressureTemperature for the Carnahan-Starling EOS.
1011

1112
### Changed (changing behavior/API/variables/...)
1213

‎singularity-eos/eos/eos_carnahan_starling.hpp‎

Lines changed: 44 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -71,9 +71,9 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
7171
return std::max(robust::SMALL(), (sie - _qq) / _Cv);
7272
}
7373
PORTABLE_INLINE_FUNCTION void CheckParams() const {
74-
PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
75-
PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
76-
PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
74+
PORTABLE_ALWAYS_REQUIRE(_Cv > 0, "Heat capacity must be strictly positive");
75+
PORTABLE_ALWAYS_REQUIRE(_gm1 > 0, "Gruneisen parameter must be positive");
76+
PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be strictly non-negative");
7777
_AZbar.CheckParams();
7878
}
7979
template <typename Indexer_t = Real *>
@@ -94,24 +94,51 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
9494
PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature(
9595
const Real press, const Real temperature, const Real guess = robust::SMALL(),
9696
Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
97-
Real real_root;
98-
auto poly = [=](Real dens) {
99-
return _Cv * temperature * _gm1 * dens * ZedFromDensity(dens);
97+
static constexpr Real rho_low = 0.;
98+
// Need to identify an appropriate density scaling for the tolerances.
99+
// Note: cannot use `_rho0` since `_rho0` is computed from this function itself.
100+
const Real xtol = 1.0e-9;
101+
const Real ytol = 1.0e-9;
102+
103+
// At absolute zero, P = 0, so assumed vacuum state (no cold curve exists).
104+
if (temperature <= std::numeric_limits<Real>::min()) {
105+
return rho_low;
106+
}
107+
108+
// Use ideal gas value when _bb is zero
109+
if (_bb <= std::numeric_limits<Real>::min()) {
110+
return robust::ratio(press, _Cv * _gm1 * temperature);
111+
}
112+
const Real rho_high = 1.0 / _bb;
113+
114+
// Setup lambda function for finding rho.
115+
// The equation has been rewritten to avoid division by zero (polynomial equation).
116+
auto f = [=](const Real x /* density */) {
117+
const Real eta = _bb * x;
118+
const Real term1 = x * (1. + eta + eta * eta - eta * eta * eta);
119+
const Real term2 =
120+
robust::ratio(press * math_utils::pow<3>(1. - eta), _Cv * _gm1 * temperature);
121+
return term1 - term2;
100122
};
101-
using RootFinding1D::findRoot;
102-
using RootFinding1D::Status;
103-
static constexpr Real xtol = 1.0e-12;
104-
static constexpr Real ytol = 1.0e-12;
105-
static constexpr Real lo_bound = robust::SMALL();
106-
const Real hi_bound = robust::ratio(1.0, _bb);
107-
auto status = findRoot(poly, press, guess, lo_bound, hi_bound, xtol, ytol, real_root);
108-
if (status != Status::SUCCESS) {
123+
124+
const RootFinding1D::RootCounts root_info;
125+
Real rho; // `rho` will be the root of the equation $f=0$.
126+
127+
/* The regula falsi method should always return a unique real root in the interval (0,
128+
* b^{-1}) since f(0) < 0 and f(b^{-1}) > 0. It can be shown that the derivative of f
129+
* (as a function of \eta) is always positive since \eta \in (0,1) hence the root is
130+
* unique. */
131+
auto status =
132+
regula_falsi(f, 0.0 /*target*/, guess, rho_low /*left bracket*/,
133+
rho_high /*right bracket*/, xtol, ytol, rho, &root_info, true);
134+
135+
if (status != RootFinding1D::Status::SUCCESS) {
109136
// Root finder failed even though the solution was bracketed... this is an error
110137
EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in "
111-
"Carnahan-Starling EoS (root finder util) ***\n");
112-
real_root = -1.0;
138+
"Carnahan-Starling EoS (root finder `regula_falsi`) ***\n");
139+
rho = -1.0; // guarantees zero density is returned on convergence failure.
113140
}
114-
return std::max(robust::SMALL(), real_root);
141+
return std::max(rho_low, rho);
115142
}
116143
template <typename Indexer_t = Real *>
117144
PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(

‎test/test_eos_carnahan_starling.cpp‎

Lines changed: 37 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,13 @@
2323
#endif
2424

2525
#include <singularity-eos/base/constants.hpp>
26+
#include <singularity-eos/base/robust_utils.hpp>
2627
#include <singularity-eos/eos/eos.hpp>
2728
#include <test/eos_unit_test_helpers.hpp>
2829

2930
using singularity::CarnahanStarling;
3031
using singularity::IdealGas;
32+
namespace robust = singularity::robust;
3133
using EOS = singularity::Variant<IdealGas, CarnahanStarling>;
3234

3335
SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
@@ -39,6 +41,8 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
3941
// Create the EOS
4042
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq);
4143
EOS eos = host_eos.GetOnDevice();
44+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
45+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
4246
GIVEN("Densities and energies") {
4347
constexpr int num = 4;
4448
#ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -121,7 +125,7 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
121125
#endif // PORTABILITY_STRATEGY_KOKKOS
122126
THEN("The returned P(rho, e) should be equal to the true value") {
123127
array_compare(num, density, energy, h_pressure, pressure_true, "Density",
124-
"Energy");
128+
"Energy", 1e-12 * P_0);
125129
}
126130
}
127131

@@ -133,7 +137,7 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
133137
#endif // PORTABILITY_STRATEGY_KOKKOS
134138
THEN("The returned B_S(rho, e) should be equal to the true value") {
135139
array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density",
136-
"Energy");
140+
"Energy", 1e-12 * bmod_0);
137141
}
138142
}
139143

@@ -145,7 +149,7 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
145149
#endif // PORTABILITY_STRATEGY_KOKKOS
146150
THEN("The returned B_S(rho, e) should be equal to the true value") {
147151
array_compare(num, density, energy, h_temperature, temperature_true, "Density",
148-
"Energy");
152+
"Energy", 1e-12 * T_0);
149153
}
150154
}
151155

@@ -158,7 +162,7 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
158162
#endif // PORTABILITY_STRATEGY_KOKKOS
159163
THEN("The returned Gamma(rho, e) should be equal to the true value") {
160164
array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density",
161-
"Energy");
165+
"Energy", 1e-12 * dpde_0 / rho_0);
162166
}
163167
}
164168

@@ -189,6 +193,8 @@ SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") {
189193
// Create the EOS
190194
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq, qp, T0, P0);
191195
EOS eos = host_eos.GetOnDevice();
196+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
197+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
192198
GIVEN("Densities and energies") {
193199
constexpr int num = 4;
194200
#ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -271,7 +277,7 @@ SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") {
271277
#endif // PORTABILITY_STRATEGY_KOKKOS
272278
THEN("The returned P(rho, e) should be equal to the true value") {
273279
array_compare(num, density, energy, h_pressure, pressure_true, "Density",
274-
"Energy");
280+
"Energy", 1e-12 * P_0);
275281
}
276282
}
277283

@@ -283,7 +289,7 @@ SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") {
283289
#endif // PORTABILITY_STRATEGY_KOKKOS
284290
THEN("The returned B_S(rho, e) should be equal to the true value") {
285291
array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density",
286-
"Energy");
292+
"Energy", 1e-12 * bmod_0);
287293
}
288294
}
289295

@@ -295,7 +301,7 @@ SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") {
295301
#endif // PORTABILITY_STRATEGY_KOKKOS
296302
THEN("The returned B_S(rho, e) should be equal to the true value") {
297303
array_compare(num, density, energy, h_temperature, temperature_true, "Density",
298-
"Energy");
304+
"Energy", 1e-12 * T_0);
299305
}
300306
}
301307

@@ -308,7 +314,7 @@ SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") {
308314
#endif // PORTABILITY_STRATEGY_KOKKOS
309315
THEN("The returned Gamma(rho, e) should be equal to the true value") {
310316
array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density",
311-
"Energy");
317+
"Energy", 1e-12 * dpde_0 / rho_0);
312318
}
313319
}
314320
}
@@ -325,6 +331,10 @@ SCENARIO("(C-S EoS) Isentropic Bulk Modulus Analytic vs. FD",
325331
// Create the EOS
326332
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq);
327333
EOS eos = host_eos.GetOnDevice();
334+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
335+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
336+
const Real c_0 = std::sqrt(bmod_0 / rho_0);
337+
328338
GIVEN("Density and energy") {
329339
constexpr Real density = 1.1833012071291069e-03;
330340
constexpr Real energy = 2.1407169999999998e+09;
@@ -337,7 +347,7 @@ SCENARIO("(C-S EoS) Isentropic Bulk Modulus Analytic vs. FD",
337347
INFO("Density: " << density << " Energy: " << energy
338348
<< " Sound speed: " << sound_speed
339349
<< " cm/s True sound speed: " << true_sound_speed << " cm/s");
340-
REQUIRE(isClose(sound_speed, true_sound_speed, 1e-12));
350+
REQUIRE(isClose(sound_speed, true_sound_speed, 1e-12 * c_0));
341351
}
342352
}
343353
WHEN("A finite difference approximation is used for the bulk modulus") {
@@ -359,7 +369,7 @@ SCENARIO("(C-S EoS) Isentropic Bulk Modulus Analytic vs. FD",
359369
INFO("Density: " << density << " Energy: " << energy
360370
<< " Sound speed: " << sound_speed
361371
<< " cm/s Approximate sound speed: " << ss_approx << " cm/s");
362-
REQUIRE(isClose(sound_speed, ss_approx, 1e-5));
372+
REQUIRE(isClose(sound_speed, ss_approx, 1e-5 * c_0));
363373
}
364374
}
365375
}
@@ -370,12 +380,14 @@ SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]")
370380
GIVEN("Parameters for a CarnahanStarling EOS") {
371381
constexpr Real gm1 = 0.4;
372382
constexpr Real Cv = 7180000.0;
373-
constexpr Real bb = 0;
383+
constexpr Real bb = robust::EPS();
374384
constexpr Real qq = 0;
375385
// Create the EOS
376386
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq);
377387
EOS eos = host_eos.GetOnDevice();
378388
EOS ideal_eos = singularity::IdealGas(gm1, Cv);
389+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
390+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
379391
GIVEN("Densities and energies") {
380392
constexpr int num = 1;
381393
#ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -448,7 +460,7 @@ SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]")
448460
ideal_eos.PressureFromDensityInternalEnergy(density[i], energy[i]);
449461
}
450462
array_compare(num, density, energy, h_pressure, pressure_true, "Density",
451-
"Energy");
463+
"Energy", 1e-12 * P_0);
452464
}
453465
}
454466

@@ -464,7 +476,7 @@ SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]")
464476
ideal_eos.BulkModulusFromDensityInternalEnergy(density[i], energy[i]);
465477
}
466478
array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density",
467-
"Energy");
479+
"Energy", 1e-12 * bmod_0);
468480
}
469481
}
470482

@@ -480,7 +492,7 @@ SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]")
480492
ideal_eos.TemperatureFromDensityInternalEnergy(density[i], energy[i]);
481493
}
482494
array_compare(num, density, energy, h_temperature, temperature_true, "Density",
483-
"Energy");
495+
"Energy", 1e-12 * T_0);
484496
}
485497
}
486498

@@ -497,7 +509,7 @@ SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]")
497509
ideal_eos.GruneisenParamFromDensityInternalEnergy(density[i], energy[i]);
498510
}
499511
array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density",
500-
"Energy");
512+
"Energy", 1e-12 * dpde_0 / rho_0);
501513
}
502514
}
503515
}
@@ -516,6 +528,8 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") {
516528
// Create the EOS
517529
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq, qp, T0, P0);
518530
EOS eos = host_eos.GetOnDevice();
531+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
532+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
519533
GIVEN("Densities and energies") {
520534
constexpr int num = 1;
521535
#ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -549,7 +563,7 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") {
549563
// Gold standard values for a subset of lookups
550564
// constexpr std::array<Real, num> temperature_true{4.0000000000000000e+03};
551565
constexpr std::array<Real, num> entropy_true{-2.2983150752058342e+10};
552-
566+
constexpr Real entropy_ref = 2.2983150752058342e+10;
553567
#ifdef PORTABILITY_STRATEGY_KOKKOS
554568
// Create device views for outputs and mirror those views on the host
555569
Kokkos::View<Real[num]> v_temperature("Temperature");
@@ -582,7 +596,7 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") {
582596
#endif // PORTABILITY_STRATEGY_KOKKOS
583597
THEN("The returned S(rho, e) should be equal to the true value") {
584598
array_compare(num, density, energy, h_entropy, entropy_true, "Density",
585-
"Energy");
599+
"Energy", 1e-12 * entropy_ref);
586600
}
587601
}
588602
WHEN("A S(rho, T(rho,e)) lookup is performed") {
@@ -594,7 +608,7 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") {
594608
#endif // PORTABILITY_STRATEGY_KOKKOS
595609
THEN("The returned S(rho, e) should be equal to the true value") {
596610
array_compare(num, density, energy, h_entropy, entropy_true, "Density",
597-
"Energy");
611+
"Energy", 1e-12 * entropy_ref);
598612
}
599613
}
600614
}
@@ -610,6 +624,8 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") {
610624
// Create the EOS
611625
EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq);
612626
EOS eos = host_eos.GetOnDevice();
627+
Real rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0;
628+
host_eos.ValuesAtReferenceState(rho_0, T_0, e_0, P_0, cv_0, bmod_0, dpde_0, dvdt_0);
613629
GIVEN("Pressure and temperature") {
614630
constexpr int num = 4;
615631
#ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -648,8 +664,8 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") {
648664

649665
// Gold standard values for a subset of lookups
650666
constexpr std::array<Real, num> density_true{
651-
1.1833012071291069e-03, 7.8290736890381501e-03, 8.5453943327882340e-03,
652-
8.8197619601121831e-03};
667+
1.18330186346892515e-03, 7.82907368898192423e-03, 8.54539433271512064e-03,
668+
8.81976196003179773e-03};
653669

654670
#ifdef PORTABILITY_STRATEGY_KOKKOS
655671
// Create device views for outputs and mirror those views on the host
@@ -685,7 +701,7 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") {
685701

686702
THEN("The returned rho(P, T) should be equal to the true value") {
687703
array_compare(num, pressure, temperature, h_density, density_true, "Pressure",
688-
"Temperature");
704+
"Temperature", 1e-7 * rho_0);
689705
}
690706
}
691707
}

0 commit comments

Comments
 (0)