From eeab2b70f139597ff4ec490ae79bcd3ea38b5b26 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 24 Sep 2022 12:18:04 +0200 Subject: [PATCH 01/10] [component][benzene] Implement some mandatory interfaces --- dumux/material/components/benzene.hh | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/dumux/material/components/benzene.hh b/dumux/material/components/benzene.hh index 01e2728bb7..4139d3c2ac 100644 --- a/dumux/material/components/benzene.hh +++ b/dumux/material/components/benzene.hh @@ -59,6 +59,30 @@ public: static constexpr Scalar molarMass() { return 0.07811; } + /*! + * \brief Returns true if the gas phase is assumed to be compressible + */ + static constexpr bool gasIsCompressible() + { return true; } + + /*! + * \brief Returns true if the liquid phase is assumed to be compressible + */ + static constexpr bool liquidIsCompressible() + { return false; } + + /*! + * \brief Returns true if the gas phase viscosity is constant + */ + static constexpr bool gasViscosityIsConstant() + { return true; } + + /*! + * \brief Returns true if the liquid phase viscosity is constant + */ + static constexpr bool liquidViscosityIsConstant() + { return true; } + /*! * \brief The density of benzene steam at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$. * -- GitLab From 18dea855a768bdb470d92d83767aebd2fdabfb4a Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 24 Sep 2022 13:27:02 +0200 Subject: [PATCH 02/10] [test][1p][compressible] Make sure water is liquid in the requested temperature range --- test/porousmediumflow/1p/compressible/instationary/problem.hh | 4 ++-- test/porousmediumflow/1p/compressible/stationary/problem.hh | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/porousmediumflow/1p/compressible/instationary/problem.hh b/test/porousmediumflow/1p/compressible/instationary/problem.hh index bb5f6c99c8..d904cc1b2e 100644 --- a/test/porousmediumflow/1p/compressible/instationary/problem.hh +++ b/test/porousmediumflow/1p/compressible/instationary/problem.hh @@ -55,8 +55,8 @@ public: OnePTestProblem(std::shared_ptr gridGeometry) : ParentType(gridGeometry) { - Components::TabulatedComponent>::init(272.15, 294.15, 10, - 1.0e4, 1.0e6, 200); + Components::TabulatedComponent> + ::init(273.15, 294.15, 10, 1.0e4, 1.0e6, 200); } /*! diff --git a/test/porousmediumflow/1p/compressible/stationary/problem.hh b/test/porousmediumflow/1p/compressible/stationary/problem.hh index 4560eba4d1..d0b8353493 100644 --- a/test/porousmediumflow/1p/compressible/stationary/problem.hh +++ b/test/porousmediumflow/1p/compressible/stationary/problem.hh @@ -58,7 +58,7 @@ public: : ParentType(gridGeometry) { Components::TabulatedComponent> - ::init(272.15, 294.15, 10, 1.0e4, 1.0e6, 200); + ::init(273.15, 294.15, 10, 1.0e4, 1.0e6, 200); } /*! -- GitLab From b28bc16f1a2e54f8fd9e830ccb587a9a0096d40c Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 24 Sep 2022 12:18:24 +0200 Subject: [PATCH 03/10] [doc][simpleh2o] Fix typo --- dumux/material/components/simpleh2o.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/material/components/simpleh2o.hh b/dumux/material/components/simpleh2o.hh index 418a340970..fe1bb60ba7 100644 --- a/dumux/material/components/simpleh2o.hh +++ b/dumux/material/components/simpleh2o.hh @@ -217,13 +217,13 @@ public: { return false; } /*! - * \brief Returns true if the gas phase viscostiy is constant + * \brief Returns true if the gas phase viscosity is constant */ static constexpr bool gasViscosityIsConstant() { return true; } /*! - * \brief Returns true if the liquid phase viscostiy is constant + * \brief Returns true if the liquid phase viscosity is constant */ static constexpr bool liquidViscosityIsConstant() { return true; } -- GitLab From bf8f58c1e1db4256ad96d7fcd92b54be1ca5ca26 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 17:15:50 +0200 Subject: [PATCH 04/10] [cleanup][components] Use nested namespaces --- dumux/material/components/tabulatedcomponent.hh | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh index ce347e8464..91c8eb5840 100644 --- a/dumux/material/components/tabulatedcomponent.hh +++ b/dumux/material/components/tabulatedcomponent.hh @@ -38,12 +38,13 @@ #include #include -namespace Dumux { -namespace Components { +namespace Dumux::Components { // forward declaration template class TabulatedComponent; -} // end namespace Components +} // end namespace Dumux::Components + +namespace Dumux { //! component traits for tabulated component template @@ -60,8 +61,9 @@ struct ComponentTraits, RawComponent>::value; }; +} // end namespace Dumux -namespace Components { +namespace Dumux::Components { /*! * \ingroup Components @@ -1088,8 +1090,6 @@ struct IsAqueous; template struct IsAqueous> : public IsAqueous {}; -} // end namespace Components - -} // end namespace Dumux +} // end namespace Dumux::Components #endif -- GitLab From a71ba31749f4dd7e562e482d9da2f6d2d8c2ae75 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 22 Sep 2022 01:53:00 +0200 Subject: [PATCH 05/10] [test][tabularization] Use multithreading and fix range for water In the given range we are solid and then the IAPWS is outside the applicable range. --- test/material/tabulation/test_tabulation.cc | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/test/material/tabulation/test_tabulation.cc b/test/material/tabulation/test_tabulation.cc index 8c12d7d36f..7265ab56ef 100644 --- a/test/material/tabulation/test_tabulation.cc +++ b/test/material/tabulation/test_tabulation.cc @@ -30,6 +30,8 @@ #include +#include + #include #include #include @@ -68,6 +70,9 @@ int main(int argc, char *argv[]) using namespace Dumux; using Scalar = double; + // initialize MPI+X (tabulation uses multithreading features) + Dumux::initialize(argc, argv); + // test IapwsH2O in detail { using IapwsH2O = Components::H2O; @@ -89,7 +94,7 @@ int main(int argc, char *argv[]) std::cout << "Checking tabulation\n"; const int m = nTemp*3; const int n = nPress*3; - for (int i = 0; i < m; ++i) + Dumux::parallelFor(m, [&](int i) { const Scalar T = tempMin + (tempMax - tempMin)*Scalar(i)/Scalar(m-1); checkEquality("vaporPressure", TabulatedH2O::vaporPressure(T), IapwsH2O::vaporPressure(T), eps); @@ -120,8 +125,7 @@ int main(int argc, char *argv[]) } - - } + }); } // test if other components can be tabulated @@ -135,7 +139,7 @@ int main(int argc, char *argv[]) Components::TabulatedComponent, false>::init(273, 275, 3, 1e5, 1e6, 3); Components::TabulatedComponent, false>::init(273, 275, 3, 1e5, 1e6, 3); Components::TabulatedComponent, false>::init(273, 275, 3, 1e5, 1e6, 3); - Components::TabulatedComponent>::init(273, 275, 3, 1e5, 1e6, 3); + Components::TabulatedComponent>::init(273.15, 275, 3, 1e5, 1e6, 3); Components::TabulatedComponent>::init(273, 275, 3, 1e5, 1e6, 3); Components::TabulatedComponent>::init(273, 275, 3, 1e5, 1e6, 3); Components::TabulatedComponent>::init(273, 275, 3, 1e5, 1e6, 3); -- GitLab From 358ef2b2936e8e2e8a52484d52f32f922cf41c71 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 17:51:58 +0200 Subject: [PATCH 06/10] [fix][component] Only use lazy tabularization in serial mode --- .../material/components/tabulatedcomponent.hh | 558 ++++++++++++++---- 1 file changed, 436 insertions(+), 122 deletions(-) diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh index 91c8eb5840..365af7685b 100644 --- a/dumux/material/components/tabulatedcomponent.hh +++ b/dumux/material/components/tabulatedcomponent.hh @@ -35,7 +35,11 @@ #include #include +#include + +#include #include +#include #include namespace Dumux::Components { @@ -63,12 +67,63 @@ struct ComponentTraits struct AlwaysFalse : public std::true_type {}; +}// end namespace Dumux + +namespace Dumux::Components::Detail { + +template using CompHasNoLiquidEnthalpy = decltype(C::template liquidEnthalpy(0.0, 0.0)); +template using CompHasNoLiquidDensity = decltype(C::template liquidDensity(0.0, 0.0)); +template using CompHasNoLiquidThermalCond = decltype(C::template liquidThermalConductivity(0.0, 0.0)); +template using CompHasNoLiquidHeatCapacity = decltype(C::template liquidHeatCapacity(0.0, 0.0)); +template using CompHasNoLiquidViscosity = decltype(C::template liquidViscosity(0.0, 0.0)); +template using CompHasLiquidPressure = decltype(C::liquidPressure(0.0, 0.0)); + +template using CompHasNoGasEnthalpy = decltype(C::template gasEnthalpy(0.0, 0.0)); +template using CompHasNoGasDensity = decltype(C::template gasDensity(0.0, 0.0)); +template using CompHasNoGasThermalCond = decltype(C::template gasThermalConductivity(0.0, 0.0)); +template using CompHasNoGasHeatCapacity = decltype(C::template gasHeatCapacity(0.0, 0.0)); +template using CompHasNoGasViscosity = decltype(C::template gasViscosity(0.0, 0.0)); +template using CompHasGasPressure = decltype(C::gasPressure(0.0, 0.0)); + +template constexpr inline bool hasLiquidEnthalpy() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } +template constexpr inline bool hasLiquidDensity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } +template constexpr inline bool hasLiquidThermalConductivity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } +template constexpr inline bool hasLiquidHeatCapacity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } +template constexpr inline bool hasLiquidViscosity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } +template constexpr inline bool hasLiquidPressure() +{ return Dune::Std::is_detected::value && ComponentTraits::hasLiquidState; } + +template constexpr inline bool hasGasEnthalpy() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasGasState; } +template constexpr inline bool hasGasDensity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasGasState; } +template constexpr inline bool hasGasThermalConductivity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasGasState; } +template constexpr inline bool hasGasHeatCapacity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasGasState; } +template constexpr inline bool hasGasViscosity() +{ return !Dune::Std::is_detected::value && ComponentTraits::hasGasState; } +template constexpr inline bool hasGasPressure() +{ return Dune::Std::is_detected::value && ComponentTraits::hasGasState; } + +} // end namespace Dumux::Components::Detail + namespace Dumux::Components { /*! * \ingroup Components - * \brief Tabulates all thermodynamic properties of a given - * untabulated chemical species. + * \brief Tabulates all thermodynamic properties of a given component * * At the moment, this class can only handle the sub-critical fluids * since it tabulates along the vapor pressure curve. @@ -82,6 +137,7 @@ namespace Dumux::Components { template class TabulatedComponent { + using ThisType = TabulatedComponent; public: //! export scalar type using Scalar = typename RawComponent::Scalar; @@ -108,7 +164,11 @@ public: pressMin_ = pressMin; pressMax_ = pressMax; nPress_ = nPress; - nDensity_ = nPress_; + nDensity_ = nPress; + +#ifndef NDEBUG + warningPrinted_ = false; +#endif std::cout << "-------------------------------------------------------------------------\n" << "Initializing tables for the " << RawComponent::name() @@ -123,49 +183,86 @@ public: assert(std::numeric_limits::has_quiet_NaN); const auto NaN = std::numeric_limits::quiet_NaN(); + // initialize vapor pressure array depending on useVaporPressure vaporPressure_.resize(nTemp_, NaN); - minGasDensity_.resize(nTemp_, NaN); - maxGasDensity_.resize(nTemp_, NaN); - minLiquidDensity_.resize(nTemp_, NaN); - maxLiquidDensity_.resize(nTemp_, NaN); - - const std::size_t numEntriesTp = nTemp_*nPress_; // = nTemp_*nDensity_ - gasEnthalpy_.resize(numEntriesTp, NaN); - liquidEnthalpy_.resize(numEntriesTp, NaN); - gasHeatCapacity_.resize(numEntriesTp, NaN); - liquidHeatCapacity_.resize(numEntriesTp, NaN); - gasDensity_.resize(numEntriesTp, NaN); - liquidDensity_.resize(numEntriesTp, NaN); - gasViscosity_.resize(numEntriesTp, NaN); - liquidViscosity_.resize(numEntriesTp, NaN); - gasThermalConductivity_.resize(numEntriesTp, NaN); - liquidThermalConductivity_.resize(numEntriesTp, NaN); - gasPressure_.resize(numEntriesTp, NaN); - liquidPressure_.resize(numEntriesTp, NaN); - - // reset all flags - minMaxLiquidDensityInitialized_ = false; - minMaxGasDensityInitialized_ = false; - gasEnthalpyInitialized_ = false; - liquidEnthalpyInitialized_ = false; - gasHeatCapacityInitialized_ = false; - liquidHeatCapacityInitialized_ = false; - gasDensityInitialized_ = false; - liquidDensityInitialized_ = false; - gasViscosityInitialized_ = false; - liquidViscosityInitialized_ = false; - gasThermalConductivityInitialized_ = false; - liquidThermalConductivityInitialized_ = false; - gasPressureInitialized_ = false; - liquidPressureInitialized_ = false; - - //! initialize vapor pressure array depending on useVaporPressure initVaporPressure_(); -#ifndef NDEBUG + if constexpr (ComponentTraits::hasGasState) + { + minGasDensity_.resize(nTemp_, NaN); + maxGasDensity_.resize(nTemp_, NaN); + const std::size_t numEntriesTp = nTemp_*nPress_; + gasEnthalpy_.resize(numEntriesTp, NaN); + gasHeatCapacity_.resize(numEntriesTp, NaN); + gasDensity_.resize(numEntriesTp, NaN); + gasViscosity_.resize(numEntriesTp, NaN); + gasThermalConductivity_.resize(numEntriesTp, NaN); + + if constexpr (RawComponent::gasIsCompressible()) + gasPressure_.resize(numEntriesTp, NaN); + + minMaxGasDensityInitialized_ = false; + gasEnthalpyInitialized_ = false; + gasHeatCapacityInitialized_ = false; + gasDensityInitialized_ = false; + gasViscosityInitialized_ = false; + gasThermalConductivityInitialized_ = false; + gasPressureInitialized_ = false; + } + + if constexpr (ComponentTraits::hasLiquidState) + { + minLiquidDensity_.resize(nTemp_, NaN); + maxLiquidDensity_.resize(nTemp_, NaN); + + const std::size_t numEntriesTp = nTemp_*nPress_; + liquidEnthalpy_.resize(numEntriesTp, NaN); + liquidHeatCapacity_.resize(numEntriesTp, NaN); + liquidDensity_.resize(numEntriesTp, NaN); + liquidViscosity_.resize(numEntriesTp, NaN); + liquidThermalConductivity_.resize(numEntriesTp, NaN); + + if constexpr (RawComponent::liquidIsCompressible()) + liquidPressure_.resize(numEntriesTp, NaN); + + // reset all flags + minMaxLiquidDensityInitialized_ = false; + liquidEnthalpyInitialized_ = false; + liquidHeatCapacityInitialized_ = false; + liquidDensityInitialized_ = false; + liquidViscosityInitialized_ = false; + liquidThermalConductivityInitialized_ = false; + liquidPressureInitialized_ = false; + } + + // in the multithreaded case precompute all array because lazy-loading + // leads to potential data races + if constexpr (!Multithreading::isSerial()) + { + if constexpr (ComponentTraits::hasGasState) + { + minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_(); + gasPressureInitialized_ = tabularizeGasPressure_(); + gasEnthalpyInitialized_ = tabularizeGasEnthalpy_(); + gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_(); + gasDensityInitialized_ = tabularizeGasDensity_(); + gasViscosityInitialized_ = tabularizeGasViscosity_(); + gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_(); + } + + if constexpr (ComponentTraits::hasLiquidState) + { + minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_(); + liquidPressureInitialized_ = tabularizeLiquidPressure_(); + liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_(); + liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_(); + liquidDensityInitialized_ = tabularizeLiquidDensity_(); + liquidViscosityInitialized_ = tabularizeLiquidViscosity_(); + liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_(); + } + } + initialized_ = true; - warningPrinted_ = false; -#endif } /*! @@ -244,12 +341,15 @@ public: using std::isnan; if (isnan(result)) { - if (!gasEnthalpyInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasEnth = [] (auto T, auto p) { return RawComponent::gasEnthalpy(T, p); }; - initTPArray_(gasEnth, minGasPressure_, maxGasPressure_, gasEnthalpy_); - gasEnthalpyInitialized_ = true; - return gasEnthalpy(temperature, pressure); + // lazy tabularization + if (!gasEnthalpyInitialized_) + { + gasEnthalpyInitialized_ = tabularizeGasEnthalpy_(); + if (gasEnthalpyInitialized_) + return gasEnthalpy(temperature, pressure); + } } printWarning_("gasEnthalpy", temperature, pressure); @@ -271,12 +371,15 @@ public: using std::isnan; if (isnan(result)) { - if (!liquidEnthalpyInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqEnth = [] (auto T, auto p) { return RawComponent::liquidEnthalpy(T, p); }; - initTPArray_(liqEnth, minLiquidPressure_, maxLiquidPressure_, liquidEnthalpy_); - liquidEnthalpyInitialized_ = true; - return liquidEnthalpy(temperature, pressure); + // lazy tabularization + if (!liquidEnthalpyInitialized_) + { + liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_(); + if (liquidEnthalpyInitialized_) + return liquidEnthalpy(temperature, pressure); + } } printWarning_("liquidEnthalpy", temperature, pressure); @@ -298,12 +401,15 @@ public: using std::isnan; if (isnan(result)) { - if (!gasHeatCapacityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasHC = [] (auto T, auto p) { return RawComponent::gasHeatCapacity(T, p); }; - initTPArray_(gasHC, minGasPressure_, maxGasPressure_, gasHeatCapacity_); - gasHeatCapacityInitialized_ = true; - return gasHeatCapacity(temperature, pressure); + // lazy tabularization + if (!gasHeatCapacityInitialized_) + { + gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_(); + if (gasHeatCapacityInitialized_) + return gasHeatCapacity(temperature, pressure); + } } printWarning_("gasHeatCapacity", temperature, pressure); @@ -325,12 +431,15 @@ public: using std::isnan; if (isnan(result)) { - if (!liquidHeatCapacityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqHC = [] (auto T, auto p) { return RawComponent::liquidHeatCapacity(T, p); }; - initTPArray_(liqHC, minLiquidPressure_, maxLiquidPressure_, liquidHeatCapacity_); - liquidHeatCapacityInitialized_ = true; - return liquidHeatCapacity(temperature, pressure); + // lazy tabularization + if (!liquidHeatCapacityInitialized_) + { + liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_(); + if (liquidHeatCapacityInitialized_) + return liquidHeatCapacity(temperature, pressure); + } } printWarning_("liquidHeatCapacity", temperature, pressure); @@ -369,24 +478,24 @@ public: */ static Scalar gasPressure(Scalar temperature, Scalar density) { - //! make sure the minimum/maximum densities have been computed - if (!minMaxGasDensityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasRho = [] (auto T, auto p) { return RawComponent::gasDensity(T, p); }; - initMinMaxRhoArray_(gasRho, minGasPressure_, maxGasPressure_, minGasDensity_, maxGasDensity_); - minMaxGasDensityInitialized_ = true; + if (!minMaxGasDensityInitialized_) + minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_(); } Scalar result = interpolateTRho_(gasPressure_, temperature, density, densityGasIdx_); using std::isnan; if (isnan(result)) { - if (!gasPressureInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasPFunc = [] (auto T, auto rho) { return RawComponent::gasPressure(T, rho); }; - initPressureArray_(gasPressure_, gasPFunc, minGasDensity_, maxGasDensity_); - gasPressureInitialized_ = true; - return gasPressure(temperature, density); + if (!gasPressureInitialized_) + { + gasPressureInitialized_ = tabularizeGasPressure_(); + if (gasPressureInitialized_) + return gasPressure(temperature, density); + } } printWarning_("gasPressure", temperature, density); @@ -403,24 +512,25 @@ public: */ static Scalar liquidPressure(Scalar temperature, Scalar density) { - //! make sure the minimum/maximum densities have been computed - if (!minMaxLiquidDensityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqRho = [] (auto T, auto p) { return RawComponent::liquidDensity(T, p); }; - initMinMaxRhoArray_(liqRho, minLiquidPressure_, maxLiquidPressure_, minLiquidDensity_, maxLiquidDensity_); - minMaxLiquidDensityInitialized_ = true; + if (!minMaxLiquidDensityInitialized_) + minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_(); } Scalar result = interpolateTRho_(liquidPressure_, temperature, density, densityLiquidIdx_); using std::isnan; if (isnan(result)) { - if (!liquidPressureInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqPFunc = [] (auto T, auto rho) { return RawComponent::liquidPressure(T, rho); }; - initPressureArray_(liquidPressure_, liqPFunc, minLiquidDensity_, maxLiquidDensity_); - liquidPressureInitialized_ = true; - return liquidPressure(temperature, density); + // lazy tabularization + if (!liquidPressureInitialized_) + { + liquidPressureInitialized_ = tabularizeLiquidPressure_(); + if (liquidPressureInitialized_) + return liquidPressure(temperature, density); + } } printWarning_("liquidPressure", temperature, density); @@ -462,12 +572,15 @@ public: using std::isnan; if (isnan(result)) { - if (!gasDensityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasRho = [] (auto T, auto p) { return RawComponent::gasDensity(T, p); }; - initTPArray_(gasRho, minGasPressure_, maxGasPressure_, gasDensity_); - gasDensityInitialized_ = true; - return gasDensity(temperature, pressure); + // lazy tabularization + if (!gasDensityInitialized_) + { + gasDensityInitialized_ = tabularizeGasDensity_(); + if (gasDensityInitialized_) + return gasDensity(temperature, pressure); + } } printWarning_("gasDensity", temperature, pressure); @@ -501,16 +614,15 @@ public: using std::isnan; if (isnan(result)) { - if (!liquidDensityInitialized_) + if constexpr (Multithreading::isSerial()) { - // TODO: we could get rid of the lambdas and pass the functor irectly. But, - // currently Brine is a component (and not a fluid system) expecting a - // third argument with a default, which cannot be wrapped in a function pointer. - // For this reason we have to wrap this into a lambda here. - auto liqRho = [] (auto T, auto p) { return RawComponent::liquidDensity(T, p); }; - initTPArray_(liqRho, minLiquidPressure_, maxLiquidPressure_, liquidDensity_); - liquidDensityInitialized_ = true; - return liquidDensity(temperature, pressure); + // lazy tabularization + if (!liquidDensityInitialized_) + { + liquidDensityInitialized_ = tabularizeLiquidDensity_(); + if (liquidDensityInitialized_) + return liquidDensity(temperature, pressure); + } } printWarning_("liquidDensity", temperature, pressure); @@ -544,12 +656,15 @@ public: using std::isnan; if (isnan(result)) { - if (!gasViscosityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasVisc = [] (auto T, auto p) { return RawComponent::gasViscosity(T, p); }; - initTPArray_(gasVisc, minGasPressure_, maxGasPressure_, gasViscosity_); - gasViscosityInitialized_ = true; - return gasViscosity(temperature, pressure); + // lazy tabularization + if (!gasViscosityInitialized_) + { + gasViscosityInitialized_ = tabularizeGasViscosity_(); + if (gasViscosityInitialized_) + return gasViscosity(temperature, pressure); + } } printWarning_("gasViscosity", temperature, pressure); @@ -571,12 +686,15 @@ public: using std::isnan; if (isnan(result)) { - if (!liquidViscosityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqVisc = [] (auto T, auto p) { return RawComponent::liquidViscosity(T, p); }; - initTPArray_(liqVisc, minLiquidPressure_, maxLiquidPressure_, liquidViscosity_); - liquidViscosityInitialized_ = true; - return liquidViscosity(temperature, pressure); + // lazy tabularization + if (!liquidViscosityInitialized_) + { + liquidViscosityInitialized_ = tabularizeLiquidViscosity_(); + if (liquidViscosityInitialized_) + return liquidViscosity(temperature, pressure); + } } printWarning_("liquidViscosity",temperature, pressure); @@ -598,12 +716,15 @@ public: using std::isnan; if (isnan(result)) { - if (!gasThermalConductivityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto gasTC = [] (auto T, auto p) { return RawComponent::gasThermalConductivity(T, p); }; - initTPArray_(gasTC, minGasPressure_, maxGasPressure_, gasThermalConductivity_); - gasThermalConductivityInitialized_ = true; - return gasThermalConductivity(temperature, pressure); + // lazy tabularization + if (!gasThermalConductivityInitialized_) + { + gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_(); + if (gasThermalConductivityInitialized_) + return gasThermalConductivity(temperature, pressure); + } } printWarning_("gasThermalConductivity", temperature, pressure); @@ -625,12 +746,15 @@ public: using std::isnan; if (isnan(result)) { - if (!liquidThermalConductivityInitialized_) + if constexpr (Multithreading::isSerial()) { - auto liqTC = [] (auto T, auto p) { return RawComponent::liquidThermalConductivity(T, p); }; - initTPArray_(liqTC, minLiquidPressure_, maxLiquidPressure_, liquidThermalConductivity_); - liquidThermalConductivityInitialized_ = true; - return liquidThermalConductivity(temperature, pressure); + // lazy tabularization + if (!liquidThermalConductivityInitialized_) + { + liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_(); + if (liquidThermalConductivityInitialized_) + return liquidThermalConductivity(temperature, pressure); + } } printWarning_("liquidThermalConductivity", temperature, pressure); @@ -888,7 +1012,7 @@ private: return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin); } - //! returns the minimum tabulized liquid pressure at a given temperature index + //! returns the minimum tabularized liquid pressure at a given temperature index static Scalar minLiquidPressure_(int tempIdx) { using std::max; @@ -898,7 +1022,7 @@ private: return max(pressMin_, vaporPressure_[tempIdx] / 1.1); } - //! returns the maximum tabulized liquid pressure at a given temperature index + //! returns the maximum tabularized liquid pressure at a given temperature index static Scalar maxLiquidPressure_(int tempIdx) { using std::max; @@ -908,7 +1032,7 @@ private: return max(pressMax_, vaporPressure_[tempIdx] * 1.1); } - //! returns the minimum tabulized gas pressure at a given temperature index + //! returns the minimum tabularized gas pressure at a given temperature index static Scalar minGasPressure_(int tempIdx) { using std::min; @@ -918,7 +1042,7 @@ private: return min(pressMin_, vaporPressure_[tempIdx] / 1.1 ); } - //! returns the maximum tabulized gas pressure at a given temperature index + //! returns the maximum tabularized gas pressure at a given temperature index static Scalar maxGasPressure_(int tempIdx) { using std::min; @@ -928,9 +1052,198 @@ private: return min(pressMax_, vaporPressure_[tempIdx] * 1.1); } -#ifndef NDEBUG + template + static bool tabularizeGasEnthalpy_() + { + if constexpr (Detail::hasGasEnthalpy()) + { + auto gasEnth = [] (auto T, auto p) { return RC::gasEnthalpy(T, p); }; + initTPArray_(gasEnth, minGasPressure_, maxGasPressure_, gasEnthalpy_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidEnthalpy_() + { + if constexpr (Detail::hasLiquidEnthalpy()) + { + auto liqEnth = [] (auto T, auto p) { return RC::liquidEnthalpy(T, p); }; + initTPArray_(liqEnth, minLiquidPressure_, maxLiquidPressure_, liquidEnthalpy_); + return true; + } + + return false; + } + + template + static bool tabularizeGasHeatCapacity_() + { + if constexpr (Detail::hasGasHeatCapacity()) + { + auto gasHC = [] (auto T, auto p) { return RC::gasHeatCapacity(T, p); }; + initTPArray_(gasHC, minGasPressure_, maxGasPressure_, gasHeatCapacity_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidHeatCapacity_() + { + if constexpr (Detail::hasLiquidHeatCapacity()) + { + auto liqHC = [] (auto T, auto p) { return RC::liquidHeatCapacity(T, p); }; + initTPArray_(liqHC, minLiquidPressure_, maxLiquidPressure_, liquidHeatCapacity_); + return true; + } + + return false; + } + + template + static bool tabularizeMinMaxGasDensity_() + { + if constexpr (Detail::hasGasDensity()) + { + auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); }; + initMinMaxRhoArray_(gasRho, minGasPressure_, maxGasPressure_, minGasDensity_, maxGasDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeMinMaxLiquidDensity_() + { + if constexpr (Detail::hasGasEnthalpy()) + { + auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); }; + initMinMaxRhoArray_(liqRho, minLiquidPressure_, maxLiquidPressure_, minLiquidDensity_, maxLiquidDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeGasPressure_() + { + // pressure is only defined if the gas is compressible (this is usually the case) + if constexpr (Detail::hasGasPressure() && RC::gasIsCompressible()) + { + auto gasPFunc = [] (auto T, auto rho) { return RC::gasPressure(T, rho); }; + initPressureArray_(gasPressure_, gasPFunc, minGasDensity_, maxGasDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidPressure_() + { + // pressure is only defined if the liquid is compressible (this is often not the case) + if constexpr (Detail::hasLiquidPressure() && RC::liquidIsCompressible()) + { + auto liqPFunc = [] (auto T, auto rho) { return RC::liquidPressure(T, rho); }; + initPressureArray_(liquidPressure_, liqPFunc, minLiquidDensity_, maxLiquidDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeGasDensity_() + { + if constexpr (Detail::hasGasDensity()) + { + auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); }; + initTPArray_(gasRho, minGasPressure_, maxGasPressure_, gasDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidDensity_() + { + if constexpr (Detail::hasLiquidDensity()) + { + // TODO: we could get rid of the lambdas and pass the functor directly. But, + // currently Brine is a component (and not a fluid system) expecting a + // third argument with a default, which cannot be wrapped in a function pointer. + // For this reason we have to wrap this into a lambda here. + auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); }; + initTPArray_(liqRho, minLiquidPressure_, maxLiquidPressure_, liquidDensity_); + return true; + } + + return false; + } + + template + static bool tabularizeGasViscosity_() + { + if constexpr (Detail::hasGasViscosity()) + { + auto gasVisc = [] (auto T, auto p) { return RC::gasViscosity(T, p); }; + initTPArray_(gasVisc, minGasPressure_, maxGasPressure_, gasViscosity_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidViscosity_() + { + if constexpr (Detail::hasLiquidViscosity()) + { + auto liqVisc = [] (auto T, auto p) { return RC::liquidViscosity(T, p); }; + initTPArray_(liqVisc, minLiquidPressure_, maxLiquidPressure_, liquidViscosity_); + return true; + } + + return false; + } + + template + static bool tabularizeGasThermalConductivity_() + { + if constexpr (Detail::hasGasThermalConductivity()) + { + auto gasTC = [] (auto T, auto p) { return RC::gasThermalConductivity(T, p); }; + initTPArray_(gasTC, minGasPressure_, maxGasPressure_, gasThermalConductivity_); + return true; + } + + return false; + } + + template + static bool tabularizeLiquidThermalConductivity_() + { + if constexpr (Detail::hasLiquidThermalConductivity()) + { + auto liqTC = [] (auto T, auto p) { return RC::liquidThermalConductivity(T, p); }; + initTPArray_(liqTC, minLiquidPressure_, maxLiquidPressure_, liquidThermalConductivity_); + return true; + } + + return false; + } + // specifies whether the table was initialized static bool initialized_; + +#ifndef NDEBUG // specifies whether some warning was printed static bool warningPrinted_; #endif @@ -992,10 +1305,11 @@ private: static unsigned nDensity_; }; -#ifndef NDEBUG + template bool TabulatedComponent::initialized_ = false; +#ifndef NDEBUG template bool TabulatedComponent::warningPrinted_ = false; #endif -- GitLab From 40cba4e9cf71be93b990258fd60723144cdb2898 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 17:54:46 +0200 Subject: [PATCH 07/10] [component] Exploit parallelism in tabularization --- dumux/material/components/tabulatedcomponent.hh | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh index 365af7685b..0e383459f9 100644 --- a/dumux/material/components/tabulatedcomponent.hh +++ b/dumux/material/components/tabulatedcomponent.hh @@ -40,6 +40,7 @@ #include #include #include +#include #include namespace Dumux::Components { @@ -790,11 +791,11 @@ private: static void initVaporPressure_() { // fill the temperature-pressure arrays - for (unsigned iT = 0; iT < nTemp_; ++ iT) + Dumux::parallelFor(nTemp_, [&](unsigned iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; vaporPressure_[iT] = RawComponent::vaporPressure(temperature); - } + }); } //! if !useVaporPressure, do nothing here @@ -818,7 +819,7 @@ private: template static void initTPArray_(PropFunc&& f, MinPFunc&& minP, MaxPFunc&& maxP, std::vector& values) { - for (unsigned iT = 0; iT < nTemp_; ++ iT) + Dumux::parallelFor(nTemp_, [&](unsigned iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; @@ -829,7 +830,7 @@ private: Scalar pressure = iP * (pMax - pMin)/(nPress_ - 1) + pMin; values[iT + iP*nTemp_] = f(temperature, pressure); } - } + }); } /*! @@ -854,7 +855,7 @@ private: std::vector& rhoMin, std::vector& rhoMax) { - for (unsigned iT = 0; iT < nTemp_; ++ iT) + Dumux::parallelFor(nTemp_, [&](unsigned iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; @@ -863,7 +864,7 @@ private: rhoMax[iT] = rho(temperature, maxP(iT + 1)); else rhoMax[iT] = rho(temperature, maxP(iT)); - } + }); } /*! @@ -881,7 +882,7 @@ private: const std::vector& rhoMin, const std::vector& rhoMax) { - for (unsigned iT = 0; iT < nTemp_; ++ iT) + Dumux::parallelFor(nTemp_, [&](unsigned iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; @@ -892,7 +893,7 @@ private: + rhoMin[iT]; pressure[iT + iRho*nTemp_] = p(temperature, density); } - } + }); } //! returns an interpolated value depending on temperature -- GitLab From e629211189f467dda95ec3f742fb5da5c2a73758 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 18:04:48 +0200 Subject: [PATCH 08/10] [changelog] Add entry about tabluarized component bugfix --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2499c7d10a..b63ae12313 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ Differences Between DuMux 3.6 and DuMux 3.5 ### Improvements and Enhancements +- __Components__: Fixed a bug in `TabularizedComponent` that caused data races in multithreaded applications + - __Discretization__: There is now defaults for `FluxVariablesCache` and `FluxVariablesCacheFiller` for box/cctpfa/ccmpfa/staggered that implement an empty cache (i.e. nothing is cached for the flux variables). -- GitLab From 537e9829e6e2a2afb408f3372f0a13a3e4b199c0 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 19:00:57 +0200 Subject: [PATCH 09/10] [h2on2][cleanup] Use same alias --- dumux/material/fluidsystems/h2on2.hh | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/dumux/material/fluidsystems/h2on2.hh b/dumux/material/fluidsystems/h2on2.hh index e75b3ffe12..ccbae655f0 100644 --- a/dumux/material/fluidsystems/h2on2.hh +++ b/dumux/material/fluidsystems/h2on2.hh @@ -347,11 +347,8 @@ public: std::cout << " - use N2 heat conductivity as gas mixture heat conductivity: " << std::boolalpha << Policy::useN2HeatConductivityAsGasMixtureHeatConductivity() << "\n"; std::cout << " - use ideal gas heat capacities: " << std::boolalpha << Policy::useIdealGasHeatCapacities() << std::endl; - if (H2O::isTabulated) - { - TabulatedH2O::init(tempMin, tempMax, nTemp, - pressMin, pressMax, nPress); - } + if constexpr (H2O::isTabulated) + H2O::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress); } using Base::density; -- GitLab From 6edc15f63dec9a4ec804cc271b15e6425de3f5f8 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 21 Sep 2022 19:34:29 +0200 Subject: [PATCH 10/10] [component][tabulate] Use explicit C++ style type casts --- .../material/components/tabulatedcomponent.hh | 53 +++++++++---------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh index 0e383459f9..742ef7a669 100644 --- a/dumux/material/components/tabulatedcomponent.hh +++ b/dumux/material/components/tabulatedcomponent.hh @@ -34,6 +34,7 @@ #include #include #include +#include #include @@ -791,7 +792,7 @@ private: static void initVaporPressure_() { // fill the temperature-pressure arrays - Dumux::parallelFor(nTemp_, [&](unsigned iT) + Dumux::parallelFor(nTemp_, [&](std::size_t iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; vaporPressure_[iT] = RawComponent::vaporPressure(temperature); @@ -819,13 +820,13 @@ private: template static void initTPArray_(PropFunc&& f, MinPFunc&& minP, MaxPFunc&& maxP, std::vector& values) { - Dumux::parallelFor(nTemp_, [&](unsigned iT) + Dumux::parallelFor(nTemp_, [&](std::size_t iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; Scalar pMax = maxP(iT); Scalar pMin = minP(iT); - for (unsigned iP = 0; iP < nPress_; ++ iP) + for (std::size_t iP = 0; iP < nPress_; ++ iP) { Scalar pressure = iP * (pMax - pMin)/(nPress_ - 1) + pMin; values[iT + iP*nTemp_] = f(temperature, pressure); @@ -855,7 +856,7 @@ private: std::vector& rhoMin, std::vector& rhoMax) { - Dumux::parallelFor(nTemp_, [&](unsigned iT) + Dumux::parallelFor(nTemp_, [&](std::size_t iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; @@ -882,11 +883,11 @@ private: const std::vector& rhoMin, const std::vector& rhoMax) { - Dumux::parallelFor(nTemp_, [&](unsigned iT) + Dumux::parallelFor(nTemp_, [&](std::size_t iT) { Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_; - for (unsigned iRho = 0; iRho < nDensity_; ++ iRho) + for (std::size_t iRho = 0; iRho < nDensity_; ++ iRho) { Scalar density = Scalar(iRho)/(nDensity_ - 1) * (rhoMax[iT] - rhoMin[iT]) @@ -903,7 +904,7 @@ private: if (alphaT < 0 || alphaT >= nTemp_ - 1) return std::numeric_limits::quiet_NaN(); - unsigned iT = (unsigned) alphaT; + const auto iT = static_cast(alphaT); alphaT -= iT; return values[iT ]*(1 - alphaT) + @@ -919,16 +920,15 @@ private: if (alphaT < 0 || alphaT >= nTemp_ - 1) { return std::numeric_limits::quiet_NaN(); } - using std::min; - using std::max; - unsigned iT = max(0, min(nTemp_ - 2, (int) alphaT)); + using std::clamp; + const auto iT = clamp(static_cast(alphaT), 0, nTemp_ - 2); alphaT -= iT; Scalar alphaP1 = getPIdx(p, iT); Scalar alphaP2 = getPIdx(p, iT + 1); - unsigned iP1 = max(0, min(nPress_ - 2, (int) alphaP1)); - unsigned iP2 = max(0, min(nPress_ - 2, (int) alphaP2)); + const auto iP1 = clamp(static_cast(alphaP1), 0, nPress_ - 2); + const auto iP2 = clamp(static_cast(alphaP2), 0, nPress_ - 2); alphaP1 -= iP1; alphaP2 -= iP2; @@ -954,16 +954,15 @@ private: template static Scalar interpolateTRho_(const std::vector& values, Scalar T, Scalar rho, GetRhoIdx&& rhoIdx) { - using std::min; - using std::max; + using std::clamp; Scalar alphaT = tempIdx_(T); - unsigned iT = max(0, min(nTemp_ - 2, (int) alphaT)); + const auto iT = clamp(static_cast(alphaT), 0, nTemp_ - 2); alphaT -= iT; Scalar alphaP1 = rhoIdx(rho, iT); Scalar alphaP2 = rhoIdx(rho, iT + 1); - unsigned iP1 = max(0, min(nDensity_ - 2, (int) alphaP1)); - unsigned iP2 = max(0, min(nDensity_ - 2, (int) alphaP2)); + const auto iP1 = clamp(static_cast(alphaP1), 0, nDensity_ - 2); + const auto iP2 = clamp(static_cast(alphaP2), 0, nDensity_ - 2); alphaP1 -= iP1; alphaP2 -= iP2; @@ -980,7 +979,7 @@ private: } //! returns the index of an entry in a pressure field - static Scalar pressLiquidIdx_(Scalar pressure, unsigned tempIdx) + static Scalar pressLiquidIdx_(Scalar pressure, std::size_t tempIdx) { Scalar plMin = minLiquidPressure_(tempIdx); Scalar plMax = maxLiquidPressure_(tempIdx); @@ -989,7 +988,7 @@ private: } //! returns the index of an entry in a temperature field - static Scalar pressGasIdx_(Scalar pressure, unsigned tempIdx) + static Scalar pressGasIdx_(Scalar pressure, std::size_t tempIdx) { Scalar pgMin = minGasPressure_(tempIdx); Scalar pgMax = maxGasPressure_(tempIdx); @@ -998,7 +997,7 @@ private: } //! returns the index of an entry in a density field - static Scalar densityLiquidIdx_(Scalar density, unsigned tempIdx) + static Scalar densityLiquidIdx_(Scalar density, std::size_t tempIdx) { Scalar densityMin = minLiquidDensity_[tempIdx]; Scalar densityMax = maxLiquidDensity_[tempIdx]; @@ -1006,7 +1005,7 @@ private: } //! returns the index of an entry in a density field - static Scalar densityGasIdx_(Scalar density, unsigned tempIdx) + static Scalar densityGasIdx_(Scalar density, std::size_t tempIdx) { Scalar densityMin = minGasDensity_[tempIdx]; Scalar densityMax = maxGasDensity_[tempIdx]; @@ -1295,15 +1294,15 @@ private: // temperature, pressure and density ranges static Scalar tempMin_; static Scalar tempMax_; - static unsigned nTemp_; + static std::size_t nTemp_; static Scalar pressMin_; static Scalar pressMax_; - static unsigned nPress_; + static std::size_t nPress_; static Scalar densityMin_; static Scalar densityMax_; - static unsigned nDensity_; + static std::size_t nDensity_; }; @@ -1383,19 +1382,19 @@ typename RawComponent::Scalar TabulatedComponent template typename RawComponent::Scalar TabulatedComponent::tempMax_; template -unsigned TabulatedComponent::nTemp_; +std::size_t TabulatedComponent::nTemp_; template typename RawComponent::Scalar TabulatedComponent::pressMin_; template typename RawComponent::Scalar TabulatedComponent::pressMax_; template -unsigned TabulatedComponent::nPress_; +std::size_t TabulatedComponent::nPress_; template typename RawComponent::Scalar TabulatedComponent::densityMin_; template typename RawComponent::Scalar TabulatedComponent::densityMax_; template -unsigned TabulatedComponent::nDensity_; +std::size_t TabulatedComponent::nDensity_; // forward declaration template -- GitLab