3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/3pwateroil/volumevariables.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
26#ifndef DUMUX_3P2CNI_VOLUME_VARIABLES_HH
27#define DUMUX_3P2CNI_VOLUME_VARIABLES_HH
28
29#include <vector>
30#include <iostream>
31
32#include <dumux/common/math.hh>
35
39
42
44
46
47namespace Dumux {
48
49namespace Detail {
50// helper struct and function detecting if the fluid matrix interaction features a adsorptionModel() function
51#ifndef DOXYGEN // hide from doxygen
52template <class FluidMatrixInteraction>
53using AdsorptionModelDetector = decltype(std::declval<FluidMatrixInteraction>().adsorptionModel());
54#endif // DOXYGEN
55
56template<class FluidMatrixInteraction>
57static constexpr bool hasAdsorptionModel()
58{ return Dune::Std::is_detected<AdsorptionModelDetector, FluidMatrixInteraction>::value; }
59
60}
61
67template <class Traits>
70, public EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >
71{
74 using Scalar = typename Traits::PrimaryVariables::value_type;
75 using ModelTraits = typename Traits::ModelTraits;
76 using FS = typename Traits::FluidSystem;
77 static constexpr int numFluidComps = ParentType::numFluidComponents();
78
79 enum {
81
82 wCompIdx = FS::wCompIdx,
83 nCompIdx = FS::nCompIdx,
84
85 wPhaseIdx = FS::wPhaseIdx,
86 gPhaseIdx = FS::gPhaseIdx,
87 nPhaseIdx = FS::nPhaseIdx,
88
89 switch1Idx = ModelTraits::Indices::switch1Idx,
90 switch2Idx = ModelTraits::Indices::switch2Idx,
91 pressureIdx = ModelTraits::Indices::pressureIdx
92 };
93
94 // present phases
95 enum {
96 threePhases = ModelTraits::Indices::threePhases,
97 wPhaseOnly = ModelTraits::Indices::wPhaseOnly,
98 gnPhaseOnly = ModelTraits::Indices::gnPhaseOnly,
99 wnPhaseOnly = ModelTraits::Indices::wnPhaseOnly,
100 gPhaseOnly = ModelTraits::Indices::gPhaseOnly,
101 wgPhaseOnly = ModelTraits::Indices::wgPhaseOnly
102 };
103
104 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
105 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
106
107public:
109 using FluidState = typename Traits::FluidState;
111 using FluidSystem = typename Traits::FluidSystem;
113 using Indices = typename ModelTraits::Indices;
115 using SolidState = typename Traits::SolidState;
117 using SolidSystem = typename Traits::SolidSystem;
121 static constexpr bool onlyGasPhaseCanDisappear()
122 { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); }
123
132 template<class ElemSol, class Problem, class Element, class Scv>
133 void update(const ElemSol &elemSol,
134 const Problem &problem,
135 const Element &element,
136 const Scv& scv)
137 {
138 ParentType::update(elemSol, problem, element, scv);
139 const auto& priVars = elemSol[scv.localDofIndex()];
140 const auto phasePresence = priVars.state();
141
142 if constexpr (!onlyGasPhaseCanDisappear())
143 {
144 /* first the saturations */
145 if (phasePresence == threePhases)
146 {
147 sw_ = priVars[switch1Idx];
148 sn_ = priVars[switch2Idx];
149 sg_ = 1. - sw_ - sn_;
150 }
151 else if (phasePresence == wPhaseOnly)
152 {
153 sw_ = 1.;
154 sn_ = 0.;
155 sg_ = 0.;
156 }
157 else if (phasePresence == gnPhaseOnly)
158 {
159 sw_ = 0.;
160 sn_ = priVars[switch1Idx];
161 sg_ = 1. - sn_;
162 }
163 else if (phasePresence == wnPhaseOnly)
164 {
165 sn_ = priVars[switch2Idx];
166 sw_ = 1. - sn_;
167 sg_ = 0.;
168 }
169 else if (phasePresence == gPhaseOnly)
170 {
171 sw_ = 0.;
172 sn_ = 0.;
173 sg_ = 1.;
174 }
175 else if (phasePresence == wgPhaseOnly)
176 {
177 sw_ = priVars[switch1Idx];
178 sn_ = 0.;
179 sg_ = 1. - sw_;
180 }
181 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
182
183 fluidState_.setSaturation(wPhaseIdx, sw_);
184 fluidState_.setSaturation(gPhaseIdx, sg_);
185 fluidState_.setSaturation(nPhaseIdx, sn_);
186
187 // old material law interface is deprecated: Replace this by
188 // const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
189 // after the release of 3.3, when the deprecated interface is no longer supported
190 const auto fluidMatrixInteraction = Deprecated::makePcKrSw<3>(Scalar{}, problem.spatialParams(), element, scv, elemSol);
191
192 // calculate capillary pressures
193 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
194 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
195 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
196
197 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
198 const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
199
200 /* now the pressures */
201 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
202 {
203 pg_ = priVars[pressureIdx];
204 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
205 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
206 }
207 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
208 {
209 pw_ = priVars[pressureIdx];
210 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
211 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
212 }
213 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
214
215 fluidState_.setPressure(wPhaseIdx, pw_);
216 fluidState_.setPressure(gPhaseIdx, pg_);
217 fluidState_.setPressure(nPhaseIdx, pn_);
218
219 /* now the temperature */
220 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
221 {
222 temp_ = priVars[switch1Idx];
223 }
224 else if (phasePresence == threePhases)
225 {
226 // temp from inverse pwsat and pnsat which have to sum up to pg
227 Scalar temp = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); // initial guess
228 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
229 {
230 fluidState_.setTemperature(phaseIdx, temp);
231 }
232 solidState_.setTemperature(temp);
233 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
234 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
235
236 using std::abs;
237 while(abs(defect) > 0.01) // simply a small number chosen ...
238 {
239 Scalar deltaT = 1.e-8 * temp;
240 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
241 {
242 fluidState_.setTemperature(phaseIdx, temp+deltaT);
243 }
244 solidState_.setTemperature(temp+deltaT);
245 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
246 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
247
248 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
249 {
250 fluidState_.setTemperature(phaseIdx, temp-deltaT);
251 }
252 solidState_.setTemperature(temp-deltaT);
253 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
254 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
255
256 temp = temp - defect * 2. * deltaT / (fUp - fDown);
257
258 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
259 {
260 fluidState_.setTemperature(phaseIdx, temp);
261 }
262 solidState_.setTemperature(temp);
263 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
264 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
265 }
266 temp_ = temp;
267 }
268 else if (phasePresence == wgPhaseOnly)
269 {
270 // temp from inverse pwsat
271 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
272 }
273 else if (phasePresence == gnPhaseOnly)
274 {
275 // temp from inverse pnsat
276 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
277 }
278 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
279
280 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
281 {
282 fluidState_.setTemperature(phaseIdx, temp_);
283 }
284 solidState_.setTemperature(temp_);
285
286 // now comes the tricky part: calculate phase composition
287 if (phasePresence == threePhases) {
288
289 // all phases are present, phase compositions are a
290 // result of the the gas <-> liquid equilibrium.
291 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
292 Scalar partPressNAPL = pg_ - partPressH2O;
293
294 Scalar xgn = partPressNAPL/pg_;
295 Scalar xgw = partPressH2O/pg_;
296
297 // Henry
298 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
299 Scalar xww = 1.-xwn;
300
301 // Not yet filled with real numbers for the NAPL phase
302 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
303 Scalar xnn = 1.-xnw;
304
305 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
306 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
307 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
308 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
309 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
310 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
311
312 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
313 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
314 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
315 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
316 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
317 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
318
319 fluidState_.setDensity(wPhaseIdx, rhoW);
320 fluidState_.setDensity(gPhaseIdx, rhoG);
321 fluidState_.setDensity(nPhaseIdx, rhoN);
322 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
323 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
324 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
325 }
326 else if (phasePresence == wPhaseOnly) {
327 // only the water phase is present, water phase composition is
328 // stored explicitly.
329
330 // extract mole fractions in the water phase
331 Scalar xwn = priVars[switch2Idx];
332 Scalar xww = 1 - xwn;
333
334 // write water mole fractions in the fluid state
335 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
336 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
337
338 // note that the gas phase is actually not existing!
339 // thus, this is used as phase switch criterion
340 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
341 Scalar xgw = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
342
343
344 // note that the NAPL phase is actually not existing!
345 // thus, this is used as phase switch criterion
346 // maybe solubility would be better than this approach via Henry
347 Scalar xnn = xwn * FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
348 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
349
350 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
351 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
352 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
353 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
354
355 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
356 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
357 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
358 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
359 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
360 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
361
362 fluidState_.setDensity(wPhaseIdx, rhoW);
363 fluidState_.setDensity(gPhaseIdx, rhoG);
364 fluidState_.setDensity(nPhaseIdx, rhoN);
365 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
366 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
367 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
368 }
369 else if (phasePresence == gnPhaseOnly) {
370
371 // only gas and NAPL phases are present
372
373 Scalar xnw = priVars[switch2Idx];
374 Scalar xnn = 1.-xnw;
375 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
376 Scalar xgw = 1.-xgn;
377
378 // note that the water phase is actually not present
379 // the values are used as switching criteria
380 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
381
382 // write mole fractions in the fluid state
383 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
384 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
385 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
386 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
387 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
388
389 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
390 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
391 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
392 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
393 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
394 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
395
396 fluidState_.setDensity(wPhaseIdx, rhoW);
397 fluidState_.setDensity(gPhaseIdx, rhoG);
398 fluidState_.setDensity(nPhaseIdx, rhoN);
399 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
400 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
401 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
402
403 }
404 else if (phasePresence == wnPhaseOnly) {
405 // water and NAPL are present, phase compositions are a
406 // mole fractions of non-existing gas phase are used as switching criteria
407 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
408 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
409
410 Scalar xgn = partPressNAPL/pg_;
411 Scalar xgw = partPressH2O/pg_;
412
413 // Henry
414 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
415 Scalar xww = 1.-xwn;
416
417 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
418 Scalar xnn = 1.-xnw;
419
420 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
421 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
422 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
423 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
424 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
425 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
426
427 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
428 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
429 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
430 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
431 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
432 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
433
434 fluidState_.setDensity(wPhaseIdx, rhoW);
435 fluidState_.setDensity(gPhaseIdx, rhoG);
436 fluidState_.setDensity(nPhaseIdx, rhoN);
437 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
438 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
439 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
440 }
441 else if (phasePresence == gPhaseOnly) {
442 // only the gas phase is present, gas phase composition is
443 // stored explicitly here below.
444
445 const Scalar xgn = priVars[switch2Idx];
446 Scalar xgw = 1 - xgn;
447
448 // write mole fractions in the fluid state
449 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
450 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
451
452 // note that the water and NAPL phase is actually not present
453 // the values are used as switching criteria
454 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
455 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
456
457 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
458 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
459
460 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
461 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
462 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
463 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
464 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
465 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
466
467 fluidState_.setDensity(wPhaseIdx, rhoW);
468 fluidState_.setDensity(gPhaseIdx, rhoG);
469 fluidState_.setDensity(nPhaseIdx, rhoN);
470 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
471 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
472 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
473 }
474 else if (phasePresence == wgPhaseOnly) {
475 // only water and gas phases are present
476 const Scalar xgn = priVars[switch2Idx];
477 Scalar xgw = 1 - xgn;
478
479 // write mole fractions in the fluid state
480 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
481 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
482
483
484 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
485 Scalar xww = 1.-xwn;
486
487 // write mole fractions in the fluid state
488 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
489 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
490
491 // note that the NAPL phase is actually not existing!
492 // thus, this is used as phase switch criterion
493 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
494
495 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
496
497 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
498 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
499 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
500 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
501 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
502 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
503
504 fluidState_.setDensity(wPhaseIdx, rhoW);
505 fluidState_.setDensity(gPhaseIdx, rhoG);
506 fluidState_.setDensity(nPhaseIdx, rhoN);
507 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
508 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
509 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
510 }
511 else
512 assert(false); // unhandled phase state
513 } // end of if(!UseSimpleModel), i.e. the more complex version with six phase states
514
515 else // use the simpler model with only two phase states
516 {
517 /* first the saturations */
518 if (phasePresence == threePhases)
519 {
520 sw_ = priVars[switch1Idx];
521 sn_ = priVars[switch2Idx];
522 sg_ = 1. - sw_ - sn_;
523 }
524 else if (phasePresence == wnPhaseOnly)
525 {
526 sn_ = priVars[switch2Idx];
527 sw_ = 1. - sn_;
528 sg_ = 0.;
529 }
530 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
531
532 fluidState_.setSaturation(wPhaseIdx, sw_);
533 fluidState_.setSaturation(gPhaseIdx, sg_);
534 fluidState_.setSaturation(nPhaseIdx, sn_);
535
536 // old material law interface is deprecated: Replace this by
537 // const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
538 // after the release of 3.3, when the deprecated interface is no longer supported
539 const auto fluidMatrixInteraction = Deprecated::makePcKrSw<3>(Scalar{}, problem.spatialParams(), element, scv, elemSol);
540
541 // calculate capillary pressures
542 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
543 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
544 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
545
546 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
547 const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
548
549 /* now the pressures */
550 if (phasePresence == threePhases)
551 {
552 pg_ = priVars[pressureIdx];
553 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
554 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
555 }
556 else if (phasePresence == wnPhaseOnly)
557 {
558 pw_ = priVars[pressureIdx];
559 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
560 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
561 }
562 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
563
564 fluidState_.setPressure(wPhaseIdx, pw_);
565 fluidState_.setPressure(gPhaseIdx, pg_);
566 fluidState_.setPressure(nPhaseIdx, pn_);
567
568 /* now the temperature */
569 if (phasePresence == wnPhaseOnly)
570 {
571 temp_ = priVars[switch1Idx];
572 }
573 else if (phasePresence == threePhases)
574 {
575 if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
576 {
577 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
578 temp_ = tempOnlyWater;
579 }
580 if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
581 {
582 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
583 temp_ = tempOnlyNAPL;
584 }
585 else
586 {
587 // temp from inverse pwsat and pnsat which have to sum up to pg
588 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
589 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
590 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
591 {
592 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
593 }
594 solidState_.setTemperature(tempOnlyWater);
595 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
596 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
597
598 Scalar temp = tempOnlyWater; // initial guess
599 int counter = 0;
600 using std::abs;
601 while(abs(defect) > 0.01) // simply a small number chosen ...
602 {
603 Scalar deltaT = 1.e-6; // fixed number, but T should always be in the order of a few hundred Kelvin
604 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
605 {
606 fluidState_.setTemperature(phaseIdx, temp+deltaT);
607 }
608 solidState_.setTemperature(temp+deltaT);
609 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
610 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
611
612 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
613 {
614 fluidState_.setTemperature(phaseIdx, temp-deltaT);
615 }
616 solidState_.setTemperature(temp-deltaT);
617 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
618 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
619
620 temp = temp - defect * 2. * deltaT / (fUp - fDown);
621
622 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
623 {
624 fluidState_.setTemperature(phaseIdx, temp);
625 }
626 solidState_.setTemperature(temp);
627 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
628 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
629 counter +=1;
630 if (counter>10) break;
631 }
632 if ((sw_>1.e-10)&&(sw_<0.01))
633 temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
634 if ((sn_>1.e-10)&&(sn_<0.01))
635 temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
636 temp_ = temp;
637 }
638 }
639 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
640
641 for (int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
642 fluidState_.setTemperature(phaseIdx, temp_);
643
644 solidState_.setTemperature(temp_);
645
646 // now comes the tricky part: calculate phase composition
647 if (phasePresence == threePhases) {
648
649 // all phases are present, phase compositions are a
650 // result of the the gas <-> liquid equilibrium.
651 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
652 Scalar partPressNAPL = pg_ - partPressH2O;
653 // regularized evaporation for small liquid phase saturations
654 // avoids negative saturations of liquid phases
655 if (sw_<0.02) partPressH2O *= sw_/0.02;
656 if (partPressH2O < 0.) partPressH2O = 0;
657 if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
658 if (partPressNAPL < 0.) partPressNAPL = 0;
659
660 Scalar xgn = partPressNAPL/pg_;
661 Scalar xgw = partPressH2O/pg_;
662
663 // Immiscible liquid phases, mole fractions are just dummy values
664 Scalar xwn = 0;
665 Scalar xww = 1.-xwn;
666
667 // Not yet filled with real numbers for the NAPL phase
668 Scalar xnw = 0;
669 Scalar xnn = 1.-xnw;
670
671 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
672 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
673 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
674 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
675 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
676 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
677
678 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
679 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
680 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
681 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
682 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
683 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
684
685 fluidState_.setDensity(wPhaseIdx, rhoW);
686 fluidState_.setDensity(gPhaseIdx, rhoG);
687 fluidState_.setDensity(nPhaseIdx, rhoN);
688 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
689 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
690 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
691 }
692 else if (phasePresence == wnPhaseOnly) {
693 // mole fractions of non-existing gas phase are used as switching criteria
694 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
695 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
696
697 Scalar xgn = partPressNAPL/pg_;
698 Scalar xgw = partPressH2O/pg_;
699
700 // immiscible liquid phases, mole fractions are just dummy values
701 Scalar xwn = 0;
702 Scalar xww = 1.-xwn;
703
704 Scalar xnw = 0;
705 Scalar xnn = 1.-xnw;
706
707 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
708 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
709 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
710 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
711 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
712 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
713
714 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
715 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
716 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
717 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
718 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
719 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
720
721 fluidState_.setDensity(wPhaseIdx, rhoW);
722 fluidState_.setDensity(gPhaseIdx, rhoG);
723 fluidState_.setDensity(nPhaseIdx, rhoN);
724 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
725 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
726 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
727 }
728 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
729 }
730
731 // old material law interface is deprecated: Replace this by
732 // const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
733 // after the release of 3.3, when the deprecated interface is no longer supported
734 const auto fluidMatrixInteraction = Deprecated::makePcKrSw<3>(Scalar{}, problem.spatialParams(), element, scv, elemSol);
735
736 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
737 {
738 // Mobilities
739 const Scalar mu =
741 phaseIdx);
742 fluidState_.setViscosity(phaseIdx,mu);
743
744 const Scalar kr = fluidMatrixInteraction.kr(phaseIdx,
745 fluidState_.saturation(wPhaseIdx),
746 fluidState_.saturation(nPhaseIdx));
747 mobility_[phaseIdx] = kr / mu;
748 }
749
750 // material dependent parameters for NAPL adsorption (only if law is provided)
751 if constexpr (Detail::hasAdsorptionModel<std::decay_t<decltype(fluidMatrixInteraction)>>())
752 bulkDensTimesAdsorpCoeff_ = fluidMatrixInteraction.adsorptionModel().bulkDensTimesAdsorpCoeff();
753
754 // porosity
755 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
756 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
757
758 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
759 {
760 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
761 };
762
763 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
764
765 // permeability
766 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
767
768 fluidState_.setTemperature(temp_);
769 // the enthalpies (internal energies are directly calculated in the fluidstate
770 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
771 {
772 Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
773 fluidState_.setEnthalpy(phaseIdx, h);
774 }
775
776 EnergyVolVars::updateEffectiveThermalConductivity();
777 }
778
782 const FluidState &fluidState() const
783 { return fluidState_; }
784
788 const SolidState &solidState() const
789 { return solidState_; }
790
796 Scalar averageMolarMass(int phaseIdx) const
797 { return fluidState_.averageMolarMass(phaseIdx); }
798
805 Scalar saturation(const int phaseIdx) const
806 { return fluidState_.saturation(phaseIdx); }
807
815 Scalar massFraction(const int phaseIdx, const int compIdx) const
816 { return fluidState_.massFraction(phaseIdx, compIdx); }
817
825 Scalar moleFraction(const int phaseIdx, const int compIdx) const
826 { return fluidState_.moleFraction(phaseIdx, compIdx); }
827
834 Scalar density(const int phaseIdx) const
835 { return fluidState_.density(phaseIdx); }
836
843 Scalar molarDensity(const int phaseIdx) const
844 { return fluidState_.molarDensity(phaseIdx); }
845
852 Scalar pressure(const int phaseIdx) const
853 { return fluidState_.pressure(phaseIdx); }
854
862 Scalar temperature() const
863 { return fluidState_.temperature(/*phaseIdx=*/0); }
864
871 Scalar mobility(const int phaseIdx) const
872 { return mobility_[phaseIdx]; }
873
874 Scalar viscosity(const int phaseIdx) const
875 { return fluidState_.viscosity(phaseIdx); }
876
880 Scalar capillaryPressure() const
881 { return fluidState_.capillaryPressure(); }
882
886 Scalar porosity() const
887 { return solidState_.porosity(); }
888
892 Scalar permeability() const
893 { return permeability_; }
894
895 /*
896 * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
897 */
898 Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
899 {
900 if (phaseIdx != nPhaseIdx)
901 return FluidSystem::diffusionCoefficient(fluidState_, phaseIdx);
902 else
903 return 1.e-10;
904 }
905
909 Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
910 { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
911
916 {
917 if (bulkDensTimesAdsorpCoeff_)
918 return bulkDensTimesAdsorpCoeff_.value();
919 else
920 DUNE_THROW(Dune::NotImplemented, "Your spatialParams do not provide an adsorption model");
921 }
922
929 Scalar internalEnergy(int phaseIdx) const
930 { return fluidState_.internalEnergy(phaseIdx); };
931
938 Scalar enthalpy(int phaseIdx) const
939 { return fluidState_.enthalpy(phaseIdx); };
940
941protected:
944
945private:
946 Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
947
948 Scalar moleFrac_[numPs][numFluidComps];
949 Scalar massFrac_[numPs][numFluidComps];
950
951 Scalar permeability_;
952 Scalar mobility_[numPs];
953 OptionalScalar<Scalar> bulkDensTimesAdsorpCoeff_;
954
956 DiffusionCoefficients effectiveDiffCoeff_;
957
958};
959} // end namespace Dumux
960
961#endif
Some exceptions thrown in DuMux
Helpers for deprecation.
Define some often used mathematical functions.
A wrapper that can either contain a valid Scalar or NaN.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
A central place for various physical constants occuring in some equations.
void updateSolidVolumeFractions(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv, SolidState &solidState, const int solidVolFracOffset)
update the solid volume fractions (inert and reacitve) and set them in the solidstate
Definition: updatesolidvolumefractions.hh:36
Definition: adapt.hh:29
decltype(std::declval< FluidMatrixInteraction >().adsorptionModel()) AdsorptionModelDetector
Definition: porousmediumflow/3p3c/volumevariables.hh:47
static constexpr bool hasAdsorptionModel()
Definition: porousmediumflow/3p3c/volumevariables.hh:50
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
T value() const
Definition: optionalscalar.hh:48
The primary variable switch controlling the phase presence state variable.
Definition: 3pwateroil/primaryvariableswitch.hh:39
Contains the quantities which are are constant within a finite volume in the three-phase,...
Definition: porousmediumflow/3pwateroil/volumevariables.hh:71
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:782
SolidState solidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:943
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:898
Scalar molarDensity(const int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:843
Scalar massFraction(const int phaseIdx, const int compIdx) const
Returns the mass fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:815
typename Traits::FluidState FluidState
The type of the object returned by the fluidState() method.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:109
Scalar permeability() const
Returns the permeability within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:892
typename Traits::FluidSystem FluidSystem
The type of the fluid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:111
Scalar viscosity(const int phaseIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:874
Scalar bulkDensTimesAdsorpCoeff() const
Returns the adsorption information.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:915
typename Traits::SolidState SolidState
Export type of solid state.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:115
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:886
Scalar enthalpy(int phaseIdx) const
Returns the total enthalpy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:938
Scalar mobility(const int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:871
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:880
Scalar moleFraction(const int phaseIdx, const int compIdx) const
Returns the mole fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:825
Scalar saturation(const int phaseIdx) const
Returns the effective saturation of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:805
Scalar density(const int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:834
typename Traits::SolidSystem SolidSystem
Export type of solid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:117
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:133
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the effective diffusion coefficients for a phase in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:909
static constexpr bool onlyGasPhaseCanDisappear()
State if only the gas phase is allowed to disappear.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:121
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:862
Scalar internalEnergy(int phaseIdx) const
Returns the total internal energy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:929
typename ModelTraits::Indices Indices
Export the indices.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:113
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:788
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:796
FluidState fluidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:939
Scalar pressure(const int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:852
Definition: porousmediumflow/nonisothermal/volumevariables.hh:75
The isothermal base class.
Definition: porousmediumflow/volumevariables.hh:40
static constexpr int numFluidComponents()
Return number of components considered by the model.
Definition: porousmediumflow/volumevariables.hh:52
const PrimaryVariables & priVars() const
Returns the vector of primary variables.
Definition: porousmediumflow/volumevariables.hh:76
static constexpr int numFluidPhases()
Return number of phases considered by the model.
Definition: porousmediumflow/volumevariables.hh:50
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/volumevariables.hh:64
Base class for the model specific class which provides access to all volume averaged quantities.
Base class for the model specific class which provides access to all volume averaged quantities.
The primary variable switch for the extended Richards model.