From 1b7bfa8df2fd94ea154f03489294a64e26818c50 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 23 Jun 2023 15:51:36 -0600 Subject: [PATCH 01/45] builds with chris changes but not yet hierarchical grids --- sesame2spiner/io_eospac.hpp | 4 +- sesame2spiner/test.cpp | 4 +- singularity-eos/eos/eos_spiner.hpp | 48 +++++++------ singularity-eos/eos/eos_stellar_collapse.hpp | 76 ++++++++++---------- test/compare_to_eospac.cpp | 5 +- test/profile_eos.cpp | 4 +- test/profile_eospac.cpp | 4 +- utils/ports-of-call | 2 +- utils/spiner | 2 +- 9 files changed, 77 insertions(+), 72 deletions(-) diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 82518dbf17..9ff15bbd41 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -37,8 +37,8 @@ #include using EospacWrapper::Verbosity; -using Spiner::DataBox; -using Spiner::RegularGrid1D; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; // For logarithmic interpolation, quantities may be negative. // If they are, use offset to ensure negative values make sense. diff --git a/sesame2spiner/test.cpp b/sesame2spiner/test.cpp index d3f6f8f1ee..c7bd1bd766 100644 --- a/sesame2spiner/test.cpp +++ b/sesame2spiner/test.cpp @@ -39,8 +39,8 @@ #include "io_eospac.hpp" #include "parser.hpp" -using Spiner::DataBox; -using Spiner::RegularGrid1D; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; constexpr int matid = 5030; // dry air constexpr int matid2 = 4272; // stainless steel diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 1af778e455..8be31d8844 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -65,6 +65,8 @@ using namespace eos_base; */ class SpinerEOSDependsRhoT : public EosBase { public: + using DataBox = Spiner::DataBox; + // A weakly typed index map for lambdas struct Lambda { enum Index { lRho = 0, lT = 1 }; @@ -252,12 +254,12 @@ class SpinerEOSDependsRhoT : public EosBase { thermalqs::density | thermalqs::temperature; // static constexpr const char _eos_type[] {"SpinerEOSDependsRhoT"}; static constexpr const int numDataBoxes_ = 12; - Spiner::DataBox P_, sie_, bMod_, dPdRho_, dPdE_, dTdRho_, dTdE_, dEdRho_, dEdT_; - Spiner::DataBox PMax_, sielTMax_, dEdTMax_, gm1Max_; - Spiner::DataBox lTColdCrit_; - Spiner::DataBox PCold_, sieCold_, bModCold_; - Spiner::DataBox dPdRhoCold_, dPdECold_, dTdRhoCold_, dTdECold_, dEdTCold_; - Spiner::DataBox rho_at_pmin_; + DataBox P_, sie_, bMod_, dPdRho_, dPdE_, dTdRho_, dTdE_, dEdRho_, dEdT_; + DataBox PMax_, sielTMax_, dEdTMax_, gm1Max_; + DataBox lTColdCrit_; + DataBox PCold_, sieCold_, bModCold_; + DataBox dPdRhoCold_, dPdECold_, dTdRhoCold_, dTdECold_, dEdTCold_; + DataBox rho_at_pmin_; int numRho_, numT_; Real lRhoMin_, lRhoMax_, rhoMax_; Real lRhoMinSearch_; @@ -298,8 +300,9 @@ class SpinerEOSDependsRhoT : public EosBase { */ class SpinerEOSDependsRhoSie : public EosBase { public: + using DataBox = Spiner::DataBox; struct SP5Tables { - Spiner::DataBox P, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho; + DataBox P, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho; }; // Generic functions provided by the base class. These contain // e.g. the vector overloads that use the scalar versions declared @@ -437,23 +440,23 @@ class SpinerEOSDependsRhoSie : public EosBase { return FastMath::pow10(lx) - offset; } PORTABLE_INLINE_FUNCTION - Real interpRhoT_(const Real rho, const Real T, const Spiner::DataBox &db, + Real interpRhoT_(const Real rho, const Real T, const DataBox &db, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION - Real interpRhoSie_(const Real rho, const Real sie, const Spiner::DataBox &db, + Real interpRhoSie_(const Real rho, const Real sie, const DataBox &db, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real lRhoFromPlT_(const Real P, const Real lT, Real *lambda) const; - Spiner::DataBox sie_; // depends on (rho,T) - Spiner::DataBox T_; // depends on (rho, sie) + DataBox sie_; // depends on (rho,T) + DataBox T_; // depends on (rho, sie) SP5Tables dependsRhoT_; SP5Tables dependsRhoSie_; int numRho_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoMin_, lRhoMax_, rhoMax_; - Spiner::DataBox PlRhoMax_, dPdRhoMax_; + DataBox PlRhoMax_, dPdRhoMax_; Real lRhoOffset_, lTOffset_, lEOffset_; // offsets must be non-negative @@ -479,15 +482,16 @@ class SpinerEOSDependsRhoSie : public EosBase { // replace lambdas with callable namespace callable_interp { +using DataBox = Spiner::DataBox; class l_interp { private: - const Spiner::DataBox &field; + const DataBox &field; const Real fixed; public: PORTABLE_INLINE_FUNCTION - l_interp(const Spiner::DataBox &field_, const Real fixed_) + l_interp(const DataBox &field_, const Real fixed_) : field{field_}, fixed{fixed_} {} PORTABLE_INLINE_FUNCTION Real operator()(const Real x) const { @@ -497,12 +501,12 @@ class l_interp { class r_interp { private: - const Spiner::DataBox &field; + const DataBox &field; const Real fixed; public: PORTABLE_INLINE_FUNCTION - r_interp(const Spiner::DataBox &field_, const Real fixed_) + r_interp(const DataBox &field_, const Real fixed_) : field{field_}, fixed{fixed_} {} PORTABLE_INLINE_FUNCTION Real operator()(const Real x) const { @@ -512,12 +516,12 @@ class r_interp { class prod_interp_1d { private: - const Spiner::DataBox &field1, field2; + const DataBox &field1, field2; const Real r; public: PORTABLE_INLINE_FUNCTION - prod_interp_1d(const Spiner::DataBox &field1_, const Spiner::DataBox &field2_, + prod_interp_1d(const DataBox &field1_, const DataBox &field2_, const Real r_) : field1{field1_}, field2{field2_}, r{r_} {} PORTABLE_INLINE_FUNCTION Real operator()(const Real x) const { @@ -527,11 +531,11 @@ class prod_interp_1d { class interp { private: - const Spiner::DataBox &field; + const DataBox &field; public: PORTABLE_INLINE_FUNCTION - interp(const Spiner::DataBox &field_) : field(field_) {} + interp(const DataBox &field_) : field(field_) {} PORTABLE_INLINE_FUNCTION Real operator()(const Real x) const { return field.interpToReal(x); } @@ -1839,7 +1843,7 @@ void SpinerEOSDependsRhoSie::ValuesAtReferenceState(Real &rho, Real &temp, Real PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::interpRhoT_(const Real rho, const Real T, - const Spiner::DataBox &db, Real *lambda) const { + const DataBox &db, Real *lambda) const { const Real lRho = toLog_(rho, lRhoOffset_); const Real lT = toLog_(T, lTOffset_); if (lambda != nullptr) { @@ -1850,7 +1854,7 @@ Real SpinerEOSDependsRhoSie::interpRhoT_(const Real rho, const Real T, PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoSie::interpRhoSie_(const Real rho, const Real sie, - const Spiner::DataBox &db, + const DataBox &db, Real *lambda) const { const Real lRho = toLog_(rho, lRhoOffset_); const Real lE = toLog_(sie, lEOffset_); diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index cb5d0340be..1fab1c55ca 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -58,6 +58,8 @@ using namespace eos_base; // and introduce extrapolation as needed. class StellarCollapse : public EosBase { public: + using DataBox = Spiner::DataBox; + // A weakly typed index map for lambdas struct Lambda { enum Index { Ye = 0, lT = 1 }; @@ -190,18 +192,32 @@ class StellarCollapse : public EosBase { static std::string EosPyType() { return EosType(); } private: + class LogT { + public: + PORTABLE_INLINE_FUNCTION + LogT(const DataBox &field, const Real Ye, const Real lRho) + : field_(field), Ye_(Ye), lRho_(lRho) {} + PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { + return field_.interpToReal(Ye_, lT, lRho_); + } + + private: + const DataBox &field_; + const Real Ye_, lRho_; + }; + inline void LoadFromSP5File_(const std::string &filename); inline void LoadFromStellarCollapseFile_(const std::string &filename); inline int readSCInt_(const hid_t &file_id, const std::string &name); inline void readBounds_(const hid_t &file_id, const std::string &name, int size, Real &lo, Real &hi); inline void readSCDset_(const hid_t &file_id, const std::string &name, - Spiner::DataBox &db); + DataBox &db); - inline void medianFilter_(Spiner::DataBox &db); - inline void medianFilter_(const Spiner::DataBox &in, Spiner::DataBox &out); + inline void medianFilter_(DataBox &db); + inline void medianFilter_(const DataBox &in, DataBox &out); inline void fillMedianBuffer_(Real buffer[], int width, int iY, int iT, int irho, - const Spiner::DataBox &tab) const; + const DataBox &tab) const; inline Real findMedian_(Real buffer[], int size) const; inline void computeBulkModulus_(); inline void computeColdAndHotCurves_(); @@ -283,19 +299,19 @@ class StellarCollapse : public EosBase { thermalqs::density | thermalqs::temperature; // Dependent variables - Spiner::DataBox lP_, lE_, dPdRho_, dPdE_, dEdT_, lBMod_; - Spiner::DataBox entropy_; // kb/baryon - Spiner::DataBox Xa_; // mass fraction of alpha particles - Spiner::DataBox Xh_; // mass fraction of heavy ions - Spiner::DataBox Xn_; // mass fraction of neutrons - Spiner::DataBox Xp_; // mass fraction of protons - Spiner::DataBox Abar_; // Average atomic mass - Spiner::DataBox Zbar_; // Average atomic number - // Spiner::DataBox gamma_; // polytropic index. dlog(P)/dlog(rho). + DataBox lP_, lE_, dPdRho_, dPdE_, dEdT_, lBMod_; + DataBox entropy_; // kb/baryon + DataBox Xa_; // mass fraction of alpha particles + DataBox Xh_; // mass fraction of heavy ions + DataBox Xn_; // mass fraction of neutrons + DataBox Xp_; // mass fraction of protons + DataBox Abar_; // Average atomic mass + DataBox Zbar_; // Average atomic number + // DataBox gamma_; // polytropic index. dlog(P)/dlog(rho). // dTdRho_, dTdE_, dEdRho_, dEdT_; // Bounds of dependent variables. Needed for root finding. - Spiner::DataBox eCold_, eHot_; + DataBox eCold_, eHot_; // Independent variable bounds int numRho_, numT_, numYe_; @@ -338,24 +354,6 @@ class StellarCollapse : public EosBase { // Implementation details below // ====================================================================== -namespace callable_interp { - -class LogT { - public: - PORTABLE_INLINE_FUNCTION - LogT(const Spiner::DataBox &field, const Real Ye, const Real lRho) - : field_(field), Ye_(Ye), lRho_(lRho) {} - PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { - return field_.interpToReal(Ye_, lT, lRho_); - } - - private: - const Spiner::DataBox &field_; - const Real Ye_, lRho_; -}; - -} // namespace callable_interp - // For some reason, the linker doesn't like this being a member field // of StellarCollapse. So we'll make it a global variable. constexpr char METADATA_NAME[] = "Metadata"; @@ -817,7 +815,7 @@ inline void StellarCollapse::readBounds_(const hid_t &file_id, const std::string * https://forum.hdfgroup.org/t/is-this-a-bug-in-hdf5-1-8-6/2211 */ inline void StellarCollapse::readSCDset_(const hid_t &file_id, const std::string &name, - Spiner::DataBox &db) { + DataBox &db) { herr_t exists = H5LTfind_dataset(file_id, name.c_str()); if (!exists) { std::string msg = "Tried to read dataset " + name + " but it doesn't exist\n"; @@ -836,15 +834,15 @@ inline void StellarCollapse::readSCDset_(const hid_t &file_id, const std::string db.setRange(0, lRhoMin_, lRhoMax_, numRho_); } -inline void StellarCollapse::medianFilter_(Spiner::DataBox &db) { - Spiner::DataBox tmp; +inline void StellarCollapse::medianFilter_(DataBox &db) { + DataBox tmp; tmp.copy(db); medianFilter_(tmp, db); free(tmp); } -inline void StellarCollapse::medianFilter_(const Spiner::DataBox &in, - Spiner::DataBox &out) { +inline void StellarCollapse::medianFilter_(const DataBox &in, + DataBox &out) { Real buffer[MF_S]; // filter, overwriting as needed for (int iY = MF_W; iY < numYe_ - MF_W; ++iY) { @@ -863,7 +861,7 @@ inline void StellarCollapse::medianFilter_(const Spiner::DataBox &in, inline void StellarCollapse::fillMedianBuffer_(Real buffer[], int width, int iY, int iT, int irho, - const Spiner::DataBox &tab) const { + const DataBox &tab) const { int i = 0; for (int iWy = -width; iWy <= width; iWy++) { for (int iWt = -width; iWt <= width; iWt++) { @@ -989,7 +987,7 @@ Real StellarCollapse::lTFromlRhoSie_(const Real lRho, const Real sie, } // Get log(sie) Real lE = e2le_(sie); - const callable_interp::LogT lEFunc(lE_, Ye, lRho); + const LogT lEFunc(lE_, Ye, lRho); status = regula_falsi(lEFunc, lE, lTGuess, lTMin_, lTMax_, ROOT_THRESH, ROOT_THRESH, lT, pcounts); if (status != RootFinding1D::Status::SUCCESS) { diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 03c86c8089..bcd1f5685b 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -47,10 +47,13 @@ #include -using namespace Spiner; +//using namespace Spiner; using namespace singularity; using namespace EospacWrapper; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; + using duration = std::chrono::microseconds; constexpr char diffFileName[] = "diffs.sp5"; constexpr int NTIMES = 10; diff --git a/test/profile_eos.cpp b/test/profile_eos.cpp index dfae6b17e4..fec104438f 100644 --- a/test/profile_eos.cpp +++ b/test/profile_eos.cpp @@ -40,8 +40,8 @@ using namespace singularity; using duration = std::chrono::microseconds; using dvec = std::vector; using ivec = std::vector; -using Spiner::DataBox; -using Spiner::RegularGrid1D; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; diff --git a/test/profile_eospac.cpp b/test/profile_eospac.cpp index 07f732667d..654d0b02c3 100644 --- a/test/profile_eospac.cpp +++ b/test/profile_eospac.cpp @@ -42,8 +42,8 @@ using namespace singularity; using duration = std::chrono::microseconds; using dvec = std::vector; using ivec = std::vector; -using Spiner::DataBox; -using Spiner::RegularGrid1D; +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; diff --git a/utils/ports-of-call b/utils/ports-of-call index 1b64df8189..1ae17ad8b6 160000 --- a/utils/ports-of-call +++ b/utils/ports-of-call @@ -1 +1 @@ -Subproject commit 1b64df818933bbe2f755e32e59a57735886fa123 +Subproject commit 1ae17ad8b6d51b14e46a13d329a09f76c7e8aa29 diff --git a/utils/spiner b/utils/spiner index cba97e8cd9..c76462824c 160000 --- a/utils/spiner +++ b/utils/spiner @@ -1 +1 @@ -Subproject commit cba97e8cd90af34d77499cc72830a31f33a20d1b +Subproject commit c76462824ca56e83ce6257f05da3454af042ac0c From 1b4081e5bab56ba073d6e653f9260d81e2dcf467 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 4 Aug 2023 14:53:33 -0600 Subject: [PATCH 02/45] update to spiner main with piecewise grids --- utils/spiner | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/spiner b/utils/spiner index cba97e8cd9..db0c8e6230 160000 --- a/utils/spiner +++ b/utils/spiner @@ -1 +1 @@ -Subproject commit cba97e8cd90af34d77499cc72830a31f33a20d1b +Subproject commit db0c8e623023037650e31e97557e5e6c37bc0892 From e85a8fbdac7528dfb765aea8d3e95df9b2adbf62 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 4 Aug 2023 15:26:40 -0600 Subject: [PATCH 03/45] fix so it compiles --- singularity-eos/eos/eos_spiner.hpp | 4 ---- singularity-eos/eos/eos_stellar_collapse.hpp | 2 -- 2 files changed, 6 deletions(-) diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 128bdaf720..bc433dadf7 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -65,8 +65,6 @@ using namespace eos_base; we use log-linear extrapolation. */ class SpinerEOSDependsRhoT : public EosBase { - using DataBox = Spiner::DataBox; - public: using DataBox = Spiner::DataBox; @@ -302,8 +300,6 @@ class SpinerEOSDependsRhoT : public EosBase { mitigated by Ye and (1-Ye) to control how important each term is. */ class SpinerEOSDependsRhoSie : public EosBase { - using DataBox = Spiner::DataBox; - public: using DataBox = Spiner::DataBox; struct SP5Tables { diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 617fa55117..b6d4cb1675 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -58,8 +58,6 @@ using namespace eos_base; // is linear extrapolation in log-log space. We should reconsider this // and introduce extrapolation as needed. class StellarCollapse : public EosBase { - using DataBox = Spiner::DataBox; - public: using DataBox = Spiner::DataBox; From 792e4f7b103ab9e9e96d8475957d52a4a7c791a6 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 4 Aug 2023 15:39:40 -0600 Subject: [PATCH 04/45] NGRIDS --- singularity-eos/eos/eos_spiner.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index bc433dadf7..1c94a3355b 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -66,7 +66,9 @@ using namespace eos_base; */ class SpinerEOSDependsRhoT : public EosBase { public: - using DataBox = Spiner::DataBox; + static constexpr int NGRIDS = 3; + using Grid_t = Spiner::PiecewiseGrid1D; + using DataBox = Spiner::DataBox; // A weakly typed index map for lambdas struct Lambda { @@ -301,7 +303,8 @@ class SpinerEOSDependsRhoT : public EosBase { */ class SpinerEOSDependsRhoSie : public EosBase { public: - using DataBox = Spiner::DataBox; + using Grid_t = SpinerEOSDependsRhoT::Grid_t; + using DataBox = SpinerEOSDependsRhoT::DataBox; struct SP5Tables { DataBox P, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho; }; @@ -483,9 +486,7 @@ class SpinerEOSDependsRhoSie : public EosBase { // replace lambdas with callable namespace callable_interp { -using DataBox = Spiner::DataBox; - -using DataBox = Spiner::DataBox; +using DataBox = SpinerEOSDependsRhoT::DataBox; class l_interp { private: const DataBox &field; @@ -736,8 +737,7 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, // fill in minimum pressure as a function of temperature rho_at_pmin_.resize(numT_); - // rho_at_pmin_.setRange(0, P_.range(0)); - rho_at_pmin_.setRange(0, lTMin_, lTMax_, numT_); + rho_at_pmin_.setRange(0, P_.range(0)); for (int i = 0; i < numT_; i++) { Real pmin = std::numeric_limits::max(); int jmax = -1; From ad8d9aee6b174d64a5b57f2e252aa9ebfc97b9fa Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 4 Aug 2023 15:49:19 -0600 Subject: [PATCH 05/45] slowly modifying bounds object to use piecewise grids --- sesame2spiner/io_eospac.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 86b780b415..8e430a1010 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -37,8 +37,10 @@ #include using EospacWrapper::Verbosity; -using DataBox = Spiner::DataBox; +constexpr int NGRIDS = 3; using RegularGrid1D = Spiner::RegularGrid1D; +using Grid_t = Spiner::PiecewiseGrid1D; +using DataBox = Spiner::DataBox; // For logarithmic interpolation, quantities may be negative. // If they are, use offset to ensure negative values make sense. @@ -47,7 +49,7 @@ class Bounds { Bounds() {} Bounds(Real min, Real max, int N, Real offset) - : grid(RegularGrid1D(min, max, N)), offset(offset) {} + : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)}), offset(offset) {} Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) @@ -90,7 +92,7 @@ class Bounds { } } - grid = RegularGrid1D(min, max, N); + grid = Grid_t(min, max, N); } inline Real log2lin(Real xl) const { return singularity::FastMath::pow10(xl) - offset; } @@ -105,7 +107,7 @@ class Bounds { } public: - RegularGrid1D grid; + Grid_t grid; Real offset; }; From d3b5ad1160d8fcec78885d968b1c9bbf081563ec Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 29 Jul 2024 15:47:30 -0600 Subject: [PATCH 06/45] move bounds into singularity-eos. rewrite it for piecewise grids and begin testing. --- sesame2spiner/generate_files.cpp | 6 ++ sesame2spiner/io_eospac.hpp | 79 +-------------- singularity-eos/CMakeLists.txt | 1 + singularity-eos/base/table_bounds.hpp | 138 ++++++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/test_bounds.cpp | 49 +++++++++ 6 files changed, 198 insertions(+), 76 deletions(-) create mode 100644 singularity-eos/base/table_bounds.hpp create mode 100644 test/test_bounds.cpp diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 7fea0e706b..6c74f82121 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -289,6 +289,12 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real rhoAnchor = metadata.normalDensity; Real TAnchor = 298.15; + // Piecewise grids stuff + Real rho_fine_center = rhoAnchor; + Real rho_fine_diameter = 0.5; + Real T_fine_center = TAnchor; + Real T_fine_diameter = 1.0; + // Forces density and temperature to be in a region where an offset // is not needed. This improves resolution at low densities and // temperatures. diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 8e430a1010..6ecb454d1e 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 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 @@ -17,10 +17,6 @@ #ifndef _SESAME2SPINER_IO_EOSPAC_HPP_ #define _SESAME2SPINER_IO_EOSPAC_HPP_ -#include -#include -#include -#include #include #include // eospac API @@ -31,86 +27,17 @@ #include #include +#include #include -#include #include using EospacWrapper::Verbosity; constexpr int NGRIDS = 3; -using RegularGrid1D = Spiner::RegularGrid1D; +using Bounds = singularity::Bounds; using Grid_t = Spiner::PiecewiseGrid1D; using DataBox = Spiner::DataBox; -// For logarithmic interpolation, quantities may be negative. -// If they are, use offset to ensure negative values make sense. -class Bounds { - public: - Bounds() {} - - Bounds(Real min, Real max, int N, Real offset) - : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)}), offset(offset) {} - - Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, - Real anchor_point = std::numeric_limits::signaling_NaN()) - : offset(0) { - if (convertToLog) { - // Log scales can't handle negative numbers or exactly zero. To - // deal with that, we offset. - constexpr Real epsilon = std::numeric_limits::epsilon(); - const Real min_offset = 10 * std::abs(epsilon); - // 1.1 so that the y-intercept isn't identically zero - // when min < 0. - // min_offset to handle the case where min=0 - if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - - min += offset; - max += offset; - - min = singularity::FastMath::log10(std::abs(min)); - max = singularity::FastMath::log10(std::abs(max)); - Real delta = max - min; - min += 0.5 * shrinkRange * delta; - max -= 0.5 * shrinkRange * delta; - - if (!(std::isnan(anchor_point))) { - anchor_point += offset; - anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - } - } - - if (!(std::isnan(anchor_point))) { - if (min < anchor_point && anchor_point < max) { - Real dxguess = (max - min) / ((Real)N - 1); - int Nmax = static_cast((max - min) / dxguess); - int Nanchor = static_cast((anchor_point - min) / dxguess); - Real dx = (anchor_point - min) / static_cast(Nanchor + 1); - int Nmax_new = static_cast((max - min) / dx); - Nmax_new = std::max(Nmax, Nmax_new); - max = dx * (Nmax_new) + min; - N = Nmax_new; - } - } - - grid = Grid_t(min, max, N); - } - - inline Real log2lin(Real xl) const { return singularity::FastMath::pow10(xl) - offset; } - inline Real i2lin(int i) const { return log2lin(grid.x(i)); } - - friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { - os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" - << " + " << b.offset << ", " - << "[N,dx] = [" << b.grid.nPoints() << ", " << b.grid.dx() << "]" - << "\n"; - return os; - } - - public: - Grid_t grid; - Real offset; -}; - void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds, DataBox &P, DataBox &T, DataBox &bMods, DataBox &dPdRho, DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, DataBox &dEdRho, diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 399702f989..f4c0ea2d63 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -27,6 +27,7 @@ register_headers( base/fast-math/logs.hpp base/robust_utils.hpp base/root-finding-1d/root_finding.hpp + base/table_bounds.hpp base/variadic_utils.hpp base/math_utils.hpp base/constants.hpp diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp new file mode 100644 index 0000000000..14c36a1a70 --- /dev/null +++ b/singularity-eos/base/table_bounds.hpp @@ -0,0 +1,138 @@ +//------------------------------------------------------------------------------ +// © 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. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ +#define SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ +#ifdef SINGULARITY_USE_SPINER + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace singularity { + +// For logarithmic interpolation, quantities may be negative. +// If they are, use offset to ensure negative values make sense. +template +class Bounds { + public: + using RegularGrid1D = Spiner::RegularGrid1D; + using Grid_t = Spiner::PiecewiseGrid1D; + + Bounds() {} + + Bounds(Real min, Real max, int N, Real offset) + : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), + offset(offset), piecewise(false) {} + + Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, + Real anchor_point = std::numeric_limits::signaling_NaN()) + : offset(0), piecewise(true) { + if (convertToLog) { + convertBoundsToLog_(min, max, shrinkRange); + if (!(std::isnan(anchor_point))) { + anchor_point += offset; + anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); + } + } + if (!(std::isnan(anchor_point))) { + adjustForAnchor_(min, max, N, anchor_point); + } + grid = Grid_t(std::vector{RegularGrid1D(min, max, N)}); + } + + Bounds(Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, + int N_fine, Real N_factor, Real shrinkRange = 0) + : offset(0), piecewise(true) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + anchor_point += offset; + anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); + + int N_coarse = ((N_factor >= 1) && (N_fine > N_factor)) + ? static_cast(std::ceil(N_fine / N_factor)) + : N_fine; + + Real mid_min = anchor_point - 0.5 * log_fine_diameter; + Real mid_max = anchor_point + 0.5 * log_fine_diameter; + adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); + + RegularGrid1D grid_middle(mid_min, mid_max, N_fine); + RegularGrid1D grid_lower(global_min, mid_min, N_coarse); + RegularGrid1D grid_upper(mid_max, mid_max, N_coarse); + + grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); + } + + inline Real log2lin(Real xl) const { return singularity::FastMath::pow10(xl) - offset; } + inline Real i2lin(int i) const { return log2lin(grid.x(i)); } + + friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { + os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" + << " + " << b.offset << ", " + << "[N,dx] = [" << b.grid.nPoints() << ", " << b.grid.dx() << "]" + << "\n"; + return os; + } + + private: + void convertBoundsToLog_(Real &min, Real &max, Real shrinkRange = 0) { + // Log scales can't handle negative numbers or exactly zero. To + // deal with that, we offset. + constexpr Real epsilon = std::numeric_limits::epsilon(); + const Real min_offset = 10 * std::abs(epsilon); + // 1.1 so that the y-intercept isn't identically zero + // when min < 0. + // min_offset to handle the case where min=0 + if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; + + min += offset; + max += offset; + + min = singularity::FastMath::log10(std::abs(min)); + max = singularity::FastMath::log10(std::abs(max)); + Real delta = max - min; + min += 0.5 * shrinkRange * delta; + max -= 0.5 * shrinkRange * delta; + } + + void adjustForAnchor_(Real min, Real &max, int &N, Real anchor_point) { + if (min < anchor_point && anchor_point < max) { + Real dxguess = (max - min) / ((Real)N - 1); + int Nmax = static_cast((max - min) / dxguess); + int Nanchor = static_cast((anchor_point - min) / dxguess); + Real dx = (anchor_point - min) / static_cast(Nanchor + 1); + int Nmax_new = static_cast((max - min) / dx); + Nmax_new = std::max(Nmax, Nmax_new); + max = dx * (Nmax_new) + min; + N = Nmax_new; + } + } + + public: + Grid_t grid; + Real offset = 0; + bool piecewise = false; +}; + +} // namespace singularity + +#endif // SINGULARITY_USE_SPINER +#endif // SINGULARITY_EOS_BASE_TABLE_BOUNDS_HPP_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4c654975ec..c0a127feb6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -34,6 +34,7 @@ add_executable( test_eos_vector.cpp test_math_utils.cpp test_variadic_utils.cpp + test_bounds.cpp ) add_executable( diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp new file mode 100644 index 0000000000..079df01399 --- /dev/null +++ b/test/test_bounds.cpp @@ -0,0 +1,49 @@ +//------------------------------------------------------------------------------ +// © 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. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +constexpr int NGRIDS = 3; +using Bounds = singularity::Bounds; +using RegularGrid1D = Spiner::RegularGrid1D; +using Grid_t = Bounds::Grid_t; + +constexpr Real rho_normal = 3; // g/cm^3 +constexpr Real T_normal = 300; // Kelvin +constexpr Real rho_min = 1e-6; +constexpr Real rho_max = 1e3; +constexpr Real T_min = 1; +constexpr Real T_max = 1e9; +constexpr int N_per_decade_fine = 10; +constexpr Real N_factor = 5; + +SCENARIO("linear bounds") { + WHEN("We compute linear bounds ") { + constexpr Real min = -2; + constexpr Real max = 5; + constexpr int N = 21; + Bounds lin_bounds(min, max, N); + THEN("The min and max and point number correct") { + REQUIRE( lin_bounds.grid.min() == min ); + REQUIRE( lin_bounds.grid.max() == max ); + REQUIRE( lin_bounds.grid.nPoints() == N ); + } + } +} From f2502d3add4f0ef0cabb179caafeeba8ef047cbe Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 30 Jul 2024 16:57:46 -0600 Subject: [PATCH 07/45] add table bounds tests --- sesame2spiner/generate_files.cpp | 12 +--- singularity-eos/base/table_bounds.hpp | 47 +++++++++----- singularity-eos/eos/eos_spiner.hpp | 2 +- singularity-eos/eos/eos_stellar_collapse.hpp | 8 +-- test/compare_to_eospac.cpp | 2 +- test/test_bounds.cpp | 68 ++++++++++++++++++-- 6 files changed, 104 insertions(+), 35 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 6c74f82121..fd25137511 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -274,13 +274,13 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params } int ppdRho = params.Get("numrho/decade", PPD_DEFAULT); - int numRhoDefault = getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); + int numRhoDefault = Bounds::getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); int ppdT = params.Get("numT/decade", PPD_DEFAULT); - int numTDefault = getNumPointsFromPPD(TMin, TMax, ppdT); + int numTDefault = Bounds::getNumPointsFromPPD(TMin, TMax, ppdT); int ppdSie = params.Get("numSie/decade", PPD_DEFAULT); - int numSieDefault = getNumPointsFromPPD(sieMin, sieMax, ppdSie); + int numSieDefault = Bounds::getNumPointsFromPPD(sieMin, sieMax, ppdSie); int numRho = params.Get("numrho", numRhoDefault); int numT = params.Get("numT", numTDefault); @@ -322,9 +322,3 @@ bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin } return true; } - -int getNumPointsFromPPD(Real min, Real max, int ppd) { - Bounds b(min, max, 3, true); - Real ndecades = b.grid.max() - b.grid.min(); - return static_cast(std::ceil(ppd * ndecades)); -} diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 14c36a1a70..91c6710dc8 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -31,7 +32,7 @@ namespace singularity { // For logarithmic interpolation, quantities may be negative. // If they are, use offset to ensure negative values make sense. -template +template class Bounds { public: using RegularGrid1D = Spiner::RegularGrid1D; @@ -60,23 +61,27 @@ class Bounds { } Bounds(Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, - int N_fine, Real N_factor, Real shrinkRange = 0) + Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) : offset(0), piecewise(true) { + const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; + convertBoundsToLog_(global_min, global_max, shrinkRange); anchor_point += offset; anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - int N_coarse = ((N_factor >= 1) && (N_fine > N_factor)) - ? static_cast(std::ceil(N_fine / N_factor)) - : N_fine; - Real mid_min = anchor_point - 0.5 * log_fine_diameter; Real mid_max = anchor_point + 0.5 * log_fine_diameter; - adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); + // add a point just to make sure we have enough points after adjusting for anchor + int N_fine = getNumPointsFromPPD(mid_min, mid_max, ppd_fine) + 1; + adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); RegularGrid1D grid_middle(mid_min, mid_max, N_fine); - RegularGrid1D grid_lower(global_min, mid_min, N_coarse); - RegularGrid1D grid_upper(mid_max, mid_max, N_coarse); + + const int N_lower = getNumPointsFromPPD(global_min, mid_min, ppd_coarse); + RegularGrid1D grid_lower(global_min, mid_min, N_lower); + + const int N_upper = getNumPointsFromPPD(mid_max, global_max, N_upper); + RegularGrid1D grid_upper(mid_max, global_max, N_upper); grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); } @@ -92,6 +97,13 @@ class Bounds { return os; } + template + static int getNumPointsFromPPD(Real min, Real max, T ppd) { + Bounds b(min, max, 3, true); + Real ndecades = b.grid.max() - b.grid.min(); + return std::max(2, static_cast(std::ceil(ppd * ndecades))); + } + private: void convertBoundsToLog_(Real &min, Real &max, Real shrinkRange = 0) { // Log scales can't handle negative numbers or exactly zero. To @@ -115,14 +127,19 @@ class Bounds { void adjustForAnchor_(Real min, Real &max, int &N, Real anchor_point) { if (min < anchor_point && anchor_point < max) { - Real dxguess = (max - min) / ((Real)N - 1); - int Nmax = static_cast((max - min) / dxguess); - int Nanchor = static_cast((anchor_point - min) / dxguess); - Real dx = (anchor_point - min) / static_cast(Nanchor + 1); + Real Nfrac = (anchor_point - min) / (max - min); + PORTABLE_REQUIRE((0 < Nfrac && Nfrac < 1), "anchor in bounds"); + + int Nanchor = static_cast(std::ceil(N * Nfrac)); + if ((Nanchor < 2) || (Nanchor >= N)) return; // not possible to shift this safely + + Real dx = (anchor_point - min) / static_cast(Nanchor - 1); int Nmax_new = static_cast((max - min) / dx); - Nmax_new = std::max(Nmax, Nmax_new); - max = dx * (Nmax_new) + min; + Real max_new = dx * (Nmax_new - 1) + min; + PORTABLE_REQUIRE(max_new <= max, "must not exceed table bounds"); + N = Nmax_new; + max = max_new; } } diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 5e43209593..8b1d57343c 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -70,7 +70,7 @@ class SpinerEOSDependsRhoT : public EosBase { static constexpr int NGRIDS = 3; using Grid_t = Spiner::PiecewiseGrid1D; using DataBox = Spiner::DataBox; - + // A weakly typed index map for lambdas struct Lambda { enum Index { lRho = 0, lT = 1 }; diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index b8aab80bb0..1bbf1dfc9a 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -221,15 +221,15 @@ class StellarCollapse : public EosBase { private: class LogT { - public: + public: PORTABLE_INLINE_FUNCTION LogT(const DataBox &field, const Real Ye, const Real lRho) - : field_(field), Ye_(Ye), lRho_(lRho) {} + : field_(field), Ye_(Ye), lRho_(lRho) {} PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { return field_.interpToReal(Ye_, lT, lRho_); } - - private: + + private: const DataBox &field_; const Real Ye_, lRho_; }; diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 447b69815c..dee58fbcf1 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -47,7 +47,7 @@ #include -//using namespace Spiner; +// using namespace Spiner; using namespace singularity; using namespace EospacWrapper; diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 079df01399..8907f076ee 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -12,7 +12,10 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +#include + #include +#include #include #ifndef CATCH_CONFIG_FAST_COMPILE @@ -32,18 +35,73 @@ constexpr Real rho_max = 1e3; constexpr Real T_min = 1; constexpr Real T_max = 1e9; constexpr int N_per_decade_fine = 10; -constexpr Real N_factor = 5; +constexpr Real N_factor = 2; + +SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { + WHEN("We compute the number of points from points per decade") { + int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); + THEN("We get the right number") { REQUIRE(np == 90); } + } +} -SCENARIO("linear bounds") { +SCENARIO("Linear bounds in the bounds object", "[Bounds]") { WHEN("We compute linear bounds ") { constexpr Real min = -2; constexpr Real max = 5; constexpr int N = 21; Bounds lin_bounds(min, max, N); THEN("The min and max and point number correct") { - REQUIRE( lin_bounds.grid.min() == min ); - REQUIRE( lin_bounds.grid.max() == max ); - REQUIRE( lin_bounds.grid.nPoints() == N ); + REQUIRE(lin_bounds.grid.min() == min); + REQUIRE(lin_bounds.grid.max() == max); + REQUIRE(lin_bounds.grid.nPoints() == N); + } + } +} + +SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { + WHEN("We compute logarithmic single-grid bounds") { + int np = Bounds::getNumPointsFromPPD(rho_min, rho_max, N_per_decade_fine); + Bounds lRhoBounds(rho_min, rho_max, np, true, 0.0, rho_normal); + THEN("The lower and upper bounds are right") { + REQUIRE(std::abs(lRhoBounds.grid.min() - singularity::FastMath::log10(rho_min)) <= + 1e-12); + REQUIRE(lRhoBounds.grid.max() <= + singularity::FastMath::log10(rho_max)); // shifted due to anchor + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(rho_normal); + int ianchor; + Spiner::weights_t w; + lRhoBounds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= 1e-12); + REQUIRE(std::abs(w[1]) <= 1e-12); + } + } + } +} + +SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { + WHEN("We compute a piecewise bounds object") { + Bounds bnds(rho_min, rho_max, rho_normal, 0.5, N_per_decade_fine, N_factor); + THEN("The bounds are right") { + Real lrmin = singularity::FastMath::log10(rho_min); + Real lrmax = singularity::FastMath::log10(rho_max); + REQUIRE(std::abs(bnds.grid.min() - lrmin) <= 1e-12); + REQUIRE(std::abs(bnds.grid.max() - lrmax) <= 1e-12); + REQUIRE(bnds.grid.nGrids() == 3); + AND_THEN( + "The total number of points is less than a uniform fine spacing woudl imply") { + printf("npoints = %ld %e\n", bnds.grid.nPoints(), + N_per_decade_fine * (lrmax - lrmin)); + REQUIRE(bnds.grid.nPoints() < N_per_decade_fine * (lrmax - lrmin)); + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(rho_normal); + int ianchor; + Spiner::weights_t w; + bnds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= 1e-12); + REQUIRE(std::abs(w[1]) <= 1e-12); + } + } } } } From 39391be98169276a9f8e5ca8bbfaa5e8cbcf3395 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 10:58:18 -0600 Subject: [PATCH 08/45] modified bounds to support 2 and 3 grid modes. --- sesame2spiner/generate_files.cpp | 54 +++++++++++++++++++++------ sesame2spiner/generate_files.hpp | 4 +- sesame2spiner/io_eospac.cpp | 2 +- singularity-eos/base/table_bounds.hpp | 35 ++++++++++++++++- test/test_bounds.cpp | 33 +++++++++++++--- 5 files changed, 108 insertions(+), 20 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index fd25137511..de2f0dece8 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -273,27 +274,30 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params << "shrinkleBounds > 0 and siemin or siemax set" << std::endl; } - int ppdRho = params.Get("numrho/decade", PPD_DEFAULT); + int ppdRho = params.Get("numrho/decade", PPD_DEFAULT_RHO); int numRhoDefault = Bounds::getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); - int ppdT = params.Get("numT/decade", PPD_DEFAULT); + int ppdT = params.Get("numT/decade", PPD_DEFAULT_T); int numTDefault = Bounds::getNumPointsFromPPD(TMin, TMax, ppdT); - int ppdSie = params.Get("numSie/decade", PPD_DEFAULT); + int ppdSie = params.Get("numSie/decade", PPD_DEFAULT_T); int numSieDefault = Bounds::getNumPointsFromPPD(sieMin, sieMax, ppdSie); int numRho = params.Get("numrho", numRhoDefault); int numT = params.Get("numT", numTDefault); int numSie = params.Get("numsie", numSieDefault); - Real rhoAnchor = metadata.normalDensity; - Real TAnchor = 298.15; + const Real rhoAnchor = metadata.normalDensity; + constexpr Real TAnchor = 298.15; // Piecewise grids stuff + bool do_piecewise_grids = params.Get("piecewise", true); + Real ppd_factor_rho = params.Get("rhoCoarseFactor", COARSE_FACTOR_DEFAULT); + Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT); + Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT); + Real TSplitPoint = params.Get("TSplitPoint", 1e4); + Real rho_fine_center = rhoAnchor; - Real rho_fine_diameter = 0.5; - Real T_fine_center = TAnchor; - Real T_fine_diameter = 1.0; // Forces density and temperature to be in a region where an offset // is not needed. This improves resolution at low densities and @@ -304,9 +308,37 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params if (rhoMin < STRICTLY_POS_MIN) rhoMin = STRICTLY_POS_MIN; if (TMin < STRICTLY_POS_MIN) TMin = STRICTLY_POS_MIN; - lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); - lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); - leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); + if (do_piecewise_grids) { + lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, + rho_fine_diameter, ppdRho, ppd_factor_rho, shrinklRhoBounds); + lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, + ppd_factor_T, shrinklTBounds); + + // compute temperature as a reasonable anchor point + constexpr int NT = 1; + constexpr EOS_INTEGER nXYPairs = 2; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER tableType[NT] = {EOS_Ut_DT}; + EOS_REAL rho[2], T[2], sie[2], dx[2], dy[2]; + { + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Ut_DT"}, Verbosity::Quiet); + EOS_INTEGER eospacEofRT = tableHandle[0]; + rho[0] = rho[1] = densityToSesame(rhoAnchor); + T[0] = temperatureToSesame(TAnchor); + T[1] = temperatureToSesame(TSplitPoint); + eosSafeInterpolate(&eospacEofRT, nXYPairs, rho, T, sie, dx, dy, "EofRT", + Verbosity::Quiet); + eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); + } + Real sieAnchor = sieFromSesame(sie[0]); + Real sieSplitPoint = sieFromSesame(sie[1]); + leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, + ppdSie, ppd_factor_sie, shrinkleBounds); + } else { + lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); + lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); + leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); + } return; } diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 90fb7d2276..6a9ffd4ab2 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -30,8 +30,10 @@ using namespace EospacWrapper; -constexpr int PPD_DEFAULT = 50; +constexpr int PPD_DEFAULT_RHO = 200; +constexpr int PPD_DEFAULT_T = 140; constexpr Real STRICTLY_POS_MIN = 1e-9; +constexpr Real COARSE_FACTOR_DEFAULT = 20; herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds, const Bounds &lTBounds, const Bounds &leBounds, diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/io_eospac.cpp index 0953029223..e074fb8fd0 100644 --- a/sesame2spiner/io_eospac.cpp +++ b/sesame2spiner/io_eospac.cpp @@ -30,7 +30,7 @@ void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds DataBox &mask, Verbosity eospacWarn) { using namespace EospacWrapper; - constexpr int NT = 3; + constexpr int NT = 4; constexpr EOS_INTEGER nXYPairs = 1; EOS_INTEGER tableHandle[NT]; EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 91c6710dc8..d86d4271e6 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -37,12 +38,17 @@ class Bounds { public: using RegularGrid1D = Spiner::RegularGrid1D; using Grid_t = Spiner::PiecewiseGrid1D; + // for tag dispatch in constructors + using OneGrid = std::integral_constant; + using TwoGrids = std::integral_constant; + using ThreeGrids = std::integral_constant; Bounds() {} Bounds(Real min, Real max, int N, Real offset) : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), offset(offset), piecewise(false) {} + Bounds(OneGrid, Real min, Real max, int N, Real offset) : Bounds(min, max, N, offset) {} Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) @@ -59,12 +65,37 @@ class Bounds { } grid = Grid_t(std::vector{RegularGrid1D(min, max, N)}); } + Bounds(OneGrid, Real min, Real max, int N, bool convertToLog = false, + Real shrinkRange = 0, + Real anchor_point = std::numeric_limits::signaling_NaN()) + : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} - Bounds(Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, + Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) : offset(0), piecewise(true) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; + convertBoundsToLog_(global_min, global_max, shrinkRange); + anchor_point += offset; + anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); + splitPoint += offset; + splitPoint = singularity::FastMath::log10(std::abs(splitPoint)); + + // add a point just to make sure we have enough points after adjusting for anchor + int N_fine = getNumPointsFromPPD(global_min, splitPoint, ppd_fine) + 1; + adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); + RegularGrid1D grid_lower(global_min, splitPoint, N_fine); + + const int N_upper = getNumPointsFromPPD(splitPoint, global_max, ppd_coarse); + RegularGrid1D grid_upper(splitPoint, global_max, N_upper); + + grid = Grid_t(std::vector{grid_lower, grid_upper}); + } + + Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, + Real log_fine_diameter, Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) { + const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; + convertBoundsToLog_(global_min, global_max, shrinkRange); anchor_point += offset; anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); @@ -80,7 +111,7 @@ class Bounds { const int N_lower = getNumPointsFromPPD(global_min, mid_min, ppd_coarse); RegularGrid1D grid_lower(global_min, mid_min, N_lower); - const int N_upper = getNumPointsFromPPD(mid_max, global_max, N_upper); + const int N_upper = getNumPointsFromPPD(mid_max, global_max, ppd_coarse); RegularGrid1D grid_upper(mid_max, global_max, N_upper); grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 8907f076ee..7a1f4049d4 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -34,6 +34,7 @@ constexpr Real rho_min = 1e-6; constexpr Real rho_max = 1e3; constexpr Real T_min = 1; constexpr Real T_max = 1e9; +constexpr Real T_split = 1e4; constexpr int N_per_decade_fine = 10; constexpr Real N_factor = 2; @@ -80,8 +81,9 @@ SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { } SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { - WHEN("We compute a piecewise bounds object") { - Bounds bnds(rho_min, rho_max, rho_normal, 0.5, N_per_decade_fine, N_factor); + WHEN("We compute a piecewise bounds object with three grids") { + Bounds bnds(Bounds::ThreeGrids(), rho_min, rho_max, rho_normal, 0.5, + N_per_decade_fine, N_factor); THEN("The bounds are right") { Real lrmin = singularity::FastMath::log10(rho_min); Real lrmax = singularity::FastMath::log10(rho_max); @@ -89,9 +91,7 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { REQUIRE(std::abs(bnds.grid.max() - lrmax) <= 1e-12); REQUIRE(bnds.grid.nGrids() == 3); AND_THEN( - "The total number of points is less than a uniform fine spacing woudl imply") { - printf("npoints = %ld %e\n", bnds.grid.nPoints(), - N_per_decade_fine * (lrmax - lrmin)); + "The total number of points is less than a uniform fine spacing would imply") { REQUIRE(bnds.grid.nPoints() < N_per_decade_fine * (lrmax - lrmin)); AND_THEN("The anchor is on the mesh") { Real lanchor = singularity::FastMath::log10(rho_normal); @@ -104,4 +104,27 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { } } } + WHEN("We compute a piecewise bounds object with two grids") { + Bounds bnds(Bounds::TwoGrids(), T_min, T_max, T_normal, T_split, N_per_decade_fine, + N_factor); + THEN("The bounds are right") { + Real ltmin = singularity::FastMath::log10(T_min); + Real ltmax = singularity::FastMath::log10(T_max); + REQUIRE(std::abs(bnds.grid.min() - ltmin) <= 1e-12); + REQUIRE(std::abs(bnds.grid.max() - ltmax) <= 1e-12); + REQUIRE(bnds.grid.nGrids() == 2); + AND_THEN( + "The total number of points is less than a uniform fine spacing would imply") { + REQUIRE(bnds.grid.nPoints() < N_per_decade_fine * (ltmax - ltmin)); + AND_THEN("The anchor is on the mesh") { + Real lanchor = singularity::FastMath::log10(T_normal); + int ianchor; + Spiner::weights_t w; + bnds.grid.weights(lanchor, ianchor, w); + REQUIRE(std::abs(w[0] - 1) <= 1e-12); + REQUIRE(std::abs(w[1]) <= 1e-12); + } + } + } + } } From d22d3e90fa9399aa6457abae35b2a46a43b86472 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 11:17:44 -0600 Subject: [PATCH 09/45] typos in sesame2spiner --- sesame2spiner/generate_files.cpp | 1 + singularity-eos/base/table_bounds.hpp | 8 +++++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index de2f0dece8..8351ddff39 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -295,6 +295,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real ppd_factor_rho = params.Get("rhoCoarseFactor", COARSE_FACTOR_DEFAULT); Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT); Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT); + Real rho_fine_diameter = params.Get("rhoFineDiameterDecades", 0.5); Real TSplitPoint = params.Get("TSplitPoint", 1e4); Real rho_fine_center = rhoAnchor; diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index d86d4271e6..034352be5e 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -122,9 +122,11 @@ class Bounds { friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" - << " + " << b.offset << ", " - << "[N,dx] = [" << b.grid.nPoints() << ", " << b.grid.dx() << "]" - << "\n"; + << " + " << b.offset << "\n" + << "\tN = " << b.grid.nPoints() << "\n"; + for (int ig = 0; ig < b.grid.nGrids(); ++ig) { + os << "\t[ig,dx] = [" << ig << ", " << b.grid.dx(ig) << "]" << "\n"; + } return os; } From 6cde9e48dc65ed6b50f29c26ba9e89efdf2dd615 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 11:18:09 -0600 Subject: [PATCH 10/45] formatting --- singularity-eos/base/table_bounds.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 034352be5e..3e3aeb82de 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -125,7 +125,8 @@ class Bounds { << " + " << b.offset << "\n" << "\tN = " << b.grid.nPoints() << "\n"; for (int ig = 0; ig < b.grid.nGrids(); ++ig) { - os << "\t[ig,dx] = [" << ig << ", " << b.grid.dx(ig) << "]" << "\n"; + os << "\t[ig,dx] = [" << ig << ", " << b.grid.dx(ig) << "]" + << "\n"; } return os; } From b6c2589f0bed78e5b5a113c967a7b6fec3b066de Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 11:25:54 -0600 Subject: [PATCH 11/45] fix type issues in compare_to_eospac --- test/compare_to_eospac.cpp | 84 +++++++++++++++++++------------------- 1 file changed, 42 insertions(+), 42 deletions(-) diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index dee58fbcf1..959abc85b7 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -64,7 +64,7 @@ class Bounds { Bounds() {} Bounds(Real min, Real max, int N, Real offset) - : grid(RegularGrid1D(min, max, N)), offset(offset) {} + : grid(RegularGrid1D(min, max, N)), offset(offset) {} Bounds(Real min, Real max, int N) : offset(0) { constexpr Real epsilon = std::numeric_limits::epsilon(); @@ -74,7 +74,7 @@ class Bounds { max += offset; min = std::log10(std::abs(min)); max = std::log10(std::abs(max)); - grid = RegularGrid1D(min, max, N); + grid = RegularGrid1D(min, max, N); } PORTABLE_INLINE_FUNCTION Real log2lin(Real xl) const { return pow(10., xl) - offset; } @@ -89,7 +89,7 @@ class Bounds { } public: - RegularGrid1D grid; + RegularGrid1D grid; Real offset; }; @@ -149,59 +149,59 @@ int main(int argc, char *argv[]) { constexpr int NT = 3; EOS_INTEGER nXYPairs = (nFineRho) * (nFineT); - DataBox xVals(nFineRho, nFineT); - DataBox yVals(nFineRho, nFineT); - DataBox vars(nFineRho, nFineT); - DataBox dx(nFineRho, nFineT); - DataBox dy(nFineRho, nFineT); - DataBox pressEOSPAC(nFineRho, nFineT); - DataBox pressDiff_h(nFineRho, nFineT); - DataBox pressDiff_d(nFineRho, nFineT); + DataBox xVals(nFineRho, nFineT); + DataBox yVals(nFineRho, nFineT); + DataBox vars(nFineRho, nFineT); + DataBox dx(nFineRho, nFineT); + DataBox dy(nFineRho, nFineT); + DataBox pressEOSPAC(nFineRho, nFineT); + DataBox pressDiff_h(nFineRho, nFineT); + DataBox pressDiff_d(nFineRho, nFineT); #ifdef PORTABILITY_STRATEGY_KOKKOS using RView = Kokkos::View; RView pressSpinerDView("pressSpinerDevice", nFineRho * nFineT); typename RView::HostMirror pressSpinerHView = Kokkos::create_mirror_view(pressSpinerDView); - DataBox pressSpiner_d(pressSpinerDView.data(), nFineRho, nFineT); - DataBox pressSpiner_hm(pressSpinerHView.data(), nFineRho, nFineT); + DataBox pressSpiner_d(pressSpinerDView.data(), nFineRho, nFineT); + DataBox pressSpiner_hm(pressSpinerHView.data(), nFineRho, nFineT); #else - DataBox pressSpiner_d(nFineRho, nFineT); - DataBox pressSpiner_hm = pressSpiner_d.slice(2, 0, nFineRho); + DataBox pressSpiner_d(nFineRho, nFineT); + DataBox pressSpiner_hm = pressSpiner_d.slice(2, 0, nFineRho); #endif - DataBox pressSpiner_h(nFineRho, nFineT); - DataBox tempSpiner_h(nFineRho, nFineT); - DataBox tempSpinerE_h(nFineRho, nFineT); - DataBox tempEOSPAC(nFineRho, nFineT); - DataBox tempDiff_h(nFineRho, nFineT); - DataBox tempDiff_d(nFineRho, nFineT); - DataBox tempDiffE_h(nFineRho, nFineT); - DataBox tempDiffE_d(nFineRho, nFineT); + DataBox pressSpiner_h(nFineRho, nFineT); + DataBox tempSpiner_h(nFineRho, nFineT); + DataBox tempSpinerE_h(nFineRho, nFineT); + DataBox tempEOSPAC(nFineRho, nFineT); + DataBox tempDiff_h(nFineRho, nFineT); + DataBox tempDiff_d(nFineRho, nFineT); + DataBox tempDiffE_h(nFineRho, nFineT); + DataBox tempDiffE_d(nFineRho, nFineT); #ifdef PORTABILITY_STRATEGY_KOKKOS RView tempSpinerDView("tempSpinerDevice", nFineRho * nFineT); RView tempSpinerDViewE("tempSpinerSieDevice", nFineRho * nFineT); RView::HostMirror tempSpinerHView = Kokkos::create_mirror_view(tempSpinerDView); RView::HostMirror tempSpinerHViewE = Kokkos::create_mirror_view(tempSpinerDViewE); - DataBox tempSpiner_d(tempSpinerDView.data(), nFineRho, nFineT); - DataBox tempSpiner_hm(tempSpinerHView.data(), nFineRho, nFineT); - DataBox tempSpinerE_d(tempSpinerDViewE.data(), nFineRho, nFineT); - DataBox tempSpinerE_hm(tempSpinerHViewE.data(), nFineRho, nFineT); + DataBox tempSpiner_d(tempSpinerDView.data(), nFineRho, nFineT); + DataBox tempSpiner_hm(tempSpinerHView.data(), nFineRho, nFineT); + DataBox tempSpinerE_d(tempSpinerDViewE.data(), nFineRho, nFineT); + DataBox tempSpinerE_hm(tempSpinerHViewE.data(), nFineRho, nFineT); RView rhos_v("rhos", nFineRho); RView Ts_v("Ts", nFineT); RView sies_v("sies", nFineT); - DataBox rhos(rhos_v.data(), nFineRho); - DataBox Ts(Ts_v.data(), nFineT); - DataBox sies(sies_v.data(), nFineT); + DataBox rhos(rhos_v.data(), nFineRho); + DataBox Ts(Ts_v.data(), nFineT); + DataBox sies(sies_v.data(), nFineT); #else - DataBox tempSpiner_d(nFineRho, nFineT); - DataBox tempSpiner_hm = tempSpiner_d.slice(2, 0, nFineRho); - DataBox tempSpinerE_d(nFineRho, nFineT); - DataBox tempSpinerE_hm = tempSpiner_h.slice(2, 0, nFineRho); - DataBox rhos(nFineRho); - DataBox Ts(nFineT); - DataBox sies(nFineT); + DataBox tempSpiner_d(nFineRho, nFineT); + DataBox tempSpiner_hm = tempSpiner_d.slice(2, 0, nFineRho); + DataBox tempSpinerE_d(nFineRho, nFineT); + DataBox tempSpinerE_hm = tempSpiner_h.slice(2, 0, nFineRho); + DataBox rhos(nFineRho); + DataBox Ts(nFineT); + DataBox sies(nFineT); #endif std::cout << "Comparing to EOSPAC for materials..." << std::endl; @@ -242,8 +242,8 @@ int main(int argc, char *argv[]) { (Real *)PORTABLE_MALLOC(sizeof(Real) * nFineRho * nFineT * eos_host.nlambda()); Real *lambda_hp = (Real *)malloc(sizeof(Real) * nFineRho * nFineT * eos_host.nlambda()); - DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); - DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); + DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); + DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); Real rhoMin = eosE_host.rhoMin(); Real rhoMax = eosE_host.rhoMax(); @@ -464,10 +464,10 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTIMES; n++) { for (int j = 0; j < nFineRho; j++) { Real rho = lRhoBounds.i2lin(j); - DataBox lambda_hj = lambda_h.slice(j); + DataBox lambda_hj = lambda_h.slice(j); for (int i = 0; i < nFineT; i++) { Real sie = leBounds.i2lin(i); - DataBox lambda = lambda_hj.slice(i); + DataBox lambda = lambda_hj.slice(i); tempSpiner_h(j, i) = eos_host.TemperatureFromDensityInternalEnergy(rho, sie, lambda.data()); } @@ -487,7 +487,7 @@ int main(int argc, char *argv[]) { PORTABLE_LAMBDA(const int &j, const int &i) { Real rho = rhos(j); Real sie = sies(i); - DataBox lambda = lambda_d.slice(j).slice(i); + DataBox lambda = lambda_d.slice(j).slice(i); tempSpiner_d(j, i) = eos_device.TemperatureFromDensityInternalEnergy( rho, sie, lambda.data()); }); From b720356325376c7639fe4a71ca1d29f1c257f3be Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 12:14:41 -0600 Subject: [PATCH 12/45] play with bounds --- singularity-eos/base/table_bounds.hpp | 16 ++++++++++++++-- test/test_bounds.cpp | 6 +++--- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 3e3aeb82de..0dd2a3d3eb 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -131,10 +131,22 @@ class Bounds { return os; } + // This uses real logs template static int getNumPointsFromPPD(Real min, Real max, T ppd) { - Bounds b(min, max, 3, true); - Real ndecades = b.grid.max() - b.grid.min(); + constexpr Real epsilon = std::numeric_limits::epsilon(); + const Real min_offset = 10 * std::abs(epsilon); + // 1.1 so that the y-intercept isn't identically zero + // when min < 0. + // min_offset to handle the case where min=0 + Real offset = 0; + if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; + min += offset; + max += offset; + + Real lmin = std::log10(min); + Real lmax = std::log10(max); + Real ndecades = lmax - lmin; return std::max(2, static_cast(std::ceil(ppd * ndecades))); } diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 7a1f4049d4..f11fc2ec87 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -35,13 +35,13 @@ constexpr Real rho_max = 1e3; constexpr Real T_min = 1; constexpr Real T_max = 1e9; constexpr Real T_split = 1e4; -constexpr int N_per_decade_fine = 10; -constexpr Real N_factor = 2; +constexpr int N_per_decade_fine = 200; +constexpr Real N_factor = 5; SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { WHEN("We compute the number of points from points per decade") { int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); - THEN("We get the right number") { REQUIRE(np == 90); } + THEN("We get the right number") { REQUIRE(np == 1800); } } } From f80bf9f3fba3e2c46675f8e132d66d2d2da5b9bc Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 12:24:23 -0600 Subject: [PATCH 13/45] num points from PPD took a second log --- singularity-eos/base/table_bounds.hpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 0dd2a3d3eb..7e97c016d7 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -82,11 +82,11 @@ class Bounds { splitPoint = singularity::FastMath::log10(std::abs(splitPoint)); // add a point just to make sure we have enough points after adjusting for anchor - int N_fine = getNumPointsFromPPD(global_min, splitPoint, ppd_fine) + 1; + int N_fine = getNumPointsFromDensity(global_min, splitPoint, ppd_fine) + 1; adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); RegularGrid1D grid_lower(global_min, splitPoint, N_fine); - const int N_upper = getNumPointsFromPPD(splitPoint, global_max, ppd_coarse); + const int N_upper = getNumPointsFromDensity(splitPoint, global_max, ppd_coarse); RegularGrid1D grid_upper(splitPoint, global_max, N_upper); grid = Grid_t(std::vector{grid_lower, grid_upper}); @@ -104,16 +104,16 @@ class Bounds { Real mid_max = anchor_point + 0.5 * log_fine_diameter; // add a point just to make sure we have enough points after adjusting for anchor - int N_fine = getNumPointsFromPPD(mid_min, mid_max, ppd_fine) + 1; + int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); RegularGrid1D grid_middle(mid_min, mid_max, N_fine); - const int N_lower = getNumPointsFromPPD(global_min, mid_min, ppd_coarse); + const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_coarse); RegularGrid1D grid_lower(global_min, mid_min, N_lower); - const int N_upper = getNumPointsFromPPD(mid_max, global_max, ppd_coarse); + const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_coarse); RegularGrid1D grid_upper(mid_max, global_max, N_upper); - + printf("N_lower, N_mid, N_upper = %d %d %d\n", N_lower, N_fine, N_upper); grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); } @@ -146,8 +146,14 @@ class Bounds { Real lmin = std::log10(min); Real lmax = std::log10(max); - Real ndecades = lmax - lmin; - return std::max(2, static_cast(std::ceil(ppd * ndecades))); + int N = getNumPointsFromDensity(lmin, lmax, ppd); + return N; + } + template + static int getNumPointsFromDensity(Real min, Real max, T density) { + Real delta = max - min; + int N = std::max(2, static_cast(std::ceil(density * delta))); + return N; } private: From 1037cef17b1838d433157ae7cf604b29a84cff07 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 31 Jul 2024 12:25:39 -0600 Subject: [PATCH 14/45] remove print --- singularity-eos/base/table_bounds.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 7e97c016d7..355040af8d 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -113,7 +113,7 @@ class Bounds { const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_coarse); RegularGrid1D grid_upper(mid_max, global_max, N_upper); - printf("N_lower, N_mid, N_upper = %d %d %d\n", N_lower, N_fine, N_upper); + grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); } From 2f4a61e4e8c103d7f318a798379be02188820bd0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 1 Aug 2024 16:38:56 -0600 Subject: [PATCH 15/45] it all seems to work pretty nicely --- sesame2spiner/examples/aluminum.dat | 21 ++++++++++ sesame2spiner/generate_files.cpp | 58 ++++++++++++++++++--------- sesame2spiner/generate_files.hpp | 13 ++++-- sesame2spiner/io_eospac.cpp | 2 +- singularity-eos/base/table_bounds.hpp | 29 +++++++++----- test/compare_to_eospac.cpp | 49 ++++++++++++---------- test/test_bounds.cpp | 2 +- 7 files changed, 118 insertions(+), 56 deletions(-) create mode 100644 sesame2spiner/examples/aluminum.dat diff --git a/sesame2spiner/examples/aluminum.dat b/sesame2spiner/examples/aluminum.dat new file mode 100644 index 0000000000..318cfc5706 --- /dev/null +++ b/sesame2spiner/examples/aluminum.dat @@ -0,0 +1,21 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 2021-2023. 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. All rights in the program are +# reserved by Triad National Security, LLC, and the U.S. Department of +# Energy/National Nuclear Security Administration. The Government is +# granted for itself and others acting on its behalf a nonexclusive, +# paid-up, irrevocable worldwide license in this material to reproduce, +# prepare derivative works, distribute copies to the public, perform +# publicly and display publicly, and to permit others to do so. +# ====================================================================== + +matid=3720 +name=aluminum +rhomin = 1e-5 +Tmin = 1 diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 8351ddff39..f48faf65e1 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -242,7 +242,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real rhoMin = params.Get("rhomin", TinyShift(metadata.rhoMin, 1)); Real rhoMax = params.Get("rhomax", metadata.rhoMax); Real TMin = params.Get("Tmin", TinyShift(metadata.TMin, 1)); - Real TMax = params.Get("Tmax", metadata.TMax); + Real TMax = params.Get("Tmax", metadata.TMax); Real sieMin = params.Get("siemin", TinyShift(metadata.sieMin, 1)); Real sieMax = params.Get("siemax", metadata.sieMax); @@ -253,9 +253,9 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params checkValInMatBounds(matid, "sieMin", sieMin, metadata.sieMin, metadata.sieMax); checkValInMatBounds(matid, "sieMax", sieMax, metadata.sieMin, metadata.sieMax); - Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.0); - Real shrinklTBounds = params.Get("shrinklTBounds", 0.0); - Real shrinkleBounds = params.Get("shrinkleBounds", 0.0); + Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.); + Real shrinklTBounds = params.Get("shrinklTBounds", 0.); + Real shrinkleBounds = params.Get("shrinkleBounds", 0.); shrinklRhoBounds = std::min(1., std::max(shrinklRhoBounds, 0.)); shrinklTBounds = std::min(1., std::max(shrinklTBounds, 0.)); @@ -287,16 +287,20 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params int numT = params.Get("numT", numTDefault); int numSie = params.Get("numsie", numSieDefault); + const Real TAnchor = 298.15; const Real rhoAnchor = metadata.normalDensity; - constexpr Real TAnchor = 298.15; // Piecewise grids stuff - bool do_piecewise_grids = params.Get("piecewise", true); - Real ppd_factor_rho = params.Get("rhoCoarseFactor", COARSE_FACTOR_DEFAULT); - Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT); - Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT); - Real rho_fine_diameter = params.Get("rhoFineDiameterDecades", 0.5); - Real TSplitPoint = params.Get("TSplitPoint", 1e4); + bool piecewiseRho = params.Get("piecewiseRho", true); + bool piecewiseT = params.Get("piecewiseT", true); + bool piecewiseSie = params.Get("piecewiseSie", true); + + Real ppd_factor_rho_lo = params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); + Real ppd_factor_rho_hi = params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); + Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); + Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); + Real rho_fine_diameter = params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); + Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); Real rho_fine_center = rhoAnchor; @@ -306,15 +310,23 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params // Extrapolation and other resolution tricks will be explored in the // future. - if (rhoMin < STRICTLY_POS_MIN) rhoMin = STRICTLY_POS_MIN; - if (TMin < STRICTLY_POS_MIN) TMin = STRICTLY_POS_MIN; + if (rhoMin < STRICTLY_POS_MIN_RHO) rhoMin = STRICTLY_POS_MIN_RHO; + if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; - if (do_piecewise_grids) { + if (piecewiseRho) { lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, - rho_fine_diameter, ppdRho, ppd_factor_rho, shrinklRhoBounds); + rho_fine_diameter, ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, + shrinklRhoBounds); + } else { + lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); + } + if (piecewiseT) { lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, ppd_factor_T, shrinklTBounds); - + } else { + lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); + } + if (piecewiseSie) { // compute temperature as a reasonable anchor point constexpr int NT = 1; constexpr EOS_INTEGER nXYPairs = 2; @@ -331,16 +343,22 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Verbosity::Quiet); eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); } - Real sieAnchor = sieFromSesame(sie[0]); - Real sieSplitPoint = sieFromSesame(sie[1]); + Real sieAnchor = sie[0]; + Real sieSplitPoint = sie[1]; leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, ppdSie, ppd_factor_sie, shrinkleBounds); } else { - lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); - lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); } + std::cout << "lRho bounds are\n" + << lRhoBounds + << "lT bounds are\n" + << lTBounds + << "lSie bounds are \n" + << leBounds + << std::endl; + return; } diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 6a9ffd4ab2..622f81a99c 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -30,10 +30,15 @@ using namespace EospacWrapper; -constexpr int PPD_DEFAULT_RHO = 200; -constexpr int PPD_DEFAULT_T = 140; -constexpr Real STRICTLY_POS_MIN = 1e-9; -constexpr Real COARSE_FACTOR_DEFAULT = 20; +constexpr int PPD_DEFAULT_RHO = 350; +constexpr int PPD_DEFAULT_T = 100; +constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; +constexpr Real STRICTLY_POS_MIN_T = 1e-2; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; +constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; +constexpr Real RHO_FINE_DIAMETER_DEFAULT = 1.5; +constexpr Real T_SPLIT_POINT_DEFAULT = 1e4; herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds, const Bounds &lTBounds, const Bounds &leBounds, diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/io_eospac.cpp index e074fb8fd0..0953029223 100644 --- a/sesame2spiner/io_eospac.cpp +++ b/sesame2spiner/io_eospac.cpp @@ -30,7 +30,7 @@ void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds DataBox &mask, Verbosity eospacWarn) { using namespace EospacWrapper; - constexpr int NT = 4; + constexpr int NT = 3; constexpr EOS_INTEGER nXYPairs = 1; EOS_INTEGER tableHandle[NT]; EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 7e97c016d7..9984bee916 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -43,16 +43,16 @@ class Bounds { using TwoGrids = std::integral_constant; using ThreeGrids = std::integral_constant; - Bounds() {} + Bounds() : offset(0), piecewise(false), linmin_(0), linmax_(0) {} Bounds(Real min, Real max, int N, Real offset) : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), - offset(offset), piecewise(false) {} + offset(offset), piecewise(false), linmin_(min), linmax_(max) {} Bounds(OneGrid, Real min, Real max, int N, Real offset) : Bounds(min, max, N, offset) {} Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) - : offset(0), piecewise(true) { + : offset(0), piecewise(true), linmin_(min), linmax_(max) { if (convertToLog) { convertBoundsToLog_(min, max, shrinkRange); if (!(std::isnan(anchor_point))) { @@ -72,7 +72,7 @@ class Bounds { Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) - : offset(0), piecewise(true) { + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; convertBoundsToLog_(global_min, global_max, shrinkRange); @@ -93,8 +93,10 @@ class Bounds { } Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, - Real log_fine_diameter, Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) { - const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; + Real log_fine_diameter, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, + Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; + const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; convertBoundsToLog_(global_min, global_max, shrinkRange); anchor_point += offset; @@ -108,16 +110,23 @@ class Bounds { adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); RegularGrid1D grid_middle(mid_min, mid_max, N_fine); - const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_coarse); + const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_lo); RegularGrid1D grid_lower(global_min, mid_min, N_lower); - const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_coarse); + const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_hi); RegularGrid1D grid_upper(mid_max, global_max, N_upper); printf("N_lower, N_mid, N_upper = %d %d %d\n", N_lower, N_fine, N_upper); grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); } - inline Real log2lin(Real xl) const { return singularity::FastMath::pow10(xl) - offset; } + inline Real log2lin(Real xl) const { + // JMM: Need to guard this with the linear bounds passed in. The + // reason is that our fast math routines, while completely + // invertible at the level of machine epsilon, do introduce error + // at machine epsilon, which can bring us out of the interpolation + // range of eospac. + return std::min(linmax_, std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); + } inline Real i2lin(int i) const { return log2lin(grid.x(i)); } friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { @@ -199,6 +208,8 @@ class Bounds { Grid_t grid; Real offset = 0; bool piecewise = false; + private: + Real linmin_, linmax_; }; } // namespace singularity diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 959abc85b7..12759470f4 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -225,6 +225,7 @@ int main(int argc, char *argv[]) { SesameMetadata metadata; eosGetMetadata(matid, metadata, Verbosity::Debug); + std::cout << metadata << std::endl; hid_t idGroup = H5Gcreate(file, std::to_string(matid).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); @@ -245,12 +246,12 @@ int main(int argc, char *argv[]) { DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); - Real rhoMin = eosE_host.rhoMin(); - Real rhoMax = eosE_host.rhoMax(); - Real TMin = eosE_host.TMin(); - Real TMax = eosE_host.TMax(); - Real sieMin = eosE_host.sieMin(); - Real sieMax = eosE_host.sieMax(); + Real rhoMin = 1.1*std::max(metadata.rhoMin, 1e-5); + Real rhoMax = 0.9*metadata.rhoMax; + Real TMin = 1.1*std::max(metadata.TMin, 1.0); + Real TMax = 0.9*metadata.TMax; + Real sieMin = metadata.sieMin + 0.1*std::abs(metadata.sieMin); + Real sieMax = 0.9*metadata.sieMax; Bounds lRhoBounds(rhoMin, rhoMax, nFineRho); Bounds lTBounds(TMin, TMax, nFineT); @@ -355,10 +356,11 @@ int main(int argc, char *argv[]) { const Real p_true = pressureFromSesame(vars(j, i)); pressEOSPAC(j, i) = p_true; const Real diff = pressSpiner_h(j, i) - p_true; - pressDiff_h(j, i) = diff; - diffL2 += diff * diff; const Real mean_p = 0.5 * diff + p_true; - L2 += mean_p * mean_p; + pressDiff_h(j, i) = diff; + const Real reldiff = diff / (1e-10 + std::abs(mean_p)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -382,10 +384,11 @@ int main(int argc, char *argv[]) { // pressSpiner_d has been copied into pressSpiner_h const Real p_true = pressEOSPAC(j, i); const Real diff = pressSpiner_hm(j, i) - p_true; - pressDiff_d(j, i) = diff; const Real mean_p = 0.5 * diff + p_true; - diffL2 += diff * diff; - L2 += mean_p * mean_p; + pressDiff_d(j, i) = diff; + const Real reldiff = diff / (1e-10 + std::abs(mean_p)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -423,7 +426,7 @@ int main(int argc, char *argv[]) { Real rho = lRhoBounds.i2lin(j); for (int i = 0; i < nFineT; i++) { Real sie = leBounds.i2lin(i); - xVals(j, i) = rho; + xVals(j, i) = densityToSesame(rho); yVals(j, i) = sieToSesame(sie); } } @@ -558,8 +561,9 @@ int main(int argc, char *argv[]) { diff = 0; } #endif - diffL2 += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2 += reldiff * reldiff; + L2 += 1; } } diffL2 /= L2; @@ -572,8 +576,9 @@ int main(int argc, char *argv[]) { const Real diff = tempSpinerE_h(j, i) - t_true; tempDiffE_h(j, i) = diff; const Real mean_t = 0.5 * diff + t_true; - diffL2E += diff * diff; - L2 += mean_t * mean_t; + const Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2E += reldiff * reldiff; + L2 += 1; } } diffL2E /= L2; @@ -612,8 +617,9 @@ int main(int argc, char *argv[]) { diff = 0; } #endif - diffL2_d += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2_d += reldiff * reldiff; + L2 += 1; } } diffL2_d /= L2; @@ -626,8 +632,9 @@ int main(int argc, char *argv[]) { const Real diff = tempSpinerE_hm(j, i) - t_true; tempDiffE_d(j, i) = diff; const Real mean_t = 0.5 * diff + t_true; - diffL2E_d += diff * diff; - L2 += mean_t * mean_t; + Real reldiff = diff / (1e-10 + std::abs(mean_t)); + diffL2E_d += reldiff * reldiff; + L2 += 1; } } diffL2E_d /= L2; diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index f11fc2ec87..68ee3b037d 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -83,7 +83,7 @@ SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { WHEN("We compute a piecewise bounds object with three grids") { Bounds bnds(Bounds::ThreeGrids(), rho_min, rho_max, rho_normal, 0.5, - N_per_decade_fine, N_factor); + N_per_decade_fine, N_factor, N_factor); THEN("The bounds are right") { Real lrmin = singularity::FastMath::log10(rho_min); Real lrmax = singularity::FastMath::log10(rho_max); From 9c9d88fac1066e1fa89600ea601baffeba980b16 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 1 Aug 2024 16:39:46 -0600 Subject: [PATCH 16/45] formatting --- sesame2spiner/generate_files.cpp | 20 +++++++++----------- singularity-eos/base/table_bounds.hpp | 11 +++++++---- test/compare_to_eospac.cpp | 12 ++++++------ 3 files changed, 22 insertions(+), 21 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index f48faf65e1..15cd2a397e 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -242,7 +242,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real rhoMin = params.Get("rhomin", TinyShift(metadata.rhoMin, 1)); Real rhoMax = params.Get("rhomax", metadata.rhoMax); Real TMin = params.Get("Tmin", TinyShift(metadata.TMin, 1)); - Real TMax = params.Get("Tmax", metadata.TMax); + Real TMax = params.Get("Tmax", metadata.TMax); Real sieMin = params.Get("siemin", TinyShift(metadata.sieMin, 1)); Real sieMax = params.Get("siemax", metadata.sieMax); @@ -299,7 +299,8 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real ppd_factor_rho_hi = params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); - Real rho_fine_diameter = params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); + Real rho_fine_diameter = + params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); Real rho_fine_center = rhoAnchor; @@ -314,9 +315,9 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; if (piecewiseRho) { - lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, - rho_fine_diameter, ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, - shrinklRhoBounds); + lRhoBounds = + Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, + ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, shrinklRhoBounds); } else { lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); } @@ -352,12 +353,9 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params } std::cout << "lRho bounds are\n" - << lRhoBounds - << "lT bounds are\n" - << lTBounds - << "lSie bounds are \n" - << leBounds - << std::endl; + << lRhoBounds << "lT bounds are\n" + << lTBounds << "lSie bounds are \n" + << leBounds << std::endl; return; } diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 10666e924c..dbb957df58 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -52,7 +52,7 @@ class Bounds { Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) - : offset(0), piecewise(true), linmin_(min), linmax_(max) { + : offset(0), piecewise(true), linmin_(min), linmax_(max) { if (convertToLog) { convertBoundsToLog_(min, max, shrinkRange); if (!(std::isnan(anchor_point))) { @@ -72,7 +72,7 @@ class Bounds { Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) - : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; convertBoundsToLog_(global_min, global_max, shrinkRange); @@ -94,7 +94,8 @@ class Bounds { Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + Real shrinkRange = 0) + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; @@ -125,7 +126,8 @@ class Bounds { // invertible at the level of machine epsilon, do introduce error // at machine epsilon, which can bring us out of the interpolation // range of eospac. - return std::min(linmax_, std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); + return std::min(linmax_, + std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); } inline Real i2lin(int i) const { return log2lin(grid.x(i)); } @@ -208,6 +210,7 @@ class Bounds { Grid_t grid; Real offset = 0; bool piecewise = false; + private: Real linmin_, linmax_; }; diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 12759470f4..4ff052c555 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -246,12 +246,12 @@ int main(int argc, char *argv[]) { DataBox lambda_h(lambda_hp, nFineRho, nFineT, eos_host.nlambda()); DataBox lambda_d(lambda_dp, nFineRho, nFineT, eos_host.nlambda()); - Real rhoMin = 1.1*std::max(metadata.rhoMin, 1e-5); - Real rhoMax = 0.9*metadata.rhoMax; - Real TMin = 1.1*std::max(metadata.TMin, 1.0); - Real TMax = 0.9*metadata.TMax; - Real sieMin = metadata.sieMin + 0.1*std::abs(metadata.sieMin); - Real sieMax = 0.9*metadata.sieMax; + Real rhoMin = 1.1 * std::max(metadata.rhoMin, 1e-5); + Real rhoMax = 0.9 * metadata.rhoMax; + Real TMin = 1.1 * std::max(metadata.TMin, 1.0); + Real TMax = 0.9 * metadata.TMax; + Real sieMin = metadata.sieMin + 0.1 * std::abs(metadata.sieMin); + Real sieMax = 0.9 * metadata.sieMax; Bounds lRhoBounds(rhoMin, rhoMax, nFineRho); Bounds lTBounds(TMin, TMax, nFineT); From bb4037c408f90d4ac252a8dd17f96e786cafd3fa Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 1 Aug 2024 16:58:19 -0600 Subject: [PATCH 17/45] add documentation --- doc/sphinx/src/models.rst | 87 ++++++++++++++++++++++++++++++++++----- 1 file changed, 76 insertions(+), 11 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 205aa421c6..237f3e8fc0 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1471,8 +1471,7 @@ key-value pairs. For exampe the following input deck is for air: matid = 5030 # These set the number of grid points per decade - # for each variable. The default is 50 points - # per decade. + # for each variable. numrho/decade = 40 numT/decade = 40 numSie/decade = 40 @@ -1491,16 +1490,82 @@ key-value pairs. For exampe the following input deck is for air: shrinklTBounds = 0.15 shrinkleBounds = 0.5 -The only required value in an input file is the matid, in this -case 5030. All other values will be inferred from the original sesame -database if possible and if no value in the input file is -provided. Comments are prefixed with ``#``. +Comments are prefixed with ``#``. `eospac`_ uses environment variables +and files to locate files in the `sesame`_ database, and +``sesame2spiner`` uses `eospac`_. So the location of the ``sesame`` +database need not be provided by the command line. For how to specify +`sesame`_ file locations, see the `eospac`_ manual. -`eospac`_ uses environment variables and files to locate files in the -`sesame`_ database, and ``sesame2spiner`` uses `eospac`_. So the -location of the ``sesame`` database need not be provided by the -command line. For how to specify `sesame`_ file locations, see the -`eospac`_ manual. +Piecewise Spiner Grids +```````````````````````` + +``SpinerEOS`` models and ``sesame2spiner`` also support grids with +different resolutions in different parts of the table. We call these +**piecewise** grids. This can be enabled or disabled with the +parameters + +.. code-block:: + + piecewiseRho = true + piecewiseT = true + piecewiseSie = true + +These options may be true or false. The default is true. When +piecewise grids are active, the density grid gets split into three +pieces, a piece of diameter ``rhoFineDiameterDecades`` (in log space) +around the material's normal density and a piece above and below. The +``numrho/decade`` parameter sets the number of points per decade in +this central refined region. The regions at lower and higher density +have ``rhoCoarseFactorLo`` and ``rhoCoarseFactorHi`` fewer points per +decade respectively compared to the finer region. + +The temperature grid has two regions, a more finely spaced region at +low temperatures and a more finely spaced region at high +temperatures. The regions are spearated by a temperature +``TSplitPoint``. The default is :math:`10^4`. The energy grid follows +the temperature grid, with the energy split point corresponding to the +temperature split point. The coarser high-temperature temperature and +energy grids are coarsened by a factor of ``TCoarseFactor`` and +``sieCoarseFactor`` respectively. + +Thus the input block for piecewise grid might look like this: + +.. code-block:: + + # defaults are true + piecewiseRho = true + piecewiseT = true + piecewiseSie = true + + # the fine resolution for rho + numrho/decade = 350 + # width of the fine region for rho + rhoFineDiameterDecades = 1.5 + # the lower density region is 3x less refined + rhoCoarseFactorLo = 3 + # the higher density region is 5x less refined + rhoCoarseFactorHi = 5 + + # the fine resolution for T + numT/decade = 100 + # the point demarking the coarse and fine regions in temperature + TSplitPoint = 1e4 + # it's usually wise to to not let + # temperature get too small in log space if you do this + Tmin = 1 + # The coarser region (above the split point) is 50 percent less refined + TCoarseFactor = 1.5 + + # energy has the split point sie(rhonormal, TSplitPoint) + # but we may still specify the resolution + numSie/decade = 100 + sieCoarseFactor = 1.5 + +.. note:: + + For all grid types, the only required value in an input file is the + matid. All other values will be inferred from the original sesame + database if possible and if no value in the input file is provided. SAP Polynomial EOS `````````````````` From 42fab4451fd6f759e5ff6a75ad1087c0f4d59530 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 1 Aug 2024 17:18:36 -0600 Subject: [PATCH 18/45] CHANGELOG --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9340c409c0..f969865a7d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,8 +3,10 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Piecewise grids for Spiner EOS. ### Fixed (Repair bugs, etc) +- [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Includes a fix for extrapolation of specific internal energy in SpinerEOS. ### Changed (changing behavior/API/variables/...) From ba5c89937b23d09e538773aea37869ff869b4db0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 1 Aug 2024 17:19:41 -0600 Subject: [PATCH 19/45] CC --- sesame2spiner/generate_files.cpp | 2 +- sesame2spiner/generate_files.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 15cd2a397e..50b9e1be9d 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 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 diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 622f81a99c..914ec62ed8 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 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 From 79ce8c9f4cba80d825e5966f50ea17092140e7ff Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 12:20:21 -0600 Subject: [PATCH 20/45] add robustness for refined region out of table bounds or normal density ill defined --- sesame2spiner/generate_files.cpp | 8 ++++- singularity-eos/base/table_bounds.hpp | 43 +++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index f48faf65e1..103a0c7a7d 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -288,7 +288,13 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params int numSie = params.Get("numsie", numSieDefault); const Real TAnchor = 298.15; - const Real rhoAnchor = metadata.normalDensity; + Real rhoAnchor = metadata.normalDensity; + if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { + std::cerr << "WARNING [" << matid << "] " + << "normal density ill defined. Setting it to a sensible default." + << std::endl; + rhoAnchor = 1; + } // Piecewise grids stuff bool piecewiseRho = params.Get("piecewiseRho", true); diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 10666e924c..431af1a1a0 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -81,6 +81,27 @@ class Bounds { splitPoint += offset; splitPoint = singularity::FastMath::log10(std::abs(splitPoint)); + if (splitPoint <= global_min) { + PORTABLE_ALWAYS_WARN("Split point less than global minimum. Adjusting."); + Real eps = 0.1*std::abs(global_min); + splitPoint = global_min + eps; + } + if (splitPoint >= global_max) { + PORTABLE_ALWAYS_WARN("Split point greater than global maximum. Adjusting."); + Real eps = 0.1*std::abs(global_max); + splitPoint = global_max - eps; + } + if (anchor_point <= global_min) { + PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); + Real eps = 0.1*std::abs(global_min); + anchor_point = global_min + eps; + } + if (anchor_point >= splitPoint) { + PORTABLE_ALWAYS_WARN("Anchor point greater than split point. Adjusting."); + Real eps = 0.1*std::abs(splitPoint); + anchor_point = splitPoint - eps; + } + // add a point just to make sure we have enough points after adjusting for anchor int N_fine = getNumPointsFromDensity(global_min, splitPoint, ppd_fine) + 1; adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); @@ -102,9 +123,31 @@ class Bounds { anchor_point += offset; anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); + if (anchor_point <= global_min) { + PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); + Real eps = 0.1*std::abs(global_min); + anchor_point = global_min + eps; + } + if (anchor_point >= global_max) { + PORTABLE_ALWAYS_WARN("Anchor point greater than global maximum. Adjusting."); + Real eps = 0.1*std::abs(global_max); + anchor_point = global_max - eps; + } + Real mid_min = anchor_point - 0.5 * log_fine_diameter; Real mid_max = anchor_point + 0.5 * log_fine_diameter; + if (mid_min <= global_min) { + PORTABLE_ALWAYS_WARN("Table bounds refined minimum lower than global minimum. Adjusting."); + Real delta = std::abs(anchor_point - global_min); + mid_min = anchor_point - 0.9*delta; + } + if (mid_max >= global_max) { + PORTABLE_ALWAYS_WARN("Table bounds refined maximum greater than global maximum. Adjusting."); + Real delta = std::abs(global_max - anchor_point); + mid_max = anchor_point + 0.9*delta; + } + // add a point just to make sure we have enough points after adjusting for anchor int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); From 2ab925566e81fda43d63e40e3e7477cf16c0dd79 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 12:21:17 -0600 Subject: [PATCH 21/45] add helium, which is another ill-behaved one --- sesame2spiner/examples/helium.dat | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 sesame2spiner/examples/helium.dat diff --git a/sesame2spiner/examples/helium.dat b/sesame2spiner/examples/helium.dat new file mode 100644 index 0000000000..2668c9eceb --- /dev/null +++ b/sesame2spiner/examples/helium.dat @@ -0,0 +1,18 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 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. All rights in the program are +# reserved by Triad National Security, LLC, and the U.S. Department of +# Energy/National Nuclear Security Administration. The Government is +# granted for itself and others acting on its behalf a nonexclusive, +# paid-up, irrevocable worldwide license in this material to reproduce, +# prepare derivative works, distribute copies to the public, perform +# publicly and display publicly, and to permit others to do so. +# ====================================================================== +matid = 5762 +name = helium From 1aba45f758e1f3b18594c1434fbe08813df6db02 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 12:21:51 -0600 Subject: [PATCH 22/45] format --- singularity-eos/base/table_bounds.hpp | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 0e67520956..1a8f62a248 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -83,22 +83,22 @@ class Bounds { if (splitPoint <= global_min) { PORTABLE_ALWAYS_WARN("Split point less than global minimum. Adjusting."); - Real eps = 0.1*std::abs(global_min); + Real eps = 0.1 * std::abs(global_min); splitPoint = global_min + eps; } if (splitPoint >= global_max) { PORTABLE_ALWAYS_WARN("Split point greater than global maximum. Adjusting."); - Real eps = 0.1*std::abs(global_max); + Real eps = 0.1 * std::abs(global_max); splitPoint = global_max - eps; } if (anchor_point <= global_min) { PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); - Real eps = 0.1*std::abs(global_min); + Real eps = 0.1 * std::abs(global_min); anchor_point = global_min + eps; } if (anchor_point >= splitPoint) { PORTABLE_ALWAYS_WARN("Anchor point greater than split point. Adjusting."); - Real eps = 0.1*std::abs(splitPoint); + Real eps = 0.1 * std::abs(splitPoint); anchor_point = splitPoint - eps; } @@ -126,12 +126,12 @@ class Bounds { if (anchor_point <= global_min) { PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); - Real eps = 0.1*std::abs(global_min); + Real eps = 0.1 * std::abs(global_min); anchor_point = global_min + eps; } if (anchor_point >= global_max) { PORTABLE_ALWAYS_WARN("Anchor point greater than global maximum. Adjusting."); - Real eps = 0.1*std::abs(global_max); + Real eps = 0.1 * std::abs(global_max); anchor_point = global_max - eps; } @@ -139,14 +139,16 @@ class Bounds { Real mid_max = anchor_point + 0.5 * log_fine_diameter; if (mid_min <= global_min) { - PORTABLE_ALWAYS_WARN("Table bounds refined minimum lower than global minimum. Adjusting."); + PORTABLE_ALWAYS_WARN( + "Table bounds refined minimum lower than global minimum. Adjusting."); Real delta = std::abs(anchor_point - global_min); - mid_min = anchor_point - 0.9*delta; + mid_min = anchor_point - 0.9 * delta; } if (mid_max >= global_max) { - PORTABLE_ALWAYS_WARN("Table bounds refined maximum greater than global maximum. Adjusting."); + PORTABLE_ALWAYS_WARN( + "Table bounds refined maximum greater than global maximum. Adjusting."); Real delta = std::abs(global_max - anchor_point); - mid_max = anchor_point + 0.9*delta; + mid_max = anchor_point + 0.9 * delta; } // add a point just to make sure we have enough points after adjusting for anchor From 1dbb8fe695e01c426be99ce9766dd8837788a6f0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 12:51:10 -0600 Subject: [PATCH 23/45] chad comments part 1 --- doc/sphinx/src/models.rst | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 237f3e8fc0..8b7243e87e 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1499,16 +1499,16 @@ database need not be provided by the command line. For how to specify Piecewise Spiner Grids ```````````````````````` -``SpinerEOS`` models and ``sesame2spiner`` also support grids with -different resolutions in different parts of the table. We call these -**piecewise** grids. This can be enabled or disabled with the -parameters +``sesame2spiner`` also supports grids with different resolutions in +different parts of the table. We call these **piecewise** grids. By +default grids are now piecewise. Piecewise grids can be disabled with .. code-block:: - piecewiseRho = true - piecewiseT = true - piecewiseSie = true + # defaults are true + piecewiseRho = false + piecewiseT = false + piecewiseSie = dalse These options may be true or false. The default is true. When piecewise grids are active, the density grid gets split into three @@ -1520,24 +1520,24 @@ have ``rhoCoarseFactorLo`` and ``rhoCoarseFactorHi`` fewer points per decade respectively compared to the finer region. The temperature grid has two regions, a more finely spaced region at -low temperatures and a more finely spaced region at high +low temperatures and a less finely spaced region at high temperatures. The regions are spearated by a temperature -``TSplitPoint``. The default is :math:`10^4`. The energy grid follows -the temperature grid, with the energy split point corresponding to the -temperature split point. The coarser high-temperature temperature and -energy grids are coarsened by a factor of ``TCoarseFactor`` and -``sieCoarseFactor`` respectively. +``TSplitPoint``. The default is :math:`10^4` Kelvin. The energy grid +follows the temperature grid, with the energy split point +corresponding to the temperature split point. The coarser +high-temperature temperature and energy grids are coarsened by a +factor of ``TCoarseFactor`` and ``sieCoarseFactor`` respectively. Thus the input block for piecewise grid might look like this: .. code-block:: - # defaults are true + # Below, all right-hand-sides are set to their default values. piecewiseRho = true piecewiseT = true piecewiseSie = true - # the fine resolution for rho + # the fine resolution for rho. numrho/decade = 350 # width of the fine region for rho rhoFineDiameterDecades = 1.5 From 1374d7f7a48247e3702872a6121e15d6ec536404 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 13:23:06 -0600 Subject: [PATCH 24/45] add ability to set rhoFineMin, rhoFineMax --- sesame2spiner/generate_files.cpp | 22 +++- singularity-eos/base/table_bounds.hpp | 159 ++++++++++++++------------ 2 files changed, 106 insertions(+), 75 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index aa48d51377..608ace5e48 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -311,6 +311,16 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real rho_fine_center = rhoAnchor; + // These override the rho center/diameter settings + Real rho_fine_min = params.Get("rhoFineMin", -1); + Real rho_fine_max = params.Get("rhoFineMax", -1); + if (rho_fine_min * rho_fine_max < 0) { + std::cerr << "WARNING [" << matid << "]: " + << "Either rhoFineMin or rhoFineMax is set while the other is still unset." + << " Both must be set to be sensible. Ignoring." << std::endl; + rho_fine_min = rho_fine_max = -1; + } + // Forces density and temperature to be in a region where an offset // is not needed. This improves resolution at low densities and // temperatures. @@ -321,9 +331,15 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; if (piecewiseRho) { - lRhoBounds = - Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, - ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, shrinklRhoBounds); + if (rho_fine_min > 0) { + lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, + rho_fine_min, rho_fine_max, ppdRho, ppd_factor_rho_lo, + ppd_factor_rho_hi, shrinklRhoBounds, true); + } else { + lRhoBounds = + Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, + ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, shrinklRhoBounds); + } } else { lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); } diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 1a8f62a248..e9f4245da8 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -71,37 +71,19 @@ class Bounds { : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, - Real ppd_fine, Real ppd_factor, Real shrinkRange = 0) + Real ppd_fine, Real ppd_factor, Real shrinkRange = 0, bool convertToLog = true) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; - convertBoundsToLog_(global_min, global_max, shrinkRange); - anchor_point += offset; - anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - splitPoint += offset; - splitPoint = singularity::FastMath::log10(std::abs(splitPoint)); - - if (splitPoint <= global_min) { - PORTABLE_ALWAYS_WARN("Split point less than global minimum. Adjusting."); - Real eps = 0.1 * std::abs(global_min); - splitPoint = global_min + eps; - } - if (splitPoint >= global_max) { - PORTABLE_ALWAYS_WARN("Split point greater than global maximum. Adjusting."); - Real eps = 0.1 * std::abs(global_max); - splitPoint = global_max - eps; - } - if (anchor_point <= global_min) { - PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); - Real eps = 0.1 * std::abs(global_min); - anchor_point = global_min + eps; - } - if (anchor_point >= splitPoint) { - PORTABLE_ALWAYS_WARN("Anchor point greater than split point. Adjusting."); - Real eps = 0.1 * std::abs(splitPoint); - anchor_point = splitPoint - eps; + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + toLog_(anchor_point, offset); + toLog_(splitPoint, offset); } + checkInterval_(splitPoint, global_min, global_max, "Split point"); + checkInterval_(anchor_point, global_min, splitPoint, "Anchor point"); + // add a point just to make sure we have enough points after adjusting for anchor int N_fine = getNumPointsFromDensity(global_min, splitPoint, ppd_fine) + 1; adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); @@ -113,56 +95,39 @@ class Bounds { grid = Grid_t(std::vector{grid_lower, grid_upper}); } + Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real fine_min, + Real fine_max, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, + Real shrinkRange = 0, bool convertToLog = true) + : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { + + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + toLog_(anchor_point, offset); + toLog_(fine_min, offset); + toLog_(fine_max, offset); + } + checkInterval_(anchor_point, global_min, global_max, "Anchor point"); + + grid = gridFromIntervals_(ThreeGrids, global_min, global_max, anchor_point, fine_min, + fine_max, ppd_fine, pdd_factor_lo, ppd_factor_hi); + } + Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - Real shrinkRange = 0) + Real shrinkRange = 0, bool convertToLog = true) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { - const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; - const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; - convertBoundsToLog_(global_min, global_max, shrinkRange); - anchor_point += offset; - anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - - if (anchor_point <= global_min) { - PORTABLE_ALWAYS_WARN("Anchor point less than global minimum. Adjusting."); - Real eps = 0.1 * std::abs(global_min); - anchor_point = global_min + eps; - } - if (anchor_point >= global_max) { - PORTABLE_ALWAYS_WARN("Anchor point greater than global maximum. Adjusting."); - Real eps = 0.1 * std::abs(global_max); - anchor_point = global_max - eps; + if (convertToLog) { + convertBoundsToLog_(global_min, global_max, shrinkRange); + toLog_(anchor_point, offset); } + checkInterval_(anchor_point, global_min, global_max, "Anchor point"); Real mid_min = anchor_point - 0.5 * log_fine_diameter; Real mid_max = anchor_point + 0.5 * log_fine_diameter; - if (mid_min <= global_min) { - PORTABLE_ALWAYS_WARN( - "Table bounds refined minimum lower than global minimum. Adjusting."); - Real delta = std::abs(anchor_point - global_min); - mid_min = anchor_point - 0.9 * delta; - } - if (mid_max >= global_max) { - PORTABLE_ALWAYS_WARN( - "Table bounds refined maximum greater than global maximum. Adjusting."); - Real delta = std::abs(global_max - anchor_point); - mid_max = anchor_point + 0.9 * delta; - } - - // add a point just to make sure we have enough points after adjusting for anchor - int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; - adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); - RegularGrid1D grid_middle(mid_min, mid_max, N_fine); - - const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_lo); - RegularGrid1D grid_lower(global_min, mid_min, N_lower); - - const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_hi); - RegularGrid1D grid_upper(mid_max, global_max, N_upper); - - grid = Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); + grid = gridFromIntervals_(ThreeGrids, global_min, global_max, anchor_point, mid_min, + mid_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); } inline Real log2lin(Real xl) const { @@ -213,6 +178,11 @@ class Bounds { } private: + Real toLog_(Real val, Real offset) { + val += offset; + return singularity::FastMath::log10(std::abs(val)); + } + void convertBoundsToLog_(Real &min, Real &max, Real shrinkRange = 0) { // Log scales can't handle negative numbers or exactly zero. To // deal with that, we offset. @@ -223,17 +193,15 @@ class Bounds { // min_offset to handle the case where min=0 if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - min += offset; - max += offset; + toLog_(min, offset); + toLog_(max, offset); - min = singularity::FastMath::log10(std::abs(min)); - max = singularity::FastMath::log10(std::abs(max)); Real delta = max - min; min += 0.5 * shrinkRange * delta; max -= 0.5 * shrinkRange * delta; } - void adjustForAnchor_(Real min, Real &max, int &N, Real anchor_point) { + static void adjustForAnchor_(Real min, Real &max, int &N, Real anchor_point) { if (min < anchor_point && anchor_point < max) { Real Nfrac = (anchor_point - min) / (max - min); PORTABLE_REQUIRE((0 < Nfrac && Nfrac < 1), "anchor in bounds"); @@ -251,6 +219,53 @@ class Bounds { } } + static void checkInterval_(Real &p, Real min, Real max, const std::string &name) { + if (p <= min) { + PORTABLE_ALWAYS_WARN(name + " less than minimum. Adjusting."); + Real eps = 0.1 * std::abs(min); + p = min + eps; + } + if (p >= max) { + PORTABLE_ALWAYS_WARN(name + " greater than maximum. Adjusting."); + Real eps = 0.1 * std::abs(max); + p = max - eps; + } + } + + static Grid_t gridFromIntervals_(ThreeGRids, Real global_min, Real global_max, + Real anchor_point, Real mid_min, Real mid_max, + Real ppd_fine, Real ppd_factor_lo, + Real ppd_factor_hi) { + const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; + const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; + + if (mid_min <= global_min) { + PORTABLE_ALWAYS_WARN( + "Table bounds refined minimum lower than global minimum. Adjusting."); + Real delta = std::abs(anchor_point - global_min); + mid_min = anchor_point - 0.9 * delta; + } + if (mid_max >= global_max) { + PORTABLE_ALWAYS_WARN( + "Table bounds refined maximum greater than global maximum. Adjusting."); + Real delta = std::abs(global_max - anchor_point); + mid_max = anchor_point + 0.9 * delta; + } + + // add a point just to make sure we have enough points after adjusting for anchor + int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; + adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); + RegularGrid1D grid_middle(mid_min, mid_max, N_fine); + + const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_lo); + RegularGrid1D grid_lower(global_min, mid_min, N_lower); + + const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_hi); + RegularGrid1D grid_upper(mid_max, global_max, N_upper); + + return Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); + } + public: Grid_t grid; Real offset = 0; From 058033a06eace16e455d903607bb97f1a764ebc7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 14:02:13 -0600 Subject: [PATCH 25/45] fix bugs --- sesame2spiner/examples/steel.dat | 16 ++++++++++---- sesame2spiner/generate_files.cpp | 8 +++---- singularity-eos/base/table_bounds.hpp | 30 +++++++++++++-------------- test/test_bounds.cpp | 4 ++-- 4 files changed, 33 insertions(+), 25 deletions(-) diff --git a/sesame2spiner/examples/steel.dat b/sesame2spiner/examples/steel.dat index ea8d8f0c2b..77aa6301de 100644 --- a/sesame2spiner/examples/steel.dat +++ b/sesame2spiner/examples/steel.dat @@ -17,10 +17,18 @@ matid = 4272 name = steel -rhomin = 1e-2 +rhomin = 1e-5 Tmin = 1 + +#piecewiseRho = false +#piecewiseT = false +#piecewiseSie = false +#numrho = 1334 +#numT = 738 +#numsie = 634 + # These shrink logarithm of bounds # by a fraction of the total interval <= 1 -shrinklRhoBounds = 0.15, -shrinklTBounds = 0.15, -shrinkleBounds = 0.5 +#shrinklRhoBounds = 0.15, +#shrinklTBounds = 0.15, +#shrinkleBounds = 0.5 diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 608ace5e48..c497334651 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -334,18 +334,18 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params if (rho_fine_min > 0) { lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_min, rho_fine_max, ppdRho, ppd_factor_rho_lo, - ppd_factor_rho_hi, shrinklRhoBounds, true); + ppd_factor_rho_hi, true, shrinklRhoBounds); } else { lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, - ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, shrinklRhoBounds); + ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, true, shrinklRhoBounds); } } else { lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); } if (piecewiseT) { lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, - ppd_factor_T, shrinklTBounds); + ppd_factor_T, true, shrinklTBounds); } else { lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); } @@ -369,7 +369,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Real sieAnchor = sie[0]; Real sieSplitPoint = sie[1]; leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, - ppdSie, ppd_factor_sie, shrinkleBounds); + ppdSie, ppd_factor_sie, true, shrinkleBounds); } else { leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); } diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index e9f4245da8..168ab90304 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -71,14 +71,14 @@ class Bounds { : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, - Real ppd_fine, Real ppd_factor, Real shrinkRange = 0, bool convertToLog = true) + Real ppd_fine, Real ppd_factor, bool convertToLog, Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; if (convertToLog) { convertBoundsToLog_(global_min, global_max, shrinkRange); - toLog_(anchor_point, offset); - toLog_(splitPoint, offset); + anchor_point = toLog_(anchor_point, offset); + splitPoint = toLog_(splitPoint, offset); } checkInterval_(splitPoint, global_min, global_max, "Split point"); @@ -97,36 +97,36 @@ class Bounds { Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real fine_min, Real fine_max, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - Real shrinkRange = 0, bool convertToLog = true) + bool convertToLog, Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { if (convertToLog) { convertBoundsToLog_(global_min, global_max, shrinkRange); - toLog_(anchor_point, offset); - toLog_(fine_min, offset); - toLog_(fine_max, offset); + anchor_point = toLog_(anchor_point, offset); + fine_min = toLog_(fine_min, offset); + fine_max = toLog_(fine_max, offset); } checkInterval_(anchor_point, global_min, global_max, "Anchor point"); - grid = gridFromIntervals_(ThreeGrids, global_min, global_max, anchor_point, fine_min, - fine_max, ppd_fine, pdd_factor_lo, ppd_factor_hi); + grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, fine_min, + fine_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); } Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real log_fine_diameter, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - Real shrinkRange = 0, bool convertToLog = true) + bool convertToLog, Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { if (convertToLog) { convertBoundsToLog_(global_min, global_max, shrinkRange); - toLog_(anchor_point, offset); + anchor_point = toLog_(anchor_point, offset); } checkInterval_(anchor_point, global_min, global_max, "Anchor point"); Real mid_min = anchor_point - 0.5 * log_fine_diameter; Real mid_max = anchor_point + 0.5 * log_fine_diameter; - grid = gridFromIntervals_(ThreeGrids, global_min, global_max, anchor_point, mid_min, + grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, mid_min, mid_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); } @@ -193,8 +193,8 @@ class Bounds { // min_offset to handle the case where min=0 if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - toLog_(min, offset); - toLog_(max, offset); + min = toLog_(min, offset); + max = toLog_(max, offset); Real delta = max - min; min += 0.5 * shrinkRange * delta; @@ -232,7 +232,7 @@ class Bounds { } } - static Grid_t gridFromIntervals_(ThreeGRids, Real global_min, Real global_max, + static Grid_t gridFromIntervals_(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real mid_min, Real mid_max, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi) { diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 68ee3b037d..5d0e2a594f 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -83,7 +83,7 @@ SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { WHEN("We compute a piecewise bounds object with three grids") { Bounds bnds(Bounds::ThreeGrids(), rho_min, rho_max, rho_normal, 0.5, - N_per_decade_fine, N_factor, N_factor); + N_per_decade_fine, N_factor, N_factor, true); THEN("The bounds are right") { Real lrmin = singularity::FastMath::log10(rho_min); Real lrmax = singularity::FastMath::log10(rho_max); @@ -106,7 +106,7 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { } WHEN("We compute a piecewise bounds object with two grids") { Bounds bnds(Bounds::TwoGrids(), T_min, T_max, T_normal, T_split, N_per_decade_fine, - N_factor); + N_factor, true); THEN("The bounds are right") { Real ltmin = singularity::FastMath::log10(T_min); Real ltmax = singularity::FastMath::log10(T_max); From eb6bf666d410c6dada66cc95667cb65bcc558c44 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 14:02:48 -0600 Subject: [PATCH 26/45] formating --- singularity-eos/base/table_bounds.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 168ab90304..3c6d0e98ee 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -108,8 +108,8 @@ class Bounds { } checkInterval_(anchor_point, global_min, global_max, "Anchor point"); - grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, fine_min, - fine_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); + grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, + fine_min, fine_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); } Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, From ac0db6e80c7cf2634bdb009f4b34c96ee3b52294 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 14:20:47 -0600 Subject: [PATCH 27/45] add diagram for mesh. fix some build issues with the docs. --- doc/sphinx/src/building.rst | 8 ++++---- doc/sphinx/src/contributing.rst | 14 +++++++------- doc/sphinx/src/integration.rst | 4 ++-- doc/sphinx/src/models.rst | 31 ++++++++++++++++++++++++++----- doc/sphinx/src/phase_diagram.png | Bin 0 -> 35401 bytes doc/sphinx/src/phase_diagram.tex | 31 +++++++++++++++++++++++++++++++ doc/sphinx/src/sphinx-doc.rst | 2 +- doc/sphinx/src/using-eos.rst | 8 ++++++-- 8 files changed, 77 insertions(+), 21 deletions(-) create mode 100644 doc/sphinx/src/phase_diagram.png create mode 100644 doc/sphinx/src/phase_diagram.tex diff --git a/doc/sphinx/src/building.rst b/doc/sphinx/src/building.rst index a25a60087a..a716f931ef 100644 --- a/doc/sphinx/src/building.rst +++ b/doc/sphinx/src/building.rst @@ -643,11 +643,11 @@ packages. This can save lots of time! Note, however, that external packages are loosely constrained and may not be correctly configured for the requested package. -*NB*: By default, Spack will try to download the package source from the -repository associated with the package. This behavior can be overrided -with Spack *mirrors* , but that is beyond the scope of this doc. +.. note:: -.. code:: bash + By default, Spack will try to download the package source from the + repository associated with the package. This behavior can be overrided + with Spack *mirrors* , but that is beyond the scope of this doc. Now, we can use Spack similarly to ``module load``, diff --git a/doc/sphinx/src/contributing.rst b/doc/sphinx/src/contributing.rst index 1344cffe37..3dbd81cac8 100644 --- a/doc/sphinx/src/contributing.rst +++ b/doc/sphinx/src/contributing.rst @@ -256,19 +256,19 @@ instantiated when ``singularity-eos`` is compiled. Therefore to exercise all code paths, it is best to create an ``EOS`` type instantiated as -.. code-block:: c++ +.. code-block:: cpp - #include - using EOS = singularity::Variant;``. - EOS my_eos = MyNewEOS(parameter1, parameter2, ...) + #include + using EOS = singularity::Variant; + EOS my_eos = MyNewEOS(parameter1, parameter2, ...) in order to properly test the functionality of a new EOS. Simply using the new class as the type such as -.. code-block:: c++ +.. code-block:: cpp - #include - auto my_eos = my_new_eos(parameter1, parameter2, ...) + #include + auto my_eos = my_new_eos(parameter1, parameter2, ...) won't ensure that the new EOS is working correctly in singularity with the static polymorphism of the ``EOS`` type. diff --git a/doc/sphinx/src/integration.rst b/doc/sphinx/src/integration.rst index 495b85a07d..cd7c3a7b09 100644 --- a/doc/sphinx/src/integration.rst +++ b/doc/sphinx/src/integration.rst @@ -8,7 +8,7 @@ or setting ``CMAKE_PREFIX_PATH`` manually, a CMake project can integrate Singularity-EOS via ``find_package(singularity-eos)`` and using the provided targets to link to it. -.. code:: cmake +.. code:: find_package(singularity-eos) ... @@ -40,7 +40,7 @@ dependencies, but only the compiled Fortran bindings provided by Singularity-EOS. To avoid unnecessary dependency checks by CMake, a Fortran project would integrate Singularity-EOS as follows: -.. code:: cmake +.. code:: find_package(singularity-eos COMPONENTS Library) ... diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 8b7243e87e..63ac375c1a 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1364,10 +1364,12 @@ the energy offset of the products EOS is given by Practically, this means :math:`e_0` should be positive for any energetic material. To provide the energy offset to the Davis Products EOS, `the energy shift -modifier`_ should be used. Note that the convention there +modifier`_ should be used. Note that the convention there is that the shift is positive, so :math:`-e_0` should be provided to the shift modifier. +.. _the energy shift modifier: modifiers shifted EOS + The constructor for the Davis Products EOS is .. code-block:: cpp @@ -1511,14 +1513,23 @@ default grids are now piecewise. Piecewise grids can be disabled with piecewiseSie = dalse These options may be true or false. The default is true. When -piecewise grids are active, the density grid gets split into three -pieces, a piece of diameter ``rhoFineDiameterDecades`` (in log space) -around the material's normal density and a piece above and below. The +piecewise grids are active, the density-temperature (or +density-energy) grid is built as a Cartesian product grid of grids of +non-uniform resolutions. The density grid gets split into three +pieces, a region ``[rhoMin, rhoFineMin]``, a region ``[rhoFineMin, +rhoFineMax]``, and a region ``[rhoFineMin, rhoMax]``. The ``numrho/decade`` parameter sets the number of points per decade in -this central refined region. The regions at lower and higher density +the central refined region. The regions at lower and higher density have ``rhoCoarseFactorLo`` and ``rhoCoarseFactorHi`` fewer points per decade respectively compared to the finer region. +Typically the fine region should be roughly centered around the normal +density for a material, which is usually a challenging region to +capture. If you neglect to set ``rhoFineMin`` and ``rhoFineMax``, +``sesame2spiner`` will set the central refined region to be a region +of diameter ``rhoFineDiameterDecades`` (in log space) around the +material's normal density. + The temperature grid has two regions, a more finely spaced region at low temperatures and a less finely spaced region at high temperatures. The regions are spearated by a temperature @@ -1528,6 +1539,16 @@ corresponding to the temperature split point. The coarser high-temperature temperature and energy grids are coarsened by a factor of ``TCoarseFactor`` and ``sieCoarseFactor`` respectively. +A diagram of a density-temperature grid is shown below. The region +with temperatures below ``TSplitPoint`` is refined in temperature. The +region between ``rhoFineMin`` and ``rhoFineMax`` is refined in +density. + +.. image:: phase_diagram.png + :width: 400 + :alt: An example piecewise density-temperature grid. + + Thus the input block for piecewise grid might look like this: .. code-block:: diff --git a/doc/sphinx/src/phase_diagram.png b/doc/sphinx/src/phase_diagram.png new file mode 100644 index 0000000000000000000000000000000000000000..469e91c5aa63e5da7c8a29e61c22e68464e0f357 GIT binary patch literal 35401 zcmeFZby$|$zAlV~$wWdCP*9W>P*U0el$MqjkQPuv8dMAf5$W!d7LhI$l@tU7q?K-z zls@-6_g-h6bW} ze%VQZUmolyB8vaR-7xbwM%YoW%0rZ!o11}wft!0izgcXmuWZFvNBRBx_jm8!ef;KC zhFYfjtLqQENQv%~(}#xmG~KmO z8*uu#8xc|PpLcI9+S$l{FD+TMBtDi(a(?e?zP>yoqc!<8q59X)>YlgeDtU&IuCA^W z6cl82wY3*q7GFh0MR70mme{4qMHrf!rz@wowzeuOD-RA1mRD4Gu8xFNKPC0O79-?v zVV;P{iEd}mhnAK+!%seBWMo`U=bfFLZfI&+b-y+L^XsQzQ~bB;XWow=f5P)-rl%M6 zvo-Tg7Dm6+Jg4TlH~l?ZyNKVs_0eg<(#~v44HOm@aye1=Rl2X&hp zVh-8c)YL?5&YGDeFK5yieu{yxm9a(dl9ki-esOmfSZaiE>VAGs9aJOYxw#ZJA1@om z>D;KEt;I*!ZpIhXk}tVPN&QC}!p6>DxpCvh#fuj&?NqU~eN$0UK^{d!bW!?}xZLWm zg$0$&_P$bw1|zz@r3F45k*_w3=b>JC=kLv;V0YB)JLIh}~p^$M4z z#8S<Xt6n(0gWSVf|I@TTDF+moCn;@cZ`^TTyXQhwamWPMswWa$NMKlb*Pdv0%? zKYyM$w<%71du?f|`Uy!}3ZrDrixYcE>5d;g%I~@5X6%!mKJiK^DFbUhuk~pHmzh+2K@ivo6sgAt9IMcM$@BW|fD#lT&e0lF`s-zlCtCUq9Y6yfLbI(NknI z#iUniFP*LKd+o)_cVQtRpKEZsn>TNs=MsZ4Qc+RC9{=-Ix=I}~z292(@XvKwbmr)m z*p57Z{@kjkaDqw9W8;SNOCtgQ(r2emonmA>cUSZ;#GIgn@$vCYja;8=DKa7RF(p!d zGz{}N;PKRNS(+-&hjHre9(@vhqK=qLsvRnod{#t6KOP^HQ&Li@_jD@%EhacPn1q(U zxTM67R=|4vYXVHGU6*1*z`k=xM_2b?-VPOed)^vGBBJIiH(U95*sjeZ@Tpx6EGa3m z9-QvVkBN!d-dJcxWSwj-hfS%ZDOCLaOber^w|YS1$0Q`Qwm8w2Z}P?B+pCtAq^lI< zqct{L>dhnBPkHdCFvz?4{^7luZW20?3yh2- z5eGRL7%o$l9LU_!R_?s;f=-l;A}|!17GblwvP_!W%KhkA}~@H93e`(Yn7}GeucYP!Q&mqE%qw%ykFRM)Mh0bHlx! zLMvgr(MCd6dGa<&#Q9Br{sNrIq0eYbrIc6>aqn^K{fVy$SJ*XjzKPDXXNLLvD?UDO zJT;ccF;=f2Ufl!rOu# z8%};wNr!9$Py0#b-<#12ktuLnwH?tDJi*MI#wZ&~nY{lS4?>d;oOPzVz|Ku=;zWHz z!xVLB*1-LAyRoL`^2eRI`Zg_zpOBP8(nMVrzg8skuKIX;zq*k(_}0A5Tb_uBGGU|V z93wL`e19e>j7xuZXei%pb>3rhZ6us_v^ma`5G%#EiV??%xnK@Xb%BhCkjQnNGcz;t za&oF{h%q4r1qE=?pdc04`hnK$i*yekK9u;_(ZOG%0qf&n`1$i^^Xpri!s}XZ&3-p( ze@8^{_TD%08I|+r^1;(e)<%1;V#43vvLyr zR^6&;3NZ(cupxg{*VoHL)YR8cHj8f`7qq=KQ0b2C5jXs}Pq)XSKugqdwijX4bK#2s zk7-k5T;J-o;mkDaZuf*Le4`$4WXCPxmlqYBaoi(0nX>wm0BQ)X21^ zprb;4SlD?ZiB$zgCYd|0qeWPM{~nyfbei9cY}W%F4$0dFSvK zVI11(g}1wMbj!L?o`#3tx^*i#B_+Lm^2d)K+C|o@fBvvMGwm;T)_40YX-)sGyIY8n zk&%~|*VXc_qvMLA*lIH3XQ^=aqTjLehRB!{%nCaBCi4Bp!?Uwn$iZ5gnwnZ#S~qSS z>%GIwAQP&b*$y{E^y8}9u-zdmjx5g0`_^S?vh;b1cCmHU85A^(jCae9Z-vl_CSTsu zl{eHCcIS8ETPyY}PincUMEFD*F1xIRe)3Jx;36YqW-f72o9zqIDmY-;VAFf}>Q!nn zokOKGnnenu&0+MOl$&MJ%X485b{Y3x=qk>u^>JChH_@kEA374hOD~uv6;?#)#+qBP z+-;d#f$Vei@L>e&PDEM!6UwH>Rg$UdH9%TQ$qyy_n&bnxiebFCxiczE)ABX3N6@9rMGr$)x|N@r4Y zdq;}m_H;p72z!HRob>hUqiso78AZy!wgz4^kq=2v|GDTm5HfD&yy<|&Me=`%zHsxt zldEi!l%(c2*hs#*&n}jx1NmK-5GP@@+-6@-5-QQwR-C5>(hrJlk=8)(@W{x&P6#c+ zd$Fg!M>$-c%;Xgnb$&{nZf$D|DY$c(Xj>+Xlbk$owfJ^d)*r*$K}^pgBIIJ$zDY1M zv%YhZ|D4!T`(W2`GNK3RcQ*bE`QtIf6^DjX9i}=E%gPXbc+1UAF)Fqusw3gV$FMzAdeAWjuFX;$xG7z30Xm z=7!DnM0L|IN9p4V-t8sj7T_vlRvd1Oc$lL_)kJ*gC*5ANAmG_sk)UfPBE)x%H-%Y9BrvSW6F)8bD3vy?ghyTs`HyM7fASzAHDv z<-$)E+4vR|oV{^!sK$4eZKK4lH9JZ*izUvD!E=*Q-;-mkxgpA?zW6SV@3Y? z7#$Fxgk%AS9^v=0ZO}6_Gg}S|7H(iu%LpE#)h@iv%HX;@-Ex#!vRMC@$SZyQ_0P23 zLCgh8O5JQ~Z|UeF+FVyU3N6$6N_Z9~W_q;9W_}Nq=0rwxHg6fa>w@xA?PoGhaG#l)7;oubo+&Ul*UDDU zY%%L8Y_HtjQjs+q(CaD6Rd=MjE@QH0KNG?pp;Bnw-ra3sWW*D<_Qxzc>ungPv6%ax z>r+4CMV+}9CYoZj$gCtjs$^c1Y3)c)+Fb6*VvP~D6yM}sxNG(bA8#JWFQn!A=XIen z<4lF?@!E^$Rd@5rA)E8?@oiRZuUD1Vd};qe=xY#cO>Zb=s%JNT`I3~BWOw^>a&PP0 z@UVp0=Y=Fy*8TflD#RdM?mrkL&|m8^G254w<}l^+`0@B>KCLRf96k4pru2JOuq-CI zb2^pT+23+F3M^$~gaVC(KTQj;eDA7@5msC>`}#_$BSTd5Q}va87xt-kgX-ar`%4@$ zFUDz@H*O6}=g`Ind)X3SPfxZp{r&4)uLu{H^GrdS3=_*Wt-{9W*5bS~Og>efqirh# zVNXa>)G`~>M(U-TV*^)W3L8S}DqOiulVqu0(tb~Ma1cUkd}u0Zf6G;5+|bru0g!k z=ER8;M~)yMW_CwM{faj4I6F2WAL%{&LWzPR@Act6{j*h@Ob(OBr2LL@40q}?GD+2k zG(;j+nk-JVoZjv~>mVneVsK^C#-yX5e5giJA~5Jv8;_mO8?_xf0%iTrC@QTfH4JH& z+ozw^Wn&l|I{n}Q`QejYWqJ*I+-qqEQ=ZBhXSTMDESgVzKW}L%`SHVLH}{3^T0{gL z!4hX4OBvfY5|Ywx<-wdW&3e18-wvm&(^_Ad^W0opLbp_C**TkC)GuT|ejV6@g++!` z(*FJX1<`y*3);~LSa1CKJ%keNOZfTI{*varG%ZsFwePSilDClKiQ86YdeU0Vm=vOG z$?lKITZd0HFo*rhA=9g{6`u?*e$Mz-g-7}Mc%JoM<1w^y}_jro9Hd3h|G=$h3|AJP0ot-Cm-u(O);BVqbV`3)Wy3Kt&Ah3{SiE-Pe^D(-JgiLuRW&|b-|ucFv)!_m{X9G(?MCc zGTVpBqJgkB&J+G4DIV@6b~SerZKEVN*GeigEjurN*Gi#D_i#^F8KbWf;8#b+*U%5? z6=IQA45Jf0BdRMR%Cj(#aHX{BZLIOxune_N{d!qh$+~45rZ1jT+SoWRshRiQ$kA;W z<7j@}NwiTU=W}mqbaZsr-oK*7J-137raVw&UU~W#`o4)aMe!4;OzF(cQ5K;S%>w?0 zH$jE#cn*b*0JCB5D=)1Q!G_apOUESH~_~YFJB+$))?X9%I zGkd(R&Gp-kDU9|l2~lx+2L-8Snf0nw&HIE$)CHVYls|MR*s^opb1qRbvE|yVmuJIV zRm01&9a>&cWtFFbf za(n-}wu9#G&Xeah+4VhEQ>v;w@7%fL_V*VM8~ko&_NuCCyU?oVkO@GxE!58*8^3D? zJh$EJgIRHa+ly_A=oaH_mX2JXPfn>Dsu5IFl;w#|w7kaa{N}BBP*>XhK@_~M%jr9$ zf>XsymJ*YzZUm&2-j|ULTC{JeGRvu&OqDloI<8s5ooOB*b$Y-pNG(4_u5<3Yv_z~( zfb+s|m_C9{bw>Kzw^;@39FOcp^HP$p*}qLrm63XK2#%u9{$UC+!wb#%2Y$Ak7PLKMM;kmcR@_NM$(&(2ewbXUm%i(B z;I0G1<4T%k>L^q#EhoM!tu5A-PtEiSY8J}juR>%>9m6#9KMv>@S_S2srl_*sbqJZO z3VUPBc1l*}<~f z>d&t$g7@UC2ST5ap!7sh$|U#lAk`TrDIJ|#f9~A15B}sk^8LKDY!Ig|a%}RIy<`Cn zQ@3jVqUODZmZ2eRq+%uc(yo$w@?zUMZtkC&6{lxV@Sf+P$W}KrJj1VEY9HV)9?EU} zrek)55qr5RU3M>It}@7WWO&G*5Wd|tCfh0W^_{ofv4#va%;DJU2P%9P}h}@-CC*g4W+5eEmv)az2zd$R90$I^Rh~d z{ke$ZU~5w{J2d{}>iFrr{2;Z=Tg`MOZkieDAz}KV=_;Q-JvO2~dG5ixM}r)EY8#s= zE)J}Bs@^=7$_(i#g5$Aia1JvvGfH2)|Gm3jmg||f8E;_lneEDltG|u*{}n3bXL~G- z>#Q8&z~rd#=3ut2_Vo0H4gli>kX*z@j%8>3*Nuf8ZbX0GjF*Uik{&)oY7s8;XN%(A z@eB7KjFIvE7g#z=gwpPd28QJ)09x&+-77%-pqNy|zxd2vV7oU#2;grC0Pf#N{JVt* zM+bfo5&sjD|Ih#Me?av2-|hZifBt{_ zZM?WM^^xF{9sl|t|G9Vg-eV*~cc?4=+ZU7ydU^353yC*-b*%A?(~N)B+=Ouz5z%$G z>(&3VrT^T+PvO*(M3gikGQ|}w#d|h%`o2%i&C%1)XtDTO00{@idEv>9)3Kh#?iIzG zbS$Lbxpp%h<2#p;)zP`6tC$nA5dHxL0 zvlR>Ub7wg?Bqg%dnHd>H-Tqv!Vi6G$>313;d13I*2^RB1aCV-#$`OdSLI|X=P=u*RQ{?;3x4Dm3!pXs63pN24U)Py?X;>#x7S zV`S)--Fw0FEnVfR^VdYlgI-~`Z{Lmvs2LNZtE)TUw(uo9{MbGAo2uW=?Qr@T45(f) zPK-erO(D8Ns&*3*w?B8bHdaeg7|G7_G;Bn5eEY^*)gm0{O+;CV#OLJd z3h<>2(4XY7@F$0;=jXx}KGWLZM2v<4kDp`+3=apEP>pv46w;Zgp%`(DiYiSmOTWtF zf%T=tf~OXx#QK>Gz@ym;w@xlcsJcUEM z(Ei@NO8W`&GUCVV{QS!)VR~%Al6<$nzvP!=T-gX3N+134-U~H_B}gaK5pAul-}!bk zz59yfY>HyVVe;FyQ-$lwtf<}S^Q~YkI7Y-7rx`L$d8`^QWZ99W$>Ix)v!S*&{+jIg zDszp>7oe4#I>qN~hp*i`@L+!}aQjVE3H;W*#oHe3N7&WR{*Z{jzEcID%K=OC($Swk zHg=u432>5ZV`)F}-F*Z$N;|nMw&;tzQM~@cXstYW0oJ&nf zahUCO9&LapPt(9H)=1qNV?KK4jem=e!ek-jYCV3#_cEhy+$2%SHt9%rH zUzKN-PKhlnpOuvr9M{^T{S;Z**>MS~il9428>x-+tE;Q2&ffT%Z>rdrX;I*U%Dz4m zU=>(FDlR=$J)t{4-hs)xa^(sjN^E&ou71_e>Zjs%?F!h0gLy%}!}L8O<__Lm6raED zCQ|L;yg05zmy0(B9yL<&2#W&AVcI`>)t@dUCB?j;dvS4bI@f|X5?mQ&&>+fBXG)FB zmmhiq>AUJVXoGfup=x`}b;pYSDeIYZl{A!Xomraspld+2l-)xH|EZEH*J8#iI9nfZ z`YMqAm3I$zl>+2{22`!Es7Sx?wkBPscGA_-?+tZ`qlH*8_=36Fi0w?CtFZVjEu-4GsnfJQm`84e$A0 z%h1q}n1n{Pz@i%E;0jnw^-Mn>pS88AT*+jdm9yuL?GzXhYtOG9?1Jf7yyw~y z*M|jvsU$2-sq$dG=BwdYPg{|8mc18cWDbbG__Tvhj#2H-onK&20GR40Tz>5B?G3=T zF-oBM^Jgi6YI@I&Ja`a%Cn7TP$&)8Q6z5|c`_JE~MMDOZIYT|$1)K5u^=nDrqkLvx zFRn~}d`1BdhU&x#Q8*32G*0cpX@pZia8Uv_>glR2EiKsQ9M$wduqAL$IeB@-SdkKd zm*6nQ3-2^V3!+10QAvHCn`;B_9&3yUawJ(mzu~#HI)=Rza+qW}ao+bDqS)L(6+J%4 zb879Nx12qf#lg9jkz;30KoIg-p`VjeIvX=!P-3{^qUeXs2cTmGR1KvOBcx!8s< zN7b|Gijb3V8k(o%E~s#&Np+9hid%e1no?~r>r@0m6epVgJ49L0ho4Wem+C1AY)4# z1YS08P#iiVpCjxNKu1jw#lX>LUQc+al`RH7-aA^%T`R&998%Tha+JTn|D_YktMenC zU_cdOgp9c!Aw(ne4^(+>V+o~VumWTd0)vdV#-RZFT2mBOjs*FH^R?+O@5<8T@sl$8 zl>zqq%@?ey9psS284eDs)@12?(`M`Tj1X=DpCh*MtNQELuK?tmV?|xz(b1j<5O4zn zTY;J41YW&*1(vLC>kgnd)CU zJ2jkLuy_Rh6QVAvMXx|)zn~Wvn`;h1fTE;Sg`=CAn)XQ)XjPat(~+%C=Qo=)$7W|f z!pEDTVzA9Jt0VO9t+`C_I?J!|4|3dmv`K=H>9e!5LqgQh)LVOI6xs~**VJ4UU<4_T zNMT}Pg4kWcq5g~y=xqI#G8*F#A3oT!3z1Rt-@ak(tL8#vYi}QY&fwwJsXVh52GCTx zoANM6JDmbU9A1}YIPsxlCmEvnty;Gqj?;Pm`T3M#p0H3mK`@nqt|UP(W@UOi^*t-V ziS?Si&d-;HASA&ed@E`}QAVa7&M9#Vf`-}IS%i7ekf<>JWXKScsj`mAyrR9oQs^CTyXsuDi!ZqaE(i**JNanEwu|QLab## z(V;PRSZb&7OO=b@_S#Fze?u3ZHr1Y)wsE+>%F`3JjE_8fw(;Y~Gh$+To=v3aG4{)) zr>BG8MHBFLD-(PX@A>mBj2R*q1O+GJUcc3{x_vv#`x}TM28Itm#90wpL(Xefn?N5^@bK?_@0%dV82C2R=WNN!P$q1nEd`$`2_{-+xvSf zT#LL_b96Mfm3_A~PG7k2{oA*1A3iX!u#9)rYfdEa7QM5d>fjh@_zvqn?n%YN%^ghR zi)3DT;lV_e$Ts|d$E5U)o%l}#rr1}Ym=GvUgR!Yr_LSmv}eCQ zpQT&IkDBMhqy1|X2=qc0-&j%t=8(<+ec@}#^mSn&-MO3ZHdYtfEeg0(K;M95d6UzX zC6fX!TM5PI`^jKSUL=3Y0AxFqtvV^ntQu;m3NabtC#%7Lr1tJEgJt3;Cy2E)?5Dpg zWwrywVHqF|uTa!BmI3U)WK*rK6(J!dHCpx#x$c6naFKbNlt2oshaqsY`siLTun0Wl z0q2cA4_~PMJye6tFFB0gOQ)z{=gNP1iG$?3|3ej~+c5&byZq7LO8e>aVKlo{soSyEfUGkQ$B#ASe!ZWR#kGiG-;r zfaGT*ioMiNWoDndJl=A;Mi9hjaApb;NwnwlyIV??e_p(JQNKS8F$n?i2pI$uHa2r8 zCsUG>$4v-nA!Ki~=gl*34cMHa1j&t-SY$3OD`O`mr=UQP>_z~cwt?oQA7R_3w>aVX zzJw`e4p?qY4YQDt5ONC0taVy2)RyN`K+oc!)j3z$hf+xCM80=*X*#P-WbrdGjUpbc z&AllofKklyc;4AuR-0eBlPqZcdIxw%?EW)vpnr!BA?T*gGY=tQ`25^LasE^}sSN_YhV(9=ISZLJTk`fZb-(S8$ z`lCdtpH@S0Yik}sR+|VGKu3Ky5hw!<$v_x{~ z5bq?N`S;fbbWcOcZwPL0gpI`^d=~)dpZrFkv1eLSe0_Zz5U3s4M^Nh~JV_t^y|O|@ zw)x~8F`)`OZd0Tld24NLjZKk{<_q9HJ+Nnw4^p9ktC+86N7`)} z@Ryrh9{TzNuV0_*JC?T=Ec_%eC@8-wdE>%_M>~^T_#G%*P!F&>+RG|9v=EC&2tlL( zXR;}WX2qsFK|K>$8k(dEV2XtD|Aer5*K4#Qle27WTF_6;;$ZyzzMw!e=uMcgvfAF{ zMovw86#lS^@#M*itC#9mqGx~nh<9yZ&n*vhU^vIaj*`ljoUoUu%1~)8{rY(YuF~gz z%Q*~-Ph-i9dln>Ovh3Ax^C@w zu%DPVKB+PsZpDLC?8NP%2*O=n$2I;q+pm@VzPseC9m_zOG^_V zTI-&}$Bq?4yaS)nkw9~g4ijJP7E3s3Sbk=7+Rk0O^5b4Zf48+Bhk|Srtqi)$8;*rC z^78E!E-3kMe4w$9K8;x)uA?ySnVz0Lu*I38mcAw5=lqdovuGTNbz^IN4$GvapjcS$ zvGG54J`t-2ifmCJc+CkSExG6m0SlzhAIaJ> zHYZC5P0!3=L(R{*k2S^Mn<#*+j_Ou<=%$BtNGmFCgUKEM-jL#;hP#AsE|quhQ;qeJSP*vHbHtNL5J$P zSXNPy<8SKXZ{6KmtFA&dRV@onVdDimn!RS_tX}3oQqX~At_*lqq51%L3?Sw%vc283 z{kndq{`bn2epe(4AAf&wyfA#d%wZ~ZBl6xvt7PXR^-jqzh@XGz?Y8P>(fOUfkOQFl z!uI`lkJ41LwH)A+NR*ug3)k5WL3SD<-U;^OJB$IAG^sBafi{&OH=KPYDY2T_*VCh2 z^a%B>T==)wC&EX*TH1n>TRwHh2Kg zANyQL5;$~A-=3_#FYY`qQDyGSdEnUYog?PUk%7~QBxFdx4o^7GvI4^17Z9_5c-vn^9=asMN&6Yj=?haMKrsqi7RKtp9zDsM1=$dftVHti zymN@9f)DzaBW~s1*c_Ip@|B|GkRnTkIXX~;Skgm>+;sO2Z6C)?f_#{pk&~12%oz!R zYJKqS3;NQCx11(e#2MbO#80!7s`QO$@RjN&H3~|RTbhsDt2ri3LD}R`VU^P|Mss5bOmyOo-9PK&(+HGl4 z&Bn&Y_|e&v^HbFa8)}9jXbY)?}eg{Vs&+O4z0I$0BO8=qn!Fmm2K0Y))q-Nq1VP3 z3YfvcRv)6X7NW~D4A||)ArADHGJAJJ=GY?tjf5z1V51}6IcWX|Jzr*h9``z; zTC>p;z%#+@8VZgLQEviZ^T+|;*l$*oi{YohSJ$S;UwiR0QD{4uOJ8q-I31A_RctPx z7i;&Il{d!dzdqUTlmKXeMf5J}$@Bh4$pHI?;f0`HmP$0*<Xu&}T{ z4WEp-BrB=yl*t4(sSEGmn6u$e6dv) zrlTM$>k746IIkRWs;v@hv33^pq)C@5-I zwW~gKtgx$;ra(k!UVZeup5g#PQ2nmR&7GeHxsC<*^2@PUujjrF4Mr)@+tP+6(koSb ztT_&SlLI>KWkmXfiDSpYE2w>J>6iA>yZxA)EH$Vm z_R;`ciTEahieNB<%}FX=o;~6&Se~5mG?O zh9-ml*s)jat4L*pfT8zrZup^V>#NEbMzReyO&>+}l#043emIcFVx z7nimF0BMcP&7FaFp?4g_qIc}r0S)38*VjLL?f(xwB+4sTiM_#L3JRd{NALnrclM3R zprVzNd;a2udX54rf+uhpu2ZK^pFV&7F3QXjcb}bFU`$}BD0uM&@G)esRyIExJbax# zeZ~35$cS;4b=JZ4W5ls>6e8k{yO}5#=%Z?DYtcM^`t*rl*x1}X6B!poMf+-M#@>7F z#h?ZnF6W7qumf)r!P0noZUZ7iSilg-w|Gg5zjcx~Ait6?s{`nJ_KROrrmbd<0?SFiqq}O{#->$8jzjTR+eCWXu{TBmEPGaZ& z3j6Q5MEvN3{VZ(HiOl{TLUC6Z&A&q^hD!av5{lq>?-oDTdl!7J=G;F?%zv_tfA7qH z{_&!c!S0JmROc8Q&ma2VzTG z#b)4S^e;0p8Cu0GEq>UuIz;;Hd2pCsU`rC@P^(KC)FQ_zG8ML8laXJUk9Q)CQ{PS?DEDh;AC`)QVE&&a7+QEs(oHdMGqV0?qUrDbrI6B1fK zT5i)9-*UOi9I1HEJ$~A-wfV-SpnT0j{pQ+?!0Luw*Xv7>d`3?RN5^8qOQV!I@|h%E zU)-+#jWM?!DolpCiZH>E4vvX$tP8&%?6O^a<>dTKhLNE)OVebLwOGpURtp1YhxBJP z{*d0-E={UDF#k^xu`w6*K`X*@Uw*%F*RDZEmch`1jlnQ+MZ3E=uY4`z4C6EJ+kI?( zrN*g@k5f|IHkZWpa!XcK6sb7KxwxnW=hL5=Hp^(ylz9^tujKnYIFw@{E?e^8vlmCn zN*om*K0Ksv)h1=pxup?hV9jaOEfUhoOCxya{^k{^7UGYc6QD93+c2H_dMaN%{nD|YLPFwpQfhFt=WT}N#K&)*CA4Oj{#h{!X=(ZRj7aJ{D4dJR?+97(*^_7OTrHCbE zVMm|L%xIwyJF5=MTeK=XlN}+o*CMzjLpd5EFHPi`FHdI`{;Yl?oo_PVBwWJT)+l$DgH?;*$jt$*DaMOt zUBya}?*}obC^KpnpNlJ$=_|$1TyK@)PqOfr^dXb!p2V$fB&M}uip!6OFDM!s-qlfE zc0Q1I*k>&9JKuC9)xA#2=J~UfbaYQtgd74LrT`cKPRJEd&t_w=UFR2n@@bz{xKDa{ zS-$`e``7V3dm?H&3NLVGoLrzZJVn_onP>T(kE=F^HFiXkWqF!RK%n|rZaD*+n%K;n zD_0~XX!zYiL!_j9idMv@=qNz)6J5V$-TSi4k%XDek)hD+_cw_Y?Hw^#fU{k?dUiZJ z*jp{=5a)lLGFH5~(}!oMefo4hL8kF@kdYC{V0P20nxB-^^K-Xz-rmc8ej>`bMSd(y z{|%k16Y%^u&UPKaZ+WZr%G+AcoIgL6-rDs(Mwo%2$Xe8>opSy&qaYE{R)p7HsD-7a z1E86*H(wp46?9e4%DS~^HqkF?$$#)5yf8vgWwPUlhK5r5tpTz92LkCmd3ECb{M_cB zy|w0`;;PRx8XO7H+%rNDGo9AcgCAL2!>@J9m~irsA7|-vt;{uKs4*}USk!H-&Tp-A zn+N_`s}E$HW;ngUm$2EKGcp*rvEz#eYr1Dl{?PM{JwG$wpVxT^-F&7x1H-9PgEJuI zS~Ei4T5o<6SzoS;zeG-XCWP%-Rn@mwzo&PNKY3i2bhR~k|B=U1mq=g^QVGe7Wa*kb z`Ezq3H#9-WtfghB^Od(gy(~eCAnhZEmVEnfam^L5BoFr(8l1AJWYzbScc!%D-N$0G z?!_-h^y@=y?O~O+<8)r{b`xxjku+@hzs{atD{}t(?0J)S|Nh_y0so=c@gX*yQ^E8vGM@+golgC`GM|4hEB~7^pZ`BA<3#?E!T&$b%E>88|Il}_(6ABb zpVeC*nsasO;({sA|JSDQ%kAcI&VQ$jAtz_qb^1Tb9ZHhVf$Xs4JQyvD$18*XO)2+J zJNtF{<`3n6=i$e8AmaV^#b@HNr?T1Vn5^8hcW>1`$Cu|rd+nSL?P7cg5m>MU)K%z0 zh0wL4RtB7o<^v!&FkNT>j7&|ht`dq&r-YQ$REPxZAdvw+B;>JC>>)D$>(?>hU8u#s zH2$%mdV_wutEZ3GNg0DOQ* z4+57-FzfY24t82v-FW0%0CsQR@(K%=;envA{OCm2x~xn5eSLf07|VVf05l)`LZJG? zhfsi+*z~t=-=arsjuzx`T^Rw|1npQgU=I^87gkm$=;->PGH~p7Ca{9eS%gtXPyz?{ z?YpL^m<#>_c0l~saWf3DqjA6U^8=`FfId)8^i@|Uq0<9$-36l9dckFJJisXKuLFKk z3mEdB{hl`jH4{KbOoahS@**a20g3|J`jClwW;-qa?dwM~r^CQvvT>aNod*CGViZiK zf^VDd%pJgMOjJ~agGmRGKreH)|M0CaBzjW3e+pnqXtPfZ^*&5zz7gsLE7%=7sj*Qjc5A|{+hS>3}&CSgaRvJGOfGp4g0?Wcdt*6C7F6Bx9!_+cE z-&g^?)XMH?Xs~>De+Qt5Hj9G#v{)X~7xQ~sD(>PhfaKGPdu#wy3}%JIY_5Fag-lAM zpk1m|Qe9;ZL`ArvpgmWApsh_M33{za!o)qi6H1ymkF^PaHSxE%!TL$w0(t?e9k4bO zgPfqb-Z+C}Xqxc?_?3k7$J=S}YZrY+V~tjO3b@p{aB(rQN<0X{hR*GdnvWl0r&ixz zpAr<*%xssElyrx%j-fSGp0=;%?6aC*3&~I!fGh;paGcNL8^Ci29heyi2ZaS79Du?x zFfmaFd_8S*8)zC(jDzm37`BFT{x)V;Ayf^CH#~Ic5CT_fR8&v=B`+ZSp@D&h=c;k> z0+@Osh~xHCFzqHLHt}4CY5=0TbHEz`m_X|~md6PUw^?jG3rIJEp^c3=kQ;pZq~6lv zqUr|~giWzOKfKB-k00OM%@iC=9|f(RtOOv4ke)*T=IhU&V>L7~BCu<_Q$e{wan<$X zECPxl;);pMQA*0-A#zT=@^MUuPESW~P8ZsZDUjWDaQJdQo$@0Hli+ulH@t#zcYsLQ z>QQ)ab8}V(oOCCt9Jq`MoZERLspXM!kUypOQ@9uz2~spz2H?BSIFHG1tRXULP`p8& z#G-;RMS_Y6;|Eit(lRp2QUMJd=MWp4-oG~hR@SB13&X;E0<0CZU)}Et@ZKr&*^2V= zC-`rF2LO1;BplP}-sQjpcc(IhVKfcCLbyvI^{j5Gsze81ueh5~xWsMl)tfL;6j!85 z1UVWO|IQu)M{GCYx&+n;=*@VcwczJh7w{@8D+A9xAA%FY+(#(ph-cd&$|JlcV6yr) z7@Rp6+tpZSPTkTo3>YlcS`Gu#PykM3qw2WZj(u`Y);rq}<`?$Km3Bu0-@k`N*HqO7 zQ-fE|ucD-)(!lJKib~u4`}dDi_2=jckL(=GcQ1t`m*s&uI29;+TIBTsjRys_;cXor zOu!1T1wdHkWn~*Gqk$|TT;Tt{eED+vF%R<4w2d?DRPq*>Rp;CNz{qV6${~qDr+OP& zC1sA_q)t6oAYF3X$+`4bLTol68HyRw)Q))H)WmuIe84-^w~0wf%yQxJ;>Lg+8Qxg; z32i>Jsa#Ec-?Ro0Jw$x}Xj6>ka}wN5*y$s%zlYKN^0QU!fDt1{k#5_ zw}M@fixU&A*@{LAV`N~)u+vuSReE;}_%Sqi~THc$LC3NboZ18@*O z`AvTNhZiY_Og0 z7o_Iq<|S5obnYeu&JsY}T-A0ZR9)5!Zi_8f5#AA*nko;Bc=^`Gd2XeDp(k+fV1i8 z>PV-dvR1QW!?njeVNadbW!C4nsAX?4Yg5(<9UB!ed47J`;Z{T*j}LXj$@Ql1A&jhN zv}Jtxk7yOl^{sDgT@ch=>9>{+8(i>kcJ8yAwQCYx5o~r}i?xe=tPP}IzzA%!wcA4w zkdK}`adCCEvtj{u92F37s#1}HA_xjkQc?2H-xDQ$3EDWk<$Uthce{Tb`AK)+mZ>Qz zEx(+}PSjjE=LI zfrk5*_r_^Kw?D7{#AR-5xG}tW{raVXjk$Rp-zW(wX_%)|Yb*rVNH7e8SA-&<8js{P zM_J$uMzVFe7xzX)#3d#tCnqJHUMPRZoSKo*4bqo8h+zL4s>H;^h1Pvq+S)xZT`0mv z7r*b>vj^m0U1Q^!L}+ct$IYzcwak$rA$d8QoyGlGuhK0dER1=fFDvUTZaVq?{l>E; zNl8i6H*jf*;MmxluX>5VW|2uli*L4dbUWPxz`d-fT8{CW z9q!!!#QV|ZpI{Ad=^d;)MR~vzc@;MAfc=e*I(c*FpGjVuejRD)^tV-eCC=#f#EN^0 zF*85)E*`Vyo|u}lbBK<9og{+kxszjkAaBV65)w`gDClS*TTJM79RH8oM6z*{E}n+X!5u&}U#f+eiwQ&5H5>Kj-lw5g{Va7hrWN@}_Ly1K{% zBpG}bw8qCzoIsR5OU5Q2!A(#bN!Sm4CIL;AW4DV8f`dMga9_Xv!h;LdLCK{ z7#qptsBV{q&ZwkBkc@C|*43+HptDi59j-F_DOB2lYbtQX5yk>~-M>Px3Bne<)!h{; zfZ@JO*d93qYvs@Ffhhi6TcwX|r8{(zFf8*t!qFt{_JLBY)xodYS+^Yge! zWJv>L1QMBI__<3>j}PojoxC*IB;s0il%0)jkm0VN`L&KKyC-)N=JxfGBQFzjo=?7b@}-@sn?IcxiH1RXMuv<(vxEd}cF5+@of8@CZlP+b zsvk|MLL=fN7D417IxqNw!3b7p!#!iO=Z}+V-nj7%5`y#Rl?E$C9aVc%?x7~pDzeti zQyH3?(r!{w!95a~@`W{dXl%UG=nc9Q5f>gE-8jM%*RJ21n`+TA@?}Km0l5;UK}bkQRv$`8 zAY{W{2lb-D!(XMRV-D!$!EngR<6!YlE-pUXLYT*Zw4#ZxcyxB^(c{N=XI61f8^;(e zRIzVGni*H$Ul{iF_9ooLfc}CQjB<7aCP{Ie;Ez_B+#!caTS;pn*dgu{fu1|fPJpTOL zI}wcS#VA17IC1kz)@`@~O8t?=z$w`89~xY21y<57m%yfr?xqyjva+@1F8cFDVBk5u z_!h;O9%>ILI#BX)<@L9<A+ z9icBnCa{$W{+P`baE8!+lfgI_Reu(Kv9CnqN}vq3)lx`|<>yY(I|nM-_NvHL*41CH>`wYRtLMU4xK zW9(>1NKU?GXvol4mK_bQGsU&fPWR`>nnn|N@#dB^ zA`?|+U?Yo1c!Yz7;A5rIxo`i+rj;HnA zJ0BoCg8#xj$4G1H+E2v zksLdB6JrYGqQvN-7(GP}di21knVFb;1PDzwiSkw8z)#^kjJQBVZK{r^Yz(IB(pAfQAcSI;>B5!Ez&!NL{@OfyMaLRD*P#U*&1jKP z9?a>T|3DbbS_1_r09W`MTXW4;R(ve1L-9ACCqF3nFgEZmAHOFHeLpJ8}68HS1L zKCnZVBp3o%goUH3tTupNAeHoVb=}?327C#xb&-M{c#u76G=n`4Yf47)wJl z){a5~O#~WrGqjNyTo~3`OiM{w8gHq7_wJ&sY$gCG1g9~9UbHCCX2DbNg^WBjJWNSX zUj+#nY|+Bfl9rYhncw)U?Hg9fq#|LlK!Szlj& z@7}$rnUMpoUadiHLOV$iMrv!9K}~)7>I*q(X@lp+CMF8_RQLhDL;v;7OAscaBK|Pb z@jpuY@^CEse_Kh3lFX!1=BbP&nd_-cg-lVVP?3}=nW;=68AHZ0g$NZgHc&)ls$`xr z^kyD2oORdl?7h!Edtc|A>-gh+uPcw|xu5&~em}$dthK(VV!FGH@!Ww<;%@dcvZ6zv21t=@Yax3ZPPHYBr&YL#rzyw)6|+KgcUV)`Mk6AwZ-xv~1vG*~?EM z5Fp&aIVBg(EEG|wsTGt6IVNzg=YWpdveNSM2;~r((KBAV_!qXskZsU*Ay0?8#{TNj z)s+?1@3OEpsH@Ampq*&fa~eqD8tomid2Ajk3Ua9VJYSu+0$9Wb`t@QQsHo`e1rHv4 zVZV{LhE*!{COFZ9>gzDka^#Ey;xl`0OHA)=cu-MLPz?U~p>O&fP7!DtrmSqU)1e~+ zM;5-9L--ovl~%^a-DuOWgKnT(2kA+}LFXnnrRLSa;oAPyXHlNDAF&;e>vw%rs2IQgO=4%tqpy~(##JB zqc9h)6TxRIJw3onIdN~-6Iif-+_Ym05;&w4DR9rUw8V4M>mVCh>6?=4dKbbbF^ob4 zcN=znG*ULncM#Sgf`X@-akIvyMLDm*z}tEL;O;F>s&SA`oqGD_yQgYaRu)i~2N8-! zMjmK)l9r$Ccsj71L1!v^2Hze?7tw|_m;jdCy<=^|5-;39z040Gn=;!Z02jUO@jq zd%>c-?oZXr5Fyqf;@-ab0lB>CQCpq&GZ_J!3zd|Vz#Ie%X~;OdX2+`q#58~F#% z7CR%i@H-MtvDSf-F9^g>wH|Lyaze?1LLJxET&RgE-w}`vgw!yAHN|_ODS$0!PDVx% z;vs+&l}qOMcytfCy1Gju#wcvijYLMi^L+`r5=7@H^vHb__wGPnfjT51!2`TA+A&k# z6Fp>6lLJ*87b(;=tqKZNZ2bNaIfpI>BsCz->FF1t_1SsQdR9gJ@HuTke`ueufn1u zE>3eZ2a+I}La`gWNWu8`_gjk~%d@l`LU;UO*Kw62c5u%N3fekNLxxYU(*Gm0JRM(* z^?;mcQtqm2SO7K|5YMZCobfQC4z}zf10X;qaG5J-a$BzHc3L*^LDVI{kRhA~@IfF^ z=rIwa&?-V9kGdC#uPKDwVPW8`uVdKoXxm2b!O)H6`B}p~8{;30Ex+H)clx!${m>OF zMrLLh7Yc}qio(sF_{S_K`R=*VY8QM^<`L}ce}DJ?bX^tO zV*J(m*n)SRVB7p(T+4ra^}k(Y)sI_jWK4;#KbZab^2B@N8{+4G|6;^17&kLK5R|n*Zz zxbMvW^-cc!Hj^c$n-zWk-n)N2fdBZVasBfCz=4VJKw8>9A+lZu6u2G0u9`fz>)hX@ zIftHi9L+d}MhHeHZON+mxEQHV2D64GBwV$!`u0Sd&teO%P4FtTA4gNJV{)Xss{S2nY!gj=CboWzU{JC)v9Hqp-bIUiY7g4GZKyHhJ0Lvaqn)e87jm z&z6`9dDzMVl$&PQ;(Jt@iq)K><;|dCQ_@SrNx$(n?%g8(^Ncy|yW@aABx!r*B0zR$o?knhv zZr#2;0i7@m8HmYk04B);0|Nw5!Fe|%MEA^@yFo!R*L%q-*OB61m6SL^E4gt|IsiLivF&BO4) zU}t|A9K6l;l#9dYsJCv#@>3=OIqcQHcnRHw+`|A2RA9ta<@FwNL+H8JZWkQbzrUfq9d-%jo&A^y z&dtq5L_`SPSkv@lTHUC`%YP(D=$TgRoCex1jLPsgpRP?qN(ysBB-g1c-4$p0CQ}gN z{+touF+MU96BSjxAH$cpm>4*h8Lr%xPhTU{ijaM4@((g2*&E=Epe$Z8*)JvK{QA}= zlH*VvlWmAHH*)zH1T&D-c6D{h$;uYbp-v&aqfwl68RO8L4c%xkX-T1uCLS-u6k{b0 zQwDl?8qXA0@{l85c$E<`9gEQcU~bO4DxkYsV?%16GND~n=unE06VpCkTDlEiS) zb_xw4hRnePr5n1Eld_hsuCF2@7);Ewva{7%d6A^90vYl2tXR2*QQ1B^QM)8H60n(} zV~RB~<>elT6iY|))j7k^X4nGq=}&>fXa{14ko_9RVJZqB2f9BIuL@oJ)s57Chty4OJ$m#Ao4^tMKRThEWvUknk%MYvDE%Lc^3}l=40Vu6z9C$)IZ}OGzaBW`G$Y^^)%$*dpld z2w0e7i9nDZ+$KE-_a8k7>Vm+Dx|$kZ7|dWxO;_$+2o%@1*H~DL!tfIDZ0H(UA9}x` zp%*A+JWt_X9h8)Wu25(;fW5Jpsnd__&%K#l{sDP5;t}fA@UKnqOX219+v;81QG@B* zij@vey8_}@;j?F#QD%@<9s1t*idULi&1o1LCoW2y%v}R=&hOPu+=JwJY^=`fSfodAQL=qi2B%CLn3^(fF}9!{WGd-lNpX9h43y%wAX zXZ^8um2+ZI?84V#d12xjEIY4U`3?*kW9Dsa;&ptM5UOwY7IzXWWJyC@#~XCogwt?JHi}S-?r^LpKy$Xk_NauDKaSu zq6b#cBGBIn%xPF#TU%HRf_ngW6rKtdAVV0Q0pbSoE111j+t9FoeK^xAi>1?SnI3~G z`v=--JU7=?r0HOmacuc1W)T-HEZ)`B*azXG(Ja_vb))ebZ+(Kvr?Q}>CjS(Nw3ibT zE9>H0S_Oc)gn1vpbO?zi=jVsJLjnStK()V_l$JIG%|r3Wrn_3&+S!wYF~xQonioYy zEHF`9^N5X-A=9y4CE7kaDIQMF3uughu>)N`Bq1>nXYb{{%0Ne_L_$cSpKC0(8^Edo z!g_m|^1A=%mCxYRaH7t4vZ4Z%qev1L@xTgz1yu3Y9N<>S$JkIa*h`(f_R0$sfk8o6 zzCP;%S*G|4=Dh}(3E~|wsH)Z70BBm*+SUjAUr;~*n8h)79%lns5SvB{k@>bOaE1kp zXTIJq?>1lmR|bW?A2mFTX#3U{m;Hkj@MSgEtKL`bP8U%^BLYm`xa{V>Rpu{o@R7b7 z08uVa_fW@$b7#+HW@j^ehQ~L!3sC!-$s5vhXbfTLOkZQf86%>Ym$-Ey)??%WJmMf0 za%lzODWazrJ-gAf=IrUypDxuFA#lF9`dv?7U%zZKflREgt%3%lnM=lKxPAw{mQy+m;@AGj>6pb^T&_bxHJ$Th8eCqe=Jez249BXtN0@_ z6LO4NYbu7l1HTL#A4+K$8Lgsyr9TJd;sD+~NxlUp796rSn?UZEilW@UOIca@x|7q> z-qQ-#SC|%w>Q+cLNT__x?#2yU*g|5G&wWT1)-l-n=#*7R$e&p6pgyon#Bds~U;igo z@cg>)xtGG++=D=V>CZL9OV8Xk+^svhOu0i;$pyyEkAVJGu1p6l9*oV-7V1%X*U)gv zz~Ea$Lx`6Y7uUd#B(FWFdtee%d=OpcDOe^12KKdJn{8u16kPhv80!h%+qTCz`ya{$ zq6d#%WWr-S3)rt@I2Kt7kNY>8rILAyc$5`@xDXPI)`|w?wfavm7KQ&dxxgbOf6Et{95SzI&!`h&_wiTOiAXwXuHvqi?ciH6|W(L2?-06R%a}Y zPnq(5%yQXa-D%HTi!cd*T)V+J{3BlG_pm z=CuHxJm6G>Gzt6)k8g4yQ~}ijP&Ps>NuZ1XH!zaHnvPFMcna?zGcz%q8G~Vz!$Wk* zuK-ve1_3xm1=*kb19gX$jZMARJD{MTK7p$Rfbb*k9+nWwr`l_WF;g7@JO zP&Qt8#mZ_G;%yWmY7mgOwLytGhl_=13l=cHCMHZYHBB#H{tELa?Bcp_|Dd8G3TPy8 zR3@cJLFl6q&G7{DpFLa5%yw8wW?b2*6qPn9qs~;5K(?FHO|Msj}q%LF1V1GO)Q_WE`6!*IK4{r0T{m<&Oa zFnk^(X%!wXG6!%0H3EPLm@{A;jRSv7umk|QAvrGqMnkbh2y`!A1a;*Jz(u%aVQ_hp z+$R|TAHW2{5-TIu3w1wvIxQz7H6(J_-`i+t)YPivS|7wjD$3eW)6#xvZiaOqj}(~- z3LoURcw^CB70$NUXTUq)g@Zq_NLN37`d&~F5~8NA?&^61SbY}7$M5>ZjEDRTg@NNq z)JR0HuWP%CQjTCMp@#2RIZh~pjmxc@k~e&bjxTN2*4O_rA?kq`&BUa8_AI~GnMiAc zb#%EvFGdlGVY9RfV^`R)JsDfc-$0?>CTDu9(O}*CRtgP+dr0wS3i{oypxP6E;!ThI zO}~XAs>gZ3DDR@?ec{1p{e`S9nnvRs=s`6{5&2<$R((?O$*#D#59fTFN!x|wgTg4= zUaW1{bnLdb;0^gEztIR-12(m_4V{aElw%Da`KJl(ecb!sG~ntq+5ZnMHEWX5|8-le zXeLVg;Ee>Ubg=lgf9sNoUhO}1lEgo4VDNrpVY!c%G5$~g?A>_GyH+^6ed*?x43-}_ z(JFbYYWM#N9pI1tQ)4};TlFrDy3uq{v`+ND_>BMfY65|vzWsk}*S#B4f#(SbP%Ztv zepM%=89+AxhhJu9ZYlBcrb_O3BYXbuhfeA@AfyJ6lar&uzt_3gs<7R$CtDQR*5Yvd0UUfn0FF_x0}he@fMgd&b=MF7HWAs& z3&9}4^sVku>4v& zI;mM%cQ20?Dn7mj6lcsQi7pLwbU^f{>-`DESf%U=#yckMi$Y>ze=kckfd``?^HfC5s1I(Jw4hMWLzQa zdDU9*_vg3oGmReZ+BGtw6`r`~$4`Tu^5^H1zHZ;WrA^9RbF{5RTbaFBK$~ASzO0P< z%kSkg4-S~9p8$soS(s%L^#%q-r~6U5=jS9&T8N}o2H6eqDh+`rqXSC?(&}FFUX(Pvy!6QOS{N)E z3$S7MJTx%Cz|Jl-7ZL=J9K@(Ve&FK2jnpm$%AGuM0`tb`o-7!P964ghpDml3o({&I z0uGyc_H0tyLfd#xbu|-9=pBl0$0NmDmkL@rBcqPX%X3FbaA>+ZI0zhyo}N()Ps@1l zbFwY_%*RU?dAmxq$6%Fts-r{PR=lfRt;h4+7+paM}_u9 zb7%Rvl}ENl););=fnai-DvAh^?!$6U@R6G;$&$ zBBOFZ5U-&H1c4lL0m3>F=?29b zsWl;i<2c;v^li$BFuJ+%lX=FU`>P#SkCYTQ74@94x3Be_p0TRgW>^-ZqcbvTR`4Q< zF6fGN^T=n%Y3{U)O6U5{w})f<*5aiH`;eHvz9JYg(6&ZMU z(3N|rUA)+uN@ZXWwcI^X28AQv=~fyA?jEtDn7bqjsmh{)u_K?Y3M zA4$Ohn>Iz2y|6Q%+N!2@(R#VtpP5-)w{2nPlIu-T&X@L zK>F-sM&nkDvst7^QIlF#3qDLlNGw}cPmJReYes6r|}kFY22Y-)a#BP=e4vj7eAZ|#xeoZ2c+LsRe^VfNyQ*jPqp zy|Tk@hu{I-nzkGf5FbPWDDIX_<=#LJtS zeu9>RS;V%3q@^Wk*YDOS(%EV8)1{wuMoXZ?@3b11QP-KmqsNaIFg<;Fh%PMGR;B2w zL}B>b%$1DC2E>%T-1PPGUG3UKohkXHrtqYH6@In1R~bZi~OR0OlsRcmX{Q z)+pJZKK%#XPW7isJ&Z3yLQVm}1_U4n(!R7C#3gIPU2IZPQjk9Z&>3RJpcA0SnI1XM zF~o-(H*nP~OvbK(%DSNSU#=p%8_GBJ*^!fLwChWl@X_^^V;YHN&W%s3v8Prwz7Kt3sb=_!asPa;#GAwFnsP6!xZKKbUeU|#qFopa z5(9_6Ce%wLKb+pcu(>)22X0P0#W`_{q3^4FW`9jw`DXY+kg6*(Qu)-=!IbYI?A0FQ zH%m$)W0Xfnb+X&tWW*#f=sa@Wdu=5`M|ZS7Hhgc+3b(k@@;4J6X`n^G>96N+9z~x7 z_e;wmZK1UkKa3^f;tXGxa+O{0yM^D7*u#*B;n|8qzmBhWWhOpVV{=w@ZDhV+0F4X{ zXTL%={twiTU`czTF?xcDy&|jy*Ui)ZKI0~hvyPb=x7^XT)+e*$d9#C59y2s@UdHxwoLAp%*(r7_ zfS=bdQJ%-7r~l-ge!aZ(Qs>&wVHHbo8diEinZ)wQDKl}AIfp)#{iH*Pib-X9jxZb7 z+>>nW%b8`HGgf~;Sh%nH`|>Z1Te%0GWJa!xcSNnMSV&F8i2H_zk5sPTUH5FwU40$P zWyei@O756PVJ&*yVYlC0Y!KbtGaJd8USQ6!byFnUs4DkDP!c$_`Y9?5{<;L8iP*6qQ- zzzZ(z#pI)>XC2-Dz!w;f0r{Q{n~^FI26|1`PiuKsL+j8&^w-+FJCzBT_YkKqVg|Nr z*DmfGA6d0+CUtxg9>#_1d}7q_=Zu!PD_Gu{!OH$1;jBrLA8p2MN~<~RGrCa=6jy(p z?|2ckNAA(ZhJeh$c_#CTp~hSzROu)uTc5If+zWe=9S24~pcr^Y?MRrDm+LgUm6w_- zq*dw3li+vQ*LQYNjAmEnq()I8s!M7rmuo{yGd#5CS%pfH_4U#N8Adw}9=Unn#wN?= z9!j7f5+{2;7vu~@;01C|)g0is3=tVj#RDgtjsDX5bkk@(3&FvLE0??>4E^>C3UWR| zafyhoQ}T0WG=o@&X?d6KrJ5qx$-0XP*S4y4OzEG~9VAQTQ>X zkceEuOPiEHhM>nhkHYgnK~b4}w@Oj5o~ZW}90Jvv?wo0H2?;=#zvsrK;Ps8U0f>iH zaD9hcPtL9)H;g7_0S^Ld-@a3F#0wHo0EZ+_Iy%5DM-Cr`Qo^>Uq8#vxRPo6)-A14* zz^cu0UM9dI(DMk8@cHLs#Gsi~%J}pNzd_d;j85O3Z0R$eIVJY;>f_?Wudb(Xzixl3 zQ*_|z*(qbW(>n7zPWk#qK2+oU>~@8F^R$4l@Y0A?%ZD02SQ0FihR3OX$qkD$oVsAT z7b(_8#zkE0#OBQtM%Eb_9~awweJ2`T)$P7b;te3BdMo+CttSoEbdiy>2UcXfFo6Mm zXi8&5wYKDFAYCY+p_v{_MycQLTs>#^a)R*&Z%erga<#bo5~(vOk#a$KxY8s5o%^T4 zVslL@;m|mPnJmU175R~~@=Us}ABB%9tde#gxqtT4{XCJS<%%!wD^}7TRLWH62?)B1 zjL@xL-n!NBnfKaqcB)U9GKu|AWal^PZSHv<#!2$1T`VojrC3CY_E%CmIW;%&?-30S zt!ix(;PV%k2zqFdAb$!7q=Ig|guh`CFC`_LRA6=@PU>y!?3Cr^e%IBt3JY;;fl@SN zAYZ`M4*1?a589}O;RbdLRPng(v5NKMal#RF7s@>pKv4Dd_4gB&c^G>UCl`C<u#XR6wAs# zQWq_WKlwdx(=8jFU}T?ZRO`uG>|M4WJZR;s$ z-Q3)O8N)bt!>li&tg(i>pkmh?^(A}Qxhi*Biv!k zKLha_PASKuRoBoUWEX@bda1KW5hcit@Y6+QsQRIxa_xegtLW9M?FBq|HWd}_?UC7A zEH4&}wxGndS@(Qtx?%iFX!*z3>1gAcuVpUPl|kEgf4#R`cz;O9o`$$XmK_~ACJHZ| zMcgL6rH_cq$)!JS@Y`nA-CfQ8>+5S9o3H7o>bnve;;N4xL#ux{X6~PJ@x?}^;aaqZ zu3cmiiB`XT*i}+O(&B}6{$4wE!%)Z@q3i7J&Wnw z^MKvbZ;m*c8XDGDc^^Mv|5ji0*fBiF`UZs?HFBG1g$TYI5;vUlPcAD?wR7LPzg5vayjknRgn}zxR?&xxZ;sKA*-}C8gSH zg>?}OqJXIg!&kIqX;&YDd<*CdV}EqoYGu_2De*LOX7#|Z<$i<>toKc}1a6YZeH9gJ zfZ$VALo!iql4hMPO$q^}d0Y1abm8J0t)uG=(l@M%#F@F^euAj6vg9EB5N3U00W~4d zH0n~x$f!Q5kd$=T`6|S@xNeWr0BJukdJ2uADhUZjxR>i=7Wck=DFp4&ab?jimI~2C ziG)F8?yrnJ@D57mq=|tr(9FW(P7YZwh5!m!Zut1C zgMc@SRtEw>5G)R>LAv8kNIYkNk7(xIw0s`@7-$j!m!Yr4G+Vp{(?hVf=?K-%6_t7l zUxI|zoBbyArC~(U!YfCBpQPFOf6jSEihytoX|Qhqh>yYviKQ$r;M|vHUoW?Ufh&-e zI^x71zXp>~EQ96BzrDlwy4gykLva?fgO^PaAZ}y<5M~1MTAP|;;Pi=cGvUc)r(=&AaTyTyREjVYM1weau_%2dr}jdz9Zj*hPzI*mb|R|mu7Z-1gfN*j zi)Voe0?i+4WRp~nyTsKF$uQ-?u`y`!!jzvQ?t+1B`H?v{Z&xFZG{-|7vLe?t1M4^(h^7WWo3b*!LQIUAXZbnS)Ih~uwqw3Ob=m3fgu*(xEwi8hcD50<3I#h#aKwdMTa;t{QNf4RaaN5R35}gmXolG zt*0rWawHCVNbDqZW6nC@t&g{4G@^xl=gLO?i1f)R zDOwq)cO$%`Q~;kHgag=d=3}|_n>7F;AQX8hdqFEI#B>CdfuOs0U*VW{)OYBO4Le%P z07b~cEfYBhxw-YME>(%_9okAiDkNIGykhA z8&n#xIz=liVHu|}I7Amx*RRzyczoP+qPbY{?$)fq0ZPXE666ejrd65h=<4wO+S}Up zUs^_L2iponKYd}-e%B?HkF~Y?ICDpxUC4QT6iv;|!4<)pllR5NzOSyf(+-GK;iiX- z^BK6;_J*V=zlW2jhbiHh;Q)4NZ0t(HB+~E;Yu&^V3B>Mj77R)YKGWAyD|0xy2$fR7syajUrHU!oeF5I$F zBAx}=Wd0vV)c?s(|8ryio!O2T{2e&!n5] (0,0) -- (9,0) node[below] {\huge $\log_{10}\rho$}; + \draw[ultra thick,->] (0,0) -- (0,9) node[left] {\huge $\log_{10}T$}; + + \draw (4.5, 2.5) node {\huge fine region}; + \draw (4.5, 7) node {\huge coarse in T}; + + \draw (1.25, 1) node {\large coarse in $\rho$}; + \draw (7.75, 1) node {\large coarse in $\rho$}; + + \draw (1.25, 7) node {\large coarse in all}; + \draw (7.75, 7) node {\large coarse in all}; +\end{tikzpicture} + +\end{document} \ No newline at end of file diff --git a/doc/sphinx/src/sphinx-doc.rst b/doc/sphinx/src/sphinx-doc.rst index d0c347139b..552e4c3c2a 100644 --- a/doc/sphinx/src/sphinx-doc.rst +++ b/doc/sphinx/src/sphinx-doc.rst @@ -27,7 +27,7 @@ Building documentation locally While you can rely on the CI to build the documentation associated with your branch, you can also very easily build sphinx documentation locally through -python. These instructions also _do not_ require admin access and are usable +python. These instructions also **do not** require admin access and are usable with shared machines or python distributions. First, ensure that you are running a modern version of python (i.e. python 3 of diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 83bdde1d35..adabdbeb33 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -404,13 +404,17 @@ A modified equation of state can be built up iteratively. To check if the equation of state currently stored in the variant can modified, you may call -.. cpp:function:: bool ModifiedInVariant() const; +.. code-block:: cpp + + bool ModifiedInVariant() const; where ``Mod`` is the type of the modifier you want to apply, for example ``ShiftedEOS``. If this function returns true, then you can apply a modifier with the function -.. cpp:function:: Variant Modify(Args &&..args) const; +.. code-block:: cpp + + Variant Modify(Args &&..args) const; where again ``Mod`` is the modifier you wish to apply, and ``args`` are the arguments to the constructor for that modifier, e.g., the From a9f23d2b944a41b9d512c2abe1a4a18a7bfecd82 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 14:29:27 -0600 Subject: [PATCH 28/45] final chad comments --- doc/sphinx/src/models.rst | 17 +++++++++++++++-- sesame2spiner/examples/aluminum.dat | 2 +- sesame2spiner/generate_files.cpp | 2 +- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 63ac375c1a..e097c846e5 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1585,8 +1585,21 @@ Thus the input block for piecewise grid might look like this: .. note:: For all grid types, the only required value in an input file is the - matid. All other values will be inferred from the original sesame - database if possible and if no value in the input file is provided. + matid. Table bounds and normal density will be inferred from teh + sesame metadata if possible and if no value in the original input + file is provided. Table densities and positions and sizes of refined + regions are not inferred from the table, but are chosen with + the default values listed in the above code block. + +.. note:: + + Both the flat and hierarchical grids attempt to align their grids so + that there is a grid point in density and temperature exactly at + room temperature and normal density. This is because normal density + and room temperature is a particularly important point in phase + space, as it is the point in phase space a piece of material sitting + on your desk would be at. This is called an *anchor* point for the + mesh. SAP Polynomial EOS `````````````````` diff --git a/sesame2spiner/examples/aluminum.dat b/sesame2spiner/examples/aluminum.dat index 318cfc5706..2d26089c5f 100644 --- a/sesame2spiner/examples/aluminum.dat +++ b/sesame2spiner/examples/aluminum.dat @@ -2,7 +2,7 @@ # Example input deck for sesame2spiner, # a tool for converting eospac to spiner # Author: Jonah Miller (jonahm@lanl.gov) -# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# © 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 diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index c497334651..57caee436b 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -288,7 +288,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params int numSie = params.Get("numsie", numSieDefault); const Real TAnchor = 298.15; - Real rhoAnchor = metadata.normalDensity; + Real rhoAnchor = params.Get("rho_fine_center", metadata.normalDensity); if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { std::cerr << "WARNING [" << matid << "] " << "normal density ill defined. Setting it to a sensible default." From 8db13ef9f4a09c87dd2174e38fed579e8c96b893 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 15:17:55 -0600 Subject: [PATCH 29/45] crazy --- sesame2spiner/generate_files.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 914ec62ed8..b865425ee0 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -32,8 +32,8 @@ using namespace EospacWrapper; constexpr int PPD_DEFAULT_RHO = 350; constexpr int PPD_DEFAULT_T = 100; -constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; -constexpr Real STRICTLY_POS_MIN_T = 1e-2; +constexpr Real STRICTLY_POS_MIN_RHO = 1e-9; +constexpr Real STRICTLY_POS_MIN_T = 1e-9; constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; From 805e16837415f7d0de374e34767198061e39c973 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 6 Aug 2024 16:06:52 -0600 Subject: [PATCH 30/45] figured out error in tests --- sesame2spiner/examples/unit_tests/copper.dat | 4 +++- sesame2spiner/generate_files.hpp | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat index 16ba63d428..8f125df4ba 100644 --- a/sesame2spiner/examples/unit_tests/copper.dat +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -17,4 +17,6 @@ matid=3337 name=copper -rhomin=0. + +piecewiseRho = false +numrho/decade = 50 diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index b865425ee0..914ec62ed8 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -32,8 +32,8 @@ using namespace EospacWrapper; constexpr int PPD_DEFAULT_RHO = 350; constexpr int PPD_DEFAULT_T = 100; -constexpr Real STRICTLY_POS_MIN_RHO = 1e-9; -constexpr Real STRICTLY_POS_MIN_T = 1e-9; +constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; +constexpr Real STRICTLY_POS_MIN_T = 1e-2; constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; From efc894c262009fcc7a689bc855b7cc50cb10c18f Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:31:58 -0600 Subject: [PATCH 31/45] Update doc/sphinx/src/models.rst Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> --- doc/sphinx/src/models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index e097c846e5..98ef1a9382 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1510,7 +1510,7 @@ default grids are now piecewise. Piecewise grids can be disabled with # defaults are true piecewiseRho = false piecewiseT = false - piecewiseSie = dalse + piecewiseSie = false These options may be true or false. The default is true. When piecewise grids are active, the density-temperature (or From dac88bc8ff6b6ac8564423c4717dd052e800cd95 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:32:10 -0600 Subject: [PATCH 32/45] Update doc/sphinx/src/models.rst Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> --- doc/sphinx/src/models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 98ef1a9382..e2db07f41b 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1585,7 +1585,7 @@ Thus the input block for piecewise grid might look like this: .. note:: For all grid types, the only required value in an input file is the - matid. Table bounds and normal density will be inferred from teh + matid. Table bounds and normal density will be inferred from the sesame metadata if possible and if no value in the original input file is provided. Table densities and positions and sizes of refined regions are not inferred from the table, but are chosen with From c103b513df6a101263a27e7ed1949050d283aa05 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:32:58 -0600 Subject: [PATCH 33/45] Update sesame2spiner/generate_files.cpp Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> --- sesame2spiner/generate_files.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 57caee436b..e2c62c98d7 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -287,7 +287,7 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params int numT = params.Get("numT", numTDefault); int numSie = params.Get("numsie", numSieDefault); - const Real TAnchor = 298.15; + constexpr Real TAnchor = 298.15; Real rhoAnchor = params.Get("rho_fine_center", metadata.normalDensity); if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { std::cerr << "WARNING [" << matid << "] " From e0425fea406b51a6405ccfa77702021def0b98fb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:35:20 -0600 Subject: [PATCH 34/45] add const to piecewise vars as suggested by @jhp-lanl --- sesame2spiner/generate_files.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index e2c62c98d7..88339ef433 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -297,9 +297,9 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params } // Piecewise grids stuff - bool piecewiseRho = params.Get("piecewiseRho", true); - bool piecewiseT = params.Get("piecewiseT", true); - bool piecewiseSie = params.Get("piecewiseSie", true); + const bool piecewiseRho = params.Get("piecewiseRho", true); + const bool piecewiseT = params.Get("piecewiseT", true); + const bool piecewiseSie = params.Get("piecewiseSie", true); Real ppd_factor_rho_lo = params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); Real ppd_factor_rho_hi = params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); From c1d2040c6b8448c4dc02a25464ae383163568966 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:42:49 -0600 Subject: [PATCH 35/45] add a bunch of const calls --- sesame2spiner/generate_files.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 88339ef433..7c608be956 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -301,15 +301,16 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params const bool piecewiseT = params.Get("piecewiseT", true); const bool piecewiseSie = params.Get("piecewiseSie", true); - Real ppd_factor_rho_lo = params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); - Real ppd_factor_rho_hi = params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); - Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); - Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); - Real rho_fine_diameter = + const Real ppd_factor_rho_lo = + params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); + const Real ppd_factor_rho_hi = + params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); + const Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real rho_fine_diameter = params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); - Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); - - Real rho_fine_center = rhoAnchor; + const Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); + const Real rho_fine_center = rhoAnchor; // These override the rho center/diameter settings Real rho_fine_min = params.Get("rhoFineMin", -1); @@ -366,8 +367,8 @@ void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params Verbosity::Quiet); eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); } - Real sieAnchor = sie[0]; - Real sieSplitPoint = sie[1]; + const Real sieAnchor = sie[0]; + const Real sieSplitPoint = sie[1]; leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, ppdSie, ppd_factor_sie, true, shrinkleBounds); } else { From f6363d34584f1c8821da570d07d7f2768df1a330 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 13:57:50 -0600 Subject: [PATCH 36/45] jhp comments --- sesame2spiner/io_eospac.hpp | 2 +- singularity-eos/base/table_bounds.hpp | 50 +++++++++++++++------------ test/test_bounds.cpp | 33 ++++++++++-------- 3 files changed, 47 insertions(+), 38 deletions(-) diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 6ecb454d1e..fc492e2472 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -34,7 +34,7 @@ using EospacWrapper::Verbosity; constexpr int NGRIDS = 3; -using Bounds = singularity::Bounds; +using Bounds = singularity::table_utils::Bounds; using Grid_t = Spiner::PiecewiseGrid1D; using DataBox = Spiner::DataBox; diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index 3c6d0e98ee..e6300a0300 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -26,11 +26,10 @@ #include #include #include -#include #include namespace singularity { - +namespace table_utils { // For logarithmic interpolation, quantities may be negative. // If they are, use offset to ensure negative values make sense. template @@ -45,12 +44,14 @@ class Bounds { Bounds() : offset(0), piecewise(false), linmin_(0), linmax_(0) {} - Bounds(Real min, Real max, int N, Real offset) + Bounds(const Real min, const Real max, const int N, const Real offset) : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), offset(offset), piecewise(false), linmin_(min), linmax_(max) {} - Bounds(OneGrid, Real min, Real max, int N, Real offset) : Bounds(min, max, N, offset) {} + Bounds(OneGrid, const Real min, const Real max, const int N, const Real offset) + : Bounds(min, max, N, offset) {} - Bounds(Real min, Real max, int N, bool convertToLog = false, Real shrinkRange = 0, + Bounds(Real min, Real max, int N, const bool convertToLog = false, + const Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) : offset(0), piecewise(true), linmin_(min), linmax_(max) { if (convertToLog) { @@ -65,13 +66,14 @@ class Bounds { } grid = Grid_t(std::vector{RegularGrid1D(min, max, N)}); } - Bounds(OneGrid, Real min, Real max, int N, bool convertToLog = false, - Real shrinkRange = 0, + Bounds(OneGrid, Real min, Real max, int N, const bool convertToLog = false, + const Real shrinkRange = 0, Real anchor_point = std::numeric_limits::signaling_NaN()) : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, - Real ppd_fine, Real ppd_factor, bool convertToLog, Real shrinkRange = 0) + const Real ppd_fine, const Real ppd_factor, const bool convertToLog, + const Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; @@ -96,8 +98,8 @@ class Bounds { } Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real fine_min, - Real fine_max, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - bool convertToLog, Real shrinkRange = 0) + Real fine_max, const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { if (convertToLog) { @@ -113,8 +115,8 @@ class Bounds { } Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, - Real log_fine_diameter, Real ppd_fine, Real ppd_factor_lo, Real ppd_factor_hi, - bool convertToLog, Real shrinkRange = 0) + Real log_fine_diameter, const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { if (convertToLog) { @@ -130,7 +132,7 @@ class Bounds { mid_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); } - inline Real log2lin(Real xl) const { + inline Real log2lin(const Real xl) const { // JMM: Need to guard this with the linear bounds passed in. The // reason is that our fast math routines, while completely // invertible at the level of machine epsilon, do introduce error @@ -139,7 +141,7 @@ class Bounds { return std::min(linmax_, std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); } - inline Real i2lin(int i) const { return log2lin(grid.x(i)); } + inline Real i2lin(const int i) const { return log2lin(grid.x(i)); } friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" @@ -154,7 +156,7 @@ class Bounds { // This uses real logs template - static int getNumPointsFromPPD(Real min, Real max, T ppd) { + static int getNumPointsFromPPD(Real min, Real max, const T ppd) { constexpr Real epsilon = std::numeric_limits::epsilon(); const Real min_offset = 10 * std::abs(epsilon); // 1.1 so that the y-intercept isn't identically zero @@ -171,19 +173,19 @@ class Bounds { return N; } template - static int getNumPointsFromDensity(Real min, Real max, T density) { + static int getNumPointsFromDensity(const Real min, const Real max, const T density) { Real delta = max - min; int N = std::max(2, static_cast(std::ceil(density * delta))); return N; } private: - Real toLog_(Real val, Real offset) { + Real toLog_(Real val, const Real offset) { val += offset; return singularity::FastMath::log10(std::abs(val)); } - void convertBoundsToLog_(Real &min, Real &max, Real shrinkRange = 0) { + void convertBoundsToLog_(Real &min, Real &max, const Real shrinkRange = 0) { // Log scales can't handle negative numbers or exactly zero. To // deal with that, we offset. constexpr Real epsilon = std::numeric_limits::epsilon(); @@ -201,7 +203,8 @@ class Bounds { max -= 0.5 * shrinkRange * delta; } - static void adjustForAnchor_(Real min, Real &max, int &N, Real anchor_point) { + static void adjustForAnchor_(const Real min, Real &max, int &N, + const Real anchor_point) { if (min < anchor_point && anchor_point < max) { Real Nfrac = (anchor_point - min) / (max - min); PORTABLE_REQUIRE((0 < Nfrac && Nfrac < 1), "anchor in bounds"); @@ -219,7 +222,8 @@ class Bounds { } } - static void checkInterval_(Real &p, Real min, Real max, const std::string &name) { + static void checkInterval_(Real &p, const Real min, const Real max, + const std::string &name) { if (p <= min) { PORTABLE_ALWAYS_WARN(name + " less than minimum. Adjusting."); Real eps = 0.1 * std::abs(min); @@ -234,8 +238,8 @@ class Bounds { static Grid_t gridFromIntervals_(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real mid_min, Real mid_max, - Real ppd_fine, Real ppd_factor_lo, - Real ppd_factor_hi) { + const Real ppd_fine, const Real ppd_factor_lo, + const Real ppd_factor_hi) { const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; @@ -274,7 +278,7 @@ class Bounds { private: Real linmin_, linmax_; }; - +} // namespace table_utils } // namespace singularity #endif // SINGULARITY_USE_SPINER diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 5d0e2a594f..44b76eef36 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -13,6 +13,7 @@ //------------------------------------------------------------------------------ #include +#include #include #include @@ -24,7 +25,7 @@ #endif constexpr int NGRIDS = 3; -using Bounds = singularity::Bounds; +using Bounds = singularity::table_utils::Bounds; using RegularGrid1D = Spiner::RegularGrid1D; using Grid_t = Bounds::Grid_t; @@ -38,10 +39,14 @@ constexpr Real T_split = 1e4; constexpr int N_per_decade_fine = 200; constexpr Real N_factor = 5; -SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { +constexpr REAL_TOL = + std::numeric_limits::epsilon * 1e3 + + SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { WHEN("We compute the number of points from points per decade") { int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); - THEN("We get the right number") { REQUIRE(np == 1800); } + constexpr int NDECADES = 9; + THEN("We get the right number") { REQUIRE(np == NDECADES * N_PER_DECADE_FINE); } } } @@ -65,7 +70,7 @@ SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { Bounds lRhoBounds(rho_min, rho_max, np, true, 0.0, rho_normal); THEN("The lower and upper bounds are right") { REQUIRE(std::abs(lRhoBounds.grid.min() - singularity::FastMath::log10(rho_min)) <= - 1e-12); + REAL_TOL); REQUIRE(lRhoBounds.grid.max() <= singularity::FastMath::log10(rho_max)); // shifted due to anchor AND_THEN("The anchor is on the mesh") { @@ -73,8 +78,8 @@ SCENARIO("Logarithmic, single-grid bounds in the bounds object", "[Bounds]") { int ianchor; Spiner::weights_t w; lRhoBounds.grid.weights(lanchor, ianchor, w); - REQUIRE(std::abs(w[0] - 1) <= 1e-12); - REQUIRE(std::abs(w[1]) <= 1e-12); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); } } } @@ -87,8 +92,8 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { THEN("The bounds are right") { Real lrmin = singularity::FastMath::log10(rho_min); Real lrmax = singularity::FastMath::log10(rho_max); - REQUIRE(std::abs(bnds.grid.min() - lrmin) <= 1e-12); - REQUIRE(std::abs(bnds.grid.max() - lrmax) <= 1e-12); + REQUIRE(std::abs(bnds.grid.min() - lrmin) <= REAL_TOL); + REQUIRE(std::abs(bnds.grid.max() - lrmax) <= REAL_TOL); REQUIRE(bnds.grid.nGrids() == 3); AND_THEN( "The total number of points is less than a uniform fine spacing would imply") { @@ -98,8 +103,8 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { int ianchor; Spiner::weights_t w; bnds.grid.weights(lanchor, ianchor, w); - REQUIRE(std::abs(w[0] - 1) <= 1e-12); - REQUIRE(std::abs(w[1]) <= 1e-12); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); } } } @@ -110,8 +115,8 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { THEN("The bounds are right") { Real ltmin = singularity::FastMath::log10(T_min); Real ltmax = singularity::FastMath::log10(T_max); - REQUIRE(std::abs(bnds.grid.min() - ltmin) <= 1e-12); - REQUIRE(std::abs(bnds.grid.max() - ltmax) <= 1e-12); + REQUIRE(std::abs(bnds.grid.min() - ltmin) <= REAL_TOL); + REQUIRE(std::abs(bnds.grid.max() - ltmax) <= REAL_TOL); REQUIRE(bnds.grid.nGrids() == 2); AND_THEN( "The total number of points is less than a uniform fine spacing would imply") { @@ -121,8 +126,8 @@ SCENARIO("Logarithmic, piecewise bounds in boudns object", "[Bounds]") { int ianchor; Spiner::weights_t w; bnds.grid.weights(lanchor, ianchor, w); - REQUIRE(std::abs(w[0] - 1) <= 1e-12); - REQUIRE(std::abs(w[1]) <= 1e-12); + REQUIRE(std::abs(w[0] - 1) <= REAL_TOL); + REQUIRE(std::abs(w[1]) <= REAL_TOL); } } } From 7ed764f51a28cd0c6cb13db3100bd6a411356e26 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sun, 11 Aug 2024 14:04:23 -0600 Subject: [PATCH 37/45] oops REAL TOL not properly declared --- test/test_bounds.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 44b76eef36..86fb664f86 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -39,7 +39,7 @@ constexpr Real T_split = 1e4; constexpr int N_per_decade_fine = 200; constexpr Real N_factor = 5; -constexpr REAL_TOL = +constexpr Real REAL_TOL = std::numeric_limits::epsilon * 1e3 SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { From fb724a62cc441f4f5fe9a52e296bf4dc87188807 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 09:57:27 -0600 Subject: [PATCH 38/45] more typos in the tests --- test/test_bounds.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 86fb664f86..0db183315b 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -39,10 +39,9 @@ constexpr Real T_split = 1e4; constexpr int N_per_decade_fine = 200; constexpr Real N_factor = 5; -constexpr Real REAL_TOL = - std::numeric_limits::epsilon * 1e3 +constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e3; - SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { +SCENARIO("Bounds can compute number of points from points per decade", "[Bounds]") { WHEN("We compute the number of points from points per decade") { int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); constexpr int NDECADES = 9; From 2c743198bca5a93977fe5b71f2a37d3b3af789a7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 10:16:47 -0600 Subject: [PATCH 39/45] typos --- test/test_bounds.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index 0db183315b..c07e82a606 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -45,7 +45,7 @@ SCENARIO("Bounds can compute number of points from points per decade", "[Bounds] WHEN("We compute the number of points from points per decade") { int np = Bounds::getNumPointsFromPPD(T_min, T_max, N_per_decade_fine); constexpr int NDECADES = 9; - THEN("We get the right number") { REQUIRE(np == NDECADES * N_PER_DECADE_FINE); } + THEN("We get the right number") { REQUIRE(np == NDECADES * N_per_decade_fine); } } } From 0dc46c4b4a26f55326b96259ca84e7214beba7a7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 11:24:13 -0600 Subject: [PATCH 40/45] update spiner so backwards compatibility with hierarchical grids retained --- utils/spiner | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/spiner b/utils/spiner index cc089ccb3a..0d0d7474f3 160000 --- a/utils/spiner +++ b/utils/spiner @@ -1 +1 @@ -Subproject commit cc089ccb3a1f57afdc21becaaab286217589cf7c +Subproject commit 0d0d7474f34cdc168c141185cc692b65ee8ad0c0 From bd485f9a7700ad9a4d8b5110d6f6a587259278ff Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 17:20:56 -0600 Subject: [PATCH 41/45] dholladay comments --- singularity-eos/base/table_bounds.hpp | 7 ++++--- singularity-eos/eos/eos_stellar_collapse.hpp | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/table_bounds.hpp index e6300a0300..80d3dbd080 100644 --- a/singularity-eos/base/table_bounds.hpp +++ b/singularity-eos/base/table_bounds.hpp @@ -57,7 +57,6 @@ class Bounds { if (convertToLog) { convertBoundsToLog_(min, max, shrinkRange); if (!(std::isnan(anchor_point))) { - anchor_point += offset; anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); } } @@ -157,7 +156,7 @@ class Bounds { // This uses real logs template static int getNumPointsFromPPD(Real min, Real max, const T ppd) { - constexpr Real epsilon = std::numeric_limits::epsilon(); + constexpr Real epsilon = std::numeric_limits::epsilon(); const Real min_offset = 10 * std::abs(epsilon); // 1.1 so that the y-intercept isn't identically zero // when min < 0. @@ -188,7 +187,7 @@ class Bounds { void convertBoundsToLog_(Real &min, Real &max, const Real shrinkRange = 0) { // Log scales can't handle negative numbers or exactly zero. To // deal with that, we offset. - constexpr Real epsilon = std::numeric_limits::epsilon(); + constexpr Real epsilon = std::numeric_limits::epsilon(); const Real min_offset = 10 * std::abs(epsilon); // 1.1 so that the y-intercept isn't identically zero // when min < 0. @@ -219,6 +218,8 @@ class Bounds { N = Nmax_new; max = max_new; + } else { + PORTABLE_ALWAYS_WARN("Anchor point out of bounds. Ignoring it."); } } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 1bbf1dfc9a..087f5f4bea 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -230,8 +230,8 @@ class StellarCollapse : public EosBase { } private: - const DataBox &field_; - const Real Ye_, lRho_; + DataBox &field_; + Real Ye_, lRho_; }; inline void LoadFromSP5File_(const std::string &filename); From 7b189118d2358a0d1eb4cc11b37e381980c7484d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 17:21:39 -0600 Subject: [PATCH 42/45] remove unnecessary line --- test/compare_to_eospac.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/compare_to_eospac.cpp b/test/compare_to_eospac.cpp index 4ff052c555..5a4bfeb51c 100644 --- a/test/compare_to_eospac.cpp +++ b/test/compare_to_eospac.cpp @@ -47,7 +47,6 @@ #include -// using namespace Spiner; using namespace singularity; using namespace EospacWrapper; From 6e36ffc5c98ffdfc0098d44c1eaeffe8925210fd Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 12 Aug 2024 17:25:57 -0600 Subject: [PATCH 43/45] change bounds file name --- sesame2spiner/io_eospac.hpp | 2 +- singularity-eos/CMakeLists.txt | 2 +- .../base/{table_bounds.hpp => spiner_table_bounds.hpp} | 0 test/test_bounds.cpp | 2 +- 4 files changed, 3 insertions(+), 3 deletions(-) rename singularity-eos/base/{table_bounds.hpp => spiner_table_bounds.hpp} (100%) diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index fc492e2472..dcbfce18eb 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include #include diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index f4c0ea2d63..0096902e09 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -27,7 +27,7 @@ register_headers( base/fast-math/logs.hpp base/robust_utils.hpp base/root-finding-1d/root_finding.hpp - base/table_bounds.hpp + base/spiner_table_bounds.hpp base/variadic_utils.hpp base/math_utils.hpp base/constants.hpp diff --git a/singularity-eos/base/table_bounds.hpp b/singularity-eos/base/spiner_table_bounds.hpp similarity index 100% rename from singularity-eos/base/table_bounds.hpp rename to singularity-eos/base/spiner_table_bounds.hpp diff --git a/test/test_bounds.cpp b/test/test_bounds.cpp index c07e82a606..34d2f7da40 100644 --- a/test/test_bounds.cpp +++ b/test/test_bounds.cpp @@ -17,7 +17,7 @@ #include #include -#include +#include #ifndef CATCH_CONFIG_FAST_COMPILE #define CATCH_CONFIG_FAST_COMPILE From 0260d0c00a3e3bf9df3a24018cbaa642e78cf82d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 13 Aug 2024 09:17:29 -0600 Subject: [PATCH 44/45] cascade of const correctness --- singularity-eos/eos/eos_stellar_collapse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 087f5f4bea..edcd60a600 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -223,7 +223,7 @@ class StellarCollapse : public EosBase { class LogT { public: PORTABLE_INLINE_FUNCTION - LogT(const DataBox &field, const Real Ye, const Real lRho) + LogT(DataBox &field, const Real Ye, const Real lRho) : field_(field), Ye_(Ye), lRho_(lRho) {} PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { return field_.interpToReal(Ye_, lT, lRho_); From b2db905d63de4f55f6347737e0309ec33d61d1c8 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 13 Aug 2024 09:26:25 -0600 Subject: [PATCH 45/45] these have to be const --- singularity-eos/eos/eos_stellar_collapse.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index edcd60a600..1bbf1dfc9a 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -223,15 +223,15 @@ class StellarCollapse : public EosBase { class LogT { public: PORTABLE_INLINE_FUNCTION - LogT(DataBox &field, const Real Ye, const Real lRho) + LogT(const DataBox &field, const Real Ye, const Real lRho) : field_(field), Ye_(Ye), lRho_(lRho) {} PORTABLE_INLINE_FUNCTION Real operator()(const Real lT) const { return field_.interpToReal(Ye_, lT, lRho_); } private: - DataBox &field_; - Real Ye_, lRho_; + const DataBox &field_; + const Real Ye_, lRho_; }; inline void LoadFromSP5File_(const std::string &filename);