diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index b5d08b7d..19057102 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -68,6 +68,12 @@ std::shared_ptr Initialize(ParameterInput *pin, const Real cfl_number = pin->GetOrAddReal("gas", "cfl", 0.8); params.Add("cfl", cfl_number); + // Floors + const Real dfloor = pin->GetOrAddReal("gas", "dfloor", 1.0e-20); + const Real siefloor = pin->GetOrAddReal("gas", "siefloor", 1.0e-20); + params.Add("dfloor", dfloor); + params.Add("siefloor", siefloor); + // Equation of state std::string eos_type = "none"; if (pin->DoesBlockExist("gas/eos/ideal") || (pin->DoesParameterExist("gas", "gamma"))) { @@ -131,10 +137,12 @@ std::shared_ptr Initialize(ParameterInput *pin, const Real ldmax = pin->GetOrAddReal(block_name, "ldmax", -3); const int nd = pin->GetOrAddInteger(block_name, "nd", 100); const int nt = pin->GetOrAddInteger(block_name, "nt", 100); + ArtemisEOS::IdealHHe eos_base(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, + save_to_file, true); + eos_base.SetFloors(siefloor, 0.0, dfloor, 0.0); + EOS eos_host = singularity::UnitSystem( - ArtemisEOS::IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, - true), - singularity::eos_units_init::LengthTimeUnitsInit(), + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); EOS eos_device = eos_host.GetOnDevice(); @@ -293,12 +301,6 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("scattering_h", scattering); params.Add("scattering_d", scattering.GetOnDevice()); - // Floors - const Real dfloor = pin->GetOrAddReal("gas", "dfloor", 1.0e-20); - const Real siefloor = pin->GetOrAddReal("gas", "siefloor", 1.0e-20); - params.Add("dfloor", dfloor); - params.Add("siefloor", siefloor); - // Dual energy switch // When internal > de_switch * total we use the total // The default turns off the switch diff --git a/src/pgen/strat.hpp b/src/pgen/strat.hpp index c23dbc26..048b2fe1 100644 --- a/src/pgen/strat.hpp +++ b/src/pgen/strat.hpp @@ -116,18 +116,17 @@ Real InitialDensity(const EOS &eos, const StratParams &pars, const Real z) { if ((pars.three_d) && (std::abs(z) > 1e-16)) { const Real dz = std::abs(z) / static_cast(pars.npoints); const Real dlnr = 1e-6; - const Real ldmin = std::log(pars.dfloor / pars.rho0); Real pres = eos.PressureFromDensityTemperature(dens, pars.temp0); Real zj = 0.0; Real ld = 0.0; - dens = std::exp(ld); + dens = pars.rho0 * std::exp(ld); for (int j = 0; j < pars.npoints; j++) { // ln(d/d0) = \int_0^z - Omega^2 z dz Real pp = eos.PressureFromDensityTemperature(dens * (1. + dlnr), pars.temp0); Real pm = eos.PressureFromDensityTemperature(dens * (1. - dlnr), pars.temp0); Real dPdrho = (pp - pm) / (dlnr * dens); ld -= 0.5 * pars.Om0 * (2 * j + 1) * SQR(dz) / (dPdrho + Fuzz()); - dens = std::exp(ld); + dens = pars.rho0 * std::exp(ld); if (dens <= pars.dfloor) return pars.dfloor; } } diff --git a/src/utils/eos/ideal_h_he.hpp b/src/utils/eos/ideal_h_he.hpp index 35819353..7199f998 100644 --- a/src/utils/eos/ideal_h_he.hpp +++ b/src/utils/eos/ideal_h_he.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2025. 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 @@ -45,19 +45,37 @@ class IdealHHe : public singularity::eos_base::EosBase { IdealHHe() = default; IdealHHe(Real X, Real Y, Real ltmin, Real ltmax, int nt, Real ldmin, Real ldmax, int nd, const std::string &save_to_file, bool use_table = true, Real dlnT = 1e-6, + Real Eoffset = 0, int max_iters = 100, const singularity::MeanAtomicProperties &AZbar = singularity::MeanAtomicProperties()) : _X(X), _Y(Y), lTmin(ltmin), lTmax(ltmax), nt(nt), lDmin(ldmin), lDmax(ldmax), - nd(nd), use_table(use_table), _dlnT(dlnT), _AZbar(AZbar) { + nd(nd), use_table(use_table), _dlnT(dlnT), Eoffset(Eoffset), ITER_MAX(max_iters), + _AZbar(AZbar) { _fp = 0.25; _fo = 1. - _fp; CheckParams(); FillTable(save_to_file); } + IdealHHe(Real X, Real Y, Real dlnT = 1e-6, int max_iters = 100, + const singularity::MeanAtomicProperties &AZbar = + singularity::MeanAtomicProperties()) + : _X(X), _Y(Y), lTmin(0.0), lTmax(0.0), nt(0), lDmin(0.0), lDmax(0.0), nd(0), + Eoffset(0), use_table(false), _dlnT(dlnT), ITER_MAX(max_iters), _AZbar(AZbar) { + _fp = 0.25; + _fo = 1. - _fp; + CheckParams(); + } IdealHHe(const std::string &filename) { // Load from file Load(filename); } + void SetFloors(const Real sie_floor_, const Real t_floor_, const Real rho_floor_, + const Real p_floor_) { + rho_floor = rho_floor_; + t_floor = t_floor_; + sie_floor = sie_floor_; + p_floor = p_floor_; + } inline void FillTable(const std::string &filename); inline void Save(const std::string &filename); inline void Load(const std::string &filename); @@ -68,18 +86,11 @@ class IdealHHe : public singularity::eos_base::EosBase { PORTABLE_ALWAYS_REQUIRE(_Y >= 0, "Y must be positive"); _AZbar.CheckParams(); } - template - PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( - const Real rho, const Real sie, - Indexer_t &&lambda = static_cast(nullptr)) const { - const Real T = TofRE(rho, sie, lambda); - return T; - } template - PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( - const Real rho, const Real temperature, - Indexer_t &&lambda = static_cast(nullptr)) const { + PORTABLE_INLINE_FUNCTION Real + EofRT(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { const auto &r = GetMassFractions(rho, temperature); const Real kT = _eV * temperature; const Real henorm = 1. + r.z2 * (r.z1 - 1.0); @@ -99,40 +110,132 @@ class IdealHHe : public singularity::eos_base::EosBase { const Real efac = singularity::robust::safe_arg_exp(-lnx); const Real EH2 = 0.5 * _X * (1. - r.y) * (1.5 + lnx * singularity::robust::ratio(efac, 1. - efac) + - _fp * dlnp + _fo * (dlno - 2. * 85.5 / temperature)); + (_fp * xp.dZ + _fo * xo.dZ) / (_fp * xp.Z + _fo * xo.Z)); Real sie = EH2 + EH + EHe + (EHpH + EHp + EHep + EHepp) / temperature; sie = std::max(_small, sie * temperature * _kb / _mp); return sie; } template - PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( - const Real rho, const Real temperature, - Indexer_t &&lambda = static_cast(nullptr)) const { + PORTABLE_INLINE_FUNCTION Real + PofRT(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { const auto &[mu, dlmut, dlmur] = MeanMass(rho, temperature); const Real P = _kb / (mu * _mp) * rho * temperature; return std::max(_small, P); } template + PORTABLE_INLINE_FUNCTION Real + CvofRT(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + + Real tp = temperature * (1.0 + _dlnT); + Real tm = temperature * (1.0 - _dlnT); + + Real ep = EofRT(rho, tp); + Real em = EofRT(rho, tm); + Real cv = singularity::robust::ratio(ep - em, 2 * _dlnT * temperature); + return std::max(_small, cv); + } + template + PORTABLE_INLINE_FUNCTION Real + G1ofRT(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + // Gamma*rho = (dP/dT) / Cv + const Real Cv = CvofRT(rho, temperature); + const Real P = PofRT(rho, temperature); + + const auto &[mu, dlmudt, dlmudr] = MeanMass(rho, temperature); + // dln(P)/dln(T) + // dln(P)/dln(rho) + const Real xt = std::max(0.0, 1. - dlmudt); + const Real xd = std::max(0.0, 1. - dlmudr); + const Real denom = Cv * rho * temperature; + const Real fac1 = singularity::robust::ratio(P, Cv * rho * temperature); + const Real G1 = fac1 * xt * xt + xd; + return singularity::robust::ratio(P, Cv * rho * temperature) * xt * xt + xd; + } + template + PORTABLE_INLINE_FUNCTION Real + BofRT(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real P = PofRT(rho, temperature); + const Real G1 = G1ofRT(rho, temperature); + return std::max(_small, G1 * P); + } + PORTABLE_INLINE_FUNCTION std::tuple FillEosofRT(const Real rho, + const Real T) { + const Real P = PofRT(rho, T); + const Real G1 = G1ofRT(rho, T); + const Real B = std::max(_small, G1 * P); + const Real Cv = CvofRT(rho, T); + return {P, B, Cv, G1}; + } + + template + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + + const Real ld = std::log10(rho); + const Real lE = std::log10(sie + Eoffset); + if ((lE <= lEmin) || (rho <= rho_floor)) return std::pow(10., lTmin); + + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lE >= lEmin) && (lE <= lEmax))) { + return std::pow(10., lT_.interpToReal(ld, lE)); + } + + // fall back to inline + const Real T = TofRE(rho, sie); + return T; + } + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + + if ((temperature <= t_floor) || (rho <= rho_floor)) return sie_floor; + + const Real ld = std::log10(rho); + const Real lT = std::log10(temperature); + + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lT >= lTmin) && (lT <= lTmax))) { + Real sie = std::pow(10., lE_.interpToReal(ld, lT)); + return sie - Eoffset; + } + + Real sie = EofRT(rho, temperature); + return sie; + } + template + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real sie = InternalEnergyFromDensityTemperature(rho, temperature); + return PressureFromDensityInternalEnergy(rho, sie); + } + template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + if ((sie <= sie_floor) || (rho <= rho_floor)) return p_floor; + const Real ld = std::log10(rho); - const Real lE = std::log10(sie); + const Real lE = std::log10(sie + Eoffset); + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lE >= lEmin) && (lE <= lEmax))) { return std::pow(10., lP_.interpToReal(ld, lE)); } + // fall back to inline const Real T = TofRE(rho, sie); - return PressureFromDensityTemperature(rho, T); + return PofRT(rho, T); } template PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - // ? kb log(T)? - // ln(P/rho^g) ? return 0.0; } template @@ -146,143 +249,75 @@ class IdealHHe : public singularity::eos_base::EosBase { const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - // NOTE(@adempsey) - // Uncomment if you want the analytic Cv - - // const auto &r = GetMassFractions(rho, temperature); - // const Real kT = _eV * temperature; - // const Real henorm = 1. + r.z2 * (r.z1 - 1.0); - // const Real dhenorm = r.z2 * r.dz1dt + (r.z1 - 1.0) * r.dz2dt; - - // Real EH = 1.5 * _X * (1. + r.x) * r.y; // H - // Real cH = EH + 1.5 * _X * ((1. + r.x) * r.dydt + r.y * r.dxdt); - - // Real EHe = 0.375 * _Y * (1. + (r.z1 + r.z1 * r.z2) / henorm); // He - // Real cHe = EHe + 0.375 * _Y * - // (((1. + r.z2) * r.dz1dt + r.z1 * r.dz2dt) / henorm - - // r.z1 * (1. + r.z2) * dhenorm / (henorm * henorm)); - - // Real EHpH = _ye * _X * r.y * 0.5; - // Real cHpH = _ye * _X * r.dydt * 0.5; - - // Real EHp = _xe * _X * r.x * r.y; - // Real cHp = _xe * _X * (r.x * r.dydt + r.y * r.dxdt); - - // Real EHep = _z1e * 0.25 * _Y * r.z1 / henorm; - // Real cHep = - // _z1e * 0.25 * _Y * (r.dz1dt / henorm - r.z1 * dhenorm / (henorm * henorm)); - - // Real EHepp = _z2e * 0.25 * _Y * (r.z1 * r.z2 / henorm); - // Real cHepp = _z2e * 0.25 * _Y * - // ((r.z1 * r.dz2dt + r.z2 * r.dz1dt) / henorm - - // r.z1 * r.z2 * dhenorm / (henorm * henorm)); - - // // molecular hydrogen - // // NOTE: yes, the coefficients are in K - // constexpr Real tiny = std::numeric_limits::min(); - // const auto &[xp, xo] = H2PartitionFunctions(temperature); - // const Real dlno = (xo.Z > tiny) ? singularity::robust::ratio(xo.dZ, xo.Z) : 0.0; - // const Real dlnp = singularity::robust::ratio(xp.dZ, xp.Z); - // const Real d2lno = (xo.Z > tiny) ? singularity::robust::ratio(xo.d2Z, xo.Z) : 0.0; - // const Real d2lnp = singularity::robust::ratio(xp.d2Z, xp.Z); - // const Real lnx = 6140. / temperature; - // const Real xv = std::min(dlno * 0.5, 85.5 / temperature); - // const Real efac = singularity::robust::safe_arg_exp(-lnx); - // const Real evf = lnx * singularity::robust::ratio(efac, 1. - efac); - // // d1= 0.25 * (58.4792 / 58.6465) + 0.997148 * (0.00855 - 2*0) - // const Real d1 = _fp * dlnp + _fo * (dlno - 2. * xv); - // Real d2 = _fp * d2lnp + _fo * d2lno - 4. * _fo * (_fo * dlno + _fp * dlnp) * xv - - // _fp * _fo * SQR(dlnp - dlno) + 4. * _fo * (1. + _fo * xv) * xv; - - // // d2 = 0.0; - // Real EH2 = 0.5 * (1.5 + evf + _fp * dlnp + _fo * (dlno - 2. * xv)); - // Real cH2 = (1.5 + 2. * d1 + d2 - d1 * d1 + - // SQR(lnx) * singularity::robust::ratio(efac, SQR(1. - efac))) * - // 0.5; - // cH2 = _X * ((1. - r.y) * cH2 - EH2 * r.dydt); - - // // const Real EH2 = 0.5 * _X * (1. - r.y) * - // // (1.5 + lnx * singularity::robust::ratio(efac, 1. - efac) + - // // _fp * dlnp + _fo * (dlno - 2. * 85.5 / temperature)); - - // EH2 *= _X * (1. - r.y); - // Real cv = cH2 + cH + cHe + (cHpH + cHp + cHep + cHepp) / temperature; - - // cv = std::max(_small, cv * _kb / _mp); - // return cv; - - Real tp = temperature * (1.0 + _dlnT); - Real tm = temperature * (1.0 - _dlnT); - - Real ep = InternalEnergyFromDensityTemperature(rho, tp); - Real em = InternalEnergyFromDensityTemperature(rho, tm); - Real cv = singularity::robust::ratio(ep - em, 2 * _dlnT * temperature); - return std::max(_small, cv); + const Real sie = InternalEnergyFromDensityTemperature(rho, temperature); + return SpecificHeatFromDensityInternalEnergy(rho, sie); } + template PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + if ((sie <= sie_floor) || (rho <= rho_floor)) return sie_floor / t_floor; + const Real ld = std::log10(rho); - const Real lE = std::log10(sie); + const Real lE = std::log10(sie + Eoffset); + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lE >= lEmin) && (lE <= lEmax))) { return Cv_.interpToReal(ld, lE); } // fall back to inline const Real T = TofRE(rho, sie); - return SpecificHeatFromDensityTemperature(rho, T); + return CvofRT(rho, T); } + template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real P = PressureFromDensityTemperature(rho, temperature); - const Real G1 = GruneisenParamFromDensityTemperature(rho, temperature); - return std::max(_small, G1 * P); + const Real sie = InternalEnergyFromDensityTemperature(rho, temperature); + return BulkModulusFromDensityInternalEnergy(rho, sie); } template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + if ((sie <= sie_floor) || (rho <= rho_floor)) return p_floor; + const Real ld = std::log10(rho); - const Real lE = std::log10(sie); + const Real lE = std::log10(sie + Eoffset); + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lE >= lEmin) && (lE <= lEmax))) { return std::pow(10., lB_.interpToReal(ld, lE)); } + // fall back to inline const Real T = TofRE(rho, sie); - return BulkModulusFromDensityTemperature(rho, T); + return BofRT(rho, T); } + template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - // Gamma*rho = (dP/dT) / Cv - const Real Cv = SpecificHeatFromDensityTemperature(rho, temperature); - const Real P = PressureFromDensityTemperature(rho, temperature); - - const auto &[mu, dlmudt, dlmudr] = MeanMass(rho, temperature); - // dln(P)/dln(T) - // dln(P)/dln(rho) - const Real xt = std::max(0.0, 1. - dlmudt); - const Real xd = std::max(0.0, 1. - dlmudr); - const Real denom = Cv * rho * temperature; - const Real fac1 = singularity::robust::ratio(P, Cv * rho * temperature); - const Real G1 = fac1 * xt * xt + xd; - return singularity::robust::ratio(P, Cv * rho * temperature) * xt * xt + xd; + const Real sie = GruneisenParamFromDensityTemperature(rho, temperature); + return GruneisenParamFromDensityInternalEnergy(rho, sie); } template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + if ((sie <= sie_floor) || (rho <= rho_floor)) return 0.0; + const Real ld = std::log10(rho); - const Real lE = std::log10(sie); + const Real lE = std::log10(sie + Eoffset); + if (use_table && ((ld >= lDmin) && (ld <= lDmax) && (lE >= lEmin) && (lE <= lEmax))) { return Gm_.interpToReal(ld, lE); } + // fall back to inline const Real T = TofRE(rho, sie); - return GruneisenParamFromDensityTemperature(rho, T); + return G1ofRT(rho, T); } template PORTABLE_INLINE_FUNCTION void @@ -302,22 +337,7 @@ class IdealHHe : public singularity::eos_base::EosBase { SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) static constexpr unsigned long PreferredInput() { return _preferred_input; } - PORTABLE_INLINE_FUNCTION void PrintParams() const { - printf("Ideal H-He Parameters:\nX = %g\nY = %g\n", _X, _Y); - _AZbar.PrintParams(); - } - // template - // DensityEnergyFromPressureTemperature(const Real press, const Real temp, - // Indexer_t &&lambda, Real &rho, Real &sie) - // const { - // sie = std::max( - // _qq, singularity::robust::ratio(press + (_gm1 + 1.0) * _Pinf, press + _Pinf) - // * - // _Cv * temp + - // _qq); - // rho = std::max(singularity::robust::SMALL(), - // singularity::robust::ratio(press + _Pinf, _gm1 * _Cv * temp)); - // } + PORTABLE_INLINE_FUNCTION void PrintParams() const { _AZbar.PrintParams(); } inline void Finalize(); static std::string EosType() { return std::string("IdealHHe"); } static std::string EosPyType() { return EosType(); } @@ -327,11 +347,16 @@ class IdealHHe : public singularity::eos_base::EosBase { singularity::DEFAULT_SHMEM_STNGS); private: - bool use_table; + bool use_table = true; Real _X, _Y, _fp, _fo; Real lTmin, lTmax, lDmin, lDmax, _dlnT, lEmin, lEmax; + Real Eoffset = 0.0; + Real sie_floor = 1e-10; + Real t_floor = 1e-10; + Real rho_floor = 1e-20; + Real p_floor = 1e-20; int nd, nt; - DataBox lP_, lB_, lT_, Cv_, Gm_; + DataBox lP_, lB_, lT_, lE_, Cv_, Gm_; Real _small = 1e-15; Real _na = 6.02214129e23; Real _hbar = 1.0546e-27; // cm^2 g/s @@ -345,11 +370,12 @@ class IdealHHe : public singularity::eos_base::EosBase { Real _xe = 13.598433 / _eV; Real _z1e = 24.587387 / _eV; Real _z2e = 54.417760 / _eV; + int ITER_MAX = 100; singularity::MeanAtomicProperties _AZbar; static constexpr const unsigned long _preferred_input = singularity::thermalqs::density | singularity::thermalqs::temperature; -#define DBLIST &lP_, &lB_, &lT_, &Cv_, &Gm_ +#define DBLIST &lP_, &lB_, &lT_, &lE_, &Cv_, &Gm_ auto GetDataBoxPointers_() const { return std::vector{DBLIST}; } auto GetDataBoxPointers_() { return std::vector{DBLIST}; } #undef DBLIST @@ -360,7 +386,6 @@ class IdealHHe : public singularity::eos_base::EosBase { PORTABLE_INLINE_FUNCTION Real root_solve(const Real x, const Real a, const Real b, const Real c) const { // Newton-Raphson on the polynomial (b + c*y)*y - a*(1-y) - constexpr int ITER_MAX = 100; constexpr Real tol = 1e-15; // return x; Real f = (a + b + c * x) * x - a; @@ -415,7 +440,7 @@ class IdealHHe : public singularity::eos_base::EosBase { // x Real a = f1 / _X * f2e * singularity::robust::safe_arg_exp(-_xe / (T)); Real dlat = _xe / (T) + 1.5; - res.x = quadratic_root(a, 0., 1.0); + res.x = (_X == 0.0) ? 0.0 : quadratic_root(a, 0., 1.0); if ((std::abs(a) > _small) && (res.x > 0.0) && (res.x < 1.0)) { res.dxdt = singularity::robust::ratio(dlat * a * (1. - res.x), 2 * res.x + a); res.dxdr = singularity::robust::ratio(-a * (1. - res.x), (2 * res.x + a)); @@ -425,7 +450,7 @@ class IdealHHe : public singularity::eos_base::EosBase { Real efac_ = singularity::robust::safe_arg_exp(-_ye / (T)); Real fac1_ = f1 * f2p; a = 0.5 * f1 / _X * f2p * singularity::robust::safe_arg_exp(-_ye / (T)); - res.y = quadratic_root(a, 0.0, 1.0); + res.y = (_X == 0.0) ? 0.0 : quadratic_root(a, 0.0, 1.0); dlat = _ye / (T) + 1.5; if ((std::abs(a) > _small) || (res.y > 0.0) && (res.y < 1.0)) { res.dydt = singularity::robust::ratio(dlat * a * (1. - res.y), 2 * res.y + a); @@ -457,46 +482,57 @@ class IdealHHe : public singularity::eos_base::EosBase { return res; } - PORTABLE_INLINE_FUNCTION std::tuple - H2PartitionFunctions(const Real T) const { + PORTABLE_INLINE_FUNCTION H2Partition + H2PartitionFunction_(const Real T, std::array jstart) const { + const int max_iters = ITER_MAX * 200; const Real x = 85.5 / T; - // para - H2Partition para{0.0}; - for (int j = 0; j < 1000; j += 2) { - const int jj = j * (j + 1); - const Real dr = (2 * j + 1) * singularity::robust::safe_arg_exp(-jj * x); - para.Z += dr; - if (j >= 2) { - para.dZ += jj * dr; - para.d2Z += jj * (jj * x - 2.0) * dr; - if (singularity::robust::ratio(dr * std::abs(jj * (jj * x - 2.0)), - std::abs(para.d2Z)) < 1e-12) { - break; - } - } - } - para.dZ *= x; - para.d2Z *= x; - - // ortho - H2Partition ortho{0.0}; - for (int j = 1; j < 1000; j += 2) { - const int jj = j * (j + 1); - const Real dr = (2 * j + 1) * singularity::robust::safe_arg_exp(-jj * x); - ortho.Z += dr; - ortho.dZ += jj * dr; - // if (jj * x > 2.0) { - ortho.d2Z += jj * (jj * x - 2.0) * dr; - /// } - if (singularity::robust::ratio(dr * std::abs(jj * (jj * x - 2.0)), - std::abs(ortho.d2Z)) < 3e-16) { - break; - } + H2Partition part{0.0}; + + Real j = jstart[0]; + Real dr = 0.0; + Real dr_p = 0.0; + do { + part.Z += dr; + dr_p = dr; + const Real jj = j * (j + 1.); + dr = (2. * j + 1.) * singularity::robust::safe_arg_exp(-jj * x); + j += 2.0; + } while ((dr > dr_p || dr > 3e-16 * part.Z) && (j < max_iters)); + + j = jstart[1]; + dr = 0.0; + dr_p = 0.0; + do { + part.dZ += dr; + dr_p = dr; + const Real jj = j * (j + 1.); + dr = jj * (2. * j + 1.) * singularity::robust::safe_arg_exp(-jj * x); + j += 2.0; + } while ((dr > dr_p || dr > 3e-16 * part.dZ) && (j < max_iters)); + + j = jstart[2]; + dr = 0.0; + dr_p = 0.0; + while (j < max_iters) { + part.d2Z += dr; + dr_p = dr; + const Real jj = j * (j + 1.); + dr = (2. * j + 1.) * jj * (jj * x - 2.0) * + singularity::robust::safe_arg_exp(-jj * x); + j += 2.0; + if (dr < 0.0 || dr > dr_p || dr > 3e-16 * std::abs(part.d2Z)) continue; + break; } - ortho.dZ *= x; - ortho.d2Z *= x; + part.dZ *= x; + part.d2Z *= x; + return part; + } + PORTABLE_INLINE_FUNCTION std::tuple + H2PartitionFunctions(const Real T) const { + auto para = H2PartitionFunction_(T, {0.0, 2.0, 2.0}); + auto ortho = H2PartitionFunction_(T, {1.0, 1.0, 1.0}); return {para, ortho}; } PORTABLE_INLINE_FUNCTION std::tuple MeanMass(const Real rho, @@ -516,75 +552,26 @@ class IdealHHe : public singularity::eos_base::EosBase { Real dlmut = (get_mu(_X, _Y, rtp) - get_mu(_X, _Y, rtm)) / (2. * _dlnT * mu); Real dlmur = (get_mu(_X, _Y, rdp) - get_mu(_X, _Y, rdm)) / (2. * _dlnT * mu); - // NOTE(@adempsey) - // Uncomment for the analytic derivatives - - // const auto &r = GetMassFractions(rho, T); - // Real imu = - // 0.25 * (2. * _X * (1. + r.y * (1. + 2. * r.x)) + _Y * (1. + r.z1 * (1. + - // r.z2))); - // const Real mu = singularity::robust::ratio(1.0, imu); - // const Real dlmut = -0.25 * mu * - // (2 * _X * (r.dydt * (1.0 + 2 * r.x) + 2 * r.dxdt * r.y) + - // _Y * (r.dz1dt * (1.0 + r.z2) + r.dz2dt * r.z1)); - // const Real dlmur = -0.25 * mu * - // (2 * _X * (r.dydr * (1.0 + 2 * r.x) + 2 * r.dxdr * r.y) + - // _Y * (r.dz1dr * (1.0 + r.z2) + r.dz2dr * r.z1)); return {mu, dlmut, dlmur}; } template PORTABLE_INLINE_FUNCTION Real TofRE(const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - // Root finding version. - - // T_k+1 = T_k - (E(T) - E0)/cv(T) - // Use ideal gas as guess - // kb/(2.*mp) * T = E - Real T = sie * _mp / _kb; - int iter = 0; - Real sie_new = 0.0; - bool conv = false; - Real lower = std::pow(10., lTmin) * .01; - Real upper = std::pow(10., lTmax) * 100; - Real E_low = InternalEnergyFromDensityTemperature(rho, lower); - ; - Real E_high = InternalEnergyFromDensityTemperature(rho, upper); - Real dE_low = E_low - sie; - Real dE_high = E_high - sie; - - if (dE_low * dE_high > 0.0) { - printf("Initial sie %lg not bracketed at density %lg by temperature bounds [%lg, " - "%lg], [%lg, %lg]\n", - sie, rho, lower, upper, E_low, E_high); - } - - for (iter = 0; iter < 100; iter++) { - T = std::sqrt(lower * upper); - sie_new = InternalEnergyFromDensityTemperature(rho, T); - const Real dE = sie_new - sie; - if (std::abs(dE) <= 1e-8 * sie) { - conv = true; - break; - } - if (dE * dE_low < 0.0) { - upper = T; - dE_high = dE; - } else if (dE * dE_high < 0.0) { - lower = T; - dE_low = dE; - } else { - printf("Failed to converge %lg %lg %lg %lg\n", rho, sie, T, - InternalEnergyFromDensityTemperature(rho, T, lambda)); - break; - } - } - if (!conv) { - printf("Failed to converge %lg %lg %lg %lg\n", rho, sie, T, - InternalEnergyFromDensityTemperature(rho, T, lambda)); + using RootFinding1D::findRoot; + using RootFinding1D::Status; + + const Real tmin = 1e-10; + const Real tmax = 1e10; + Real temp; + const Real t_guess = 1e3; + if (findRoot([&](const Real t) { return EofRT(rho, t); }, sie, t_guess, tmin, tmax, + 1e-12 * t_guess, 1e-12, temp) != Status::SUCCESS) { + PARTHENON_DEBUG_WARN("TofRE did not converge"); + return Tiny(); } - return std::max(_small, T); + return std::max(_small, temp); } }; @@ -626,6 +613,7 @@ inline void IdealHHe::FillTable(const std::string &filename) { lT_.resize(nd, nt); lT_.setRange(0, lTmin, lTmax, nt); lT_.setRange(1, lDmin, lDmax, nd); + lE_.copyMetadata(lT_); // Determine the energy grid lEmin = std::numeric_limits::max(); lEmax = std::numeric_limits::min(); @@ -636,21 +624,17 @@ inline void IdealHHe::FillTable(const std::string &filename) { const Real T = std::pow(10., lT_.range(0).x(i)); const auto &[mu, dlmut, dlmur] = MeanMass(d, T); const auto &r = GetMassFractions(d, T); - const Real E = InternalEnergyFromDensityTemperature(d, T); - const Real P = PressureFromDensityTemperature(d, T); - const Real cv = SpecificHeatFromDensityTemperature(d, T); - const Real G = GruneisenParamFromDensityTemperature(d, T); - const Real B = BulkModulusFromDensityTemperature(d, T); - lEmin = std::min(lEmin, E); - lEmax = std::max(lEmax, E); + const Real E = EofRT(d, T); + const Real lE = std::log10(E + Eoffset); + lE_(j, i) = lE; + lEmin = std::min(lEmin, lE); + lEmax = std::max(lEmax, lE); } } if (lEmin <= 0.0 || (lEmax <= 0.0) || std::isnan(lEmin) || std::isnan(lEmax)) { PORTABLE_THROW_OR_ABORT("Failed to find positive or real energy values from given " "temperature and density grid."); } - lEmin = std::log10(lEmin); - lEmax = std::log10(lEmax); lT_.setRange(0, lEmin, lEmax, nt); lP_.copyMetadata(lT_); lB_.copyMetadata(lT_); @@ -661,13 +645,14 @@ inline void IdealHHe::FillTable(const std::string &filename) { for (int j = 0; j < nd; j++) { const Real d = std::pow(10., lT_.range(1).x(j)); for (int i = 0; i < nt; i++) { - const Real e = std::pow(10., lT_.range(0).x(i)); + const Real e = std::pow(10., lT_.range(0).x(i)) - Eoffset; const Real T = TofRE(d, e); lT_(j, i) = std::log10(T); - lP_(j, i) = std::log10(PressureFromDensityTemperature(d, T)); - lB_(j, i) = std::log10(BulkModulusFromDensityTemperature(d, T)); - Cv_(j, i) = SpecificHeatFromDensityTemperature(d, T); - Gm_(j, i) = GruneisenParamFromDensityTemperature(d, T); + const auto &[P, B, Cv, G1] = FillEosofRT(d, T); + lP_(j, i) = std::log10(P); + lB_(j, i) = std::log10(B); + Cv_(j, i) = Cv; + Gm_(j, i) = G1; } } @@ -677,7 +662,7 @@ inline void IdealHHe::FillTable(const std::string &filename) { const Real d = std::pow(10., lT_.range(1).x(j)); for (int i = 0; i < nt; i++) { const Real lE = lT_.range(0).x(i); - const Real e = std::pow(10., lT_.range(0).x(i)); + const Real e = std::pow(10., lT_.range(0).x(i)) - Eoffset; const Real T = TofRE(d, e); assert(std::abs(std::pow(10., lT_.interpToReal(ld, lE)) / T - 1) <= 1e-4); } @@ -705,6 +690,9 @@ inline void IdealHHe::Save(const std::string &filename) { status += H5LTset_attribute_double(file, METADATA_NAME, "ltmax", &lTmax, 1); status += H5LTset_attribute_double(file, METADATA_NAME, "ldmin", &lDmin, 1); status += H5LTset_attribute_double(file, METADATA_NAME, "ldmax", &lDmax, 1); + status += H5LTset_attribute_double(file, METADATA_NAME, "lemin", &lEmin, 1); + status += H5LTset_attribute_double(file, METADATA_NAME, "lemax", &lEmax, 1); + status += H5LTset_attribute_double(file, METADATA_NAME, "eoffset", &Eoffset, 1); status += H5LTset_attribute_double(file, METADATA_NAME, "dlnT", &_dlnT, 1); status += H5LTset_attribute_double(file, METADATA_NAME, "fp", &_fp, 1); status += H5LTset_attribute_double(file, METADATA_NAME, "fm", &_fo, 1); @@ -717,6 +705,7 @@ inline void IdealHHe::Save(const std::string &filename) { status += Cv_.saveHDF(file, "cv"); status += lB_.saveHDF(file, "logbulkmodulus"); status += Gm_.saveHDF(file, "grun"); + status += lE_.saveHDF(file, "logsie"); status += H5Fclose(file); if (status != H5_SUCCESS) { @@ -736,6 +725,9 @@ inline void IdealHHe::Load(const std::string &filename) { status += H5LTget_attribute_double(file, METADATA_NAME, "ltmax", &lTmax); status += H5LTget_attribute_double(file, METADATA_NAME, "ldmin", &lDmin); status += H5LTget_attribute_double(file, METADATA_NAME, "ldmax", &lDmax); + status += H5LTget_attribute_double(file, METADATA_NAME, "lemin", &lEmin); + status += H5LTget_attribute_double(file, METADATA_NAME, "lemax", &lEmax); + status += H5LTget_attribute_double(file, METADATA_NAME, "eoffset", &Eoffset); status += H5LTget_attribute_double(file, METADATA_NAME, "dlnT", &_dlnT); status += H5LTget_attribute_double(file, METADATA_NAME, "fp", &_fp); status += H5LTget_attribute_double(file, METADATA_NAME, "fm", &_fo); @@ -748,6 +740,7 @@ inline void IdealHHe::Load(const std::string &filename) { status += Cv_.loadHDF(file, "cv"); status += lB_.loadHDF(file, "logbulkmodulus"); status += Gm_.loadHDF(file, "grun"); + status += lE_.loadHDF(file, "logsie"); status += H5Fclose(file); if (status != H5_SUCCESS) { EOS_ERROR("[IdealHHe::Save]: There was a problem with HDF5\n"); diff --git a/tst/unit/test_ideal_hhe.cpp b/tst/unit/test_ideal_hhe.cpp index 58afaed2..33a187aa 100644 --- a/tst/unit/test_ideal_hhe.cpp +++ b/tst/unit/test_ideal_hhe.cpp @@ -26,25 +26,33 @@ using singularity::IdealGas; #ifdef SPINER_USE_HDF TEST_CASE("IdealHHe EOS basic construction and properties", "[IdealHHe][EOS]") { - // Use solar composition: X = 0.7, Y = 0.28 (Z = 0.02) const Real X = 0.7; - const Real Y = 0.28; - - // Temperature range: 1 K to 1e6 K (ltmin=0, ltmax=6 in log10) - const Real ltmin = 0.0; - const Real ltmax = 6.0; - const int nt = 50; - - // Density range: 1e-15 to 1e-3 g/cc (ldmin=-15, ldmax=-3 in log10) - const Real ldmin = -15.0; - const Real ldmax = -3.0; - const int nd = 50; - - // Create EOS without table lookup for testing - const bool use_table = false; - const std::string save_to_file = ""; - - IdealHHe eos(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table); + const Real Y = 0.3; + + IdealHHe eos(X, Y); + eos.SetFloors(0.0, 0.0, 0.0, 0.0); + + SECTION("2D Grid") { + const Real lTmin = -2; + const Real lTmax = 6.; + const int nt = 10; + const int nd = 4; + const Real lDmin = -25.0; + const Real lDmax = 5.; + const Real dlT = (lTmax - lTmin) / static_cast(nt); + const Real dlD = (lDmax - lDmin) / static_cast(nd); + for (int i = 0; i <= nd; i++) { + const Real rho = std::pow(10., lDmin + i * dlD); + for (int j = 0; j <= nt; j++) { + const Real T = std::pow(10., lTmin + j * dlT); + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T); + const Real cv = eos.SpecificHeatFromDensityTemperature(rho, T); + // invert + const Real T_ = eos.TemperatureFromDensityInternalEnergy(rho, sie); + REQUIRE(T_ == Approx(T).epsilon(1e-10)); + } + } + } SECTION("Temperature recovery from density and internal energy") { const Real rho = 1.0e-10; // g/cc @@ -57,7 +65,7 @@ TEST_CASE("IdealHHe EOS basic construction and properties", "[IdealHHe][EOS]") { const Real T_recovered = eos.TemperatureFromDensityInternalEnergy(rho, sie); // Temperature should be recovered to within 1% - REQUIRE(T_recovered == Approx(T_input).epsilon(0.01)); + REQUIRE(T_recovered == Approx(T_input).epsilon(1e-15)); } SECTION("Pressure calculation from density and temperature") { @@ -97,19 +105,65 @@ TEST_CASE("IdealHHe EOS basic construction and properties", "[IdealHHe][EOS]") { } } +TEST_CASE("IdealHHE EOS with X = 0") { + IdealHHe eos(0.0, 1.0); + eos.SetFloors(0.0, 0.0, 0.0, 0.0); + + SECTION("2D Grid") { + const Real lTmin = -2; + const Real lTmax = 6.; + const int nt = 10; + const int nd = 4; + const Real lDmin = -25.0; + const Real lDmax = 5.; + const Real dlT = (lTmax - lTmin) / static_cast(nt); + const Real dlD = (lDmax - lDmin) / static_cast(nd); + for (int i = 0; i <= nd; i++) { + const Real rho = std::pow(10., lDmin + i * dlD); + for (int j = 0; j <= nt; j++) { + const Real T = std::pow(10., lTmin + j * dlT); + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T); + const Real cv = eos.SpecificHeatFromDensityTemperature(rho, T); + // invert + const Real T_ = eos.TemperatureFromDensityInternalEnergy(rho, sie); + REQUIRE(T_ == Approx(T).epsilon(1e-12)); + } + } + } +} + +TEST_CASE("IdealHHE EOS with Y = 0") { + IdealHHe eos(1.0, 0.0); + eos.SetFloors(0.0, 0.0, 0.0, 0.0); + + SECTION("2D Grid") { + const Real lTmin = -2; + const Real lTmax = 6.; + const int nt = 10; + const int nd = 4; + const Real lDmin = -25.0; + const Real lDmax = 5.; + const Real dlT = (lTmax - lTmin) / static_cast(nt); + const Real dlD = (lDmax - lDmin) / static_cast(nd); + for (int i = 0; i <= nd; i++) { + const Real rho = std::pow(10., lDmin + i * dlD); + for (int j = 0; j <= nt; j++) { + const Real T = std::pow(10., lTmin + j * dlT); + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T); + const Real cv = eos.SpecificHeatFromDensityTemperature(rho, T); + // invert + const Real T_ = eos.TemperatureFromDensityInternalEnergy(rho, sie); + REQUIRE(T_ == Approx(T).epsilon(1e-9)); + } + } + } +} + TEST_CASE("IdealHHe EOS with different compositions", "[IdealHHe][EOS]") { - const Real ltmin = 0.0; - const Real ltmax = 6.0; - const int nt = 50; - const Real ldmin = -15.0; - const Real ldmax = -3.0; - const int nd = 50; - const bool use_table = false; - const std::string save_to_file = ""; SECTION("Hydrogen-rich composition (X=0.9, Y=0.1)") { - IdealHHe eos_H_rich(0.9, 0.1, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, - use_table); + IdealHHe eos_H_rich(0.9, 0.1); + eos_H_rich.SetFloors(0.0, 0.0, 0.0, 0.0); const Real rho = 1.0e-10; const Real T = 1000.0; @@ -122,8 +176,8 @@ TEST_CASE("IdealHHe EOS with different compositions", "[IdealHHe][EOS]") { } SECTION("Helium-rich composition (X=0.3, Y=0.7)") { - IdealHHe eos_He_rich(0.3, 0.7, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, - use_table); + IdealHHe eos_He_rich(0.3, 0.7); + eos_He_rich.SetFloors(0.0, 0.0, 0.0, 0.0); const Real rho = 1.0e-10; const Real T = 1000.0; @@ -140,16 +194,9 @@ TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") // Use solar composition const Real X = 0.7; const Real Y = 0.28; - const Real ltmin = 0.0; - const Real ltmax = 6.0; - const int nt = 50; - const Real ldmin = -15.0; - const Real ldmax = -3.0; - const int nd = 50; - const bool use_table = false; - const std::string save_to_file = ""; - IdealHHe eos_base(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table); + IdealHHe eos_base(X, Y); + eos_base.SetFloors(0.0, 0.0, 0.0, 0.0); SECTION("Unit conversion with CGS units") { // Define unit system: time [s], mass [g], length [cm], temperature [K] @@ -160,9 +207,8 @@ TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") // Create IdealHHe inline as temporary for UnitSystem auto eos = singularity::UnitSystem( - IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table), - singularity::eos_units_init::LengthTimeUnitsInit(), time_cgs, mass_cgs, - length_cgs, temp_cgs); + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), time_cgs, + mass_cgs, length_cgs, temp_cgs); const Real rho = 1.0e-10; // g/cc const Real T = 1000.0; // K @@ -185,11 +231,12 @@ TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") const Real length_scale = 1.496e13; const Real temp_scale = 1.0; + IdealHHe eos_base(X, Y); + eos_base.SetFloors(0.0, 0., 0.0, 0.0); // Create IdealHHe inline as temporary for UnitSystem auto eos = singularity::UnitSystem( - IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table), - singularity::eos_units_init::LengthTimeUnitsInit(), time_scale, mass_scale, - length_scale, temp_scale); + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), + time_scale, mass_scale, length_scale, temp_scale); // Density in code units (Msun/AU^3) // Convert 1e-10 g/cc to Msun/AU^3: 1e-10 * (AU^3/Msun) @@ -213,9 +260,11 @@ TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") SECTION("Consistency between wrapped and unwrapped EOS in CGS units") { // Wrap with CGS units (scale factors = 1) + IdealHHe eos_base(X, Y); + eos_base.SetFloors(0.0, 0., 0.0, 0.0); auto eos_wrapped = singularity::UnitSystem( - IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table), - singularity::eos_units_init::LengthTimeUnitsInit(), 1.0, 1.0, 1.0, 1.0); + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), 1.0, 1.0, + 1.0, 1.0); const Real rho = 1.0e-10; // g/cc const Real T = 1000.0; // K @@ -240,10 +289,11 @@ TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") const Real temp_scale = 1.0; // Create IdealHHe inline as temporary for UnitSystem + IdealHHe eos_base(X, Y); + eos_base.SetFloors(0.0, 0., 0.0, 0.0); auto eos = singularity::UnitSystem( - IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table), - singularity::eos_units_init::LengthTimeUnitsInit(), time_scale, mass_scale, - length_scale, temp_scale); + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), + time_scale, mass_scale, length_scale, temp_scale); // Density in code units const Real rho_cgs = 1.0e-10;