From a01f63fdb182b04180f700d083f6e1389bd251c3 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@math.uio.no>
Date: Fri, 21 Oct 2022 21:41:44 +0000
Subject: [PATCH] Merge branch 'fix/tabulated-component-multithreading' into
'master'
Refactor tabulated component for multithreading
Closes #1189
See merge request dumux-repositories/dumux!3321
(cherry picked from commit f7a90b13cf4351c41811fa955f33b9e107d9394c)
8ace8f1d [component] Pass const l-value references to parallelFor
2b34ccec [tabulatedcomp] Cleanup
b03f466d [tabulatedcomp] Use const& instead of && to receive member function
fd7767d3 [tabcomp] Unify checks
947309a4 [tabulatedcomponent] Capture by value in lambda functions
ae5b06cb [tabulatedcomponent] avoid explicit if/else
b6f92d58 [tabulatedcomponent] use epsilon for float comparison
3c910454 [tabulatedcomponent] Use singleton for data
b03fa10f [tabulatedcomponent][bugfix] Clamp index to allowed range
d7b8d6ac [h2o] Use Brent's method for fluid pressure
---
dumux/material/components/h2o.hh | 40 +-
.../material/components/tabulatedcomponent.hh | 1860 ++++++++---------
2 files changed, 883 insertions(+), 1017 deletions(-)
diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh
index 3d0b6cf21c..2325b05afc 100644
--- a/dumux/material/components/h2o.hh
+++ b/dumux/material/components/h2o.hh
@@ -27,6 +27,7 @@
#include <cmath>
#include <cassert>
+#include <dumux/nonlinear/findscalarroot.hh>
#include <dumux/material/idealgas.hh>
#include <dumux/common/exceptions.hh>
@@ -677,29 +678,24 @@ public:
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{
- // We use the newton method for this. For the initial value we
- // assume the pressure to be 10% higher than the vapor
- // pressure
- Scalar pressure = 1.1*vaporPressure(temperature);
- Scalar eps = pressure*1e-7;
-
- Scalar deltaP = pressure*2;
-
- using std::abs;
- for (int i = 0; i < 5 && abs(pressure*1e-9) < abs(deltaP); ++i) {
- Scalar f = liquidDensity(temperature, pressure) - density;
-
- Scalar df_dp;
- df_dp = liquidDensity(temperature, pressure + eps);
- df_dp -= liquidDensity(temperature, pressure - eps);
- df_dp /= 2*eps;
-
- deltaP = - f/df_dp;
-
- pressure += deltaP;
+ // We use Brent's method for this, with a pressure range including the surroundings of the
+ // vapor pressure line
+ Scalar minPressure = vaporPressure(temperature)/1.11;
+ Scalar maxPressure = 100e6;
+ const auto residualFunction = [&] (const Scalar pressure) {
+ return liquidDensity(temperature, pressure) - density;
+ };
+ try
+ {
+ return findScalarRootBrent(minPressure, maxPressure, residualFunction);
+ }
+ catch (const NumericalProblem& e)
+ {
+ DUNE_THROW(NumericalProblem,
+ "searched for pressure(T=" << temperature << ",rho=" << density
+ <<") in [" << minPressure << ", " << maxPressure << "]: "
+ << e.what());
}
-
- return pressure;
}
/*!
diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh
index 742ef7a669..f69d4cc54e 100644
--- a/dumux/material/components/tabulatedcomponent.hh
+++ b/dumux/material/components/tabulatedcomponent.hh
@@ -119,46 +119,28 @@ template<class C> constexpr inline bool hasGasViscosity()
template<class C> constexpr inline bool hasGasPressure()
{ return Dune::Std::is_detected<CompHasGasPressure, C>::value && ComponentTraits<C>::hasGasState; }
-} // end namespace Dumux::Components::Detail
-
-namespace Dumux::Components {
-
-/*!
- * \ingroup Components
- * \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.
- *
- * \tparam Scalar The type used for scalar values
- * \tparam RawComponent The component which ought to be tabulated
- * \tparam useVaporPressure If set to true, the min/max pressure
- * values for gas&liquid phase will be set
- * depending on the vapor pressure.
- */
-template <class RawComponent, bool useVaporPressure=true>
-class TabulatedComponent
+template<class RawComponent, bool useVaporPressure = true>
+class TabulatedComponentTable
{
- using ThisType = TabulatedComponent<RawComponent, useVaporPressure>;
-public:
- //! export scalar type
using Scalar = typename RawComponent::Scalar;
+ friend class TabulatedComponent<RawComponent, useVaporPressure>;
- //! state that we are tabulated
- static constexpr bool isTabulated = true;
+ struct GasPolicy
+ {
+ Scalar minP(std::size_t iT) const { return table.minGasPressure(iT); }
+ Scalar maxP(std::size_t iT) const { return table.maxGasPressure(iT); }
+ const TabulatedComponentTable<RawComponent, useVaporPressure>& table;
+ };
- /*!
- * \brief Initialize the tables.
- *
- * \param tempMin The minimum of the temperature range in \f$\mathrm{[K]}\f$
- * \param tempMax The maximum of the temperature range in \f$\mathrm{[K]}\f$
- * \param nTemp The number of entries/steps within the temperature range
- * \param pressMin The minimum of the pressure range in \f$\mathrm{[Pa]}\f$
- * \param pressMax The maximum of the pressure range in \f$\mathrm{[Pa]}\f$
- * \param nPress The number of entries/steps within the pressure range
- */
- static void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
- Scalar pressMin, Scalar pressMax, std::size_t nPress)
+ struct LiquidPolicy
+ {
+ Scalar minP(std::size_t iT) const { return table.minLiquidPressure(iT); }
+ Scalar maxP(std::size_t iT) const { return table.maxLiquidPressure(iT); }
+ const TabulatedComponentTable<RawComponent, useVaporPressure>& table;
+ };
+public:
+ void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
+ Scalar pressMin, Scalar pressMax, std::size_t nPress)
{
tempMin_ = tempMin;
tempMax_ = tempMax;
@@ -168,28 +150,15 @@ public:
nPress_ = nPress;
nDensity_ = nPress;
-#ifndef NDEBUG
- warningPrinted_ = false;
-#endif
-
- std::cout << "-------------------------------------------------------------------------\n"
- << "Initializing tables for the " << RawComponent::name()
- << " fluid properties (" << nTemp*nPress << " entries).\n"
- << "Temperature -> min: " << std::scientific << std::setprecision(3)
- << tempMin << ", max: " << tempMax << ", n: " << nTemp << '\n'
- << "Pressure -> min: " << std::scientific << std::setprecision(3)
- << pressMin << ", max: " << pressMax << ", n: " << nPress << '\n'
- << "-------------------------------------------------------------------------" << std::endl;
-
// resize & initialize the arrays with NaN
assert(std::numeric_limits<Scalar>::has_quiet_NaN);
const auto NaN = std::numeric_limits<Scalar>::quiet_NaN();
// initialize vapor pressure array depending on useVaporPressure
vaporPressure_.resize(nTemp_, NaN);
- initVaporPressure_();
+ tabularizeVaporPressure_();
- if constexpr (ComponentTraits<ThisType>::hasGasState)
+ if constexpr (ComponentTraits<RawComponent>::hasGasState)
{
minGasDensity_.resize(nTemp_, NaN);
maxGasDensity_.resize(nTemp_, NaN);
@@ -203,16 +172,16 @@ public:
if constexpr (RawComponent::gasIsCompressible())
gasPressure_.resize(numEntriesTp, NaN);
- minMaxGasDensityInitialized_ = false;
- gasEnthalpyInitialized_ = false;
- gasHeatCapacityInitialized_ = false;
- gasDensityInitialized_ = false;
- gasViscosityInitialized_ = false;
- gasThermalConductivityInitialized_ = false;
- gasPressureInitialized_ = false;
+ minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_();
+ gasPressureInitialized_ = tabularizeGasPressure_();
+ gasEnthalpyInitialized_ = tabularizeGasEnthalpy_();
+ gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_();
+ gasDensityInitialized_ = tabularizeGasDensity_();
+ gasViscosityInitialized_ = tabularizeGasViscosity_();
+ gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_();
}
- if constexpr (ComponentTraits<ThisType>::hasLiquidState)
+ if constexpr (ComponentTraits<RawComponent>::hasLiquidState)
{
minLiquidDensity_.resize(nTemp_, NaN);
maxLiquidDensity_.resize(nTemp_, NaN);
@@ -227,1017 +196,1051 @@ public:
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;
+ minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_();
+ liquidPressureInitialized_ = tabularizeLiquidPressure_();
+ liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_();
+ liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_();
+ liquidDensityInitialized_ = tabularizeLiquidDensity_();
+ liquidViscosityInitialized_ = tabularizeLiquidViscosity_();
+ liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_();
}
+ }
- // in the multithreaded case precompute all array because lazy-loading
- // leads to potential data races
- if constexpr (!Multithreading::isSerial())
- {
- if constexpr (ComponentTraits<ThisType>::hasGasState)
- {
- minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_();
- gasPressureInitialized_ = tabularizeGasPressure_();
- gasEnthalpyInitialized_ = tabularizeGasEnthalpy_();
- gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_();
- gasDensityInitialized_ = tabularizeGasDensity_();
- gasViscosityInitialized_ = tabularizeGasViscosity_();
- gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_();
- }
-
- if constexpr (ComponentTraits<ThisType>::hasLiquidState)
- {
- minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_();
- liquidPressureInitialized_ = tabularizeLiquidPressure_();
- liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_();
- liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_();
- liquidDensityInitialized_ = tabularizeLiquidDensity_();
- liquidViscosityInitialized_ = tabularizeLiquidViscosity_();
- liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_();
- }
- }
+ //! returns the index of an entry in a temperature field
+ inline Scalar tempIdx(Scalar temperature) const
+ {
+ return (nTemp_ - 1)*(temperature - tempMin_)/(tempMax_ - tempMin_);
+ }
- initialized_ = true;
+ //! returns the index of an entry in a pressure field
+ inline Scalar pressLiquidIdx(Scalar pressure, std::size_t tempIdx) const
+ {
+ const Scalar plMin = minLiquidPressure(tempIdx);
+ const Scalar plMax = maxLiquidPressure(tempIdx);
+ return (nPress_ - 1)*(pressure - plMin)/(plMax - plMin);
}
- /*!
- * \brief A human readable name for the component.
- */
- static std::string name()
- { return RawComponent::name(); }
+ //! returns the index of an entry in a temperature field
+ inline Scalar pressGasIdx(Scalar pressure, std::size_t tempIdx) const
+ {
+ const Scalar pgMin = minGasPressure(tempIdx);
+ const Scalar pgMax = maxGasPressure(tempIdx);
+ return (nPress_ - 1)*(pressure - pgMin)/(pgMax - pgMin);
+ }
- /*!
- * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of the component.
- */
- static constexpr Scalar molarMass()
- { return RawComponent::molarMass(); }
+ //! returns the index of an entry in a density field
+ inline Scalar densityLiquidIdx(Scalar density, std::size_t tempIdx) const
+ {
+ const Scalar densityMin = minLiquidDensity_[tempIdx];
+ const Scalar densityMax = maxLiquidDensity_[tempIdx];
+ return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
+ }
- /*!
- * \brief Returns the critical temperature in \f$\mathrm{[K]}\f$ of the component.
- */
- static Scalar criticalTemperature()
- { return RawComponent::criticalTemperature(); }
+ //! returns the index of an entry in a density field
+ inline Scalar densityGasIdx(Scalar density, std::size_t tempIdx) const
+ {
+ const Scalar densityMin = minGasDensity_[tempIdx];
+ const Scalar densityMax = maxGasDensity_[tempIdx];
+ return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
+ }
- /*!
- * \brief Returns the critical pressure in \f$\mathrm{[Pa]}\f$ of the component.
- */
- static Scalar criticalPressure()
- { return RawComponent::criticalPressure(); }
+ //! returns the minimum tabularized liquid pressure at a given temperature index
+ inline Scalar minLiquidPressure(int tempIdx) const
+ {
+ using std::max;
+ if (!useVaporPressure)
+ return pressMin_;
+ else
+ return max(pressMin_, vaporPressure_[tempIdx] / 1.1);
+ }
- /*!
- * \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
- */
- static Scalar tripleTemperature()
- { return RawComponent::tripleTemperature(); }
+ //! returns the maximum tabularized liquid pressure at a given temperature index
+ inline Scalar maxLiquidPressure(int tempIdx) const
+ {
+ using std::max;
+ if (!useVaporPressure)
+ return pressMax_;
+ else
+ return max(pressMax_, vaporPressure_[tempIdx] * 1.1);
+ }
- /*!
- * \brief Returns the pressure in \f$\mathrm{[Pa]}\f$ at the component's triple point.
- */
- static Scalar triplePressure()
- { return RawComponent::triplePressure(); }
+ //! returns the minimum tabularized gas pressure at a given temperature index
+ inline Scalar minGasPressure(int tempIdx) const
+ {
+ using std::min;
+ if (!useVaporPressure)
+ return pressMin_;
+ else
+ return min(pressMin_, vaporPressure_[tempIdx] / 1.1 );
+ }
- /*!
- * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
- * temperature.
- *
- * \param T temperature of component
- */
- static Scalar vaporPressure(Scalar T)
+ //! returns the maximum tabularized gas pressure at a given temperature index
+ inline Scalar maxGasPressure(int tempIdx) const
{
- using std::isnan;
- Scalar result = interpolateT_(vaporPressure_, T);
- if (isnan(result))
- return RawComponent::vaporPressure(T);
- return result;
+ using std::min;
+ if (!useVaporPressure)
+ return pressMax_;
+ else
+ return min(pressMax_, vaporPressure_[tempIdx] * 1.1);
}
- /*!
- * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
- * temperature.
- *
- * The method is only called by the sequential flash, so tabulating is omitted.
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar vaporTemperature(Scalar pressure)
+ //! returns the minimum tabularized liquid density at a given temperature index
+ inline Scalar minLiquidDensity(int tempIdx) const
+ { return minLiquidDensity_[tempIdx]; }
+
+ //! returns the maximum tabularized liquid density at a given temperature index
+ inline Scalar maxLiquidDensity(int tempIdx) const
+ { return maxLiquidDensity_[tempIdx]; }
+
+ //! returns the minimum tabularized gas density at a given temperature index
+ inline Scalar minGasDensity(int tempIdx) const
+ { return minGasDensity_[tempIdx]; }
+
+ //! returns the maximum tabularized gas density at a given temperature index
+ inline Scalar maxGasDensity(int tempIdx) const
+ { return maxGasDensity_[tempIdx]; }
+
+ inline std::size_t nTemp() const { return nTemp_; }
+ inline std::size_t nPress() const { return nPress_; }
+ inline std::size_t nDensity() const { return nDensity_; }
+
+ inline Scalar tempMax() const { return tempMax_; }
+ inline Scalar tempMin() const { return tempMin_; }
+
+private:
+
+ //! initializes vapor pressure if useVaporPressure = true
+ template< bool useVP = useVaporPressure, std::enable_if_t<useVP, int> = 0 >
+ void tabularizeVaporPressure_()
{
- return RawComponent::vaporTemperature(pressure);
+ // fill the temperature-pressure arrays
+ Dumux::parallelFor(nTemp_, [=](std::size_t iT)
+ {
+ Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
+ vaporPressure_[iT] = RawComponent::vaporPressure(temperature);
+ });
}
+ //! if !useVaporPressure, do nothing here
+ template< bool useVP = useVaporPressure, std::enable_if_t<!useVP, int> = 0 >
+ void tabularizeVaporPressure_() {}
+
/*!
- * \brief Specific enthalpy of the gas \f$\mathrm{[J/kg]}\f$.
+ * \brief Initializes property values as function of temperature and pressure.
*
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \tparam PropFunc Function to evaluate the property prop(T, p)
+ * \tparam MinPFunc Function to evaluate the minimum pressure for a
+ * temperature index (depends on useVaporPressure)
+ * \tparam MaxPFunc Function to evaluate the maximum pressure for a
+ * temperature index (depends on useVaporPressure)
+ *
+ * \param f property function
+ * \param minP function to evaluate minimum pressure for temp idx
+ * \param maxP function to evaluate maximum pressure for temp idx
+ * \param values container to store property values
*/
- static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
+ template<class PropFunc, class Policy>
+ void tabularizeTPArray_(const PropFunc& f, Policy policy, std::vector<typename RawComponent::Scalar>& values) const
{
- Scalar result = interpolateTP_(gasEnthalpy_, temperature, pressure,
- pressGasIdx_, minGasPressure_, maxGasPressure_);
- using std::isnan;
- if (isnan(result))
+ Dumux::parallelFor(nTemp_, [=,&values](std::size_t iT)
{
- if constexpr (Multithreading::isSerial())
+ Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
+
+ Scalar pMax = policy.maxP(iT);
+ Scalar pMin = policy.minP(iT);
+ for (std::size_t iP = 0; iP < nPress_; ++ iP)
{
- // lazy tabularization
- if (!gasEnthalpyInitialized_)
- {
- gasEnthalpyInitialized_ = tabularizeGasEnthalpy_();
- if (gasEnthalpyInitialized_)
- return gasEnthalpy(temperature, pressure);
- }
+ Scalar pressure = iP * (pMax - pMin)/(nPress_ - 1) + pMin;
+ values[iT + iP*nTemp_] = f(temperature, pressure);
}
-
- printWarning_("gasEnthalpy", temperature, pressure);
- return RawComponent::gasEnthalpy(temperature, pressure);
- }
- return result;
+ });
}
/*!
- * \brief Specific enthalpy of the liquid \f$\mathrm{[J/kg]}\f$.
+ * \brief Initializes the minimum/maximum densities on the temperature range.
*
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \tparam RhoFunc Function to evaluate the density rho(T, p)
+ * \tparam MinPFunc Function to evaluate the minimum pressure for a
+ * temperature index (depends on useVaporPressure)
+ * \tparam MaxPFunc Function to evaluate the maximum pressure for a
+ * temperature index (depends on useVaporPressure)
+ *
+ * \param rho density function
+ * \param minP function to evaluate minimum pressure for temp idx
+ * \param maxP function to evaluate maximum pressure for temp idx
+ * \param rhoMin container to store minimum density values
+ * \param rhoMax container to store maximum density values
*/
- static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
+ template<class RhoFunc, class Policy>
+ void tabularizeMinMaxRhoArray_(const RhoFunc& rho, Policy policy,
+ std::vector<typename RawComponent::Scalar>& rhoMin,
+ std::vector<typename RawComponent::Scalar>& rhoMax) const
{
- Scalar result = interpolateTP_(liquidEnthalpy_, temperature, pressure,
- pressLiquidIdx_, minLiquidPressure_, maxLiquidPressure_);
- using std::isnan;
- if (isnan(result))
+ Dumux::parallelFor(nTemp_, [=,&rhoMin,&rhoMax](std::size_t iT)
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidEnthalpyInitialized_)
- {
- liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_();
- if (liquidEnthalpyInitialized_)
- return liquidEnthalpy(temperature, pressure);
- }
- }
+ Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
- printWarning_("liquidEnthalpy", temperature, pressure);
- return RawComponent::liquidEnthalpy(temperature, pressure);
- }
- return result;
+ using std::min;
+ rhoMin[iT] = rho(temperature, policy.minP(iT));
+ rhoMax[iT] = rho(temperature, policy.maxP(min(iT + 1, nTemp_ - 1)));
+ });
}
/*!
- * \brief Specific isobaric heat capacity of the gas \f$\mathrm{[J/(kg*K)]}\f$.
+ * \brief Initializes pressure arrays as function of temperature and density.
*
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \tparam PFunc Function to evaluate the pressure p(T, rho)
+ *
+ * \param pressure container to store pressure values
+ * \param p pressure function p(T, rho)
+ * \param rhoMin container with minimum density values
+ * \param rhoMax container with maximum density values
*/
- static const Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
+ template<class PFunc>
+ void tabularizePressureArray_(std::vector<typename RawComponent::Scalar>& pressure,
+ const PFunc& p,
+ const std::vector<typename RawComponent::Scalar>& rhoMin,
+ const std::vector<typename RawComponent::Scalar>& rhoMax) const
{
- Scalar result = interpolateTP_(gasHeatCapacity_, temperature, pressure,
- pressGasIdx_, minGasPressure_, maxGasPressure_);
- using std::isnan;
- if (isnan(result))
+ Dumux::parallelFor(nTemp_, [=,&pressure,&rhoMin,&rhoMax](std::size_t iT)
{
- if constexpr (Multithreading::isSerial())
+ Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
+
+ for (std::size_t iRho = 0; iRho < nDensity_; ++ iRho)
{
- // lazy tabularization
- if (!gasHeatCapacityInitialized_)
- {
- gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_();
- if (gasHeatCapacityInitialized_)
- return gasHeatCapacity(temperature, pressure);
- }
+ Scalar density = Scalar(iRho)/(nDensity_ - 1)
+ * (rhoMax[iT] - rhoMin[iT])
+ + rhoMin[iT];
+ pressure[iT + iRho*nTemp_] = p(temperature, density);
}
-
- printWarning_("gasHeatCapacity", temperature, pressure);
- return RawComponent::gasHeatCapacity(temperature, pressure);
- }
- return result;
+ });
}
- /*!
- * \brief Specific isobaric heat capacity of the liquid \f$\mathrm{[J/(kg*K)]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static const Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeGasEnthalpy_()
{
- Scalar result = interpolateTP_(liquidHeatCapacity_, temperature, pressure,
- pressLiquidIdx_, minLiquidPressure_, maxLiquidPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasGasEnthalpy<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidHeatCapacityInitialized_)
- {
- liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_();
- if (liquidHeatCapacityInitialized_)
- return liquidHeatCapacity(temperature, pressure);
- }
- }
-
- printWarning_("liquidHeatCapacity", temperature, pressure);
- return RawComponent::liquidHeatCapacity(temperature, pressure);
+ const auto gasEnth = [] (auto T, auto p) { return RC::gasEnthalpy(T, p); };
+ tabularizeTPArray_(gasEnth, GasPolicy{ *this }, gasEnthalpy_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief Specific internal energy of the gas \f$\mathrm{[J/kg]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeLiquidEnthalpy_()
{
- return gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure);
+ if constexpr (Detail::hasLiquidEnthalpy<RC>())
+ {
+ const auto liqEnth = [] (auto T, auto p) { return RC::liquidEnthalpy(T, p); };
+ tabularizeTPArray_(liqEnth, LiquidPolicy{ *this }, liquidEnthalpy_);
+ return true;
+ }
+
+ return false;
}
- /*!
- * \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeGasHeatCapacity_()
{
- return liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
+ if constexpr (Detail::hasGasHeatCapacity<RC>())
+ {
+ const auto gasHC = [] (auto T, auto p) { return RC::gasHeatCapacity(T, p); };
+ tabularizeTPArray_(gasHC, GasPolicy{ *this }, gasHeatCapacity_);
+ return true;
+ }
+
+ return false;
}
- /*!
- * \brief The pressure of gas in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
- */
- static Scalar gasPressure(Scalar temperature, Scalar density)
+ template<class RC = RawComponent>
+ bool tabularizeLiquidHeatCapacity_()
{
- if constexpr (Multithreading::isSerial())
+ if constexpr (Detail::hasLiquidHeatCapacity<RC>())
{
- if (!minMaxGasDensityInitialized_)
- minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_();
+ const auto liqHC = [] (auto T, auto p) { return RC::liquidHeatCapacity(T, p); };
+ tabularizeTPArray_(liqHC, LiquidPolicy{ *this }, liquidHeatCapacity_);
+ return true;
}
- Scalar result = interpolateTRho_(gasPressure_, temperature, density, densityGasIdx_);
- using std::isnan;
- if (isnan(result))
- {
- if constexpr (Multithreading::isSerial())
- {
- if (!gasPressureInitialized_)
- {
- gasPressureInitialized_ = tabularizeGasPressure_();
- if (gasPressureInitialized_)
- return gasPressure(temperature, density);
- }
- }
+ return false;
+ }
- printWarning_("gasPressure", temperature, density);
- return RawComponent::gasPressure(temperature, density);
+ template<class RC = RawComponent>
+ bool tabularizeMinMaxGasDensity_()
+ {
+ if constexpr (Detail::hasGasDensity<RC>())
+ {
+ const auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
+ tabularizeMinMaxRhoArray_(gasRho, GasPolicy{ *this }, minGasDensity_, maxGasDensity_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief The pressure of liquid in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
- */
- static Scalar liquidPressure(Scalar temperature, Scalar density)
+ template<class RC = RawComponent>
+ bool tabularizeMinMaxLiquidDensity_()
{
- if constexpr (Multithreading::isSerial())
+ if constexpr (Detail::hasGasEnthalpy<RC>())
{
- if (!minMaxLiquidDensityInitialized_)
- minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_();
+ const auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); };
+ tabularizeMinMaxRhoArray_(liqRho, LiquidPolicy{ *this }, minLiquidDensity_, maxLiquidDensity_);
+ return true;
}
- Scalar result = interpolateTRho_(liquidPressure_, temperature, density, densityLiquidIdx_);
- using std::isnan;
- if (isnan(result))
- {
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidPressureInitialized_)
- {
- liquidPressureInitialized_ = tabularizeLiquidPressure_();
- if (liquidPressureInitialized_)
- return liquidPressure(temperature, density);
- }
- }
+ return false;
+ }
- printWarning_("liquidPressure", temperature, density);
- return RawComponent::liquidPressure(temperature, density);
+ template<class RC = RawComponent>
+ bool tabularizeGasPressure_()
+ {
+ // pressure is only defined if the gas is compressible (this is usually the case)
+ if constexpr (Detail::hasGasPressure<RC>() && RC::gasIsCompressible())
+ {
+ const auto gasPFunc = [] (auto T, auto rho) { return RC::gasPressure(T, rho); };
+ tabularizePressureArray_(gasPressure_, gasPFunc, minGasDensity_, maxGasDensity_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief Returns true if the gas phase is assumed to be compressible
- */
- static constexpr bool gasIsCompressible()
- { return RawComponent::gasIsCompressible(); }
+ template<class RC = RawComponent>
+ bool tabularizeLiquidPressure_()
+ {
+ // pressure is only defined if the liquid is compressible (this is often not the case)
+ if constexpr (Detail::hasLiquidPressure<RC>() && RC::liquidIsCompressible())
+ {
+ const auto liqPFunc = [] (auto T, auto rho) { return RC::liquidPressure(T, rho); };
+ tabularizePressureArray_(liquidPressure_, liqPFunc, minLiquidDensity_, maxLiquidDensity_);
+ return true;
+ }
- /*!
- * \brief Returns true if the liquid phase is assumed to be compressible
- */
- static constexpr bool liquidIsCompressible()
- { return RawComponent::liquidIsCompressible(); }
+ return false;
+ }
- /*!
- * \brief Returns true if the gas phase is assumed to be ideal
- */
- static constexpr bool gasIsIdeal()
- { return RawComponent::gasIsIdeal(); }
+ template<class RC = RawComponent>
+ bool tabularizeGasDensity_()
+ {
+ if constexpr (Detail::hasGasDensity<RC>())
+ {
+ const auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
+ tabularizeTPArray_(gasRho, GasPolicy{ *this }, gasDensity_);
+ return true;
+ }
+ return false;
+ }
- /*!
- * \brief The density of gas at a given pressure and temperature
- * \f$\mathrm{[kg/m^3]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar gasDensity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeLiquidDensity_()
{
- Scalar result = interpolateTP_(gasDensity_, temperature, pressure,
- pressGasIdx_, minGasPressure_, maxGasPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasLiquidDensity<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!gasDensityInitialized_)
- {
- gasDensityInitialized_ = tabularizeGasDensity_();
- if (gasDensityInitialized_)
- return gasDensity(temperature, pressure);
- }
- }
-
- printWarning_("gasDensity", temperature, pressure);
- return RawComponent::gasDensity(temperature, pressure);
+ // 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.
+ const auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); };
+ tabularizeTPArray_(liqRho, LiquidPolicy{ *this }, liquidDensity_);
+ return true;
}
- return result;
- }
- /*!
- * \brief The molar density of gas in \f$\mathrm{[mol/m^3]}\f$
- * at a given pressure and temperature.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- *
- */
- static Scalar gasMolarDensity(Scalar temperature, Scalar pressure)
- { return gasDensity(temperature, pressure)/molarMass(); }
+ return false;
+ }
- /*!
- * \brief The density of liquid at a given pressure and
- * temperature \f$\mathrm{[kg/m^3]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeGasViscosity_()
{
- Scalar result = interpolateTP_(liquidDensity_, temperature, pressure,
- pressLiquidIdx_, minLiquidPressure_, maxLiquidPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasGasViscosity<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidDensityInitialized_)
- {
- liquidDensityInitialized_ = tabularizeLiquidDensity_();
- if (liquidDensityInitialized_)
- return liquidDensity(temperature, pressure);
- }
- }
-
- printWarning_("liquidDensity", temperature, pressure);
- return RawComponent::liquidDensity(temperature, pressure);
+ const auto gasVisc = [] (auto T, auto p) { return RC::gasViscosity(T, p); };
+ tabularizeTPArray_(gasVisc, GasPolicy{ *this }, gasViscosity_);
+ return true;
}
- return result;
+ return false;
}
- /*!
- * \brief The molar density of liquid in \f$\mathrm{[mol/m^3]}\f$
- * at a given pressure and temperature.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- *
- */
- static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
- { return liquidDensity(temperature, pressure)/molarMass(); }
-
- /*!
- * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of gas.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar gasViscosity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeLiquidViscosity_()
{
- Scalar result = interpolateTP_(gasViscosity_, temperature, pressure,
- pressGasIdx_, minGasPressure_, maxGasPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasLiquidViscosity<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!gasViscosityInitialized_)
- {
- gasViscosityInitialized_ = tabularizeGasViscosity_();
- if (gasViscosityInitialized_)
- return gasViscosity(temperature, pressure);
- }
- }
-
- printWarning_("gasViscosity", temperature, pressure);
- return RawComponent::gasViscosity(temperature, pressure);
+ const auto liqVisc = [] (auto T, auto p) { return RC::liquidViscosity(T, p); };
+ tabularizeTPArray_(liqVisc, LiquidPolicy{ *this }, liquidViscosity_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of liquid.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeGasThermalConductivity_()
{
- Scalar result = interpolateTP_(liquidViscosity_, temperature, pressure,
- pressLiquidIdx_, minLiquidPressure_, maxLiquidPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasGasThermalConductivity<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidViscosityInitialized_)
- {
- liquidViscosityInitialized_ = tabularizeLiquidViscosity_();
- if (liquidViscosityInitialized_)
- return liquidViscosity(temperature, pressure);
- }
- }
-
- printWarning_("liquidViscosity",temperature, pressure);
- return RawComponent::liquidViscosity(temperature, pressure);
+ const auto gasTC = [] (auto T, auto p) { return RC::gasThermalConductivity(T, p); };
+ tabularizeTPArray_(gasTC, GasPolicy{ *this }, gasThermalConductivity_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief The thermal conductivity of gaseous water \f$\mathrm{[W/(m*K)]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
+ template<class RC = RawComponent>
+ bool tabularizeLiquidThermalConductivity_()
{
- Scalar result = interpolateTP_(gasThermalConductivity_, temperature, pressure,
- pressGasIdx_, minGasPressure_, maxGasPressure_);
- using std::isnan;
- if (isnan(result))
+ if constexpr (Detail::hasLiquidThermalConductivity<RC>())
{
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!gasThermalConductivityInitialized_)
- {
- gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_();
- if (gasThermalConductivityInitialized_)
- return gasThermalConductivity(temperature, pressure);
- }
- }
-
- printWarning_("gasThermalConductivity", temperature, pressure);
- return RawComponent::gasThermalConductivity(temperature, pressure);
+ const auto liqTC = [] (auto T, auto p) { return RC::liquidThermalConductivity(T, p); };
+ tabularizeTPArray_(liqTC, LiquidPolicy{ *this }, liquidThermalConductivity_);
+ return true;
}
- return result;
+
+ return false;
}
- /*!
- * \brief The thermal conductivity of liquid water \f$\mathrm{[W/(m*K)]}\f$.
- *
- * \param temperature temperature of component in \f$\mathrm{[K]}\f$
- * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
- */
- static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
- {
- Scalar result = interpolateTP_(liquidThermalConductivity_, temperature, pressure,
- pressLiquidIdx_, minLiquidPressure_, maxLiquidPressure_);
- using std::isnan;
- if (isnan(result))
- {
- if constexpr (Multithreading::isSerial())
- {
- // lazy tabularization
- if (!liquidThermalConductivityInitialized_)
- {
- liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_();
- if (liquidThermalConductivityInitialized_)
- return liquidThermalConductivity(temperature, pressure);
- }
- }
+private:
+ // 1D fields with the temperature as degree of freedom
+ std::vector<Scalar> vaporPressure_;
- printWarning_("liquidThermalConductivity", temperature, pressure);
- return RawComponent::liquidThermalConductivity(temperature, pressure);
- }
- return result;
- }
+ std::vector<Scalar> minLiquidDensity_;
+ std::vector<Scalar> maxLiquidDensity_;
+ bool minMaxLiquidDensityInitialized_;
+ std::vector<Scalar> minGasDensity_;
+ std::vector<Scalar> maxGasDensity_;
+ bool minMaxGasDensityInitialized_;
-private:
- //! prints a warning if the result is not in range or the table has not been initialized
- static void printWarning_(const std::string& quantity, Scalar arg1, Scalar arg2)
+ // 2D fields with the temperature and pressure as degrees of freedom
+ std::vector<Scalar> gasEnthalpy_;
+ std::vector<Scalar> liquidEnthalpy_;
+ bool gasEnthalpyInitialized_;
+ bool liquidEnthalpyInitialized_;
+
+ std::vector<Scalar> gasHeatCapacity_;
+ std::vector<Scalar> liquidHeatCapacity_;
+ bool gasHeatCapacityInitialized_;
+ bool liquidHeatCapacityInitialized_;
+
+ std::vector<Scalar> gasDensity_;
+ std::vector<Scalar> liquidDensity_;
+ bool gasDensityInitialized_;
+ bool liquidDensityInitialized_;
+
+ std::vector<Scalar> gasViscosity_;
+ std::vector<Scalar> liquidViscosity_;
+ bool gasViscosityInitialized_;
+ bool liquidViscosityInitialized_;
+
+ std::vector<Scalar> gasThermalConductivity_;
+ std::vector<Scalar> liquidThermalConductivity_;
+ bool gasThermalConductivityInitialized_;
+ bool liquidThermalConductivityInitialized_;
+
+ // 2D fields with the temperature and density as degrees of freedom
+ std::vector<Scalar> gasPressure_;
+ std::vector<Scalar> liquidPressure_;
+ bool gasPressureInitialized_;
+ bool liquidPressureInitialized_;
+
+ // temperature, pressure and density ranges
+ Scalar tempMin_;
+ Scalar tempMax_;
+ std::size_t nTemp_;
+
+ Scalar pressMin_;
+ Scalar pressMax_;
+ std::size_t nPress_;
+
+ Scalar densityMin_;
+ Scalar densityMax_;
+ std::size_t nDensity_;
+};
+
+} // end namespace Dumux::Components::Detail
+
+namespace Dumux::Components {
+
+/*!
+ * \ingroup Components
+ * \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.
+ *
+ * \tparam Scalar The type used for scalar values
+ * \tparam RawComponent The component which ought to be tabulated
+ * \tparam useVaporPressure If set to true, the min/max pressure
+ * values for gas&liquid phase will be set
+ * depending on the vapor pressure.
+ */
+template <class RawComponent, bool useVaporPressure=true>
+class TabulatedComponent
+{
+ using ThisType = TabulatedComponent<RawComponent, useVaporPressure>;
+ using Table = Detail::TabulatedComponentTable<RawComponent, useVaporPressure>;
+
+ struct InterpolatePolicy
{
-#ifndef NDEBUG
- if (warningPrinted_)
- return;
+ using Scalar = typename RawComponent::Scalar;
- if (!initialized_)
- std::cerr << "Warning: tabulated component '" << name()
- << "' has not been initialized. "
- << "Call FluidSystem::init() to use the tabulation in order to reduce runtime. \n";
- else
- std::cerr << "Warning: "<<quantity<<"(T="<<arg1<<", p="<<arg2<<") of component '"<<name()
- << "' is outside tabulation range: ("<< tempMin_<<"<=T<="<<tempMax_<<"), ("
- << pressMin_<<"<=p<=" <<pressMax_<<"). "
- << "Forwarded to FluidSystem for direct evaluation of "<<quantity<<". \n";
- warningPrinted_ = true;
-#endif
- }
+ const Table& table;
- //! initializes vapor pressure if useVaporPressure = true
- template< bool useVP = useVaporPressure, std::enable_if_t<useVP, int> = 0 >
- static void initVaporPressure_()
+ Scalar tempIdx(Scalar T) const { return table.tempIdx(T); }
+ Scalar tempMin() const { return table.tempMin(); }
+ Scalar tempMax() const { return table.tempMax(); }
+ std::size_t nTemp() const { return table.nTemp(); }
+ std::size_t nPress() const { return table.nPress(); }
+ std::size_t nDensity() const { return table.nDensity(); }
+ };
+
+ struct InterpolateGasPolicy : public InterpolatePolicy
{
- // fill the temperature-pressure arrays
- Dumux::parallelFor(nTemp_, [&](std::size_t iT)
- {
- Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
- vaporPressure_[iT] = RawComponent::vaporPressure(temperature);
- });
- }
+ using Scalar = typename RawComponent::Scalar;
+ Scalar pressIdx(Scalar p, std::size_t tempIdx) const { return this->table.pressGasIdx(p, tempIdx); }
+ Scalar rhoIdx(Scalar rho, std::size_t tempIdx) const { return this->table.densityGasIdx(rho, tempIdx); }
+ Scalar minP(int tempIdx) const { return this->table.minGasPressure(tempIdx); }
+ Scalar maxP(int tempIdx) const { return this->table.maxGasPressure(tempIdx); }
+ Scalar minRho(int tempIdx) const { return this->table.minGasDensity(tempIdx); }
+ Scalar maxRho(int tempIdx) const { return this->table.maxGasDensity(tempIdx); }
+ };
+
+ struct InterpolateLiquidPolicy : public InterpolatePolicy
+ {
+ using Scalar = typename RawComponent::Scalar;
+ Scalar pressIdx(Scalar p, std::size_t tempIdx) const { return this->table.pressLiquidIdx(p, tempIdx); }
+ Scalar rhoIdx(Scalar rho, std::size_t tempIdx) const { return this->table.densityLiquidIdx(rho, tempIdx); }
+ Scalar minP(int tempIdx) const { return this->table.minLiquidPressure(tempIdx); }
+ Scalar maxP(int tempIdx) const { return this->table.maxLiquidPressure(tempIdx); }
+ Scalar minRho(int tempIdx) const { return this->table.minLiquidDensity(tempIdx); }
+ Scalar maxRho(int tempIdx) const { return this->table.maxLiquidDensity(tempIdx); }
+ };
+public:
+ //! export scalar type
+ using Scalar = typename RawComponent::Scalar;
- //! if !useVaporPressure, do nothing here
- template< bool useVP = useVaporPressure, std::enable_if_t<!useVP, int> = 0 >
- static void initVaporPressure_() {}
+ //! state that we are tabulated
+ static constexpr bool isTabulated = true;
/*!
- * \brief Initializes property values as function of temperature and pressure.
- *
- * \tparam PropFunc Function to evaluate the property prop(T, p)
- * \tparam MinPFunc Function to evaluate the minimum pressure for a
- * temperature index (depends on useVaporPressure)
- * \tparam MaxPFunc Function to evaluate the maximum pressure for a
- * temperature index (depends on useVaporPressure)
+ * \brief Initialize the tables.
*
- * \param f property function
- * \param minP function to evaluate minimum pressure for temp idx
- * \param maxP function to evaluate maximum pressure for temp idx
- * \param values container to store property values
+ * \param tempMin The minimum of the temperature range in \f$\mathrm{[K]}\f$
+ * \param tempMax The maximum of the temperature range in \f$\mathrm{[K]}\f$
+ * \param nTemp The number of entries/steps within the temperature range
+ * \param pressMin The minimum of the pressure range in \f$\mathrm{[Pa]}\f$
+ * \param pressMax The maximum of the pressure range in \f$\mathrm{[Pa]}\f$
+ * \param nPress The number of entries/steps within the pressure range
*/
- template<class PropFunc, class MinPFunc, class MaxPFunc>
- static void initTPArray_(PropFunc&& f, MinPFunc&& minP, MaxPFunc&& maxP, std::vector<typename RawComponent::Scalar>& values)
+ static void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
+ Scalar pressMin, Scalar pressMax, std::size_t nPress)
{
- Dumux::parallelFor(nTemp_, [&](std::size_t iT)
- {
- Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
+#ifndef NDEBUG
+ warningPrinted_ = false;
+#endif
+ std::cout << "-------------------------------------------------------------------------\n"
+ << "Initializing tables for the " << RawComponent::name()
+ << " fluid properties (" << nTemp*nPress << " entries).\n"
+ << "Temperature -> min: " << std::scientific << std::setprecision(3)
+ << tempMin << ", max: " << tempMax << ", n: " << nTemp << '\n'
+ << "Pressure -> min: " << std::scientific << std::setprecision(3)
+ << pressMin << ", max: " << pressMax << ", n: " << nPress << '\n'
+ << "-------------------------------------------------------------------------" << std::endl;
- Scalar pMax = maxP(iT);
- Scalar pMin = minP(iT);
- for (std::size_t iP = 0; iP < nPress_; ++ iP)
- {
- Scalar pressure = iP * (pMax - pMin)/(nPress_ - 1) + pMin;
- values[iT + iP*nTemp_] = f(temperature, pressure);
- }
- });
+ table().init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
+
+ initialized_ = true;
}
/*!
- * \brief Initializes the minimum/maximum densities on the temperature range.
- *
- * \tparam RhoFunc Function to evaluate the density rho(T, p)
- * \tparam MinPFunc Function to evaluate the minimum pressure for a
- * temperature index (depends on useVaporPressure)
- * \tparam MaxPFunc Function to evaluate the maximum pressure for a
- * temperature index (depends on useVaporPressure)
- *
- * \param rho density function
- * \param minP function to evaluate minimum pressure for temp idx
- * \param maxP function to evaluate maximum pressure for temp idx
- * \param rhoMin container to store minimum density values
- * \param rhoMax container to store maximum density values
+ * \brief A human readable name for the component.
*/
- template<class RhoFunc, class MinPFunc, class MaxPFunc>
- static void initMinMaxRhoArray_(RhoFunc&& rho,
- MinPFunc&& minP,
- MaxPFunc&& maxP,
- std::vector<typename RawComponent::Scalar>& rhoMin,
- std::vector<typename RawComponent::Scalar>& rhoMax)
- {
- Dumux::parallelFor(nTemp_, [&](std::size_t iT)
- {
- Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
-
- rhoMin[iT] = rho(temperature, minP(iT));
- if (iT < nTemp_ - 1)
- rhoMax[iT] = rho(temperature, maxP(iT + 1));
- else
- rhoMax[iT] = rho(temperature, maxP(iT));
- });
- }
+ static std::string name()
+ { return RawComponent::name(); }
/*!
- * \brief Initializes pressure arrays as function of temperature and density.
- *
- * \tparam PFunc Function to evaluate the pressure p(T, rho)
- *
- * \param pressure container to store pressure values
- * \param p pressure function p(T, rho)
- * \param rhoMin container with minimum density values
- * \param rhoMax container with maximum density values
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of the component.
*/
- template<class PFunc>
- static void initPressureArray_(std::vector<typename RawComponent::Scalar>& pressure, PFunc&& p,
- const std::vector<typename RawComponent::Scalar>& rhoMin,
- const std::vector<typename RawComponent::Scalar>& rhoMax)
- {
- Dumux::parallelFor(nTemp_, [&](std::size_t iT)
- {
- Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
+ static constexpr Scalar molarMass()
+ { return RawComponent::molarMass(); }
- for (std::size_t iRho = 0; iRho < nDensity_; ++ iRho)
- {
- Scalar density = Scalar(iRho)/(nDensity_ - 1)
- * (rhoMax[iT] - rhoMin[iT])
- + rhoMin[iT];
- pressure[iT + iRho*nTemp_] = p(temperature, density);
- }
- });
- }
+ /*!
+ * \brief Returns the critical temperature in \f$\mathrm{[K]}\f$ of the component.
+ */
+ static Scalar criticalTemperature()
+ { return RawComponent::criticalTemperature(); }
- //! returns an interpolated value depending on temperature
- static Scalar interpolateT_(const std::vector<typename RawComponent::Scalar>& values, Scalar T)
- {
- Scalar alphaT = tempIdx_(T);
- if (alphaT < 0 || alphaT >= nTemp_ - 1)
- return std::numeric_limits<Scalar>::quiet_NaN();
+ /*!
+ * \brief Returns the critical pressure in \f$\mathrm{[Pa]}\f$ of the component.
+ */
+ static Scalar criticalPressure()
+ { return RawComponent::criticalPressure(); }
- const auto iT = static_cast<int>(alphaT);
- alphaT -= iT;
+ /*!
+ * \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
+ */
+ static Scalar tripleTemperature()
+ { return RawComponent::tripleTemperature(); }
- return values[iT ]*(1 - alphaT) +
- values[iT + 1]*( alphaT);
- }
+ /*!
+ * \brief Returns the pressure in \f$\mathrm{[Pa]}\f$ at the component's triple point.
+ */
+ static Scalar triplePressure()
+ { return RawComponent::triplePressure(); }
- //! returns an interpolated value depending on temperature and pressure
- template<class GetPIdx, class MinPFunc, class MaxPFunc>
- static Scalar interpolateTP_(const std::vector<typename RawComponent::Scalar>& values, Scalar T, Scalar p,
- GetPIdx&& getPIdx, MinPFunc&& minP, MaxPFunc&& maxP)
+ /*!
+ * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
+ * temperature.
+ *
+ * \param T temperature of component
+ */
+ static Scalar vaporPressure(Scalar T)
{
- Scalar alphaT = tempIdx_(T);
- if (alphaT < 0 || alphaT >= nTemp_ - 1) {
- return std::numeric_limits<Scalar>::quiet_NaN();
- }
- using std::clamp;
- const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp_ - 2);
- alphaT -= iT;
-
- Scalar alphaP1 = getPIdx(p, iT);
- Scalar alphaP2 = getPIdx(p, iT + 1);
-
- const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nPress_ - 2);
- const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nPress_ - 2);
- alphaP1 -= iP1;
- alphaP2 -= iP2;
-
-#if 0 && !defined NDEBUG
- if(!(0 <= alphaT && alphaT <= 1.0))
- DUNE_THROW(NumericalProblem, "Temperature out of range: "
- << "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
- if(!(0 <= alphaP1 && alphaP1 <= 1.0))
- DUNE_THROW(NumericalProblem, "First pressure out of range: "
- << "p=" << p << " range: [" << minP(tempIdx_(T)) << ", " << maxP(tempIdx_(T)) << "]");
- if(!(0 <= alphaP2 && alphaP2 <= 1.0))
- DUNE_THROW(NumericalProblem, "Second pressure out of range: "
- << "p=" << p << " range: [" << minP(tempIdx_(T) + 1) << ", " << maxP(tempIdx_(T) + 1) << "]");
-#endif
-
- return values[(iT ) + (iP1 )*nTemp_]*(1 - alphaT)*(1 - alphaP1) +
- values[(iT ) + (iP1 + 1)*nTemp_]*(1 - alphaT)*( alphaP1) +
- values[(iT + 1) + (iP2 )*nTemp_]*( alphaT)*(1 - alphaP2) +
- values[(iT + 1) + (iP2 + 1)*nTemp_]*( alphaT)*( alphaP2);
+ using std::isnan;
+ Scalar result = interpolateT_(table().vaporPressure_, T, InterpolatePolicy{{ table() }});
+ if (isnan(result))
+ return RawComponent::vaporPressure(T);
+ return result;
}
- //! returns an interpolated value for gas depending on temperature and density
- template<class GetRhoIdx>
- static Scalar interpolateTRho_(const std::vector<typename RawComponent::Scalar>& values, Scalar T, Scalar rho, GetRhoIdx&& rhoIdx)
+ /*!
+ * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
+ * temperature.
+ *
+ * The method is only called by the sequential flash, so tabulating is omitted.
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar vaporTemperature(Scalar pressure)
{
- using std::clamp;
- Scalar alphaT = tempIdx_(T);
- const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp_ - 2);
- alphaT -= iT;
-
- Scalar alphaP1 = rhoIdx(rho, iT);
- Scalar alphaP2 = rhoIdx(rho, iT + 1);
- const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nDensity_ - 2);
- const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nDensity_ - 2);
- alphaP1 -= iP1;
- alphaP2 -= iP2;
-
- return values[(iT ) + (iP1 )*nTemp_]*(1 - alphaT)*(1 - alphaP1) +
- values[(iT ) + (iP1 + 1)*nTemp_]*(1 - alphaT)*( alphaP1) +
- values[(iT + 1) + (iP2 )*nTemp_]*( alphaT)*(1 - alphaP2) +
- values[(iT + 1) + (iP2 + 1)*nTemp_]*( alphaT)*( alphaP2);
+ return RawComponent::vaporTemperature(pressure);
}
- //! returns the index of an entry in a temperature field
- static Scalar tempIdx_(Scalar temperature)
+ /*!
+ * \brief Specific enthalpy of the gas \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
{
- return (nTemp_ - 1)*(temperature - tempMin_)/(tempMax_ - tempMin_);
+ Scalar result = interpolateTP_(table().gasEnthalpy_, temperature, pressure, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTP_("gasEnthalpy", temperature, pressure, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasEnthalpy(temperature, pressure);
+ }
+ return result;
}
- //! returns the index of an entry in a pressure field
- static Scalar pressLiquidIdx_(Scalar pressure, std::size_t tempIdx)
+ /*!
+ * \brief Specific enthalpy of the liquid \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
{
- Scalar plMin = minLiquidPressure_(tempIdx);
- Scalar plMax = maxLiquidPressure_(tempIdx);
-
- return (nPress_ - 1)*(pressure - plMin)/(plMax - plMin);
+ Scalar result = interpolateTP_(table().liquidEnthalpy_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTP_("liquidEnthalpy", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidEnthalpy(temperature, pressure);
+ }
+ return result;
}
- //! returns the index of an entry in a temperature field
- static Scalar pressGasIdx_(Scalar pressure, std::size_t tempIdx)
+ /*!
+ * \brief Specific isobaric heat capacity of the gas \f$\mathrm{[J/(kg*K)]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
{
- Scalar pgMin = minGasPressure_(tempIdx);
- Scalar pgMax = maxGasPressure_(tempIdx);
-
- return (nPress_ - 1)*(pressure - pgMin)/(pgMax - pgMin);
+ Scalar result = interpolateTP_(table().gasHeatCapacity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTP_("gasHeatCapacity", temperature, pressure, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasHeatCapacity(temperature, pressure);
+ }
+ return result;
}
- //! returns the index of an entry in a density field
- static Scalar densityLiquidIdx_(Scalar density, std::size_t tempIdx)
+ /*!
+ * \brief Specific isobaric heat capacity of the liquid \f$\mathrm{[J/(kg*K)]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
{
- Scalar densityMin = minLiquidDensity_[tempIdx];
- Scalar densityMax = maxLiquidDensity_[tempIdx];
- return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
+ Scalar result = interpolateTP_(table().liquidHeatCapacity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTP_("liquidHeatCapacity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidHeatCapacity(temperature, pressure);
+ }
+ return result;
}
- //! returns the index of an entry in a density field
- static Scalar densityGasIdx_(Scalar density, std::size_t tempIdx)
+ /*!
+ * \brief Specific internal energy of the gas \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
{
- Scalar densityMin = minGasDensity_[tempIdx];
- Scalar densityMax = maxGasDensity_[tempIdx];
- return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
+ return gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure);
}
- //! returns the minimum tabularized liquid pressure at a given temperature index
- static Scalar minLiquidPressure_(int tempIdx)
+ /*!
+ * \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
{
- using std::max;
- if (!useVaporPressure)
- return pressMin_;
- else
- return max(pressMin_, vaporPressure_[tempIdx] / 1.1);
+ return liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
}
- //! returns the maximum tabularized liquid pressure at a given temperature index
- static Scalar maxLiquidPressure_(int tempIdx)
+ /*!
+ * \brief The pressure of gas in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
+ */
+ static Scalar gasPressure(Scalar temperature, Scalar density)
{
- using std::max;
- if (!useVaporPressure)
- return pressMax_;
- else
- return max(pressMax_, vaporPressure_[tempIdx] * 1.1);
+ Scalar result = interpolateTRho_(table().gasPressure_, temperature, density, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTRho_("gasPressure", temperature, density, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasPressure(temperature, density);
+ }
+ return result;
}
- //! returns the minimum tabularized gas pressure at a given temperature index
- static Scalar minGasPressure_(int tempIdx)
+ /*!
+ * \brief The pressure of liquid in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
+ */
+ static Scalar liquidPressure(Scalar temperature, Scalar density)
{
- using std::min;
- if (!useVaporPressure)
- return pressMin_;
- else
- return min(pressMin_, vaporPressure_[tempIdx] / 1.1 );
+ Scalar result = interpolateTRho_(table().liquidPressure_, temperature, density, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
+ {
+ printWarningTRho_("liquidPressure", temperature, density, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidPressure(temperature, density);
+ }
+ return result;
}
- //! returns the maximum tabularized gas pressure at a given temperature index
- static Scalar maxGasPressure_(int tempIdx)
- {
- using std::min;
- if (!useVaporPressure)
- return pressMax_;
- else
- return min(pressMax_, vaporPressure_[tempIdx] * 1.1);
- }
+ /*!
+ * \brief Returns true if the gas phase is assumed to be compressible
+ */
+ static constexpr bool gasIsCompressible()
+ { return RawComponent::gasIsCompressible(); }
- template<class RC = RawComponent>
- static bool tabularizeGasEnthalpy_()
+ /*!
+ * \brief Returns true if the liquid phase is assumed to be compressible
+ */
+ static constexpr bool liquidIsCompressible()
+ { return RawComponent::liquidIsCompressible(); }
+
+ /*!
+ * \brief Returns true if the gas phase is assumed to be ideal
+ */
+ static constexpr bool gasIsIdeal()
+ { return RawComponent::gasIsIdeal(); }
+
+
+ /*!
+ * \brief The density of gas at a given pressure and temperature
+ * \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasGasEnthalpy<RC>())
+ Scalar result = interpolateTP_(table().gasDensity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto gasEnth = [] (auto T, auto p) { return RC::gasEnthalpy(T, p); };
- initTPArray_(gasEnth, minGasPressure_, maxGasPressure_, gasEnthalpy_);
- return true;
+ printWarningTP_("gasDensity", temperature, pressure, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasDensity(temperature, pressure);
}
-
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeLiquidEnthalpy_()
+ /*!
+ * \brief The molar density of gas in \f$\mathrm{[mol/m^3]}\f$
+ * at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar gasMolarDensity(Scalar temperature, Scalar pressure)
+ { return gasDensity(temperature, pressure)/molarMass(); }
+
+ /*!
+ * \brief The density of liquid at a given pressure and
+ * temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasLiquidEnthalpy<RC>())
+ Scalar result = interpolateTP_(table().liquidDensity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto liqEnth = [] (auto T, auto p) { return RC::liquidEnthalpy(T, p); };
- initTPArray_(liqEnth, minLiquidPressure_, maxLiquidPressure_, liquidEnthalpy_);
- return true;
+ printWarningTP_("liquidDensity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidDensity(temperature, pressure);
}
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeGasHeatCapacity_()
+ /*!
+ * \brief The molar density of liquid in \f$\mathrm{[mol/m^3]}\f$
+ * at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
+ { return liquidDensity(temperature, pressure)/molarMass(); }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of gas.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasViscosity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasGasHeatCapacity<RC>())
+ Scalar result = interpolateTP_(table().gasViscosity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto gasHC = [] (auto T, auto p) { return RC::gasHeatCapacity(T, p); };
- initTPArray_(gasHC, minGasPressure_, maxGasPressure_, gasHeatCapacity_);
- return true;
+ printWarningTP_("gasViscosity", temperature, pressure, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasViscosity(temperature, pressure);
}
-
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeLiquidHeatCapacity_()
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of liquid.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasLiquidHeatCapacity<RC>())
+ Scalar result = interpolateTP_(table().liquidViscosity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto liqHC = [] (auto T, auto p) { return RC::liquidHeatCapacity(T, p); };
- initTPArray_(liqHC, minLiquidPressure_, maxLiquidPressure_, liquidHeatCapacity_);
- return true;
+ printWarningTP_("liquidViscosity",temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidViscosity(temperature, pressure);
}
-
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeMinMaxGasDensity_()
+ /*!
+ * \brief The thermal conductivity of gaseous water \f$\mathrm{[W/(m*K)]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasGasDensity<RC>())
+ Scalar result = interpolateTP_(table().gasThermalConductivity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
- initMinMaxRhoArray_(gasRho, minGasPressure_, maxGasPressure_, minGasDensity_, maxGasDensity_);
- return true;
+ printWarningTP_("gasThermalConductivity", temperature, pressure, InterpolateGasPolicy{{ table() }});
+ return RawComponent::gasThermalConductivity(temperature, pressure);
}
-
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeMinMaxLiquidDensity_()
+ /*!
+ * \brief The thermal conductivity of liquid water \f$\mathrm{[W/(m*K)]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
{
- if constexpr (Detail::hasGasEnthalpy<RC>())
+ Scalar result = interpolateTP_(table().liquidThermalConductivity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ using std::isnan;
+ if (isnan(result))
{
- auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); };
- initMinMaxRhoArray_(liqRho, minLiquidPressure_, maxLiquidPressure_, minLiquidDensity_, maxLiquidDensity_);
- return true;
+ printWarningTP_("liquidThermalConductivity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
+ return RawComponent::liquidThermalConductivity(temperature, pressure);
}
-
- return false;
+ return result;
}
- template<class RC = RawComponent>
- static bool tabularizeGasPressure_()
+
+private:
+ //! prints a warning if the result is not in range or the table has not been initialized
+ template<class Policy>
+ static void printWarningTP_(const std::string& quantity, Scalar T, Scalar p, Policy policy)
{
- // pressure is only defined if the gas is compressible (this is usually the case)
- if constexpr (Detail::hasGasPressure<RC>() && RC::gasIsCompressible())
- {
- auto gasPFunc = [] (auto T, auto rho) { return RC::gasPressure(T, rho); };
- initPressureArray_(gasPressure_, gasPFunc, minGasDensity_, maxGasDensity_);
- return true;
- }
+#ifndef NDEBUG
+ if (warningPrinted_)
+ return;
- return false;
+ if (!initialized_)
+ std::cerr << "Warning: tabulated component '" << name()
+ << "' has not been initialized. "
+ << "Call FluidSystem::init() to use the tabulation in order to reduce runtime. \n";
+ else
+ std::cerr << "Warning: "<<quantity<<"(T="<<T<<", p="<<p<<") of component '"<<name()
+ << "' is outside tabulation range: ("<< policy.tempMin() <<"<=T<="<< policy.tempMax() <<"), ("
+ << policy.minP(0) <<"<=p<=" << policy.maxP(policy.nTemp()-1) <<"). "
+ << "Forwarded to FluidSystem for direct evaluation of "<<quantity<<". \n";
+ warningPrinted_ = true;
+#endif
}
- template<class RC = RawComponent>
- static bool tabularizeLiquidPressure_()
+ //! prints a warning if the result is not in range or the table has not been initialized
+ template<class Policy>
+ static void printWarningTRho_(const std::string& quantity, Scalar T, Scalar rho, Policy policy)
{
- // pressure is only defined if the liquid is compressible (this is often not the case)
- if constexpr (Detail::hasLiquidPressure<RC>() && RC::liquidIsCompressible())
+#ifndef NDEBUG
+ if (warningPrinted_)
+ return;
+
+ if (!initialized_)
+ std::cerr << "Warning: tabulated component '" << name()
+ << "' has not been initialized. "
+ << "Call FluidSystem::init() to use the tabulation in order to reduce runtime. \n";
+ else
{
- auto liqPFunc = [] (auto T, auto rho) { return RC::liquidPressure(T, rho); };
- initPressureArray_(liquidPressure_, liqPFunc, minLiquidDensity_, maxLiquidDensity_);
- return true;
+ const auto [densityMin, densityMax] = [&]
+ {
+ const Scalar alphaT = policy.tempIdx(T);
+ using std::clamp;
+ const auto iT = clamp<int>(static_cast<int>(alphaT), 0, policy.nTemp() - 2);
+ return std::make_pair( policy.minRho(iT), policy.maxRho(iT) );
+ }();
+
+ std::cerr << "Warning: "<<quantity<<"(T="<<T<<", density="<<rho<<") of component '"<<name()
+ << "' is outside tabulation range: ("<< policy.tempMin() <<"<=T<="<< policy.tempMax() <<"), ("
+ << densityMin <<"<=density<=" << densityMax <<"). "
+ << "Forwarded to FluidSystem for direct evaluation of "<<quantity<<". \n";
}
-
- return false;
+ warningPrinted_ = true;
+#endif
}
- template<class RC = RawComponent>
- static bool tabularizeGasDensity_()
+ //! returns an interpolated value depending on temperature
+ template<class Policy>
+ static Scalar interpolateT_(const std::vector<Scalar>& values, Scalar T, Policy policy)
{
- if constexpr (Detail::hasGasDensity<RC>())
- {
- auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
- initTPArray_(gasRho, minGasPressure_, maxGasPressure_, gasDensity_);
- return true;
- }
+ Scalar alphaT = policy.tempIdx(T);
+ const auto nTemp = policy.nTemp();
+ if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
+ return std::numeric_limits<Scalar>::quiet_NaN();
- return false;
+ using std::clamp;
+ const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
+ alphaT -= iT;
+
+ return values[iT ]*(1 - alphaT) +
+ values[iT + 1]*( alphaT);
}
- template<class RC = RawComponent>
- static bool tabularizeLiquidDensity_()
+ //! returns an interpolated value depending on temperature and pressure
+ template<class Policy>
+ static Scalar interpolateTP_(const std::vector<Scalar>& values, const Scalar T, const Scalar p, Policy policy)
{
- if constexpr (Detail::hasLiquidDensity<RC>())
- {
- // 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;
- }
+ const auto nTemp = policy.nTemp();
+ const auto nPress = policy.nPress();
- return false;
- }
+ Scalar alphaT = policy.tempIdx(T);
+ if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
+ return std::numeric_limits<Scalar>::quiet_NaN();
- template<class RC = RawComponent>
- static bool tabularizeGasViscosity_()
- {
- if constexpr (Detail::hasGasViscosity<RC>())
- {
- auto gasVisc = [] (auto T, auto p) { return RC::gasViscosity(T, p); };
- initTPArray_(gasVisc, minGasPressure_, maxGasPressure_, gasViscosity_);
- return true;
- }
+ using std::clamp;
+ const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
+ alphaT -= iT;
- return false;
- }
+ Scalar alphaP1 = policy.pressIdx(p, iT);
+ Scalar alphaP2 = policy.pressIdx(p, iT + 1);
- template<class RC = RawComponent>
- static bool tabularizeLiquidViscosity_()
- {
- if constexpr (Detail::hasLiquidViscosity<RC>())
- {
- auto liqVisc = [] (auto T, auto p) { return RC::liquidViscosity(T, p); };
- initTPArray_(liqVisc, minLiquidPressure_, maxLiquidPressure_, liquidViscosity_);
- return true;
- }
+ const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nPress - 2);
+ const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nPress - 2);
+ alphaP1 -= iP1;
+ alphaP2 -= iP2;
- return false;
+#if 0 && !defined NDEBUG
+ const auto tempMin = policy.tempMin();
+ const auto tempMin = policy.tempMax();
+ if(!(0 <= alphaT && alphaT <= 1.0))
+ DUNE_THROW(NumericalProblem, "Temperature out of range: "
+ << "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
+ if(!(0 <= alphaP1 && alphaP1 <= 1.0))
+ DUNE_THROW(NumericalProblem, "First pressure out of range: "
+ << "p=" << p << " range: [" << policy.minP(policy.tempIdx(T)) << ", " << policy.maxP(policy.tempIdx(T)) << "]");
+ if(!(0 <= alphaP2 && alphaP2 <= 1.0))
+ DUNE_THROW(NumericalProblem, "Second pressure out of range: "
+ << "p=" << p << " range: [" << policy.minP(policy.tempIdx(T) + 1) << ", " << policy.maxP(policy.tempIdx(T) + 1) << "]");
+#endif
+ return values[(iT ) + (iP1 )*nTemp]*(1 - alphaT)*(1 - alphaP1) +
+ values[(iT ) + (iP1 + 1)*nTemp]*(1 - alphaT)*( alphaP1) +
+ values[(iT + 1) + (iP2 )*nTemp]*( alphaT)*(1 - alphaP2) +
+ values[(iT + 1) + (iP2 + 1)*nTemp]*( alphaT)*( alphaP2);
}
- template<class RC = RawComponent>
- static bool tabularizeGasThermalConductivity_()
+ //! returns an interpolated value for gas depending on temperature and density
+ template<class Policy>
+ static Scalar interpolateTRho_(const std::vector<Scalar>& values, const Scalar T, const Scalar rho, Policy policy)
{
- if constexpr (Detail::hasGasThermalConductivity<RC>())
- {
- auto gasTC = [] (auto T, auto p) { return RC::gasThermalConductivity(T, p); };
- initTPArray_(gasTC, minGasPressure_, maxGasPressure_, gasThermalConductivity_);
- return true;
- }
+ const auto nTemp = policy.nTemp();
+ const auto nDensity = policy.nDensity();
- return false;
- }
+ using std::clamp;
+ Scalar alphaT = policy.tempIdx(T);
+ if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
+ return std::numeric_limits<Scalar>::quiet_NaN();
- template<class RC = RawComponent>
- static bool tabularizeLiquidThermalConductivity_()
- {
- if constexpr (Detail::hasLiquidThermalConductivity<RC>())
- {
- auto liqTC = [] (auto T, auto p) { return RC::liquidThermalConductivity(T, p); };
- initTPArray_(liqTC, minLiquidPressure_, maxLiquidPressure_, liquidThermalConductivity_);
- return true;
- }
+ const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
+ alphaT -= iT;
- return false;
+ Scalar alphaP1 = policy.rhoIdx(rho, iT);
+ Scalar alphaP2 = policy.rhoIdx(rho, iT + 1);
+ const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nDensity - 2);
+ const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nDensity - 2);
+ alphaP1 -= iP1;
+ alphaP2 -= iP2;
+
+ return values[(iT ) + (iP1 )*nTemp]*(1 - alphaT)*(1 - alphaP1) +
+ values[(iT ) + (iP1 + 1)*nTemp]*(1 - alphaT)*( alphaP1) +
+ values[(iT + 1) + (iP2 )*nTemp]*( alphaT)*(1 - alphaP2) +
+ values[(iT + 1) + (iP2 + 1)*nTemp]*( alphaT)*( alphaP2);
}
// specifies whether the table was initialized
@@ -1248,64 +1251,13 @@ private:
static bool warningPrinted_;
#endif
- // 1D fields with the temperature as degree of freedom
- static std::vector<typename RawComponent::Scalar> vaporPressure_;
-
- static std::vector<typename RawComponent::Scalar> minLiquidDensity_;
- static std::vector<typename RawComponent::Scalar> maxLiquidDensity_;
- static bool minMaxLiquidDensityInitialized_;
-
- static std::vector<typename RawComponent::Scalar> minGasDensity_;
- static std::vector<typename RawComponent::Scalar> maxGasDensity_;
- static bool minMaxGasDensityInitialized_;
-
- // 2D fields with the temperature and pressure as degrees of freedom
- static std::vector<typename RawComponent::Scalar> gasEnthalpy_;
- static std::vector<typename RawComponent::Scalar> liquidEnthalpy_;
- static bool gasEnthalpyInitialized_;
- static bool liquidEnthalpyInitialized_;
-
- static std::vector<typename RawComponent::Scalar> gasHeatCapacity_;
- static std::vector<typename RawComponent::Scalar> liquidHeatCapacity_;
- static bool gasHeatCapacityInitialized_;
- static bool liquidHeatCapacityInitialized_;
-
- static std::vector<typename RawComponent::Scalar> gasDensity_;
- static std::vector<typename RawComponent::Scalar> liquidDensity_;
- static bool gasDensityInitialized_;
- static bool liquidDensityInitialized_;
-
- static std::vector<typename RawComponent::Scalar> gasViscosity_;
- static std::vector<typename RawComponent::Scalar> liquidViscosity_;
- static bool gasViscosityInitialized_;
- static bool liquidViscosityInitialized_;
-
- static std::vector<typename RawComponent::Scalar> gasThermalConductivity_;
- static std::vector<typename RawComponent::Scalar> liquidThermalConductivity_;
- static bool gasThermalConductivityInitialized_;
- static bool liquidThermalConductivityInitialized_;
-
- // 2D fields with the temperature and density as degrees of freedom
- static std::vector<typename RawComponent::Scalar> gasPressure_;
- static std::vector<typename RawComponent::Scalar> liquidPressure_;
- static bool gasPressureInitialized_;
- static bool liquidPressureInitialized_;
-
- // temperature, pressure and density ranges
- static Scalar tempMin_;
- static Scalar tempMax_;
- static std::size_t nTemp_;
-
- static Scalar pressMin_;
- static Scalar pressMax_;
- static std::size_t nPress_;
-
- static Scalar densityMin_;
- static Scalar densityMax_;
- static std::size_t nDensity_;
+ static Table& table()
+ {
+ static Table t;
+ return t;
+ }
};
-
template <class RawComponent, bool useVaporPressure>
bool TabulatedComponent<RawComponent, useVaporPressure>::initialized_ = false;
@@ -1314,88 +1266,6 @@ template <class RawComponent, bool useVaporPressure>
bool TabulatedComponent<RawComponent, useVaporPressure>::warningPrinted_ = false;
#endif
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::minMaxLiquidDensityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::minMaxGasDensityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasEnthalpyInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidEnthalpyInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasHeatCapacityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidHeatCapacityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasDensityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidDensityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasViscosityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidViscosityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasThermalConductivityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidThermalConductivityInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::gasPressureInitialized_ = false;
-template <class RawComponent, bool useVaporPressure>
-bool TabulatedComponent<RawComponent, useVaporPressure>::liquidPressureInitialized_ = false;
-
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::vaporPressure_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::minLiquidDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::maxLiquidDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::minGasDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::maxGasDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasEnthalpy_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidEnthalpy_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasHeatCapacity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidHeatCapacity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidDensity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasViscosity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidViscosity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasThermalConductivity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidThermalConductivity_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::gasPressure_;
-template <class RawComponent, bool useVaporPressure>
-std::vector<typename RawComponent::Scalar> TabulatedComponent<RawComponent, useVaporPressure>::liquidPressure_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::tempMin_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::tempMax_;
-template <class RawComponent, bool useVaporPressure>
-std::size_t TabulatedComponent<RawComponent, useVaporPressure>::nTemp_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::pressMin_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::pressMax_;
-template <class RawComponent, bool useVaporPressure>
-std::size_t TabulatedComponent<RawComponent, useVaporPressure>::nPress_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::densityMin_;
-template <class RawComponent, bool useVaporPressure>
-typename RawComponent::Scalar TabulatedComponent<RawComponent, useVaporPressure>::densityMax_;
-template <class RawComponent, bool useVaporPressure>
-std::size_t TabulatedComponent<RawComponent, useVaporPressure>::nDensity_;
-
// forward declaration
template <class Component>
struct IsAqueous;
--
GitLab