GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/components/tabulatedcomponent.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 257 329 78.1%
Functions: 190 516 36.8%
Branches: 585 929 63.0%

Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Components
10 * \brief Tabulates all thermodynamic properties of a given
11 * untabulated chemical species.
12 *
13 * At the moment, this class can only handle the sub-critical fluids
14 * since it tabulates along the vapor pressure curve.
15 */
16 #ifndef DUMUX_TABULATED_COMPONENT_HH
17 #define DUMUX_TABULATED_COMPONENT_HH
18
19 #include <cmath>
20 #include <limits>
21 #include <cassert>
22 #include <vector>
23 #include <iostream>
24 #include <iomanip>
25 #include <algorithm>
26
27 #include <dune/common/std/type_traits.hh>
28
29 #include <dumux/common/typetraits/typetraits.hh>
30 #include <dumux/common/exceptions.hh>
31 #include <dumux/parallel/multithreading.hh>
32 #include <dumux/parallel/parallel_for.hh>
33 #include <dumux/material/components/componenttraits.hh>
34
35 namespace Dumux::Components {
36 // forward declaration
37 template<class RawComponent, bool useVaporPressure>
38 class TabulatedComponent;
39 } // end namespace Dumux::Components
40
41 namespace Dumux {
42
43 //! component traits for tabulated component
44 template<class RawComponent, bool useVaporPressure>
45 struct ComponentTraits<Components::TabulatedComponent<RawComponent, useVaporPressure>>
46 {
47 using Scalar = typename RawComponent::Scalar;
48
49 //! if the component implements a solid state
50 static constexpr bool hasSolidState = std::is_base_of<Components::Solid<Scalar, RawComponent>, RawComponent>::value;
51
52 //! if the component implements a liquid state
53 static constexpr bool hasLiquidState = std::is_base_of<Components::Liquid<Scalar, RawComponent>, RawComponent>::value;
54
55 //! if the component implements a gaseous state
56 static constexpr bool hasGasState = std::is_base_of<Components::Gas<Scalar, RawComponent>, RawComponent>::value;
57 };
58 } // end namespace Dumux
59
60 namespace Dumux::Components::Detail {
61 struct DisableStaticAssert {};
62 } // end namespace Dumux::Components::Detail
63
64 namespace Dumux {
65 template<> struct AlwaysFalse<Components::Detail::DisableStaticAssert> : public std::true_type {};
66 }// end namespace Dumux
67
68 namespace Dumux::Components::Detail {
69
70 template<class C> using CompHasNoLiquidEnthalpy = decltype(C::template liquidEnthalpy<DisableStaticAssert>(0.0, 0.0));
71 template<class C> using CompHasNoLiquidDensity = decltype(C::template liquidDensity<DisableStaticAssert>(0.0, 0.0));
72 template<class C> using CompHasNoLiquidThermalCond = decltype(C::template liquidThermalConductivity<DisableStaticAssert>(0.0, 0.0));
73 template<class C> using CompHasNoLiquidHeatCapacity = decltype(C::template liquidHeatCapacity<DisableStaticAssert>(0.0, 0.0));
74 template<class C> using CompHasNoLiquidViscosity = decltype(C::template liquidViscosity<DisableStaticAssert>(0.0, 0.0));
75 template<class C> using CompHasLiquidPressure = decltype(C::liquidPressure(0.0, 0.0));
76
77 template<class C> using CompHasNoGasEnthalpy = decltype(C::template gasEnthalpy<DisableStaticAssert>(0.0, 0.0));
78 template<class C> using CompHasNoGasDensity = decltype(C::template gasDensity<DisableStaticAssert>(0.0, 0.0));
79 template<class C> using CompHasNoGasThermalCond = decltype(C::template gasThermalConductivity<DisableStaticAssert>(0.0, 0.0));
80 template<class C> using CompHasNoGasHeatCapacity = decltype(C::template gasHeatCapacity<DisableStaticAssert>(0.0, 0.0));
81 template<class C> using CompHasNoGasViscosity = decltype(C::template gasViscosity<DisableStaticAssert>(0.0, 0.0));
82 template<class C> using CompHasGasPressure = decltype(C::gasPressure(0.0, 0.0));
83
84 template<class C> constexpr inline bool hasLiquidEnthalpy()
85 { return !Dune::Std::is_detected<CompHasNoLiquidEnthalpy, C>::value && ComponentTraits<C>::hasLiquidState; }
86 template<class C> constexpr inline bool hasLiquidDensity()
87 { return !Dune::Std::is_detected<CompHasNoLiquidDensity, C>::value && ComponentTraits<C>::hasLiquidState; }
88 template<class C> constexpr inline bool hasLiquidThermalConductivity()
89 { return !Dune::Std::is_detected<CompHasNoLiquidThermalCond, C>::value && ComponentTraits<C>::hasLiquidState; }
90 template<class C> constexpr inline bool hasLiquidHeatCapacity()
91 { return !Dune::Std::is_detected<CompHasNoLiquidHeatCapacity, C>::value && ComponentTraits<C>::hasLiquidState; }
92 template<class C> constexpr inline bool hasLiquidViscosity()
93 { return !Dune::Std::is_detected<CompHasNoLiquidViscosity, C>::value && ComponentTraits<C>::hasLiquidState; }
94 template<class C> constexpr inline bool hasLiquidPressure()
95 { return Dune::Std::is_detected<CompHasLiquidPressure, C>::value && ComponentTraits<C>::hasLiquidState; }
96
97 template<class C> constexpr inline bool hasGasEnthalpy()
98 { return !Dune::Std::is_detected<CompHasNoGasEnthalpy, C>::value && ComponentTraits<C>::hasGasState; }
99 template<class C> constexpr inline bool hasGasDensity()
100 { return !Dune::Std::is_detected<CompHasNoGasDensity, C>::value && ComponentTraits<C>::hasGasState; }
101 template<class C> constexpr inline bool hasGasThermalConductivity()
102 { return !Dune::Std::is_detected<CompHasNoGasThermalCond, C>::value && ComponentTraits<C>::hasGasState; }
103 template<class C> constexpr inline bool hasGasHeatCapacity()
104 { return !Dune::Std::is_detected<CompHasNoGasHeatCapacity, C>::value && ComponentTraits<C>::hasGasState; }
105 template<class C> constexpr inline bool hasGasViscosity()
106 { return !Dune::Std::is_detected<CompHasNoGasViscosity, C>::value && ComponentTraits<C>::hasGasState; }
107 template<class C> constexpr inline bool hasGasPressure()
108 { return Dune::Std::is_detected<CompHasGasPressure, C>::value && ComponentTraits<C>::hasGasState; }
109
110 template<class RawComponent, bool useVaporPressure = true>
111 class TabulatedComponentTable
112 {
113 using Scalar = typename RawComponent::Scalar;
114 friend class TabulatedComponent<RawComponent, useVaporPressure>;
115
116 struct GasPolicy
117 {
118 172086 Scalar minP(std::size_t iT) const { return table.minGasPressure(iT); }
119 172086 Scalar maxP(std::size_t iT) const { return table.maxGasPressure(iT); }
120 const TabulatedComponentTable<RawComponent, useVaporPressure>& table;
121 };
122
123 struct LiquidPolicy
124 {
125 171504 Scalar minP(std::size_t iT) const { return table.minLiquidPressure(iT); }
126 171504 Scalar maxP(std::size_t iT) const { return table.maxLiquidPressure(iT); }
127 const TabulatedComponentTable<RawComponent, useVaporPressure>& table;
128 };
129 public:
130 201 void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
131 Scalar pressMin, Scalar pressMax, std::size_t nPress)
132 {
133 208 tempMin_ = tempMin;
134 208 tempMax_ = tempMax;
135 208 nTemp_ = nTemp;
136 208 pressMin_ = pressMin;
137 208 pressMax_ = pressMax;
138 208 nPress_ = nPress;
139 208 nDensity_ = nPress;
140
141 // resize & initialize the arrays with NaN
142 assert(std::numeric_limits<Scalar>::has_quiet_NaN);
143 208 const auto NaN = std::numeric_limits<Scalar>::quiet_NaN();
144
145 // initialize vapor pressure array depending on useVaporPressure
146 208 vaporPressure_.resize(nTemp_, NaN);
147 398 tabularizeVaporPressure_();
148
149 if constexpr (ComponentTraits<RawComponent>::hasGasState)
150 {
151 201 minGasDensity_.resize(nTemp_, NaN);
152 201 maxGasDensity_.resize(nTemp_, NaN);
153 201 const std::size_t numEntriesTp = nTemp_*nPress_;
154 201 gasEnthalpy_.resize(numEntriesTp, NaN);
155 201 gasHeatCapacity_.resize(numEntriesTp, NaN);
156 201 gasDensity_.resize(numEntriesTp, NaN);
157 201 gasViscosity_.resize(numEntriesTp, NaN);
158 201 gasThermalConductivity_.resize(numEntriesTp, NaN);
159
160 if constexpr (RawComponent::gasIsCompressible())
161 201 gasPressure_.resize(numEntriesTp, NaN);
162
163 402 minMaxGasDensityInitialized_ = tabularizeMinMaxGasDensity_();
164 392 gasPressureInitialized_ = tabularizeGasPressure_();
165 398 gasEnthalpyInitialized_ = tabularizeGasEnthalpy_();
166 392 gasHeatCapacityInitialized_ = tabularizeGasHeatCapacity_();
167 402 gasDensityInitialized_ = tabularizeGasDensity_();
168 398 gasViscosityInitialized_ = tabularizeGasViscosity_();
169 390 gasThermalConductivityInitialized_ = tabularizeGasThermalConductivity_();
170 }
171
172 if constexpr (ComponentTraits<RawComponent>::hasLiquidState)
173 {
174 190 minLiquidDensity_.resize(nTemp_, NaN);
175 190 maxLiquidDensity_.resize(nTemp_, NaN);
176
177 190 const std::size_t numEntriesTp = nTemp_*nPress_;
178 190 liquidEnthalpy_.resize(numEntriesTp, NaN);
179 190 liquidHeatCapacity_.resize(numEntriesTp, NaN);
180 190 liquidDensity_.resize(numEntriesTp, NaN);
181 190 liquidViscosity_.resize(numEntriesTp, NaN);
182 190 liquidThermalConductivity_.resize(numEntriesTp, NaN);
183
184 if constexpr (RawComponent::liquidIsCompressible())
185 180 liquidPressure_.resize(numEntriesTp, NaN);
186
187 376 minMaxLiquidDensityInitialized_ = tabularizeMinMaxLiquidDensity_();
188 368 liquidPressureInitialized_ = tabularizeLiquidPressure_();
189 376 liquidEnthalpyInitialized_ = tabularizeLiquidEnthalpy_();
190 376 liquidHeatCapacityInitialized_ = tabularizeLiquidHeatCapacity_();
191 380 liquidDensityInitialized_ = tabularizeLiquidDensity_();
192 380 liquidViscosityInitialized_ = tabularizeLiquidViscosity_();
193 376 liquidThermalConductivityInitialized_ = tabularizeLiquidThermalConductivity_();
194 }
195 201 }
196
197 //! returns the index of an entry in a temperature field
198 inline Scalar tempIdx(Scalar temperature) const
199 {
200 1232339122 return (nTemp_ - 1)*(temperature - tempMin_)/(tempMax_ - tempMin_);
201 }
202
203 //! returns the index of an entry in a pressure field
204 1770286592 inline Scalar pressLiquidIdx(Scalar pressure, std::size_t tempIdx) const
205 {
206 3554572648 const Scalar plMin = minLiquidPressure(tempIdx);
207 3554572648 const Scalar plMax = maxLiquidPressure(tempIdx);
208 1777286324 return (nPress_ - 1)*(pressure - plMin)/(plMax - plMin);
209 }
210
211 //! returns the index of an entry in a temperature field
212 591696468 inline Scalar pressGasIdx(Scalar pressure, std::size_t tempIdx) const
213 {
214 1200922640 const Scalar pgMin = minGasPressure(tempIdx);
215 1200922640 const Scalar pgMax = maxGasPressure(tempIdx);
216 600461320 return (nPress_ - 1)*(pressure - pgMin)/(pgMax - pgMin);
217 }
218
219 //! returns the index of an entry in a density field
220 inline Scalar densityLiquidIdx(Scalar density, std::size_t tempIdx) const
221 {
222 4941532 const Scalar densityMin = minLiquidDensity_[tempIdx];
223
4/4
✓ Branch 0 taken 1228521 times.
✓ Branch 1 taken 6862 times.
✓ Branch 2 taken 1228521 times.
✓ Branch 3 taken 6862 times.
2470766 const Scalar densityMax = maxLiquidDensity_[tempIdx];
224 2470766 return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
225 }
226
227 //! returns the index of an entry in a density field
228 inline Scalar densityGasIdx(Scalar density, std::size_t tempIdx) const
229 {
230 820092 const Scalar densityMin = minGasDensity_[tempIdx];
231
2/4
✓ Branch 0 taken 205023 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 205023 times.
✗ Branch 3 not taken.
410046 const Scalar densityMax = maxGasDensity_[tempIdx];
232 410046 return (nDensity_ - 1) * (density - densityMin)/(densityMax - densityMin);
233 }
234
235 //! returns the minimum tabularized liquid pressure at a given temperature index
236 inline Scalar minLiquidPressure(int tempIdx) const
237 {
238 using std::max;
239 if (!useVaporPressure)
240 return pressMin_;
241 else
242
19/34
✓ Branch 0 taken 152979022 times.
✓ Branch 1 taken 1850688 times.
✓ Branch 2 taken 152979022 times.
✓ Branch 3 taken 1850688 times.
✓ Branch 4 taken 7701674 times.
✓ Branch 5 taken 1332 times.
✓ Branch 6 taken 378354116 times.
✓ Branch 7 taken 1224748936 times.
✓ Branch 8 taken 370661532 times.
✓ Branch 9 taken 1224748904 times.
✓ Branch 10 taken 3009 times.
✓ Branch 11 taken 76 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✓ Branch 26 taken 10412965 times.
✓ Branch 27 taken 1940865 times.
✓ Branch 28 taken 10412965 times.
✓ Branch 29 taken 1940865 times.
✓ Branch 30 taken 692 times.
✓ Branch 31 taken 111 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
3540587465 return max(pressMin_, vaporPressure_[tempIdx] / 1.1);
243 }
244
245 //! returns the maximum tabularized liquid pressure at a given temperature index
246 inline Scalar maxLiquidPressure(int tempIdx) const
247 {
248 using std::max;
249 if (!useVaporPressure)
250 return pressMax_;
251 else
252
31/58
✗ Branch 0 not taken.
✓ Branch 1 taken 154829710 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 43 times.
✓ Branch 5 taken 7702963 times.
✓ Branch 6 taken 43 times.
✓ Branch 7 taken 10347 times.
✓ Branch 8 taken 43 times.
✓ Branch 9 taken 10950 times.
✓ Branch 10 taken 43 times.
✓ Branch 11 taken 13432 times.
✓ Branch 12 taken 43 times.
✓ Branch 13 taken 13432 times.
✓ Branch 14 taken 45 times.
✓ Branch 15 taken 13430 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1595413521 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3085 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2494 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 3 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 3 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 3 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 3 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 3 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 12353830 times.
✗ Branch 42 not taken.
✓ Branch 43 taken 803 times.
✗ Branch 44 not taken.
✓ Branch 45 taken 803 times.
✗ Branch 46 not taken.
✓ Branch 47 taken 803 times.
✗ Branch 48 not taken.
✓ Branch 49 taken 803 times.
✗ Branch 50 not taken.
✓ Branch 51 taken 803 times.
✗ Branch 52 not taken.
✓ Branch 53 taken 803 times.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
1770372296 return max(pressMax_, vaporPressure_[tempIdx] * 1.1);
253 }
254
255 //! returns the minimum tabularized gas pressure at a given temperature index
256 inline Scalar minGasPressure(int tempIdx) const
257 {
258 using std::min;
259 if (!useVaporPressure)
260 return pressMin_;
261 else
262
22/54
✗ Branch 0 not taken.
✓ Branch 1 taken 39503696 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 39503696 times.
✓ Branch 4 taken 1336 times.
✓ Branch 5 taken 8130 times.
✓ Branch 6 taken 119843074 times.
✓ Branch 7 taken 430299468 times.
✓ Branch 8 taken 119843114 times.
✓ Branch 9 taken 430301334 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2103 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 3 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 3 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 3 times.
✓ Branch 42 taken 3 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 3 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 145 times.
✓ Branch 47 taken 2050085 times.
✓ Branch 48 taken 145 times.
✓ Branch 49 taken 2050085 times.
✓ Branch 50 taken 111 times.
✓ Branch 51 taken 692 times.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
1183407235 return min(pressMin_, vaporPressure_[tempIdx] / 1.1 );
263 }
264
265 //! returns the maximum tabularized gas pressure at a given temperature index
266 inline Scalar maxGasPressure(int tempIdx) const
267 {
268 using std::min;
269 if (!useVaporPressure)
270 return pressMax_;
271 else
272
47/96
✓ Branch 0 taken 39503696 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9424 times.
✓ Branch 5 taken 43 times.
✓ Branch 6 taken 11329 times.
✓ Branch 7 taken 43 times.
✓ Branch 8 taken 11629 times.
✓ Branch 9 taken 43 times.
✓ Branch 10 taken 13432 times.
✓ Branch 11 taken 43 times.
✓ Branch 12 taken 13432 times.
✓ Branch 13 taken 43 times.
✓ Branch 14 taken 13430 times.
✓ Branch 15 taken 45 times.
✓ Branch 16 taken 549947895 times.
✓ Branch 17 taken 198656 times.
✓ Branch 18 taken 2103 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1803 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 3 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 3 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 3 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 3 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 3 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 3 times.
✗ Branch 47 not taken.
✓ Branch 48 taken 3 times.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✓ Branch 51 taken 3 times.
✗ Branch 52 not taken.
✓ Branch 53 taken 3 times.
✗ Branch 54 not taken.
✓ Branch 55 taken 3 times.
✗ Branch 56 not taken.
✓ Branch 57 taken 3 times.
✗ Branch 58 not taken.
✓ Branch 59 taken 3 times.
✗ Branch 60 not taken.
✓ Branch 61 taken 3 times.
✓ Branch 62 taken 3 times.
✗ Branch 63 not taken.
✓ Branch 64 taken 3 times.
✗ Branch 65 not taken.
✓ Branch 66 taken 3 times.
✗ Branch 67 not taken.
✓ Branch 68 taken 3 times.
✗ Branch 69 not taken.
✓ Branch 70 taken 3 times.
✗ Branch 71 not taken.
✓ Branch 72 taken 3 times.
✗ Branch 73 not taken.
✓ Branch 74 taken 3 times.
✗ Branch 75 not taken.
✓ Branch 76 taken 3 times.
✗ Branch 77 not taken.
✓ Branch 78 taken 2050230 times.
✗ Branch 79 not taken.
✓ Branch 80 taken 803 times.
✗ Branch 81 not taken.
✓ Branch 82 taken 803 times.
✗ Branch 83 not taken.
✓ Branch 84 taken 803 times.
✗ Branch 85 not taken.
✓ Branch 86 taken 803 times.
✗ Branch 87 not taken.
✓ Branch 88 taken 803 times.
✗ Branch 89 not taken.
✓ Branch 90 taken 803 times.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
591782206 return min(pressMax_, vaporPressure_[tempIdx] * 1.1);
273 }
274
275 //! returns the minimum tabularized liquid density at a given temperature index
276 inline Scalar minLiquidDensity(int tempIdx) const
277 { return minLiquidDensity_[tempIdx]; }
278
279 //! returns the maximum tabularized liquid density at a given temperature index
280 inline Scalar maxLiquidDensity(int tempIdx) const
281 { return maxLiquidDensity_[tempIdx]; }
282
283 //! returns the minimum tabularized gas density at a given temperature index
284 inline Scalar minGasDensity(int tempIdx) const
285 { return minGasDensity_[tempIdx]; }
286
287 //! returns the maximum tabularized gas density at a given temperature index
288 inline Scalar maxGasDensity(int tempIdx) const
289 { return maxGasDensity_[tempIdx]; }
290
291 inline std::size_t nTemp() const { return nTemp_; }
292 inline std::size_t nPress() const { return nPress_; }
293 inline std::size_t nDensity() const { return nDensity_; }
294
295 inline Scalar tempMax() const { return tempMax_; }
296 inline Scalar tempMin() const { return tempMin_; }
297
298 private:
299
300 //! initializes vapor pressure if useVaporPressure = true
301 template< bool useVP = useVaporPressure, std::enable_if_t<useVP, int> = 0 >
302 void tabularizeVaporPressure_()
303 {
304 // fill the temperature-pressure arrays
305 29008 Dumux::parallelFor(nTemp_, [=](std::size_t iT)
306 {
307 28592 Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
308 14287 vaporPressure_[iT] = RawComponent::vaporPressure(temperature);
309 });
310 }
311
312 //! if !useVaporPressure, do nothing here
313 template< bool useVP = useVaporPressure, std::enable_if_t<!useVP, int> = 0 >
314 void tabularizeVaporPressure_() {}
315
316 /*!
317 * \brief Initializes property values as function of temperature and pressure.
318 *
319 * \tparam PropFunc Function to evaluate the property prop(T, p)
320 * \tparam MinPFunc Function to evaluate the minimum pressure for a
321 * temperature index (depends on useVaporPressure)
322 * \tparam MaxPFunc Function to evaluate the maximum pressure for a
323 * temperature index (depends on useVaporPressure)
324 *
325 * \param f property function
326 * \param minP function to evaluate minimum pressure for temp idx
327 * \param maxP function to evaluate maximum pressure for temp idx
328 * \param values container to store property values
329 */
330 template<class PropFunc, class Policy>
331 void tabularizeTPArray_(const PropFunc& f, Policy policy, std::vector<typename RawComponent::Scalar>& values) const
332 {
333 56783556 Dumux::parallelFor(nTemp_, [=,&values](std::size_t iT)
334 {
335 286752 Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
336
337
43/66
✓ Branch 0 taken 43 times.
✓ Branch 1 taken 13435 times.
✓ Branch 2 taken 43 times.
✓ Branch 3 taken 13435 times.
✓ Branch 4 taken 43 times.
✓ Branch 5 taken 13435 times.
✓ Branch 6 taken 46 times.
✓ Branch 7 taken 13432 times.
✓ Branch 8 taken 46 times.
✓ Branch 9 taken 13432 times.
✓ Branch 10 taken 13435 times.
✓ Branch 11 taken 43 times.
✓ Branch 12 taken 13435 times.
✓ Branch 13 taken 43 times.
✓ Branch 14 taken 13432 times.
✓ Branch 15 taken 46 times.
✓ Branch 16 taken 13435 times.
✓ Branch 17 taken 43 times.
✓ Branch 18 taken 13435 times.
✓ Branch 19 taken 43 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 3 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✓ Branch 32 taken 3 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 3 times.
✓ Branch 40 taken 3 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 3 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 3 times.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✓ Branch 47 taken 803 times.
✗ Branch 48 not taken.
✓ Branch 49 taken 803 times.
✗ Branch 50 not taken.
✓ Branch 51 taken 803 times.
✗ Branch 52 not taken.
✓ Branch 53 taken 803 times.
✗ Branch 54 not taken.
✓ Branch 55 taken 803 times.
✓ Branch 56 taken 803 times.
✗ Branch 57 not taken.
✓ Branch 58 taken 803 times.
✗ Branch 59 not taken.
✓ Branch 60 taken 803 times.
✗ Branch 61 not taken.
✓ Branch 62 taken 803 times.
✗ Branch 63 not taken.
✓ Branch 64 taken 803 times.
✗ Branch 65 not taken.
143777 Scalar pMax = policy.maxP(iT);
338
53/66
✓ Branch 0 taken 12099 times.
✓ Branch 1 taken 1379 times.
✓ Branch 2 taken 12099 times.
✓ Branch 3 taken 1379 times.
✓ Branch 4 taken 12099 times.
✓ Branch 5 taken 1379 times.
✓ Branch 6 taken 12102 times.
✓ Branch 7 taken 1376 times.
✓ Branch 8 taken 12102 times.
✓ Branch 9 taken 1376 times.
✓ Branch 10 taken 1379 times.
✓ Branch 11 taken 12099 times.
✓ Branch 12 taken 1379 times.
✓ Branch 13 taken 12099 times.
✓ Branch 14 taken 1376 times.
✓ Branch 15 taken 12102 times.
✓ Branch 16 taken 1379 times.
✓ Branch 17 taken 12099 times.
✓ Branch 18 taken 1379 times.
✓ Branch 19 taken 12099 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 3 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✓ Branch 32 taken 3 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 3 times.
✓ Branch 40 taken 3 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 3 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 3 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 692 times.
✓ Branch 47 taken 111 times.
✓ Branch 48 taken 692 times.
✓ Branch 49 taken 111 times.
✓ Branch 50 taken 692 times.
✓ Branch 51 taken 111 times.
✓ Branch 52 taken 692 times.
✓ Branch 53 taken 111 times.
✓ Branch 54 taken 692 times.
✓ Branch 55 taken 111 times.
✓ Branch 56 taken 111 times.
✓ Branch 57 taken 692 times.
✓ Branch 58 taken 111 times.
✓ Branch 59 taken 692 times.
✓ Branch 60 taken 111 times.
✓ Branch 61 taken 692 times.
✓ Branch 62 taken 111 times.
✓ Branch 63 taken 692 times.
✓ Branch 64 taken 111 times.
✓ Branch 65 taken 692 times.
143705 Scalar pMin = policy.minP(iT);
339
132/132
✓ Branch 0 taken 2665209 times.
✓ Branch 1 taken 13580 times.
✓ Branch 2 taken 2665209 times.
✓ Branch 3 taken 13580 times.
✓ Branch 4 taken 2665209 times.
✓ Branch 5 taken 13580 times.
✓ Branch 6 taken 2665209 times.
✓ Branch 7 taken 13580 times.
✓ Branch 8 taken 2645209 times.
✓ Branch 9 taken 13480 times.
✓ Branch 10 taken 2645209 times.
✓ Branch 11 taken 13480 times.
✓ Branch 12 taken 2645209 times.
✓ Branch 13 taken 13480 times.
✓ Branch 14 taken 2645209 times.
✓ Branch 15 taken 13480 times.
✓ Branch 16 taken 2645209 times.
✓ Branch 17 taken 13480 times.
✓ Branch 18 taken 2645209 times.
✓ Branch 19 taken 13480 times.
✓ Branch 20 taken 409 times.
✓ Branch 21 taken 5 times.
✓ Branch 22 taken 409 times.
✓ Branch 23 taken 5 times.
✓ Branch 24 taken 409 times.
✓ Branch 25 taken 5 times.
✓ Branch 26 taken 409 times.
✓ Branch 27 taken 5 times.
✓ Branch 28 taken 9 times.
✓ Branch 29 taken 3 times.
✓ Branch 30 taken 9 times.
✓ Branch 31 taken 3 times.
✓ Branch 32 taken 9 times.
✓ Branch 33 taken 3 times.
✓ Branch 34 taken 9 times.
✓ Branch 35 taken 3 times.
✓ Branch 36 taken 9 times.
✓ Branch 37 taken 3 times.
✓ Branch 38 taken 9 times.
✓ Branch 39 taken 3 times.
✓ Branch 40 taken 9 times.
✓ Branch 41 taken 3 times.
✓ Branch 42 taken 9 times.
✓ Branch 43 taken 3 times.
✓ Branch 44 taken 9 times.
✓ Branch 45 taken 3 times.
✓ Branch 46 taken 9 times.
✓ Branch 47 taken 3 times.
✓ Branch 48 taken 9 times.
✓ Branch 49 taken 3 times.
✓ Branch 50 taken 9 times.
✓ Branch 51 taken 3 times.
✓ Branch 52 taken 9 times.
✓ Branch 53 taken 3 times.
✓ Branch 54 taken 9 times.
✓ Branch 55 taken 3 times.
✓ Branch 56 taken 9 times.
✓ Branch 57 taken 3 times.
✓ Branch 58 taken 9 times.
✓ Branch 59 taken 3 times.
✓ Branch 60 taken 9 times.
✓ Branch 61 taken 3 times.
✓ Branch 62 taken 9 times.
✓ Branch 63 taken 3 times.
✓ Branch 64 taken 9 times.
✓ Branch 65 taken 3 times.
✓ Branch 66 taken 9 times.
✓ Branch 67 taken 3 times.
✓ Branch 68 taken 9 times.
✓ Branch 69 taken 3 times.
✓ Branch 70 taken 9 times.
✓ Branch 71 taken 3 times.
✓ Branch 72 taken 9 times.
✓ Branch 73 taken 3 times.
✓ Branch 74 taken 9 times.
✓ Branch 75 taken 3 times.
✓ Branch 76 taken 9 times.
✓ Branch 77 taken 3 times.
✓ Branch 78 taken 9 times.
✓ Branch 79 taken 3 times.
✓ Branch 80 taken 9 times.
✓ Branch 81 taken 3 times.
✓ Branch 82 taken 9 times.
✓ Branch 83 taken 3 times.
✓ Branch 84 taken 9 times.
✓ Branch 85 taken 3 times.
✓ Branch 86 taken 9 times.
✓ Branch 87 taken 3 times.
✓ Branch 88 taken 9 times.
✓ Branch 89 taken 3 times.
✓ Branch 90 taken 9 times.
✓ Branch 91 taken 3 times.
✓ Branch 92 taken 9 times.
✓ Branch 93 taken 3 times.
✓ Branch 94 taken 9 times.
✓ Branch 95 taken 3 times.
✓ Branch 96 taken 9 times.
✓ Branch 97 taken 3 times.
✓ Branch 98 taken 9 times.
✓ Branch 99 taken 3 times.
✓ Branch 100 taken 9 times.
✓ Branch 101 taken 3 times.
✓ Branch 102 taken 9 times.
✓ Branch 103 taken 3 times.
✓ Branch 104 taken 9 times.
✓ Branch 105 taken 3 times.
✓ Branch 106 taken 9 times.
✓ Branch 107 taken 3 times.
✓ Branch 108 taken 9 times.
✓ Branch 109 taken 3 times.
✓ Branch 110 taken 9 times.
✓ Branch 111 taken 3 times.
✓ Branch 112 taken 160009 times.
✓ Branch 113 taken 803 times.
✓ Branch 114 taken 160009 times.
✓ Branch 115 taken 803 times.
✓ Branch 116 taken 160009 times.
✓ Branch 117 taken 803 times.
✓ Branch 118 taken 160009 times.
✓ Branch 119 taken 803 times.
✓ Branch 120 taken 160009 times.
✓ Branch 121 taken 803 times.
✓ Branch 122 taken 160009 times.
✓ Branch 123 taken 803 times.
✓ Branch 124 taken 160009 times.
✓ Branch 125 taken 803 times.
✓ Branch 126 taken 160009 times.
✓ Branch 127 taken 803 times.
✓ Branch 128 taken 160009 times.
✓ Branch 129 taken 803 times.
✓ Branch 130 taken 160009 times.
✓ Branch 131 taken 803 times.
28277570 for (std::size_t iP = 0; iP < nPress_; ++ iP)
340 {
341 28134194 Scalar pressure = iP * (pMax - pMin)/(nPress_ - 1) + pMin;
342 56206801 values[iT + iP*nTemp_] = f(temperature, pressure);
343 }
344 });
345 }
346
347 /*!
348 * \brief Initializes the minimum/maximum densities on the temperature range.
349 *
350 * \tparam RhoFunc Function to evaluate the density rho(T, p)
351 * \tparam MinPFunc Function to evaluate the minimum pressure for a
352 * temperature index (depends on useVaporPressure)
353 * \tparam MaxPFunc Function to evaluate the maximum pressure for a
354 * temperature index (depends on useVaporPressure)
355 *
356 * \param rho density function
357 * \param minP function to evaluate minimum pressure for temp idx
358 * \param maxP function to evaluate maximum pressure for temp idx
359 * \param rhoMin container to store minimum density values
360 * \param rhoMax container to store maximum density values
361 */
362 template<class RhoFunc, class Policy>
363 void tabularizeMinMaxRhoArray_(const RhoFunc& rho, Policy policy,
364 std::vector<typename RawComponent::Scalar>& rhoMin,
365 std::vector<typename RawComponent::Scalar>& rhoMax) const
366 {
367 172316 Dumux::parallelFor(nTemp_, [=,&rhoMin,&rhoMax](std::size_t iT)
368 {
369 57408 Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
370
371 using std::min;
372
21/27
✓ Branch 0 taken 12099 times.
✓ Branch 1 taken 1379 times.
✓ Branch 2 taken 172 times.
✓ Branch 3 taken 13306 times.
✓ Branch 4 taken 1376 times.
✓ Branch 5 taken 12100 times.
✓ Branch 6 taken 174 times.
✓ Branch 7 taken 13306 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 692 times.
✓ Branch 21 taken 111 times.
✓ Branch 22 taken 2 times.
✓ Branch 23 taken 801 times.
✓ Branch 24 taken 111 times.
✓ Branch 25 taken 692 times.
✓ Branch 26 taken 2 times.
✓ Branch 27 taken 801 times.
57649 rhoMin[iT] = rho(temperature, policy.minP(iT));
373
38/47
✓ Branch 0 taken 173 times.
✓ Branch 1 taken 13305 times.
✓ Branch 2 taken 45 times.
✓ Branch 3 taken 13431 times.
✓ Branch 4 taken 174 times.
✓ Branch 5 taken 13306 times.
✓ Branch 6 taken 13430 times.
✓ Branch 7 taken 46 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✓ Branch 21 taken 1 times.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 3 times.
✓ Branch 25 taken 1 times.
✓ Branch 26 taken 2 times.
✓ Branch 28 taken 1 times.
✓ Branch 29 taken 2 times.
✓ Branch 30 taken 3 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 2 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 3 times.
✓ Branch 36 taken 1 times.
✓ Branch 37 taken 2 times.
✓ Branch 38 taken 3 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 2 times.
✓ Branch 41 taken 801 times.
✗ Branch 42 not taken.
✓ Branch 43 taken 803 times.
✓ Branch 44 taken 2 times.
✓ Branch 45 taken 801 times.
✓ Branch 46 taken 803 times.
✗ Branch 47 not taken.
29481 rhoMax[iT] = rho(temperature, policy.maxP(min(iT + 1, nTemp_ - 1)));
374 });
375 }
376
377 /*!
378 * \brief Initializes pressure arrays as function of temperature and density.
379 *
380 * \tparam PFunc Function to evaluate the pressure p(T, rho)
381 *
382 * \param pressure container to store pressure values
383 * \param p pressure function p(T, rho)
384 * \param rhoMin container with minimum density values
385 * \param rhoMax container with maximum density values
386 */
387 template<class PFunc>
388 void tabularizePressureArray_(std::vector<typename RawComponent::Scalar>& pressure,
389 const PFunc& p,
390 const std::vector<typename RawComponent::Scalar>& rhoMin,
391 const std::vector<typename RawComponent::Scalar>& rhoMax) const
392 {
393 22558010 Dumux::parallelFor(nTemp_, [=,&pressure,&rhoMin,&rhoMax](std::size_t iT)
394 {
395 57354 Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
396
397
14/14
✓ Branch 0 taken 2665209 times.
✓ Branch 1 taken 13580 times.
✓ Branch 2 taken 2645209 times.
✓ Branch 3 taken 13480 times.
✓ Branch 4 taken 409 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 160009 times.
✓ Branch 11 taken 803 times.
✓ Branch 12 taken 160009 times.
✓ Branch 13 taken 803 times.
5659540 for (std::size_t iRho = 0; iRho < nDensity_; ++ iRho)
398 {
399 11261726 Scalar density = Scalar(iRho)/(nDensity_ - 1)
400 11261726 * (rhoMax[iT] - rhoMin[iT])
401 5630863 + rhoMin[iT];
402 11241281 pressure[iT + iRho*nTemp_] = p(temperature, density);
403 }
404 });
405 }
406
407 template<class RC = RawComponent>
408 bool tabularizeGasEnthalpy_()
409 {
410 if constexpr (Detail::hasGasEnthalpy<RC>())
411 {
412 2846126 const auto gasEnth = [] (auto T, auto p) { return RC::gasEnthalpy(T, p); };
413 370 tabularizeTPArray_(gasEnth, GasPolicy{ *this }, gasEnthalpy_);
414 return true;
415 }
416
417 return false;
418 }
419
420 template<class RC = RawComponent>
421 bool tabularizeLiquidEnthalpy_()
422 {
423 if constexpr (Detail::hasLiquidEnthalpy<RC>())
424 {
425 2805254 const auto liqEnth = [] (auto T, auto p) { return RC::liquidEnthalpy(T, p); };
426 358 tabularizeTPArray_(liqEnth, LiquidPolicy{ *this }, liquidEnthalpy_);
427 return true;
428 }
429
430 return false;
431 }
432
433 template<class RC = RawComponent>
434 bool tabularizeGasHeatCapacity_()
435 {
436 if constexpr (Detail::hasGasHeatCapacity<RC>())
437 {
438 2825681 const auto gasHC = [] (auto T, auto p) { return RC::gasHeatCapacity(T, p); };
439 182 tabularizeTPArray_(gasHC, GasPolicy{ *this }, gasHeatCapacity_);
440 return true;
441 }
442
443 return false;
444 }
445
446 template<class RC = RawComponent>
447 bool tabularizeLiquidHeatCapacity_()
448 {
449 if constexpr (Detail::hasLiquidHeatCapacity<RC>())
450 {
451 2805245 const auto liqHC = [] (auto T, auto p) { return RC::liquidHeatCapacity(T, p); };
452 179 tabularizeTPArray_(liqHC, LiquidPolicy{ *this }, liquidHeatCapacity_);
453 return true;
454 }
455
456 return false;
457 }
458
459 template<class RC = RawComponent>
460 bool tabularizeMinMaxGasDensity_()
461 {
462 if constexpr (Detail::hasGasDensity<RC>())
463 {
464
28/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 1 times.
✓ Branch 19 taken 2 times.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 2 times.
✓ Branch 22 taken 1 times.
✓ Branch 23 taken 2 times.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 2 times.
✓ Branch 26 taken 1 times.
✓ Branch 27 taken 2 times.
29000 const auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
465 187 tabularizeMinMaxRhoArray_(gasRho, GasPolicy{ *this }, minGasDensity_, maxGasDensity_);
466 return true;
467 }
468
469 return false;
470 }
471
472 template<class RC = RawComponent>
473 bool tabularizeMinMaxLiquidDensity_()
474 {
475 if constexpr (Detail::hasGasEnthalpy<RC>())
476 {
477
10/10
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 2 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 2 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 2 times.
28575 const auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); };
478 179 tabularizeMinMaxRhoArray_(liqRho, LiquidPolicy{ *this }, minLiquidDensity_, maxLiquidDensity_);
479 return true;
480 }
481
482 return false;
483 }
484
485 template<class RC = RawComponent>
486 bool tabularizeGasPressure_()
487 {
488 // pressure is only defined if the gas is compressible (this is usually the case)
489 if constexpr (Detail::hasGasPressure<RC>() && RC::gasIsCompressible())
490 {
491 2846099 const auto gasPFunc = [] (auto T, auto rho) { return RC::gasPressure(T, rho); };
492 182 tabularizePressureArray_(gasPressure_, gasPFunc, minGasDensity_, maxGasDensity_);
493 return true;
494 }
495
496 return false;
497 }
498
499 template<class RC = RawComponent>
500 bool tabularizeLiquidPressure_()
501 {
502 // pressure is only defined if the liquid is compressible (this is often not the case)
503 if constexpr (Detail::hasLiquidPressure<RC>() && RC::liquidIsCompressible())
504 {
505 2805209 const auto liqPFunc = [] (auto T, auto rho) { return RC::liquidPressure(T, rho); };
506 175 tabularizePressureArray_(liquidPressure_, liqPFunc, minLiquidDensity_, maxLiquidDensity_);
507 return true;
508 }
509
510 return false;
511 }
512
513 template<class RC = RawComponent>
514 bool tabularizeGasDensity_()
515 {
516 if constexpr (Detail::hasGasDensity<RC>())
517 {
518 2846189 const auto gasRho = [] (auto T, auto p) { return RC::gasDensity(T, p); };
519 187 tabularizeTPArray_(gasRho, GasPolicy{ *this }, gasDensity_);
520 return true;
521 }
522
523 return false;
524 }
525
526 template<class RC = RawComponent>
527 bool tabularizeLiquidDensity_()
528 {
529 if constexpr (Detail::hasLiquidDensity<RC>())
530 {
531 // TODO: we could get rid of the lambdas and pass the functor directly. But,
532 // currently Brine is a component (and not a fluid system) expecting a
533 // third argument with a default, which cannot be wrapped in a function pointer.
534 // For this reason we have to wrap this into a lambda here.
535 2805263 const auto liqRho = [] (auto T, auto p) { return RC::liquidDensity(T, p); };
536 181 tabularizeTPArray_(liqRho, LiquidPolicy{ *this }, liquidDensity_);
537 return true;
538 }
539
540 return false;
541 }
542
543 template<class RC = RawComponent>
544 bool tabularizeGasViscosity_()
545 {
546 if constexpr (Detail::hasGasViscosity<RC>())
547 {
548 2825681 const auto gasVisc = [] (auto T, auto p) { return RC::gasViscosity(T, p); };
549 185 tabularizeTPArray_(gasVisc, GasPolicy{ *this }, gasViscosity_);
550 return true;
551 }
552
553 return false;
554 }
555
556 template<class RC = RawComponent>
557 bool tabularizeLiquidViscosity_()
558 {
559 if constexpr (Detail::hasLiquidViscosity<RC>())
560 {
561 2805263 const auto liqVisc = [] (auto T, auto p) { return RC::liquidViscosity(T, p); };
562 181 tabularizeTPArray_(liqVisc, LiquidPolicy{ *this }, liquidViscosity_);
563 return true;
564 }
565
566 return false;
567 }
568
569 template<class RC = RawComponent>
570 bool tabularizeGasThermalConductivity_()
571 {
572 if constexpr (Detail::hasGasThermalConductivity<RC>())
573 {
574 2805254 const auto gasTC = [] (auto T, auto p) { return RC::gasThermalConductivity(T, p); };
575 181 tabularizeTPArray_(gasTC, GasPolicy{ *this }, gasThermalConductivity_);
576 return true;
577 }
578
579 return false;
580 }
581
582 template<class RC = RawComponent>
583 bool tabularizeLiquidThermalConductivity_()
584 {
585 if constexpr (Detail::hasLiquidThermalConductivity<RC>())
586 {
587 2805245 const auto liqTC = [] (auto T, auto p) { return RC::liquidThermalConductivity(T, p); };
588 179 tabularizeTPArray_(liqTC, LiquidPolicy{ *this }, liquidThermalConductivity_);
589 return true;
590 }
591
592 return false;
593 }
594
595 private:
596 // 1D fields with the temperature as degree of freedom
597 std::vector<Scalar> vaporPressure_;
598
599 std::vector<Scalar> minLiquidDensity_;
600 std::vector<Scalar> maxLiquidDensity_;
601 bool minMaxLiquidDensityInitialized_;
602
603 std::vector<Scalar> minGasDensity_;
604 std::vector<Scalar> maxGasDensity_;
605 bool minMaxGasDensityInitialized_;
606
607 // 2D fields with the temperature and pressure as degrees of freedom
608 std::vector<Scalar> gasEnthalpy_;
609 std::vector<Scalar> liquidEnthalpy_;
610 bool gasEnthalpyInitialized_;
611 bool liquidEnthalpyInitialized_;
612
613 std::vector<Scalar> gasHeatCapacity_;
614 std::vector<Scalar> liquidHeatCapacity_;
615 bool gasHeatCapacityInitialized_;
616 bool liquidHeatCapacityInitialized_;
617
618 std::vector<Scalar> gasDensity_;
619 std::vector<Scalar> liquidDensity_;
620 bool gasDensityInitialized_;
621 bool liquidDensityInitialized_;
622
623 std::vector<Scalar> gasViscosity_;
624 std::vector<Scalar> liquidViscosity_;
625 bool gasViscosityInitialized_;
626 bool liquidViscosityInitialized_;
627
628 std::vector<Scalar> gasThermalConductivity_;
629 std::vector<Scalar> liquidThermalConductivity_;
630 bool gasThermalConductivityInitialized_;
631 bool liquidThermalConductivityInitialized_;
632
633 // 2D fields with the temperature and density as degrees of freedom
634 std::vector<Scalar> gasPressure_;
635 std::vector<Scalar> liquidPressure_;
636 bool gasPressureInitialized_;
637 bool liquidPressureInitialized_;
638
639 // temperature, pressure and density ranges
640 Scalar tempMin_;
641 Scalar tempMax_;
642 std::size_t nTemp_;
643
644 Scalar pressMin_;
645 Scalar pressMax_;
646 std::size_t nPress_;
647
648 Scalar densityMin_;
649 Scalar densityMax_;
650 std::size_t nDensity_;
651 };
652
653 } // end namespace Dumux::Components::Detail
654
655 namespace Dumux::Components {
656
657 /*!
658 * \ingroup Components
659 * \brief Tabulates all thermodynamic properties of a given component
660 *
661 * At the moment, this class can only handle the sub-critical fluids
662 * since it tabulates along the vapor pressure curve.
663 *
664 * \tparam Scalar The type used for scalar values
665 * \tparam RawComponent The component which ought to be tabulated
666 * \tparam useVaporPressure If set to true, the min/max pressure
667 * values for gas&liquid phase will be set
668 * depending on the vapor pressure.
669 */
670 template <class RawComponent, bool useVaporPressure=true>
671 class TabulatedComponent
672 {
673 using ThisType = TabulatedComponent<RawComponent, useVaporPressure>;
674 using Table = Detail::TabulatedComponentTable<RawComponent, useVaporPressure>;
675
676 struct InterpolatePolicy
677 {
678 using Scalar = typename RawComponent::Scalar;
679
680 const Table& table;
681
682 2464678244 Scalar tempIdx(Scalar T) const { return table.tempIdx(T); }
683 13 Scalar tempMin() const { return table.tempMin(); }
684 13 Scalar tempMax() const { return table.tempMax(); }
685 1232339135 std::size_t nTemp() const { return table.nTemp(); }
686 1199436654 std::size_t nPress() const { return table.nPress(); }
687 1440406 std::size_t nDensity() const { return table.nDensity(); }
688 };
689
690 struct InterpolateGasPolicy : public InterpolatePolicy
691 {
692 using Scalar = typename RawComponent::Scalar;
693 626755876 Scalar pressIdx(Scalar p, std::size_t tempIdx) const { return this->table.pressGasIdx(p, tempIdx); }
694 820092 Scalar rhoIdx(Scalar rho, std::size_t tempIdx) const { return this->table.densityGasIdx(rho, tempIdx); }
695 Scalar minP(int tempIdx) const { return this->table.minGasPressure(tempIdx); }
696 Scalar maxP(int tempIdx) const { return this->table.maxGasPressure(tempIdx); }
697 Scalar minRho(int tempIdx) const { return this->table.minGasDensity(tempIdx); }
698 Scalar maxRho(int tempIdx) const { return this->table.maxGasDensity(tempIdx); }
699 };
700
701 struct InterpolateLiquidPolicy : public InterpolatePolicy
702 {
703 using Scalar = typename RawComponent::Scalar;
704 1798285520 Scalar pressIdx(Scalar p, std::size_t tempIdx) const { return this->table.pressLiquidIdx(p, tempIdx); }
705 4941532 Scalar rhoIdx(Scalar rho, std::size_t tempIdx) const { return this->table.densityLiquidIdx(rho, tempIdx); }
706 Scalar minP(int tempIdx) const { return this->table.minLiquidPressure(tempIdx); }
707 Scalar maxP(int tempIdx) const { return this->table.maxLiquidPressure(tempIdx); }
708 Scalar minRho(int tempIdx) const { return this->table.minLiquidDensity(tempIdx); }
709 Scalar maxRho(int tempIdx) const { return this->table.maxLiquidDensity(tempIdx); }
710 };
711 public:
712 //! export scalar type
713 using Scalar = typename RawComponent::Scalar;
714
715 //! state that we are tabulated
716 static constexpr bool isTabulated = true;
717
718 /*!
719 * \brief Initialize the tables.
720 *
721 * \param tempMin The minimum of the temperature range in \f$\mathrm{[K]}\f$
722 * \param tempMax The maximum of the temperature range in \f$\mathrm{[K]}\f$
723 * \param nTemp The number of entries/steps within the temperature range
724 * \param pressMin The minimum of the pressure range in \f$\mathrm{[Pa]}\f$
725 * \param pressMax The maximum of the pressure range in \f$\mathrm{[Pa]}\f$
726 * \param nPress The number of entries/steps within the pressure range
727 */
728 215 static void init(Scalar tempMin, Scalar tempMax, std::size_t nTemp,
729 Scalar pressMin, Scalar pressMax, std::size_t nPress)
730 {
731 #ifndef NDEBUG
732 215 warningPrinted_ = false;
733 #endif
734 std::cout << "-------------------------------------------------------------------------\n"
735 << "Initializing tables for the " << RawComponent::name()
736
2/6
✓ Branch 4 taken 192 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 192 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
860 << " fluid properties (" << nTemp*nPress << " entries).\n"
737
2/4
✓ Branch 4 taken 192 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 192 times.
✗ Branch 8 not taken.
645 << "Temperature -> min: " << std::scientific << std::setprecision(3)
738
3/6
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 192 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 192 times.
✗ Branch 10 not taken.
645 << tempMin << ", max: " << tempMax << ", n: " << nTemp << '\n'
739
3/6
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 192 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 192 times.
✗ Branch 10 not taken.
430 << "Pressure -> min: " << std::scientific << std::setprecision(3)
740
3/6
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 192 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 192 times.
✗ Branch 10 not taken.
645 << pressMin << ", max: " << pressMax << ", n: " << nPress << '\n'
741
3/6
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 192 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 192 times.
430 << "-------------------------------------------------------------------------" << std::endl;
742
743 215 table().init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
744
745 215 initialized_ = true;
746 215 }
747
748 /*!
749 * \brief A human readable name for the component.
750 */
751 static std::string name()
752 399335 { return RawComponent::name(); }
753
754 /*!
755 * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of the component.
756 */
757 static constexpr Scalar molarMass()
758 { return RawComponent::molarMass(); }
759
760 /*!
761 * \brief Returns the critical temperature in \f$\mathrm{[K]}\f$ of the component.
762 */
763 static Scalar criticalTemperature()
764 { return RawComponent::criticalTemperature(); }
765
766 /*!
767 * \brief Returns the critical pressure in \f$\mathrm{[Pa]}\f$ of the component.
768 */
769 static Scalar criticalPressure()
770 { return RawComponent::criticalPressure(); }
771
772 /*!
773 * \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
774 */
775 static Scalar tripleTemperature()
776 { return RawComponent::tripleTemperature(); }
777
778 /*!
779 * \brief Returns the pressure in \f$\mathrm{[Pa]}\f$ at the component's triple point.
780 */
781 static Scalar triplePressure()
782 { return RawComponent::triplePressure(); }
783
784 /*!
785 * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
786 * temperature.
787 *
788 * \param T temperature of component
789 */
790 31462062 static Scalar vaporPressure(Scalar T)
791 {
792 using std::isnan;
793 31462062 Scalar result = interpolateT_(table().vaporPressure_, T, InterpolatePolicy{{ table() }});
794
2/2
✓ Branch 0 taken 362045 times.
✓ Branch 1 taken 31100017 times.
31462062 if (isnan(result))
795 362045 return RawComponent::vaporPressure(T);
796 return result;
797 }
798
799 /*!
800 * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
801 * temperature.
802 *
803 * The method is only called by the sequential flash, so tabulating is omitted.
804 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
805 */
806 static Scalar vaporTemperature(Scalar pressure)
807 {
808 100333 return RawComponent::vaporTemperature(pressure);
809 }
810
811 /*!
812 * \brief Specific enthalpy of the gas \f$\mathrm{[J/kg]}\f$.
813 *
814 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
815 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
816 */
817 124929034 static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
818 {
819 124929034 Scalar result = interpolateTP_(table().gasEnthalpy_, temperature, pressure, InterpolateGasPolicy{{ table() }});
820 using std::isnan;
821
2/2
✓ Branch 0 taken 205992 times.
✓ Branch 1 taken 124723042 times.
124929034 if (isnan(result))
822 {
823
4/10
✓ Branch 2 taken 205992 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 205992 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 205992 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 205992 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
205992 printWarningTP_("gasEnthalpy", temperature, pressure, InterpolateGasPolicy{{ table() }});
824 205992 return RawComponent::gasEnthalpy(temperature, pressure);
825 }
826 return result;
827 }
828
829 /*!
830 * \brief Specific enthalpy of the liquid \f$\mathrm{[J/kg]}\f$.
831 *
832 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
833 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
834 */
835 147998035 static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
836 {
837 147998035 Scalar result = interpolateTP_(table().liquidEnthalpy_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
838 using std::isnan;
839
2/2
✓ Branch 0 taken 216394 times.
✓ Branch 1 taken 147781641 times.
147998035 if (isnan(result))
840 {
841
4/10
✓ Branch 2 taken 216394 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216394 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 216394 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 216394 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
216394 printWarningTP_("liquidEnthalpy", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
842 216394 return RawComponent::liquidEnthalpy(temperature, pressure);
843 }
844 return result;
845 }
846
847 /*!
848 * \brief Specific isobaric heat capacity of the gas \f$\mathrm{[J/(kg*K)]}\f$.
849 *
850 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
851 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
852 */
853 18189597 static const Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
854 {
855 18189597 Scalar result = interpolateTP_(table().gasHeatCapacity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
856 using std::isnan;
857
2/2
✓ Branch 0 taken 301 times.
✓ Branch 1 taken 18189296 times.
18189597 if (isnan(result))
858 {
859
4/10
✓ Branch 2 taken 301 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 301 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 301 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 301 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
301 printWarningTP_("gasHeatCapacity", temperature, pressure, InterpolateGasPolicy{{ table() }});
860 301 return RawComponent::gasHeatCapacity(temperature, pressure);
861 }
862 return result;
863 }
864
865 /*!
866 * \brief Specific isobaric heat capacity of the liquid \f$\mathrm{[J/(kg*K)]}\f$.
867 *
868 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
869 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
870 */
871 2198232 static const Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
872 {
873 2198232 Scalar result = interpolateTP_(table().liquidHeatCapacity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
874 using std::isnan;
875
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2198232 times.
2198232 if (isnan(result))
876 {
877 printWarningTP_("liquidHeatCapacity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
878 return RawComponent::liquidHeatCapacity(temperature, pressure);
879 }
880 return result;
881 }
882
883 /*!
884 * \brief Specific internal energy of the gas \f$\mathrm{[J/kg]}\f$.
885 *
886 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
887 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
888 */
889 205023 static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
890 {
891 205023 return gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure);
892 }
893
894 /*!
895 * \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
896 *
897 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
898 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
899 */
900 1235383 static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
901 {
902 1235383 return liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
903 }
904
905 /*!
906 * \brief The pressure of gas in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
907 *
908 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
909 * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
910 */
911 205023 static Scalar gasPressure(Scalar temperature, Scalar density)
912 {
913 205023 Scalar result = interpolateTRho_(table().gasPressure_, temperature, density, InterpolateGasPolicy{{ table() }});
914 using std::isnan;
915
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 205023 times.
205023 if (isnan(result))
916 {
917 printWarningTRho_("gasPressure", temperature, density, InterpolateGasPolicy{{ table() }});
918 return RawComponent::gasPressure(temperature, density);
919 }
920 return result;
921 }
922
923 /*!
924 * \brief The pressure of liquid in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
925 *
926 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
927 * \param density density of component in \f$\mathrm{[kg/m^3]}\f$
928 */
929 1235383 static Scalar liquidPressure(Scalar temperature, Scalar density)
930 {
931 1235383 Scalar result = interpolateTRho_(table().liquidPressure_, temperature, density, InterpolateLiquidPolicy{{ table() }});
932 using std::isnan;
933
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1235383 times.
1235383 if (isnan(result))
934 {
935 printWarningTRho_("liquidPressure", temperature, density, InterpolateLiquidPolicy{{ table() }});
936 return RawComponent::liquidPressure(temperature, density);
937 }
938 return result;
939 }
940
941 /*!
942 * \brief Returns true if the gas phase is assumed to be compressible
943 */
944 static constexpr bool gasIsCompressible()
945 { return RawComponent::gasIsCompressible(); }
946
947 /*!
948 * \brief Returns true if the liquid phase is assumed to be compressible
949 */
950 static constexpr bool liquidIsCompressible()
951 { return RawComponent::liquidIsCompressible(); }
952
953 /*!
954 * \brief Returns true if the gas phase is assumed to be ideal
955 */
956 static constexpr bool gasIsIdeal()
957 { return RawComponent::gasIsIdeal(); }
958
959
960 /*!
961 * \brief The density of gas at a given pressure and temperature
962 * \f$\mathrm{[kg/m^3]}\f$.
963 *
964 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
965 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
966 */
967 139977410 static Scalar gasDensity(Scalar temperature, Scalar pressure)
968 {
969 139977410 Scalar result = interpolateTP_(table().gasDensity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
970 using std::isnan;
971
2/2
✓ Branch 0 taken 687215 times.
✓ Branch 1 taken 139290195 times.
139977410 if (isnan(result))
972 {
973
4/11
✓ Branch 2 taken 687215 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 687215 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 687215 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 687215 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
687215 printWarningTP_("gasDensity", temperature, pressure, InterpolateGasPolicy{{ table() }});
974 687215 return RawComponent::gasDensity(temperature, pressure);
975 }
976 return result;
977 }
978
979 /*!
980 * \brief The molar density of gas in \f$\mathrm{[mol/m^3]}\f$
981 * at a given pressure and temperature.
982 *
983 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
984 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
985 *
986 */
987 static Scalar gasMolarDensity(Scalar temperature, Scalar pressure)
988 64471898 { return gasDensity(temperature, pressure)/molarMass(); }
989
990 /*!
991 * \brief The density of liquid at a given pressure and
992 * temperature \f$\mathrm{[kg/m^3]}\f$.
993 *
994 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
995 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
996 */
997 580814648 static Scalar liquidDensity(Scalar temperature, Scalar pressure)
998 {
999 580814648 Scalar result = interpolateTP_(table().liquidDensity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1000 using std::isnan;
1001
2/2
✓ Branch 0 taken 868529 times.
✓ Branch 1 taken 579946119 times.
580814648 if (isnan(result))
1002 {
1003
4/10
✓ Branch 2 taken 868529 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 868529 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 868529 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 868529 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
868529 printWarningTP_("liquidDensity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1004 868529 return RawComponent::liquidDensity(temperature, pressure);
1005 }
1006
1007 return result;
1008 }
1009
1010 /*!
1011 * \brief The molar density of liquid in \f$\mathrm{[mol/m^3]}\f$
1012 * at a given pressure and temperature.
1013 *
1014 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
1015 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
1016 *
1017 */
1018 static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
1019 101298805 { return liquidDensity(temperature, pressure)/molarMass(); }
1020
1021 /*!
1022 * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of gas.
1023 *
1024 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
1025 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
1026 */
1027 5798153 static Scalar gasViscosity(Scalar temperature, Scalar pressure)
1028 {
1029 5798153 Scalar result = interpolateTP_(table().gasViscosity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
1030 using std::isnan;
1031
2/2
✓ Branch 0 taken 105 times.
✓ Branch 1 taken 5798048 times.
5798153 if (isnan(result))
1032 {
1033
4/10
✓ Branch 2 taken 105 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 105 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 105 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 105 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
105 printWarningTP_("gasViscosity", temperature, pressure, InterpolateGasPolicy{{ table() }});
1034 105 return RawComponent::gasViscosity(temperature, pressure);
1035 }
1036 return result;
1037 }
1038
1039 /*!
1040 * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of liquid.
1041 *
1042 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
1043 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
1044 */
1045 68181721 static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
1046 {
1047 68181721 Scalar result = interpolateTP_(table().liquidViscosity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1048 using std::isnan;
1049
2/2
✓ Branch 0 taken 352712 times.
✓ Branch 1 taken 67829009 times.
68181721 if (isnan(result))
1050 {
1051
4/10
✓ Branch 2 taken 352712 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 352712 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 352712 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 352712 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
352712 printWarningTP_("liquidViscosity",temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1052 352712 return RawComponent::liquidViscosity(temperature, pressure);
1053 }
1054 return result;
1055 }
1056
1057 /*!
1058 * \brief The thermal conductivity of gaseous water \f$\mathrm{[W/(m*K)]}\f$.
1059 *
1060 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
1061 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
1062 */
1063 16751909 static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
1064 {
1065 16751909 Scalar result = interpolateTP_(table().gasThermalConductivity_, temperature, pressure, InterpolateGasPolicy{{ table() }});
1066 using std::isnan;
1067
2/2
✓ Branch 0 taken 139404 times.
✓ Branch 1 taken 16612505 times.
16751909 if (isnan(result))
1068 {
1069
4/10
✓ Branch 2 taken 139404 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 139404 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 139404 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 139404 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
278808 printWarningTP_("gasThermalConductivity", temperature, pressure, InterpolateGasPolicy{{ table() }});
1070 139404 return RawComponent::gasThermalConductivity(temperature, pressure);
1071 }
1072 return result;
1073 }
1074
1075 /*!
1076 * \brief The thermal conductivity of liquid water \f$\mathrm{[W/(m*K)]}\f$.
1077 *
1078 * \param temperature temperature of component in \f$\mathrm{[K]}\f$
1079 * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
1080 */
1081 94597915 static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
1082 {
1083 94597915 Scalar result = interpolateTP_(table().liquidThermalConductivity_, temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1084 using std::isnan;
1085
2/2
✓ Branch 0 taken 209888 times.
✓ Branch 1 taken 94388027 times.
94597915 if (isnan(result))
1086 {
1087
4/10
✓ Branch 2 taken 209888 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 209888 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 209888 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 209888 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
419776 printWarningTP_("liquidThermalConductivity", temperature, pressure, InterpolateLiquidPolicy{{ table() }});
1088 209888 return RawComponent::liquidThermalConductivity(temperature, pressure);
1089 }
1090 return result;
1091 }
1092
1093
1094 private:
1095 //! prints a warning if the result is not in range or the table has not been initialized
1096 template<class Policy>
1097 5361078 static void printWarningTP_(const std::string& quantity, Scalar T, Scalar p, Policy policy)
1098 {
1099 #ifndef NDEBUG
1100
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 2680522 times.
5361078 if (warningPrinted_)
1101 return;
1102
1103
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 13 times.
35 if (!initialized_)
1104 std::cerr << "Warning: tabulated component '" << name()
1105 << "' has not been initialized. "
1106
3/8
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
27 << "Call FluidSystem::init() to use the tabulation in order to reduce runtime. \n";
1107 else
1108 104 std::cerr << "Warning: "<<quantity<<"(T="<<T<<", p="<<p<<") of component '"<<name()
1109
4/10
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 13 times.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 13 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
130 << "' is outside tabulation range: ("<< policy.tempMin() <<"<=T<="<< policy.tempMax() <<"), ("
1110
6/10
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 12 times.
✓ Branch 9 taken 13 times.
✗ Branch 10 not taken.
64 << policy.minP(0) <<"<=p<=" << policy.maxP(policy.nTemp()-1) <<"). "
1111
2/4
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
78 << "Forwarded to FluidSystem for direct evaluation of "<<quantity<<". \n";
1112 35 warningPrinted_ = true;
1113 #endif
1114 }
1115
1116 //! prints a warning if the result is not in range or the table has not been initialized
1117 template<class Policy>
1118 static void printWarningTRho_(const std::string& quantity, Scalar T, Scalar rho, Policy policy)
1119 {
1120 #ifndef NDEBUG
1121 if (warningPrinted_)
1122 return;
1123
1124 if (!initialized_)
1125 std::cerr << "Warning: tabulated component '" << name()
1126 << "' has not been initialized. "
1127 << "Call FluidSystem::init() to use the tabulation in order to reduce runtime. \n";
1128 else
1129 {
1130 const auto [densityMin, densityMax] = [&]
1131 {
1132 const Scalar alphaT = policy.tempIdx(T);
1133 using std::clamp;
1134 const auto iT = clamp<int>(static_cast<int>(alphaT), 0, policy.nTemp() - 2);
1135 return std::make_pair( policy.minRho(iT), policy.maxRho(iT) );
1136 }();
1137
1138 std::cerr << "Warning: "<<quantity<<"(T="<<T<<", density="<<rho<<") of component '"<<name()
1139 << "' is outside tabulation range: ("<< policy.tempMin() <<"<=T<="<< policy.tempMax() <<"), ("
1140 << densityMin <<"<=density<=" << densityMax <<"). "
1141 << "Forwarded to FluidSystem for direct evaluation of "<<quantity<<". \n";
1142 }
1143 warningPrinted_ = true;
1144 #endif
1145 }
1146
1147 //! returns an interpolated value depending on temperature
1148 template<class Policy>
1149 31462062 static Scalar interpolateT_(const std::vector<Scalar>& values, Scalar T, Policy policy)
1150 {
1151 62924124 Scalar alphaT = policy.tempIdx(T);
1152 62924124 const auto nTemp = policy.nTemp();
1153
4/4
✓ Branch 0 taken 31141666 times.
✓ Branch 1 taken 320396 times.
✓ Branch 2 taken 31100017 times.
✓ Branch 3 taken 41649 times.
31462062 if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
1154 return std::numeric_limits<Scalar>::quiet_NaN();
1155
1156 using std::clamp;
1157
1/2
✓ Branch 0 taken 31100017 times.
✗ Branch 1 not taken.
31100017 const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
1158 31100017 alphaT -= iT;
1159
1160 31100017 return values[iT ]*(1 - alphaT) +
1161 62200034 values[iT + 1]*( alphaT);
1162 }
1163
1164 //! returns an interpolated value depending on temperature and pressure
1165 template<class Policy>
1166 2225682350 static Scalar interpolateTP_(const std::vector<Scalar>& values, const Scalar T, const Scalar p, Policy policy)
1167 {
1168 4451364700 const auto nTemp = policy.nTemp();
1169 4451364700 const auto nPress = policy.nPress();
1170
1171 4451364700 Scalar alphaT = policy.tempIdx(T);
1172
4/4
✓ Branch 0 taken 1197281233 times.
✓ Branch 1 taken 2155421 times.
✓ Branch 2 taken 1196756114 times.
✓ Branch 3 taken 525119 times.
2225682350 if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
1173 return std::numeric_limits<Scalar>::quiet_NaN();
1174
1175 using std::clamp;
1176
1/2
✓ Branch 0 taken 1196756114 times.
✗ Branch 1 not taken.
2220321272 const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
1177 2220321272 alphaT -= iT;
1178
1179 4440642544 Scalar alphaP1 = policy.pressIdx(p, iT);
1180 4440642544 Scalar alphaP2 = policy.pressIdx(p, iT + 1);
1181
1182
2/2
✓ Branch 0 taken 1194079366 times.
✓ Branch 1 taken 2676748 times.
2220321272 const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nPress - 2);
1183
2/2
✓ Branch 0 taken 1188887569 times.
✓ Branch 1 taken 7868545 times.
2220321272 const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nPress - 2);
1184 2220321272 alphaP1 -= iP1;
1185 2220321272 alphaP2 -= iP2;
1186
1187 #if 0 && !defined NDEBUG
1188 const auto tempMin = policy.tempMin();
1189 const auto tempMin = policy.tempMax();
1190 if(!(0 <= alphaT && alphaT <= 1.0))
1191 DUNE_THROW(NumericalProblem, "Temperature out of range: "
1192 << "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
1193 if(!(0 <= alphaP1 && alphaP1 <= 1.0))
1194 DUNE_THROW(NumericalProblem, "First pressure out of range: "
1195 << "p=" << p << " range: [" << policy.minP(policy.tempIdx(T)) << ", " << policy.maxP(policy.tempIdx(T)) << "]");
1196 if(!(0 <= alphaP2 && alphaP2 <= 1.0))
1197 DUNE_THROW(NumericalProblem, "Second pressure out of range: "
1198 << "p=" << p << " range: [" << policy.minP(policy.tempIdx(T) + 1) << ", " << policy.maxP(policy.tempIdx(T) + 1) << "]");
1199 #endif
1200 2220321272 return values[(iT ) + (iP1 )*nTemp]*(1 - alphaT)*(1 - alphaP1) +
1201 2220321272 values[(iT ) + (iP1 + 1)*nTemp]*(1 - alphaT)*( alphaP1) +
1202 2220321272 values[(iT + 1) + (iP2 )*nTemp]*( alphaT)*(1 - alphaP2) +
1203 4440642544 values[(iT + 1) + (iP2 + 1)*nTemp]*( alphaT)*( alphaP2);
1204 }
1205
1206 //! returns an interpolated value for gas depending on temperature and density
1207 template<class Policy>
1208 2880812 static Scalar interpolateTRho_(const std::vector<Scalar>& values, const Scalar T, const Scalar rho, Policy policy)
1209 {
1210 5761624 const auto nTemp = policy.nTemp();
1211 5761624 const auto nDensity = policy.nDensity();
1212
1213 using std::clamp;
1214 5761624 Scalar alphaT = policy.tempIdx(T);
1215
2/4
✓ Branch 0 taken 1440406 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1440406 times.
✗ Branch 3 not taken.
2880812 if (alphaT < 0 - 1e-7*nTemp || alphaT >= nTemp - 1 + 1e-7*nTemp)
1216 return std::numeric_limits<Scalar>::quiet_NaN();
1217
1218
1/2
✓ Branch 0 taken 1440406 times.
✗ Branch 1 not taken.
2880812 const auto iT = clamp<int>(static_cast<int>(alphaT), 0, nTemp - 2);
1219 2880812 alphaT -= iT;
1220
1221
2/2
✓ Branch 0 taken 1433544 times.
✓ Branch 1 taken 6862 times.
2880812 Scalar alphaP1 = policy.rhoIdx(rho, iT);
1222
2/2
✓ Branch 0 taken 1433544 times.
✓ Branch 1 taken 6862 times.
2880812 Scalar alphaP2 = policy.rhoIdx(rho, iT + 1);
1223
2/2
✓ Branch 0 taken 1433544 times.
✓ Branch 1 taken 6862 times.
2880812 const auto iP1 = clamp<int>(static_cast<int>(alphaP1), 0, nDensity - 2);
1224
1/2
✓ Branch 0 taken 1440406 times.
✗ Branch 1 not taken.
2880812 const auto iP2 = clamp<int>(static_cast<int>(alphaP2), 0, nDensity - 2);
1225 2880812 alphaP1 -= iP1;
1226 2880812 alphaP2 -= iP2;
1227
1228 2880812 return values[(iT ) + (iP1 )*nTemp]*(1 - alphaT)*(1 - alphaP1) +
1229 2880812 values[(iT ) + (iP1 + 1)*nTemp]*(1 - alphaT)*( alphaP1) +
1230 2880812 values[(iT + 1) + (iP2 )*nTemp]*( alphaT)*(1 - alphaP2) +
1231 5761624 values[(iT + 1) + (iP2 + 1)*nTemp]*( alphaT)*( alphaP2);
1232 }
1233
1234 // specifies whether the table was initialized
1235 static bool initialized_;
1236
1237 #ifndef NDEBUG
1238 // specifies whether some warning was printed
1239 static bool warningPrinted_;
1240 #endif
1241
1242 2467358978 static Table& table()
1243 {
1244
3/4
✓ Branch 0 taken 155 times.
✓ Branch 1 taken 2422070202 times.
✓ Branch 3 taken 155 times.
✗ Branch 4 not taken.
2467358978 static Table t;
1245 2467358978 return t;
1246 }
1247 };
1248
1249 template <class RawComponent, bool useVaporPressure>
1250 bool TabulatedComponent<RawComponent, useVaporPressure>::initialized_ = false;
1251
1252 #ifndef NDEBUG
1253 template <class RawComponent, bool useVaporPressure>
1254 bool TabulatedComponent<RawComponent, useVaporPressure>::warningPrinted_ = false;
1255 #endif
1256
1257 // forward declaration
1258 template <class Component>
1259 struct IsAqueous;
1260
1261 // we are aqueous if the raw compont is so
1262 template <class RawComponent, bool useVaporPressure>
1263 struct IsAqueous<TabulatedComponent<RawComponent, useVaporPressure>> : public IsAqueous<RawComponent> {};
1264
1265 } // end namespace Dumux::Components
1266
1267 #endif
1268