From 5bf0debe187bcd8a130d0ad7019ca23f02e08037 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 8 Apr 2024 10:51:57 -0600 Subject: [PATCH 1/9] Add robust::ratio and asserts to PTE solvers --- singularity-eos/closure/mixed_cell_models.hpp | 249 ++++++++++-------- 1 file changed, 136 insertions(+), 113 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 365e1dd4a6..72ad2d3bfe 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -16,7 +16,9 @@ #define _SINGULARITY_EOS_CLOSURE_MIXED_CELL_MODELS_ #include +#include #include +#include #include #include @@ -52,6 +54,9 @@ constexpr Real line_search_fac = 0.5; constexpr Real vfrac_safety_fac = 0.95; constexpr Real minimum_temperature = 1.e-9; constexpr Real maximum_temperature = 1.e9; +constexpr Real temperature_limit = 1.0e15; +constexpr Real default_tguess = 300.; +constexpr Real min_dtde = 1.0e-16; } // namespace mix_params namespace mix_impl { @@ -167,8 +172,9 @@ class PTESolverBase { Real vsum = 0; for (int m = 0; m < nmat; ++m) vsum += vfrac[m]; + PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive") for (int m = 0; m < nmat; ++m) - vfrac[m] *= vfrac_total / vsum; + vfrac[m] *= robust::ratio(vfrac_total, vsum); } // Finalize restores the temperatures, energies, and pressures to unscaled values from // the internally scaled quantities used by the solvers @@ -212,6 +218,7 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { + PORTABLE_REQUIRE(vfrac[m] >= 0., "Negative volume fraction") rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -223,13 +230,13 @@ class PTESolverBase { // set volume fractions for (int m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(T); - const Real vmax = std::min(0.9 * rhobar[m] / rho_min, 1.0); + const Real vmax = std::min(0.9 * robust::ratio(rhobar[m], rho_min), 1.0); vfrac[m] = (vfrac[m] > 0.0 ? std::min(vmax, vfrac[m]) : vmax); vsum += vfrac[m]; } // Normalize vfrac for (int m = 0; m < nmat; ++m) { - vfrac[m] *= vfrac_total / vsum; + vfrac[m] *= robust::ratio(vfrac_total, vsum); } } @@ -241,12 +248,14 @@ class PTESolverBase { if (Tnorm > 0.0) { Tguess = Tnorm; } else { - Tguess = 300.0; + Tguess = mix_params::default_tguess; for (int m = 0; m < nmat; ++m) Tguess = std::max(Tguess, temp[m]); } + PORTABLE_REQUIRE(Tguess > 0., "Non-positive temperature guess for PTE") // check for sanity. basically checks that the input temperatures weren't garbage - assert(Tguess < 1.0e15); + PORTABLE_REQUIRE(Tguess < mix_params::temperature_limit, + "Very large input temperature or temperature guess"); // iteratively increase temperature guess until all rho's are above rho_at_pmin const Real Tfactor = 10.0; bool rho_fail; @@ -256,8 +265,8 @@ class PTESolverBase { rho_fail = false; for (int m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(Tguess); - rho[m] = rhobar[m] / vfrac[m]; - if (rho[m] < rho_min) { + rho[m] = robust::ratio(rhobar[m], vfrac[m]); + if (rho[m] < rho_min and rho[m] > 0.) { rho_fail = true; Tguess *= Tfactor; break; @@ -299,7 +308,7 @@ class PTESolverBase { // intialize rhobar array and final density InitRhoBarandRho(); - uscale = rho_total * sie_total; + uscale = rho_total * abs(sie_total); // guess some non-zero temperature to start const Real Tguess = GetTguess(); @@ -311,14 +320,13 @@ class PTESolverBase { temp[m] = 1.0; sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tguess, Cache[m]); // note the scaling of pressure - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tguess, sie[m], Cache[m], - false) / - uscale; + press[m] = robust::ratio(this->GetPressureFromPreferred( + eos[m], rho[m], Tguess, sie[m], Cache[m], false), uscale); } // note the scaling of the material internal energy densities for (int m = 0; m < nmat; ++m) - u[m] = sie[m] * rhobar[m] / uscale; + u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } PORTABLE_INLINE_FUNCTION @@ -339,13 +347,16 @@ class PTESolverBase { Real rhoBsum = 0.0; Real Asum = 0.0; for (int m = 0; m < nmat; m++) { - Asum += vfrac[m] * press[m] / temp[m]; - rhoBsum += rho[m] * vfrac[m] * sie[m] / temp[m]; + Asum += vfrac[m] * robust::ratio(press[m], temp[m]); + rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } + PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess") Asum *= uscale / Tnorm; rhoBsum /= Tnorm; + PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density") Tideal = uscale / rhoBsum / Tnorm; - Pideal = Tnorm * Tideal * Asum / vfrac_total / uscale; + PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum") + Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } // Compute the ideal EOS PTE solution and replace the initial guess if it has a lower @@ -374,10 +385,10 @@ class PTESolverBase { res_norm_old += res[m] * res[m]; } // check if the volume fractions are reasonable - const Real alpha = Pideal / Tideal; + const Real alpha = robust::ratio(Pideal, Tideal); for (int m = 0; m < nmat; ++m) { - vfrac[m] *= press[m] / (temp[m] * alpha); - if (rhobar[m] / vfrac[m] < eos[m].RhoPmin(Tnorm * Tideal)) { + vfrac[m] *= robust::ratio(press[m], (temp[m] * alpha)); + if (robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { // abort because this is putting this material into a bad state for (int n = m; n >= 0; n--) vfrac[n] = vtemp[n]; @@ -386,13 +397,14 @@ class PTESolverBase { } // fill in the rest of the state for (int m = 0; m < nmat; ++m) { - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); + const Real sie_m = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * Tideal, Cache[m]); - u[m] = rhobar[m] * sie_m / uscale; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tideal, sie_m, - Cache[m], false) / - uscale; + u[m] = rhobar[m] * robust::ratio(sie_m, uscale); + press[m] = robust::ratio(this->GetPressureFromPreferred( + eos[m], rho[m], Tnorm * Tideal, sie_m, Cache[m], false), + uscale); } // fill in the residual solver->Residual(); @@ -415,7 +427,7 @@ class PTESolverBase { // did work, fill in temp and energy density for (int m = 0; m < nmat; ++m) { temp[m] = Tideal; - sie[m] = uscale * u[m] / rhobar[m]; + sie[m] = uscale * robust::ratio(u[m], rhobar[m]); } } } @@ -493,7 +505,7 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU( iter++; } - const Real alpha = (u_tot - ulo) / (uhi - ulo); + const Real alpha = robust::ratio((u_tot - ulo), (uhi - ulo)); return FastMath::pow2((1.0 - alpha) * lTlo + alpha * lThi); } @@ -611,27 +623,29 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert{}; Real e_pert = eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]); - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, - Cache[m], false) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; - dedv[m] = (rhobar[m] * e_pert / uscale - u[m]) / dv; + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, Cache[m], + false), + uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); + dedv[m] = robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dv); ////////////////////////////// // perturb temperature ////////////////////////////// Real dT = Tequil * derivative_eps; e_pert = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * (Tequil + dT), Cache[m]); - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), - e_pert, Cache[m], false) / - uscale; - dpdT[m] = (p_pert - press[m]) / dT; - dedT_sum += (rhobar[m] * e_pert / uscale - u[m]) / dT; + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), e_pert, + Cache[m], false), + uscale); + dpdT[m] = robust::ratio((p_pert - press[m]), dT); + dedT_sum += robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dT); } // Fill in the Jacobian @@ -657,7 +671,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -665,14 +679,14 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; m++) { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = -0.95 * Tequil / dx[nmat]; + scale = robust::ratio(-0.95 * Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -692,15 +706,16 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { Tequil = Ttemp + scale * dx[nmat]; for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( rho[m], Tnorm * Tequil, Cache[m]); - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); temp[m] = Tequil; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, sie[m], - Cache[m], false) / - uscale; + press[m] = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, sie[m], Cache[m], + false), + uscale); } Residual(); return ResidualNorm(); @@ -784,15 +799,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // volume fractions have been potentially reset to ensure densitites are // larger than rho(Pmin(Tequil)); set the physical density to reflect // this change in volume fraction - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tequil, Cache[m]); uscale += sie[m] * rho[m]; // note the scaling of pressure press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]); } for (int m = 0; m < nmat; ++m) { - press[m] /= uscale; - u[m] = sie[m] * rhobar[m] / uscale; + press[m] = robust::ratio(pres[m], uscale); + u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } Residual(); return ResidualNorm(); @@ -838,11 +853,11 @@ class PTESolverFixedT : public mix_impl::PTESolverBase ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); - Real p_pert = - eos[m].PressureFromDensityTemperature(rho_pert, Tequil, Cache[m]) / uscale; - dpdv[m] = (p_pert - press[m]) / dv; + Real p_pert = robust::ratio( + eos[m].PressureFromDensityTemperature(rho_pert, Tequil, Cache[m]), uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); } // Fill in the Jacobian @@ -864,15 +879,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m]); } } // control how big of a step toward rho = rho(Pmin) is allowed for (int m = 0; m < nmat; m++) { const Real rho_min = eos[m].RhoPmin(Tequil); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } } // Now apply the overall scaling @@ -891,12 +906,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase } for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature(rho[m], Tequil, Cache[m]); - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]) / uscale; + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]), uscale); } Residual(); return ResidualNorm(); @@ -990,8 +1006,9 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // note the scaling of the material internal energy densities for (int m = 0; m < nmat; ++m) { - u[m] = sie[m] * rhobar[m] / uscale; - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]) / uscale; + u[m] = robust::ratio(sie[m] * rhobar[m], uscale); + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]), uscale); } Residual(); // Set the current guess for the equilibrium temperature. Note that this is already @@ -1006,7 +1023,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Real vsum = 0.0; for (int m = 0; m < nmat; ++m) { vsum += vfrac[m]; - residual[m] = Pequil / uscale - press[m]; + residual[m] = robust::ratio(Pequil, uscale) - press[m]; } residual[nmat] = vfrac_total - vsum; } @@ -1036,23 +1053,25 @@ class PTESolverFixedP : public mix_impl::PTESolverBase ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert{}; Real e_pert{}; - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, - Cache[m], true) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, Cache[m], + true), + uscale); + dpdv[m] = robust::ratio((p_pert - press[m]), dv); ////////////////////////////// // perturb temperature ////////////////////////////// Real dT = Tequil * derivative_eps; - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), - e_pert, Cache[m], true) / - uscale; - dpdT[m] = (p_pert - press[m]) / dT; + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), e_pert, + Cache[m], true), + uscale); + dpdT[m] = robust::ratio((p_pert - press[m]), dT); } // Fill in the Jacobian @@ -1074,7 +1093,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * vfrac[m] / dx[m]; + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -1082,14 +1101,14 @@ class PTESolverFixedP : public mix_impl::PTESolverBase for (int m = 0; m < nmat; m++) { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = -0.95 * Tequil / dx[nmat]; + scale = -0.95 * robust::ratio(Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -1109,13 +1128,14 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Tequil = Ttemp + scale * dx[nmat]; for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( rho[m], Tnorm * Tequil, Cache[m]); - press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tnorm * Tequil, Cache[m]) / - uscale; - sie[m] = u[m] / rhobar[m]; - u[m] /= uscale; + press[m] = robust::ratio( + eos[m].PressureFromDensityTemperature(rho[m], Tnorm * Tequil, Cache[m]), + uscale); + sie[m] = robust::ratio(u[m], rhobar[m]); + u[m] = robust::ratio(u[m], uscale); temp[m] = Tequil; } Residual(); @@ -1228,7 +1248,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { error_p += residual[m + 1] * residual[m + 1]; error_t += residual[m + nmat] * residual[m + nmat]; } - mean_t /= rho_total; + mean_t = robust::ratio(mean_t, rho_total); error_p = std::sqrt(error_p); error_t = std::sqrt(error_t); // Check for convergence @@ -1248,31 +1268,33 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { ////////////////////////////// Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; const Real vf_pert = vfrac[m] + dv; - const Real rho_pert = rhobar[m] / vf_pert; + const Real rho_pert = robust::ratio(rhobar[m], vf_pert); Real p_pert; - Real t_pert = - eos[m].TemperatureFromDensityInternalEnergy(rho_pert, sie[m], Cache[m]) / Tnorm; - p_pert = this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, sie[m], - Cache[m], false) / - uscale; - dpdv[m] = (p_pert - press[m]) / dv; - dtdv[m] = (t_pert - temp[m]) / dv; + Real t_pert = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho_pert, sie[m], Cache[m]), Tnorm); + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, sie[m], Cache[m], + false), + uscale); + dpdv[m] = robust::ratio(p_pert - press[m], dv); + dtdv[m] = robust::ratio(t_pert - temp[m], dv); ////////////////////////////// // perturb energies ////////////////////////////// const Real de = std::abs(u[m]) * derivative_eps; - Real e_pert = (u[m] + de) / rhobar[m]; - - t_pert = - eos[m].TemperatureFromDensityInternalEnergy(rho[m], uscale * e_pert, Cache[m]) / - Tnorm; - p_pert = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, - uscale * e_pert, Cache[m], false) / - uscale; - dpde[m] = (p_pert - press[m]) / de; - dtde[m] = (t_pert - temp[m]) / de; - if (std::abs(dtde[m]) < 1.e-16) { // must be on the cold curve + Real e_pert = robust::ratio(u[m] + de, rhobar[m]); + + t_pert = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho[m], uscale * e_pert, Cache[m]), + Tnorm); + p_pert = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, uscale * e_pert, + Cache[m], false), + uscale); + dpde[m] = robust::ratio(p_pert - press[m], de); + dtde[m] = robust::ratio(t_pert - temp[m], de); + if (std::abs(dtde[m]) < mix_impl::min_dtde) { // must be on the cold curve dtde[m] = derivative_eps; } } @@ -1310,20 +1332,20 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed if (scale * dx[m] < -0.1 * vfrac[m]) { - scale = -0.1 * vfrac[m] / dx[m]; + scale = -0.1 * robust::ratio(vfrac[m], dx[m]); } // try to control steps toward T = 0 const Real dt = (dtdv[m] * dx[m] + dtde[m] * dx[m + nmat]); if (scale * dt < -0.1 * temp[m]) { - scale = -0.1 * temp[m] / dt; + scale = -0.1 * robust::ratio(temp[m], dt); } const Real tt = temp[m] + scale * dt; const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * temp[m]), eos[m].RhoPmin(Tnorm * tt)); - const Real alpha_max = rhobar[m] / rho_min; + const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * (alpha_max - vfrac[m]) / dx[m]; + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // Now apply the overall scaling @@ -1344,14 +1366,15 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { } for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; - rho[m] = rhobar[m] / vfrac[m]; + rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = utemp[m] + scale * dx[nmat + m]; - sie[m] = uscale * u[m] / rhobar[m]; - temp[m] = - eos[m].TemperatureFromDensityInternalEnergy(rho[m], sie[m], Cache[m]) / Tnorm; - press[m] = this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], sie[m], - Cache[m], false) / - uscale; + sie[m] = uscale * robust::ratio(u[m], rhobar[m]); + temp[m] = robust::ratio( + eos[m].TemperatureFromDensityInternalEnergy(rho[m], sie[m], Cache[m]), Tnorm); + press[m] = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], sie[m], + Cache[m], false), + uscale); } Residual(); return ResidualNorm(); @@ -1401,7 +1424,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // backtrack Real err_mid = s.TestUpdate(0.5); if (err_mid < err && err_mid < err_old) { - scale = 0.75 + 0.5 * (err_mid - err) / (err - 2.0 * err_mid + err_old); + scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old); } else { scale = line_search_fac; } From 504cfb0e9e3fb0513017b79201383c7c351be737 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 17:14:08 -0600 Subject: [PATCH 2/9] Update CMake to make Kokkos a common library and use a PUBLIC scope to link --- CMakeLists.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b36680e5ec..c347aa8bbe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -167,11 +167,11 @@ if(SINGULARITY_BUILD_PYTHON) set(CMAKE_POSITION_INDEPENDENT_CODE ON) endif() +target_compile_features(singularity-eos_Common INTERFACE cxx_std_14) + # checks if this is our build, or we've been imported via `add_subdirectory` NB: # this should make the `option(SINGULARITY_SUBMODULE_MODE ...)` unnecessary if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) - set(CMAKE_CXX_STANDARD 14) - set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) else() message( @@ -335,9 +335,9 @@ if(SINGULARITY_USE_SPINER) endif() if(SINGULARITY_USE_KOKKOS) - singularity_enable_kokkos(singularity-eos_Interface) + singularity_enable_kokkos(singularity-eos_Common) if(SINGULARITY_USE_KOKKOSKERNELS) - singularity_enable_kokkoskernels(singularity-eos_Interface) + singularity_enable_kokkoskernels(singularity-eos_Common) endif() endif() if(SINGULARITY_USE_EIGEN) @@ -402,7 +402,7 @@ if(SINGULARITY_BUILD_CLOSURE) add_library(singularity-eos_Library) set_target_properties(singularity-eos_Library PROPERTIES OUTPUT_NAME singularity-eos) add_library(singularity-eos::singularity-eos_Library ALIAS singularity-eos_Library) - target_link_libraries(singularity-eos_Library INTERFACE singularity-eos_Common) + target_link_libraries(singularity-eos_Library PUBLIC singularity-eos_Common) target_link_libraries(singularity-eos INTERFACE singularity-eos_Library) target_sources(singularity-eos_Library PRIVATE ${eos_srcs}) message(VERBOSE "EOS Library Sources:\n\t${eos_srcs}") From bbc44da85493831cb4de8b5572a3a51bb48019f5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 17:14:52 -0600 Subject: [PATCH 3/9] Fix typos and small bugs --- singularity-eos/closure/mixed_cell_models.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 72ad2d3bfe..9050df870e 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -172,7 +172,7 @@ class PTESolverBase { Real vsum = 0; for (int m = 0; m < nmat; ++m) vsum += vfrac[m]; - PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive") + PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive"); for (int m = 0; m < nmat; ++m) vfrac[m] *= robust::ratio(vfrac_total, vsum); } @@ -218,7 +218,7 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { - PORTABLE_REQUIRE(vfrac[m] >= 0., "Negative volume fraction") + PORTABLE_REQUIRE(vfrac[m] >= 0., "Negative volume fraction"); rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -252,7 +252,7 @@ class PTESolverBase { for (int m = 0; m < nmat; ++m) Tguess = std::max(Tguess, temp[m]); } - PORTABLE_REQUIRE(Tguess > 0., "Non-positive temperature guess for PTE") + PORTABLE_REQUIRE(Tguess > 0., "Non-positive temperature guess for PTE"); // check for sanity. basically checks that the input temperatures weren't garbage PORTABLE_REQUIRE(Tguess < mix_params::temperature_limit, "Very large input temperature or temperature guess"); @@ -350,12 +350,12 @@ class PTESolverBase { Asum += vfrac[m] * robust::ratio(press[m], temp[m]); rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } - PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess") + PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess"); Asum *= uscale / Tnorm; rhoBsum /= Tnorm; - PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density") + PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density"); Tideal = uscale / rhoBsum / Tnorm; - PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum") + PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum"); Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } @@ -806,7 +806,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]); } for (int m = 0; m < nmat; ++m) { - press[m] = robust::ratio(pres[m], uscale); + press[m] = robust::ratio(press[m], uscale); u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } Residual(); @@ -1294,7 +1294,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { uscale); dpde[m] = robust::ratio(p_pert - press[m], de); dtde[m] = robust::ratio(t_pert - temp[m], de); - if (std::abs(dtde[m]) < mix_impl::min_dtde) { // must be on the cold curve + if (std::abs(dtde[m]) < mix_params::min_dtde) { // must be on the cold curve dtde[m] = derivative_eps; } } From 44aea1daf829b08fd94c3b4ec92a274c6f02cda6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 17:35:41 -0600 Subject: [PATCH 4/9] Update copyright --- CMakeLists.txt | 2 +- singularity-eos/closure/mixed_cell_models.hpp | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c347aa8bbe..a76ca42f3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ # ------------------------------------------------------------------------------# -# © 2021-2023. Triad National Security, LLC. All rights reserved. This program +# © 2021-2024. Triad National Security, LLC. All rights reserved. This program # was produced under U.S. Government contract 89233218CNA000001 for Los Alamos # National Laboratory (LANL), which is operated by Triad National Security, LLC # for the U.S. Department of Energy/National Nuclear Security Administration. diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 9050df870e..fee315afb6 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -218,7 +218,8 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { - PORTABLE_REQUIRE(vfrac[m] >= 0., "Negative volume fraction"); + PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction"); + PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density"); rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -266,7 +267,7 @@ class PTESolverBase { for (int m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(Tguess); rho[m] = robust::ratio(rhobar[m], vfrac[m]); - if (rho[m] < rho_min and rho[m] > 0.) { + if (rho[m] < rho_min) { rho_fail = true; Tguess *= Tfactor; break; From 22bb579284e056fae09411f8a250fff9838b3e22 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 17:35:59 -0600 Subject: [PATCH 5/9] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6232dfc1d7..23633b1753 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - [[PR341]](https://github.com/lanl/singularity-eos/pull/341) Short-circuit HDF5 machinery when cray-wrappers used in-tree - [[PR340]](https://github.com/lanl/singularity-eos/pull/335) Fix in-tree builds with plugin infrastructure - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS +- [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver ### Added (new features/APIs/variables/...) - [[PR334]](https://github.com/lanl/singularity-eos/pull/334) Include plugins infrastructure From e8f37416911183ddb72fbbb12fc856e4ed9d2102 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 17:38:16 -0600 Subject: [PATCH 6/9] clang format --- singularity-eos/closure/mixed_cell_models.hpp | 80 ++++++++++--------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index fee315afb6..243c7e547f 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -321,8 +321,9 @@ class PTESolverBase { temp[m] = 1.0; sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tguess, Cache[m]); // note the scaling of pressure - press[m] = robust::ratio(this->GetPressureFromPreferred( - eos[m], rho[m], Tguess, sie[m], Cache[m], false), uscale); + press[m] = robust::ratio( + this->GetPressureFromPreferred(eos[m], rho[m], Tguess, sie[m], Cache[m], false), + uscale); } // note the scaling of the material internal energy densities @@ -403,9 +404,10 @@ class PTESolverBase { const Real sie_m = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * Tideal, Cache[m]); u[m] = rhobar[m] * robust::ratio(sie_m, uscale); - press[m] = robust::ratio(this->GetPressureFromPreferred( - eos[m], rho[m], Tnorm * Tideal, sie_m, Cache[m], false), - uscale); + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tideal, + sie_m, Cache[m], false), + uscale); } // fill in the residual solver->Residual(); @@ -629,10 +631,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { Real p_pert{}; Real e_pert = eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]); - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, Cache[m], - false), - uscale); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, + e_pert, Cache[m], false), + uscale); dpdv[m] = robust::ratio((p_pert - press[m]), dv); dedv[m] = robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dv); ////////////////////////////// @@ -641,10 +643,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { Real dT = Tequil * derivative_eps; e_pert = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * (Tequil + dT), Cache[m]); - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), e_pert, - Cache[m], false), - uscale); + p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], + Tnorm * (Tequil + dT), e_pert, + Cache[m], false), + uscale); dpdT[m] = robust::ratio((p_pert - press[m]), dT); dedT_sum += robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dT); } @@ -713,10 +715,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { sie[m] = robust::ratio(u[m], rhobar[m]); u[m] = robust::ratio(u[m], uscale); temp[m] = Tequil; - press[m] = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, sie[m], Cache[m], - false), - uscale); + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil, + sie[m], Cache[m], false), + uscale); } Residual(); return ResidualNorm(); @@ -913,7 +915,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase sie[m] = robust::ratio(u[m], rhobar[m]); u[m] = robust::ratio(u[m], uscale); press[m] = robust::ratio( - eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]), uscale); + eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]), uscale); } Residual(); return ResidualNorm(); @@ -1009,7 +1011,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase for (int m = 0; m < nmat; ++m) { u[m] = robust::ratio(sie[m] * rhobar[m], uscale); press[m] = robust::ratio( - eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]), uscale); + eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]), uscale); } Residual(); // Set the current guess for the equilibrium temperature. Note that this is already @@ -1058,20 +1060,20 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Real p_pert{}; Real e_pert{}; - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, e_pert, Cache[m], - true), - uscale); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil, + e_pert, Cache[m], true), + uscale); dpdv[m] = robust::ratio((p_pert - press[m]), dv); ////////////////////////////// // perturb temperature ////////////////////////////// Real dT = Tequil * derivative_eps; - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), e_pert, - Cache[m], true), - uscale); + p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], + Tnorm * (Tequil + dT), e_pert, + Cache[m], true), + uscale); dpdT[m] = robust::ratio((p_pert - press[m]), dT); } @@ -1274,10 +1276,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { Real p_pert; Real t_pert = robust::ratio( eos[m].TemperatureFromDensityInternalEnergy(rho_pert, sie[m], Cache[m]), Tnorm); - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, sie[m], Cache[m], - false), - uscale); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * t_pert, + sie[m], Cache[m], false), + uscale); dpdv[m] = robust::ratio(p_pert - press[m], dv); dtdv[m] = robust::ratio(t_pert - temp[m], dv); ////////////////////////////// @@ -1289,10 +1291,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { t_pert = robust::ratio( eos[m].TemperatureFromDensityInternalEnergy(rho[m], uscale * e_pert, Cache[m]), Tnorm); - p_pert = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, uscale * e_pert, - Cache[m], false), - uscale); + p_pert = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * t_pert, + uscale * e_pert, Cache[m], false), + uscale); dpde[m] = robust::ratio(p_pert - press[m], de); dtde[m] = robust::ratio(t_pert - temp[m], de); if (std::abs(dtde[m]) < mix_params::min_dtde) { // must be on the cold curve @@ -1372,10 +1374,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { sie[m] = uscale * robust::ratio(u[m], rhobar[m]); temp[m] = robust::ratio( eos[m].TemperatureFromDensityInternalEnergy(rho[m], sie[m], Cache[m]), Tnorm); - press[m] = robust::ratio( - this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], sie[m], - Cache[m], false), - uscale); + press[m] = + robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * temp[m], + sie[m], Cache[m], false), + uscale); } Residual(); return ResidualNorm(); From 9fc2d0ceaad868ace0d1b1d7ea6d838187aa6848 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 9 Apr 2024 18:36:31 -0600 Subject: [PATCH 7/9] Switch C++14 feature from 'Common' to 'Interface' --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a76ca42f3e..0cd14a6e11 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -167,7 +167,7 @@ if(SINGULARITY_BUILD_PYTHON) set(CMAKE_POSITION_INDEPENDENT_CODE ON) endif() -target_compile_features(singularity-eos_Common INTERFACE cxx_std_14) +target_compile_features(singularity-eos_Interface INTERFACE cxx_std_14) # checks if this is our build, or we've been imported via `add_subdirectory` NB: # this should make the `option(SINGULARITY_SUBMODULE_MODE ...)` unnecessary From 38aaa668617f746af5d176e5414ddd928f380318 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Apr 2024 12:54:39 -0600 Subject: [PATCH 8/9] Add CMake changes and reorganize --- CHANGELOG.md | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 23633b1753..46ede2ea26 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ ### Added (new features/APIs/variables/...) - [[PR#339]](https://github.com/lanl/singularity-eos/pull/339) Added COMPONENTS to singularity-eos CMake install, allowing to select a minimal subset needed e.g. for Fortran bindings only - [[PR#336]](https://github.com/lanl/singularity-eos/pull/336) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a pressure power law expansion in eta = 1-V/V0. PowerMG. +- [[PR334]](https://github.com/lanl/singularity-eos/pull/334) Include plugins infrastructure +- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. MGUsup. +- [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release ### Fixed (Repair bugs, etc) - [[PR343]](https://github.com/lanl/singularity-eos/pull/343) Add chemical potentials to stellar collapse gold files @@ -13,11 +16,7 @@ - [[PR340]](https://github.com/lanl/singularity-eos/pull/335) Fix in-tree builds with plugin infrastructure - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - -### Added (new features/APIs/variables/...) -- [[PR334]](https://github.com/lanl/singularity-eos/pull/334) Include plugins infrastructure -- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. MGUsup. -- [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release +- [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface ### Changed (changing behavior/API/variables/...) From ba961dfd99a36e376224daeb749e2dd4b6d1de74 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Apr 2024 12:54:54 -0600 Subject: [PATCH 9/9] Add comments to clarify changes --- CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0cd14a6e11..c1cbef85c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -334,6 +334,8 @@ if(SINGULARITY_USE_SPINER) # singularity_enable_hdf5(singularity-eos_Interface) endif() endif() +# Both the interface (headers) and the library (compiled) need to link to Kokkos +# see get_sg_eos.cpp if(SINGULARITY_USE_KOKKOS) singularity_enable_kokkos(singularity-eos_Common) if(SINGULARITY_USE_KOKKOSKERNELS) @@ -402,6 +404,8 @@ if(SINGULARITY_BUILD_CLOSURE) add_library(singularity-eos_Library) set_target_properties(singularity-eos_Library PROPERTIES OUTPUT_NAME singularity-eos) add_library(singularity-eos::singularity-eos_Library ALIAS singularity-eos_Library) + # Public scope ensures explicit Kokkos dependency for library and anything + # that links to it target_link_libraries(singularity-eos_Library PUBLIC singularity-eos_Common) target_link_libraries(singularity-eos INTERFACE singularity-eos_Library) target_sources(singularity-eos_Library PRIVATE ${eos_srcs})