GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porousmediumflow/2p1c/volumevariables.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 72 76 94.7%
Functions: 6 8 75.0%
Branches: 76 124 61.3%

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 TwoPOneCModel
10 * \copydoc Dumux::TwoPOneCVolumeVariables
11 */
12
13 #ifndef DUMUX_2P1C_VOLUME_VARIABLES_HH
14 #define DUMUX_2P1C_VOLUME_VARIABLES_HH
15
16 #include <array>
17
18 #include <dune/common/exceptions.hh>
19
20 #include <dumux/porousmediumflow/volumevariables.hh>
21 #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
22 #include <dumux/porousmediumflow/2p/formulation.hh>
23 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
24
25 #include "primaryvariableswitch.hh"
26
27 namespace Dumux {
28
29 /*!
30 * \ingroup TwoPOneCModel
31 * \brief The volume variables (i.e. secondary variables) for the two-phase one-component model.
32 */
33 template <class Traits>
34 2928336 class TwoPOneCVolumeVariables
35 : public PorousMediumFlowVolumeVariables<Traits>
36 , public EnergyVolumeVariables<Traits, TwoPOneCVolumeVariables<Traits> >
37 {
38 using ParentType = PorousMediumFlowVolumeVariables<Traits>;
39 using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPOneCVolumeVariables<Traits> >;
40 using Scalar = typename Traits::PrimaryVariables::value_type;
41 using PermeabilityType = typename Traits::PermeabilityType;
42 using FS = typename Traits::FluidSystem;
43 using Idx = typename Traits::ModelTraits::Indices;
44 static constexpr int numFluidComps = ParentType::numFluidComponents();
45
46 // primary variable indices
47 enum
48 {
49 numFluidPhases = Traits::ModelTraits::numFluidPhases(),
50 switchIdx = Idx::switchIdx,
51 pressureIdx = Idx::pressureIdx
52 };
53
54 // component indices
55 enum
56 {
57 comp0Idx = FS::comp0Idx,
58 liquidPhaseIdx = FS::liquidPhaseIdx,
59 gasPhaseIdx = FS::gasPhaseIdx
60 };
61
62 // phase presence indices
63 enum
64 {
65 twoPhases = Idx::twoPhases,
66 liquidPhaseOnly = Idx::liquidPhaseOnly,
67 gasPhaseOnly = Idx::gasPhaseOnly,
68 };
69
70 // formulations
71 static constexpr auto formulation = Traits::ModelTraits::priVarFormulation();
72
73 public:
74 //! The type of the object returned by the fluidState() method
75 using FluidState = typename Traits::FluidState;
76 //! The type of the fluid system
77 using FluidSystem = typename Traits::FluidSystem;
78 //! The type of the indices
79 using Indices = typename Traits::ModelTraits::Indices;
80 //! Export type of solid state
81 using SolidState = typename Traits::SolidState;
82 //! Export type of solid system
83 using SolidSystem = typename Traits::SolidSystem;
84 //! Export the primary variable switch
85 using PrimaryVariableSwitch = TwoPOneCPrimaryVariableSwitch;
86
87 //! Return the two-phase formulation used here
88 static constexpr TwoPFormulation priVarFormulation() { return formulation; }
89
90 // check for permissive combinations
91 static_assert(Traits::ModelTraits::numFluidPhases() == 2, "NumPhases set in the model is not two!");
92 static_assert(Traits::ModelTraits::numFluidComponents() == 1, "NumComponents set in the model is not one!");
93 static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");
94
95 /*!
96 * \brief Updates all quantities for a given control volume
97 *
98 * \param elemSol A vector containing all primary variables connected to the element
99 * \param problem The object specifying the problem which ought to
100 * be simulated
101 * \param element An element which contains part of the control volume
102 * \param scv The sub-control volume
103 */
104 template<class ElemSol, class Problem, class Element, class Scv>
105 1210607 void update(const ElemSol &elemSol,
106 const Problem &problem,
107 const Element &element,
108 const Scv& scv)
109 {
110 1210607 ParentType::update(elemSol, problem, element, scv);
111
112 1210607 completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
113
114 /////////////
115 // calculate the remaining quantities
116 /////////////
117
118 1210603 const auto& spatialParams = problem.spatialParams();
119 1210603 const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
120
121 // Second instance of a parameter cache.
122 // Could be avoided if diffusion coefficients also
123 // became part of the fluid state.
124 typename FluidSystem::ParameterCache paramCache;
125 1210603 paramCache.updateAll(fluidState_);
126 1210603 const int wPhaseIdx = fluidState_.wettingPhase();
127
2/2
✓ Branch 0 taken 2421206 times.
✓ Branch 1 taken 1210603 times.
3631809 for (int phaseIdx = 0; phaseIdx < numFluidPhases; ++phaseIdx)
128 {
129 // relative permeabilities
130 Scalar kr;
131
2/2
✓ Branch 0 taken 1210603 times.
✓ Branch 1 taken 1210603 times.
2421206 if (phaseIdx == wPhaseIdx)
132 1210603 kr = fluidMatrixInteraction.krw(saturation(wPhaseIdx));
133 else // ATTENTION: krn requires the wetting phase saturation
134 // as parameter!
135 1210603 kr = fluidMatrixInteraction.krn(saturation(wPhaseIdx));
136 4842412 relativePermeability_[phaseIdx] = kr;
137 }
138
139 // porosity & permeability
140 // porosity calculation over inert volumefraction
141 1210603 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
142 1210603 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
143 1210603 permeability_ = spatialParams.permeability(element, scv, elemSol);
144 1210603 EnergyVolVars::updateEffectiveThermalConductivity();
145 1210599 }
146
147 /*!
148 * \brief Sets complete fluid state
149 *
150 * \param elemSol A vector containing all primary variables connected to the element
151 * \param problem The object specifying the problem which ought to
152 * be simulated
153 * \param element An element which contains part of the control volume
154 * \param scv The sub-control volume
155 * \param fluidState A container with the current (physical) state of the fluid
156 * \param solidState A container with the current (physical) state of the solid
157 */
158 template<class ElemSol, class Problem, class Element, class Scv>
159 1210607 void completeFluidState(const ElemSol& elemSol,
160 const Problem& problem,
161 const Element& element,
162 const Scv& scv,
163 FluidState& fluidState,
164 SolidState& solidState)
165 {
166
167 // capillary pressure parameters
168
2/2
✓ Branch 0 taken 825160 times.
✓ Branch 1 taken 385447 times.
1210607 const auto& spatialParams = problem.spatialParams();
169
2/2
✓ Branch 0 taken 825160 times.
✓ Branch 1 taken 385447 times.
1210607 const auto wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
170 1210607 fluidState.setWettingPhase(wPhaseIdx);
171
172 1210607 const auto& priVars = elemSol[scv.localDofIndex()];
173 1210607 const auto phasePresence = priVars.state();
174
175 // set the saturations
176
2/2
✓ Branch 0 taken 2643 times.
✓ Branch 1 taken 1207964 times.
1210607 if (phasePresence == twoPhases)
177 {
178 if (formulation == TwoPFormulation::p0s1)
179 {
180 fluidState.setSaturation(gasPhaseIdx, priVars[switchIdx]);
181 fluidState.setSaturation(liquidPhaseIdx, 1.0 - priVars[switchIdx]);
182 }
183 else
184 {
185 5286 fluidState.setSaturation(liquidPhaseIdx, priVars[switchIdx]);
186 5286 fluidState.setSaturation(gasPhaseIdx, 1.0 - priVars[switchIdx]);
187 }
188 }
189
2/2
✓ Branch 0 taken 1207907 times.
✓ Branch 1 taken 57 times.
1207964 else if (phasePresence == liquidPhaseOnly)
190 {
191 1207907 fluidState.setSaturation(liquidPhaseIdx, 1.0);
192 1207907 fluidState.setSaturation(gasPhaseIdx, 0.0);
193 }
194
1/2
✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
57 else if (phasePresence == gasPhaseOnly)
195 {
196 57 fluidState.setSaturation(liquidPhaseIdx, 0.0);
197 57 fluidState.setSaturation(gasPhaseIdx, 1.0);
198 }
199 else
200 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
201
202 // set pressures of the fluid phases
203 1210607 const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
204 2421214 pc_ = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
205 if (formulation == TwoPFormulation::p0s1)
206 {
207 fluidState.setPressure(liquidPhaseIdx, priVars[pressureIdx]);
208 fluidState.setPressure(gasPhaseIdx, (wPhaseIdx == liquidPhaseIdx) ? priVars[pressureIdx] + pc_
209 : priVars[pressureIdx] - pc_);
210 }
211 else
212 {
213
4/4
✓ Branch 0 taken 825160 times.
✓ Branch 1 taken 385447 times.
✓ Branch 2 taken 825160 times.
✓ Branch 3 taken 385447 times.
2421214 fluidState.setPressure(gasPhaseIdx, priVars[pressureIdx]);
214
2/2
✓ Branch 0 taken 825160 times.
✓ Branch 1 taken 385447 times.
1596054 fluidState.setPressure(liquidPhaseIdx, (wPhaseIdx == liquidPhaseIdx) ? priVars[pressureIdx] - pc_
215 770894 : priVars[pressureIdx] + pc_);
216 }
217
218 // set the temperature
219 1210607 updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
220
221 // set the densities
222
2/2
✓ Branch 0 taken 2421210 times.
✓ Branch 1 taken 1210603 times.
3631813 for (int phaseIdx = 0; phaseIdx < numFluidPhases; ++phaseIdx)
223 {
224 2421210 Scalar rho = FluidSystem::density(fluidState, phaseIdx);
225 2421206 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, phaseIdx);
226
227 2421206 fluidState.setDensity(phaseIdx, rho);
228 4842412 fluidState.setMolarDensity(phaseIdx, rhoMolar);
229 }
230
231 //get the viscosity and mobility
232
2/2
✓ Branch 0 taken 2421206 times.
✓ Branch 1 taken 1210603 times.
3631809 for (int phaseIdx = 0; phaseIdx < numFluidPhases; ++phaseIdx)
233 {
234 // Mobilities
235 const Scalar mu =
236 2421206 FluidSystem::viscosity(fluidState,
237 phaseIdx);
238 4842412 fluidState.setViscosity(phaseIdx,mu);
239 }
240
241 // the enthalpies (internal energies are directly calculated in the fluidstate
242
2/2
✓ Branch 0 taken 2421206 times.
✓ Branch 1 taken 1210603 times.
3631809 for (int phaseIdx = 0; phaseIdx < numFluidPhases; ++phaseIdx)
243 {
244 2421206 const Scalar h = FluidSystem::enthalpy(fluidState, phaseIdx);
245 4842412 fluidState.setEnthalpy(phaseIdx, h);
246 }
247 1210603 }
248
249 //! Depending on the phase state, the fluid temperature is either obtained as a primary variable from the solution vector
250 //! or calculated from the liquid's vapor pressure.
251 template<class ElemSol, class Problem, class Element, class Scv>
252 1210607 void updateTemperature(const ElemSol& elemSol,
253 const Problem& problem,
254 const Element& element,
255 const Scv& scv,
256 FluidState& fluidState,
257 SolidState& solidState)
258 {
259 1210607 const auto& priVars = elemSol[scv.localDofIndex()];
260 1210607 const auto phasePresence = priVars.state();
261
262 // get temperature
263 Scalar fluidTemperature;
264
2/2
✓ Branch 0 taken 1207964 times.
✓ Branch 1 taken 2643 times.
1210607 if (phasePresence == liquidPhaseOnly || phasePresence == gasPhaseOnly)
265 2415928 fluidTemperature = priVars[switchIdx];
266
1/2
✓ Branch 0 taken 2643 times.
✗ Branch 1 not taken.
2643 else if (phasePresence == twoPhases)
267 2643 fluidTemperature = FluidSystem::vaporTemperature(fluidState, fluidState_.wettingPhase());
268 else
269 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
270
271 // the model assumes that all fluid phases have the same temperature
272
2/2
✓ Branch 0 taken 2421214 times.
✓ Branch 1 taken 1210607 times.
3631821 for (int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
273 4842428 fluidState.setTemperature(phaseIdx, fluidTemperature);
274
275 // the solid phase could have a different temperature
276 if (Traits::ModelTraits::numEnergyEq() == 1)
277 1210607 solidState.setTemperature(fluidTemperature);
278 else
279 {
280 const Scalar solidTemperature = elemSol[scv.localDofIndex()][Traits::ModelTraits::numEq()-1];
281 solidState.setTemperature(solidTemperature);
282 }
283 1210607 }
284
285 /*!
286 * \brief Returns the fluid state for the control-volume.
287 */
288 const FluidState &fluidState() const
289 15049134 { return fluidState_; }
290
291 /*!
292 * \brief Returns the phase state for the control volume.
293 */
294 const SolidState &solidState() const
295 1210599 { return solidState_; }
296
297 /*!
298 * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase.
299 *
300 * \param phaseIdx The phase index
301 */
302 Scalar averageMolarMass(int phaseIdx) const
303 { return fluidState_.averageMolarMass(phaseIdx); }
304
305 /*!
306 * \brief Returns the effective saturation of a given phase within
307 * the control volume.
308 *
309 * \param phaseIdx The phase index
310 */
311 Scalar saturation(const int phaseIdx) const
312
24/32
✓ Branch 0 taken 1062090 times.
✓ Branch 1 taken 2886 times.
✓ Branch 2 taken 1062090 times.
✓ Branch 3 taken 2886 times.
✓ Branch 4 taken 1062090 times.
✓ Branch 5 taken 2886 times.
✓ Branch 6 taken 1062090 times.
✓ Branch 7 taken 2886 times.
✓ Branch 8 taken 1502115 times.
✓ Branch 9 taken 4209 times.
✓ Branch 10 taken 1502115 times.
✓ Branch 11 taken 4209 times.
✓ Branch 12 taken 1502115 times.
✓ Branch 13 taken 4209 times.
✓ Branch 14 taken 1502115 times.
✓ Branch 15 taken 4209 times.
✓ Branch 18 taken 11 times.
✓ Branch 19 taken 266 times.
✓ Branch 20 taken 11 times.
✓ Branch 21 taken 266 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 12 times.
✓ Branch 29 taken 254 times.
✓ Branch 30 taken 12 times.
✓ Branch 31 taken 254 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
31373372 { return fluidState_.saturation(phaseIdx); }
313
314 /*!
315 * \brief Returns the mass density of a given phase within the
316 * control volume.
317 *
318 * \param phaseIdx The phase index
319 */
320 Scalar density(const int phaseIdx) const
321 80355600 { return fluidState_.density(phaseIdx); }
322
323 /*!
324 * \brief Returns the molar density of a given phase within the
325 * control volume.
326 *
327 * \param phaseIdx The phase index
328 */
329 Scalar molarDensity(const int phaseIdx) const
330 { return fluidState_.molarDensity(phaseIdx); }
331
332 /*!
333 * \brief Returns the effective pressure of a given phase within
334 * the control volume.
335 *
336 * \param phaseIdx The phase index
337 */
338 Scalar pressure(const int phaseIdx) const
339 25466400 { return fluidState_.pressure(phaseIdx); }
340
341 /*!
342 * \brief Returns temperature inside the sub-control volume.
343 *
344 * Note that we assume thermodynamic equilibrium, i.e. the
345 * temperatures of the rock matrix and of all fluid phases are
346 * identical.
347 */
348 Scalar temperature(const int phaseIdx = 0) const
349
20/20
✓ Branch 0 taken 2059607 times.
✓ Branch 1 taken 378 times.
✓ Branch 2 taken 2059607 times.
✓ Branch 3 taken 378 times.
✓ Branch 4 taken 1062090 times.
✓ Branch 5 taken 2886 times.
✓ Branch 6 taken 1062090 times.
✓ Branch 7 taken 2886 times.
✓ Branch 8 taken 1471397 times.
✓ Branch 9 taken 4209 times.
✓ Branch 10 taken 1471397 times.
✓ Branch 11 taken 4209 times.
✓ Branch 12 taken 1502115 times.
✓ Branch 13 taken 4209 times.
✓ Branch 14 taken 1502115 times.
✓ Branch 15 taken 4209 times.
✓ Branch 16 taken 992501 times.
✓ Branch 17 taken 2508 times.
✓ Branch 18 taken 992501 times.
✓ Branch 19 taken 2508 times.
23018400 { return fluidState_.temperature(phaseIdx); }
350
351 /*!
352 * \brief Returns the effective mobility of a given phase within
353 * the control volume.
354 *
355 * \param phaseIdx The phase index
356 */
357 Scalar mobility(const int phaseIdx) const
358 {
359 61711200 return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
360 }
361
362 /*!
363 * \brief Returns the effective capillary pressure within the control volume
364 * in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
365 */
366 Scalar capillaryPressure() const
367 { return pc_; }
368
369 /*!
370 * \brief Returns the average porosity within the control volume.
371 */
372 Scalar porosity() const
373 16245198 { return solidState_.porosity(); }
374
375 /*!
376 * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
377 */
378 const PermeabilityType& permeability() const
379 7837200 { return permeability_; }
380
381 /*!
382 * \brief Returns the vapor temperature \f$T_{vap}(p_n)\f$ of the fluid within the control volume.
383 */
384 Scalar vaporTemperature() const
385
2/4
✓ Branch 1 taken 47 times.
✓ Branch 2 taken 97261 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
97308 { return FluidSystem::vaporTemperature(fluidState_, liquidPhaseIdx);}
386
387 /*!
388 * \brief Returns the wetting phase index
389 */
390 int wettingPhase() const
391 { return fluidState_.wettingPhase(); }
392
393 protected:
394 FluidState fluidState_;
395 SolidState solidState_;
396
397 private:
398 Scalar pc_; // The capillary pressure
399 PermeabilityType permeability_; // Effective permeability within the control volume
400
401 // Relative permeability within the control volume
402 std::array<Scalar, numFluidPhases> relativePermeability_;
403 };
404
405 } // end namespace Dumux
406
407 #endif
408