3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
fvpressure.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 *****************************************************************************/
24#ifndef DUMUX_FVPRESSURE2P2C_HH
25#define DUMUX_FVPRESSURE2P2C_HH
26
27// dune environent:
28#include <dune/istl/bvector.hh>
29#include <dune/istl/operators.hh>
30#include <dune/istl/solvers.hh>
31#include <dune/istl/preconditioners.hh>
32#include <dune/common/float_cmp.hh>
33
34// dumux environment
35#include <dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh>
36#include <dumux/material/constraintsolvers/compositionalflash.hh>
37#include <dumux/common/math.hh>
38#include <dumux/io/vtkmultiwriter.hh>
39#include <dumux/porousmediumflow/2p2c/sequential/properties.hh>
40
41#include <dumux/common/deprecated.hh>
42
43namespace Dumux {
73template<class TypeTag> class FVPressure2P2C
74: public FVPressureCompositional<TypeTag>
75{
76 //the model implementation
77 using Implementation = GetPropType<TypeTag, Properties::PressureModel>;
78
79 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
80 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
81 using Problem = GetPropType<TypeTag, Properties::Problem>;
82
83 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
84
85 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
86 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
87
88 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
89 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
90
91 using CellData = GetPropType<TypeTag, Properties::CellData>;
92 enum
93 {
94 dim = GridView::dimension, dimWorld = GridView::dimensionworld
95 };
96 enum
97 {
98 pw = Indices::pressureW,
99 pn = Indices::pressureN,
100 pGlobal = Indices::pressureGlobal
101 };
102 enum
103 {
104 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
105 wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
106 contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
107 };
108
115 enum
116 {
117 rhs = 1,
118 matrix = 0
119
120 };
121
122 // using declarations to abbreviate several dune classes...
123 using Element = typename GridView::Traits::template Codim<0>::Entity;
124 using Intersection = typename GridView::Intersection;
125
126 // convenience shortcuts for Vectors/Matrices
127 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
129 using PhaseVector = Dune::FieldVector<Scalar, getPropValue<TypeTag, Properties::NumPhases>()>;
130 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
131
132 // the typenames used for the stiffness matrix and solution vector
133 using Matrix = GetPropType<TypeTag, Properties::PressureCoefficientMatrix>;
134
135protected:
137 using EntryType = Dune::FieldVector<Scalar, 2>;
138
139 Problem& problem()
140 {
141 return problem_;
142 }
143 const Problem& problem() const
144 {
145 return problem_;
146 }
147
148public:
149
150 void getSource(EntryType& sourceEntry, const Element& elementI, const CellData& cellDataI, const bool first);
151
152 void getStorage(EntryType& storageEntry, const Element& elementI, const CellData& cellDataI, const bool first);
153
154 void getFlux(EntryType& entries, const Intersection& intersection, const CellData& cellDataI, const bool first);
155
156 void getFluxOnBoundary(EntryType& entries, const Intersection& intersection, const CellData& cellDataI, const bool first);
157
158 //updates secondary variables for one cell and stores in the variables object
159 void updateMaterialLawsInElement(const Element& elementI, bool postTimeStep);
160
165 FVPressure2P2C(Problem& problem) : FVPressureCompositional<TypeTag>(problem),
166 problem_(problem)
167 {
168 ErrorTermFactor_ = getParam<Scalar>("Impet.ErrorTermFactor");
169 ErrorTermLowerBound_ = getParam<Scalar>("Impet.ErrorTermLowerBound", 0.2);
170 ErrorTermUpperBound_ = getParam<Scalar>("Impet.ErrorTermUpperBound");
171
172 enableVolumeIntegral = getParam<bool>("Impet.EnableVolumeIntegral");
173 regulateBoundaryPermeability = getPropValue<TypeTag, Properties::RegulateBoundaryPermeability>();
174 if(regulateBoundaryPermeability)
175 {
176 minimalBoundaryPermeability = getParam<Scalar>("SpatialParams.MinBoundaryPermeability");
177 Dune::dinfo << " Warning: Regulating Boundary Permeability requires correct subface indices on reference elements!"
178 << std::endl;
179 }
180 }
181
182protected:
183 Problem& problem_;
184 bool enableVolumeIntegral;
185 bool regulateBoundaryPermeability;
186 Scalar minimalBoundaryPermeability;
187 Scalar ErrorTermFactor_;
188 Scalar ErrorTermLowerBound_;
189 Scalar ErrorTermUpperBound_;
191 static constexpr int pressureType = getPropValue<TypeTag, Properties::PressureFormulation>();
192private:
194 Implementation &asImp_()
195 { return *static_cast<Implementation *>(this);}
196
198 const Implementation &asImp_() const
199 { return *static_cast<const Implementation *>(this);}
200};
201
213template<class TypeTag>
214void FVPressure2P2C<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& sourceEntry,
215 const Element& elementI,
216 const CellData& cellDataI,
217 const bool first)
218{
219 sourceEntry=0.;
220 // cell volume & perimeter, assume linear map here
221 Scalar volume = elementI.geometry().volume();
222
223 // get sources from problem
224 PrimaryVariables source(NAN);
225 problem().source(source, elementI);
226
227 if(first)
228 {
229 source[contiWEqIdx] /= cellDataI.density(wPhaseIdx);
230 source[contiNEqIdx] /= cellDataI.density(nPhaseIdx);
231 }
232 else
233 {
234 // get global coordinate of cell center
235 const GlobalPosition& globalPos = elementI.geometry().center();
236 // derivatives of the fluid volume with respect to concentration of components, or pressure
237 if (!cellDataI.hasVolumeDerivatives())
238 asImp_().volumeDerivatives(globalPos, elementI);
239
240 source[contiWEqIdx] *= cellDataI.dv(wCompIdx); // dV_[i][1] = dv_dC1 = dV/dm1
241 source[contiNEqIdx] *= cellDataI.dv(nCompIdx);
242 }
243 sourceEntry[rhs] = volume * (source[contiWEqIdx] + source[contiNEqIdx]);
244
245 return;
246}
247
263template<class TypeTag>
264void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEntry,
265 const Element& elementI,
266 const CellData& cellDataI,
267 const bool first)
268{
269 storageEntry = 0.;
270 // cell index
271 int eIdxGlobalI = problem().variables().index(elementI);
272 Scalar volume = elementI.geometry().volume();
273
274 // determine maximum error to scale error-term
275 Scalar timestep_ = problem().timeManager().timeStepSize();
276
277 // compressibility term
278 if (!first && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(timestep_, 0.0, 1.0e-30))
279 {
280 Scalar compress_term = cellDataI.dv_dp() / timestep_;
281
282 storageEntry[matrix] -= compress_term*volume;
283 // cellData has data from last TS, and pressurType points to
284 // the pressure Index used as a Primary Variable
285 storageEntry[rhs] -= cellDataI.pressure(pressureType) * compress_term * volume;
286
287 using std::isnan;
288 using std::isinf;
289 if (isnan(compress_term) || isinf(compress_term))
290 DUNE_THROW(Dune::MathError, "Compressibility term leads to NAN matrix entry at index " << eIdxGlobalI);
291
292 if(!getPropValue<TypeTag, Properties::EnableCompressibility>())
293 DUNE_THROW(Dune::NotImplemented, "Compressibility is switched off???");
294 }
295
296 // Abort error damping if there will be a possibly tiny timestep compared with last one
297 // This might be the case if the episode or simulation comes to an end.
298 if( problem().timeManager().episodeWillBeFinished()
299 || problem().timeManager().willBeFinished())
300 {
301 problem().variables().cellData(eIdxGlobalI).errorCorrection() = 0.;
302 return;
303 }
304
305 // error reduction routine: volumetric error is damped and inserted to right hand side
306 // if damping is not done, the solution method gets unstable!
307 problem().variables().cellData(eIdxGlobalI).volumeError() /= timestep_;
308 Scalar erri = fabs(cellDataI.volumeError());
309 Scalar x_lo = ErrorTermLowerBound_;
310 Scalar x_mi = ErrorTermUpperBound_;
311 Scalar fac = ErrorTermFactor_;
312 Scalar lofac = 0.;
313 Scalar hifac = 1.-x_mi;
314
315 if ((erri*timestep_ > 5e-5) && (erri > x_lo * this->maxError_))
316 {
317 if (erri <= x_mi * this->maxError_)
318 storageEntry[rhs] +=
319 problem().variables().cellData(eIdxGlobalI).errorCorrection() =
320 fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/this->maxError_)
321 * cellDataI.volumeError() * volume;
322 else
323 storageEntry[rhs] +=
324 problem().variables().cellData(eIdxGlobalI).errorCorrection() =
325 fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/this->maxError_)
326 * cellDataI.volumeError() * volume;
327 }
328 else
329 problem().variables().cellData(eIdxGlobalI).errorCorrection() = 0.;
330
331 return;
332}
333
355template<class TypeTag>
356void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
357 const Intersection& intersection,
358 const CellData& cellDataI,
359 const bool first)
360{
361 entries = 0.;
362 auto elementI = intersection.inside();
363 int eIdxGlobalI = problem().variables().index(elementI);
364
365 // get global coordinate of cell center
366 const GlobalPosition& globalPos = elementI.geometry().center();
367
368 // cell volume & perimeter, assume linear map here
369 Scalar volume = elementI.geometry().volume();
370 Scalar perimeter = cellDataI.perimeter();
371//#warning perimeter hack 2D!
372// perimeter = intersection.geometry().volume()*2;
373
374 const GlobalPosition& gravity_ = problem().gravity();
375
376 // get absolute permeability
377 DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(elementI));
378
379 // get mobilities and fractional flow factors
380 Scalar fractionalWI=0, fractionalNWI=0;
381 if (first)
382 {
383 fractionalWI = cellDataI.mobility(wPhaseIdx)
384 / (cellDataI.mobility(wPhaseIdx)+ cellDataI.mobility(nPhaseIdx));
385 fractionalNWI = cellDataI.mobility(nPhaseIdx)
386 / (cellDataI.mobility(wPhaseIdx)+ cellDataI.mobility(nPhaseIdx));
387 }
388
389 // get normal vector
390 const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
391
392 // get face volume
393 Scalar faceArea = intersection.geometry().volume();
394
395 // access neighbor
396 auto neighbor = intersection.outside();
397 int eIdxGlobalJ = problem().variables().index(neighbor);
398 CellData& cellDataJ = problem().variables().cellData(eIdxGlobalJ);
399
400 // gemotry info of neighbor
401 const GlobalPosition& globalPosNeighbor = neighbor.geometry().center();
402
403 // distance vector between barycenters
404 GlobalPosition distVec = globalPosNeighbor - globalPos;
405
406 // compute distance between cell centers
407 Scalar dist = distVec.two_norm();
408
409 GlobalPosition unitDistVec(distVec);
410 unitDistVec /= dist;
411
412 DimMatrix permeabilityJ
413 = problem().spatialParams().intrinsicPermeability(neighbor);
414
415 // compute vectorized permeabilities
416 DimMatrix meanPermeability(0);
417 harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
418
419 Dune::FieldVector<Scalar, dim> permeability(0);
420 meanPermeability.mv(unitDistVec, permeability);
421
422 // get average density for gravity flux
423 Scalar rhoMeanW = 0.5 * (cellDataI.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx));
424 Scalar rhoMeanNW = 0.5 * (cellDataI.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx));
425
426 // reset potential gradients
427 Scalar potentialW = 0;
428 Scalar potentialNW = 0;
429
430 if (first) // if we are at the very first iteration we can't calculate phase potentials
431 {
432 // get fractional flow factors in neigbor
433 Scalar fractionalWJ = cellDataJ.mobility(wPhaseIdx)
434 / (cellDataJ.mobility(wPhaseIdx)+ cellDataJ.mobility(nPhaseIdx));
435 Scalar fractionalNWJ = cellDataJ.mobility(nPhaseIdx)
436 / (cellDataJ.mobility(wPhaseIdx)+ cellDataJ.mobility(nPhaseIdx));
437
438 // perform central weighting
439 Scalar lambda = (cellDataI.mobility(wPhaseIdx) + cellDataJ.mobility(wPhaseIdx)) * 0.5
440 + (cellDataI.mobility(nPhaseIdx) + cellDataJ.mobility(nPhaseIdx)) * 0.5;
441
442 entries[0] = fabs(lambda*faceArea*fabs(permeability*unitOuterNormal)/(dist));
443
444 Scalar factor = (fractionalWI + fractionalWJ) * (rhoMeanW) * 0.5
445 + (fractionalNWI + fractionalNWJ) * (rhoMeanNW) * 0.5;
446 entries[1] = factor * lambda * faceArea * fabs(unitOuterNormal*permeability)
447 * (gravity_ * unitDistVec);
448 }
449 else
450 {
451 // determine volume derivatives
452 if (!cellDataJ.hasVolumeDerivatives())
453 asImp_().volumeDerivatives(globalPosNeighbor, neighbor);
454
455 Scalar dv_dC1 = (cellDataJ.dv(wPhaseIdx)
456 + cellDataI.dv(wPhaseIdx)) / 2; // dV/dm1= dv/dC^1
457 Scalar dv_dC2 = (cellDataJ.dv(nPhaseIdx)
458 + cellDataI.dv(nPhaseIdx)) / 2;
459
460 Scalar graddv_dC1 = (cellDataJ.dv(wPhaseIdx)
461 - cellDataI.dv(wPhaseIdx)) / dist;
462 Scalar graddv_dC2 = (cellDataJ.dv(nPhaseIdx)
463 - cellDataI.dv(nPhaseIdx)) / dist;
464
465// potentialW = problem().variables().potentialWetting(eIdxGlobalI, isIndex);
466// potentialNW = problem().variables().potentialNonwetting(eIdxGlobalI, isIndex);
467//
468// densityW = (potentialW > 0.) ? densityWI : densityWJ;
469// densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
470//
471// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
472// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
473 //jochen: central weighting for gravity term
474 Scalar densityW = rhoMeanW;
475 Scalar densityNW = rhoMeanNW;
476
477 potentialW = (cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx))/dist;
478 potentialNW = (cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx))/dist;
479
480 potentialW += densityW * (unitDistVec * gravity_);
481 potentialNW += densityNW * (unitDistVec * gravity_);
482
483 // initialize convenience shortcuts
484 Scalar lambdaW(0.), lambdaN(0.);
485 Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
486 Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) line 3 analogon dV_w
487
488
489 //do the upwinding of the mobility depending on the phase potentials
490 const CellData* upwindWCellData(0);
491 const CellData* upwindNCellData(0);
492 if (potentialW > 0.)
493 upwindWCellData = &cellDataI;
494 else if (potentialW < 0.)
495 upwindWCellData = &cellDataJ;
496 else
497 {
498 if(cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
499 upwindWCellData = &cellDataI;
500 else if(cellDataJ.isUpwindCell(intersection.indexInOutside(), contiWEqIdx))
501 upwindWCellData = &cellDataJ;
502 //else
503 // upwinding is not done!
504 }
505
506 if (potentialNW > 0.)
507 upwindNCellData = &cellDataI;
508 else if (potentialNW < 0.)
509 upwindNCellData = &cellDataJ;
510 else
511 {
512 if(cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
513 upwindNCellData = &cellDataI;
514 else if(cellDataJ.isUpwindCell(intersection.indexInOutside(), contiNEqIdx))
515 upwindNCellData = &cellDataJ;
516 //else
517 // upwinding is not done!
518 }
519
520 //perform upwinding if desired
521 if(!upwindWCellData or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementI.father() == neighbor.father()))
522 {
523 if (cellDataI.wasRefined() && cellDataJ.wasRefined())
524 {
525 problem().variables().cellData(eIdxGlobalI).setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
526 cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
527 }
528
529 Scalar averagedMassFraction[2];
530 averagedMassFraction[wCompIdx]
531 = harmonicMean(cellDataI.massFraction(wPhaseIdx, wCompIdx), cellDataJ.massFraction(wPhaseIdx, wCompIdx));
532 averagedMassFraction[nCompIdx]
533 = harmonicMean(cellDataI.massFraction(wPhaseIdx, nCompIdx), cellDataJ.massFraction(wPhaseIdx, nCompIdx));
534 Scalar averageDensity = harmonicMean(cellDataI.density(wPhaseIdx), cellDataJ.density(wPhaseIdx));
535
536 //compute means
537 dV_w = dv_dC1 * averagedMassFraction[wCompIdx] + dv_dC2 * averagedMassFraction[nCompIdx];
538 dV_w *= averageDensity;
539 gV_w = graddv_dC1 * averagedMassFraction[wCompIdx] + graddv_dC2 * averagedMassFraction[nCompIdx];
540 gV_w *= averageDensity;
541 lambdaW = harmonicMean(cellDataI.mobility(wPhaseIdx), cellDataJ.mobility(wPhaseIdx));
542 }
543 else //perform upwinding
544 {
545 dV_w = (dv_dC1 * upwindWCellData->massFraction(wPhaseIdx, wCompIdx)
546 + dv_dC2 * upwindWCellData->massFraction(wPhaseIdx, nCompIdx));
547 lambdaW = upwindWCellData->mobility(wPhaseIdx);
548 gV_w = (graddv_dC1 * upwindWCellData->massFraction(wPhaseIdx, wCompIdx)
549 + graddv_dC2 * upwindWCellData->massFraction(wPhaseIdx, nCompIdx));
550 dV_w *= upwindWCellData->density(wPhaseIdx);
551 gV_w *= upwindWCellData->density(wPhaseIdx);
552 }
553
554 if(!upwindNCellData or (cellDataI.wasRefined() && cellDataJ.wasRefined()))
555 {
556 if (cellDataI.wasRefined() && cellDataJ.wasRefined())
557 {
558 problem().variables().cellData(eIdxGlobalI).setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
559 cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
560 }
561 Scalar averagedMassFraction[2];
562 averagedMassFraction[wCompIdx]
563 = harmonicMean(cellDataI.massFraction(nPhaseIdx, wCompIdx), cellDataJ.massFraction(nPhaseIdx, wCompIdx));
564 averagedMassFraction[nCompIdx]
565 = harmonicMean(cellDataI.massFraction(nPhaseIdx, nCompIdx), cellDataJ.massFraction(nPhaseIdx, nCompIdx));
566 Scalar averageDensity = harmonicMean(cellDataI.density(nPhaseIdx), cellDataJ.density(nPhaseIdx));
567
568 //compute means
569 dV_n = dv_dC1 * averagedMassFraction[wCompIdx] + dv_dC2 * averagedMassFraction[nCompIdx];
570 dV_n *= averageDensity;
571 gV_n = graddv_dC1 * averagedMassFraction[wCompIdx] + graddv_dC2 * averagedMassFraction[nCompIdx];
572 gV_n *= averageDensity;
573 lambdaN = harmonicMean(cellDataI.mobility(nPhaseIdx), cellDataJ.mobility(nPhaseIdx));
574 }
575 else
576 {
577 dV_n = (dv_dC1 * upwindNCellData->massFraction(nPhaseIdx, wCompIdx)
578 + dv_dC2 * upwindNCellData->massFraction(nPhaseIdx, nCompIdx));
579 lambdaN = upwindNCellData->mobility(nPhaseIdx);
580 gV_n = (graddv_dC1 * upwindNCellData->massFraction(nPhaseIdx, wCompIdx)
581 + graddv_dC2 * upwindNCellData->massFraction(nPhaseIdx, nCompIdx));
582 dV_n *= upwindNCellData->density(nPhaseIdx);
583 gV_n *= upwindNCellData->density(nPhaseIdx);
584 }
585
586 //calculate current matrix entry
587 entries[matrix] = faceArea * (lambdaW * dV_w + lambdaN * dV_n)
588 * fabs((unitOuterNormal*permeability)/(dist));
589 if(enableVolumeIntegral)
590 entries[matrix] -= volume * faceArea / perimeter * (lambdaW * gV_w + lambdaN * gV_n)
591 * ((unitDistVec*permeability)/(dist)); // = boundary integral - area integral
592
593 //calculate right hand side
594 entries[rhs] = faceArea * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
595 entries[rhs] *= fabs(unitOuterNormal * permeability);
596 if(enableVolumeIntegral)
597 entries[rhs] -= volume * faceArea / perimeter * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n)
598 * (unitDistVec * permeability);
599 entries[rhs] *= (gravity_ * unitDistVec);
600
601 // calculate capillary pressure gradient
602 Scalar pCGradient = (cellDataI.capillaryPressure() - cellDataJ.capillaryPressure()) / dist;
603
604 // include capillary pressure fluxes
605 switch (pressureType)
606 {
607 case pw:
608 {
609 //add capillary pressure term to right hand side
610 entries[rhs] += lambdaN * dV_n * fabs(permeability * unitOuterNormal) * pCGradient * faceArea;
611 if(enableVolumeIntegral)
612 entries[rhs]-= lambdaN * gV_n * (permeability * unitDistVec) * pCGradient * volume * faceArea / perimeter;
613 break;
614 }
615 case pn:
616 {
617 //add capillary pressure term to right hand side
618 entries[rhs] -= lambdaW * dV_w * fabs(permeability * unitOuterNormal) * pCGradient * faceArea;
619 if(enableVolumeIntegral)
620 entries[rhs]+= lambdaW * gV_w * (permeability * unitDistVec) * pCGradient * volume * faceArea / perimeter;
621 break;
622 }
623 }
624 } // end !first
625}
626
646template<class TypeTag>
647void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entries,
648 const Intersection& intersection, const CellData& cellDataI, const bool first)
649{
650 entries = 0.;
651 // get global coordinate of cell center
652 auto elementI = intersection.inside();
653 const GlobalPosition& globalPos = elementI.geometry().center();
654
655 // get normal vector
656 const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
657 // get face volume
658 Scalar faceArea = intersection.geometry().volume();
659
660 // get volume derivatives inside the cell
661 Scalar dv_dC1 = cellDataI.dv(wCompIdx);
662 Scalar dv_dC2 = cellDataI.dv(nCompIdx);
663
664 // center of face in global coordinates
665 const GlobalPosition& globalPosFace = intersection.geometry().center();
666
667 // geometrical information
668 GlobalPosition distVec(globalPosFace - globalPos);
669 Scalar dist = distVec.two_norm();
670 GlobalPosition unitDistVec(distVec);
671 unitDistVec /= dist;
672
673 //get boundary condition for boundary face center
674 BoundaryTypes bcType;
675 problem().boundaryTypes(bcType, intersection);
676
677 // prepare pressure boundary condition
678 PhaseVector pressBC(0.);
679 Scalar pcBound (0.);
680
681 // old material law interface is deprecated: Replace this by
682 // const auto& fluidMatrixInteraction = problem().spatialParams.fluidMatrixInteractionAtPos(elementI.geometry().center());
683 // after the release of 3.3, when the deprecated interface is no longer supported
684 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem().spatialParams(), elementI);
685
686 /********** Dirichlet Boundary *************/
687 if (bcType.isDirichlet(Indices::pressureEqIdx))
688 {
689 // get absolute permeability
690 DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(elementI));
691
692 if(regulateBoundaryPermeability)
693 {
694 int axis = intersection.indexInInside() / 2;
695 if(permeabilityI[axis][axis] < minimalBoundaryPermeability)
696 permeabilityI[axis][axis] = minimalBoundaryPermeability;
697 }
698 const GlobalPosition& gravity_ = problem().gravity();
699
700 //permeability vector at boundary
701 Dune::FieldVector<Scalar, dim> permeability(0);
702 permeabilityI.mv(unitDistVec, permeability);
703
704 // create a fluid state for the boundary
705 FluidState BCfluidState;
706
707 //read boundary values
708 PrimaryVariables primaryVariablesOnBoundary(NAN);
709 problem().dirichlet(primaryVariablesOnBoundary, intersection);
710
711 if (first)
712 {
713
714 Scalar fractionalWI=0, fractionalNWI=0;
715 fractionalWI = cellDataI.mobility(wPhaseIdx)
716 / (cellDataI.mobility(wPhaseIdx)+ cellDataI.mobility(nPhaseIdx));
717 fractionalNWI = cellDataI.mobility(nPhaseIdx)
718 / (cellDataI.mobility(wPhaseIdx)+ cellDataI.mobility(nPhaseIdx));
719
720 Scalar lambda = cellDataI.mobility(wPhaseIdx)+cellDataI.mobility(nPhaseIdx);
721 entries[matrix] += lambda * faceArea * fabs(permeability * unitOuterNormal) / (dist);
722 Scalar pressBoundary = primaryVariablesOnBoundary[Indices::pressureEqIdx];
723 entries[rhs] += lambda * faceArea * pressBoundary * fabs(permeability * unitOuterNormal) / (dist);
724 Scalar rightentry = (fractionalWI * cellDataI.density(wPhaseIdx)
725 + fractionalNWI * cellDataI.density(nPhaseIdx))
726 * lambda * faceArea * fabs(unitOuterNormal * permeability)
727 * ( unitDistVec * gravity_);
728 entries[rhs] -= rightentry;
729 }
730 else //not first
731 {
732 // read boundary values
733 problem().transportModel().evalBoundary(globalPosFace,
734 intersection,
735 BCfluidState,
736 pressBC);
737 pcBound = pressBC[nPhaseIdx] - pressBC[wPhaseIdx];
738
739 // determine fluid properties at the boundary
740 Scalar lambdaWBound = 0.;
741 Scalar lambdaNWBound = 0.;
742
743 Scalar densityWBound =
744 FluidSystem::density(BCfluidState, wPhaseIdx);
745 Scalar densityNWBound =
746 FluidSystem::density(BCfluidState, nPhaseIdx);
747 Scalar viscosityWBound =
748 FluidSystem::viscosity(BCfluidState, wPhaseIdx);
749 Scalar viscosityNWBound =
750 FluidSystem::viscosity(BCfluidState, nPhaseIdx);
751
752 // mobility at the boundary
753 if(getPropValue<TypeTag, Properties::BoundaryMobility>() == Indices::satDependent)
754 {
755 lambdaWBound = BCfluidState.saturation(wPhaseIdx)
756 / viscosityWBound;
757 lambdaNWBound = BCfluidState.saturation(nPhaseIdx)
758 / viscosityNWBound;
759 }
760 else if(getPropValue<TypeTag, Properties::BoundaryMobility>() == Indices::permDependent)
761 {
762 lambdaWBound
763 = fluidMatrixInteraction.krw(BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
764 lambdaNWBound
765 = fluidMatrixInteraction.krn(BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
766 }
767 // get average density
768 Scalar rhoMeanW = 0.5 * (cellDataI.density(wPhaseIdx) + densityWBound);
769 Scalar rhoMeanNW = 0.5 * (cellDataI.density(nPhaseIdx) + densityNWBound);
770
771 Scalar potentialW = 0;
772 Scalar potentialNW = 0;
773 if (!first)
774 {
775// potentialW = problem().variables().potentialWetting(eIdxGlobalI, isIndex);
776// potentialNW = problem().variables().potentialNonwetting(eIdxGlobalI, isIndex);
777//
778// // do potential upwinding according to last potGradient vs Jochen: central weighting
779// densityW = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : densityWBound;
780// densityNW = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : densityNWBound;
781//
782// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
783// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
784 Scalar densityW=rhoMeanW;
785 Scalar densityNW=rhoMeanNW;
786
787 //calculate potential gradient
788 potentialW = (cellDataI.pressure(wPhaseIdx) - pressBC[wPhaseIdx])/dist;
789 potentialNW = (cellDataI.pressure(nPhaseIdx) - pressBC[nPhaseIdx])/dist;
790
791 potentialW += densityW * (unitDistVec * gravity_);
792 potentialNW += densityNW * (unitDistVec * gravity_);
793 } //end !first
794
795 //do the upwinding of the mobility depending on the phase potentials
796 Scalar lambdaW, lambdaNW;
797 Scalar densityW(0.), densityNW(0.);
798 Scalar dV_w, dV_n; // no gV, because dv/dC assumed to be constant at the boundary
799 // => no area integral, only boundary integral
800
801 if (potentialW >= 0.)
802 {
803 densityW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialW, 0.0, 1.0e-30)) ? rhoMeanW : cellDataI.density(wPhaseIdx);
804 dV_w = (dv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
805 + dv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
806 dV_w *= densityW;
807 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialW, 0.0, 1.0e-30)) ? 0.5 * (cellDataI.mobility(wPhaseIdx) + lambdaWBound)
808 : cellDataI.mobility(wPhaseIdx);
809 }
810 else
811 {
812 densityW = densityWBound;
813 dV_w = (dv_dC1 * BCfluidState.massFraction(wPhaseIdx, wCompIdx)
814 + dv_dC2 * BCfluidState.massFraction(wPhaseIdx, nCompIdx));
815 dV_w *= densityW;
816 lambdaW = lambdaWBound;
817 }
818 if (potentialNW >= 0.)
819 {
820 densityNW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNW, 0.0, 1.0e-30)) ? rhoMeanNW : cellDataI.density(nPhaseIdx);
821 dV_n = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
822 + dv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
823 dV_n *= densityNW;
824 lambdaNW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNW, 0.0, 1.0e-30)) ? 0.5 * (cellDataI.mobility(nPhaseIdx) + lambdaNWBound)
825 : cellDataI.mobility(nPhaseIdx);
826 }
827 else
828 {
829 densityNW = densityNWBound;
830 dV_n = (dv_dC1 * BCfluidState.massFraction(nPhaseIdx, wCompIdx)
831 + dv_dC2 * BCfluidState.massFraction(nPhaseIdx, nCompIdx));
832 dV_n *= densityNW;
833 lambdaNW = lambdaNWBound;
834 }
835
836 //calculate current matrix entry
837 Scalar entry = (lambdaW * dV_w + lambdaNW * dV_n)
838 * (fabs(unitOuterNormal * permeability) / dist) * faceArea;
839
840 //calculate right hand side
841 Scalar rightEntry = (lambdaW * densityW * dV_w + lambdaNW * densityNW * dV_n)
842 * fabs(unitOuterNormal * permeability) * (gravity_ * unitDistVec) * faceArea ;
843
844 // include capillary pressure fluxes
845 // calculate capillary pressure gradient
846 Scalar pCGradient = (cellDataI.capillaryPressure() - pcBound) / dist;
847 switch (pressureType)
848 {
849 case pw:
850 {
851 //add capillary pressure term to right hand side
852 rightEntry += lambdaNW * dV_n * pCGradient * fabs(unitOuterNormal * permeability) * faceArea;
853 break;
854 }
855 case pn:
856 {
857 //add capillary pressure term to right hand side
858 rightEntry -= lambdaW * dV_w * pCGradient * fabs(unitOuterNormal * permeability) * faceArea;
859 break;
860 }
861 }
862
863
864 // set diagonal entry and right hand side entry
865 entries[matrix] += entry;
866 entries[rhs] += entry * primaryVariablesOnBoundary[Indices::pressureEqIdx];
867 entries[rhs] -= rightEntry;
868 } //end of if(first) ... else{...
869 } // end dirichlet
870
871 /**********************************
872 * set neumann boundary condition
873 **********************************/
874 else if(bcType.isNeumann(Indices::pressureEqIdx))
875 {
876 PrimaryVariables J(NAN);
877 problem().neumann(J, intersection);
878 if (first)
879 {
880 J[contiWEqIdx] /= cellDataI.density(wPhaseIdx);
881 J[contiNEqIdx] /= cellDataI.density(nPhaseIdx);
882 }
883 else
884 {
885 J[contiWEqIdx] *= dv_dC1;
886 J[contiNEqIdx] *= dv_dC2;
887 }
888
889 entries[rhs] -= (J[contiWEqIdx] + J[contiNEqIdx]) * faceArea;
890 }
891 else
892 DUNE_THROW(Dune::NotImplemented, "Boundary Condition neither Dirichlet nor Neumann!");
893
894 return;
895}
896
907template<class TypeTag>
908void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element, bool postTimeStep)
909{
910 // get global coordinate of cell center
911 GlobalPosition globalPos = element.geometry().center();
912
913 // cell Index and cell data
914 int eIdxGlobal = problem().variables().index(element);
915 CellData& cellData = problem().variables().cellData(eIdxGlobal);
916
917 // acess the fluid state and prepare for manipulation
918 FluidState& fluidState = cellData.manipulateFluidState();
919
920 Scalar temperature_ = problem().temperatureAtPos(globalPos);
921
922 // reset to calculate new timeStep if we are at the end of a time step
923 if(postTimeStep)
924 cellData.reset();
925
926 // get the overall mass of first component Z0 = C^0 / (C^0+C^1) [-]
927 Scalar Z0 = cellData.massConcentration(wCompIdx)
928 / (cellData.massConcentration(wCompIdx)
929 + cellData.massConcentration(nCompIdx));
930
931 // make sure only physical quantities enter flash calculation
932 if(Z0 < 0. || Z0 > 1.)
933 {
934 Dune::dgrave << "Feed mass fraction unphysical: Z0 = " << Z0
935 << " at global Idx " << eIdxGlobal
936 << " , because totalConcentration(wCompIdx) = "
937 << cellData.totalConcentration(wCompIdx)
938 << " and totalConcentration(nCompIdx) = "
939 << cellData.totalConcentration(nCompIdx)<< std::endl;
940 if(Z0 < 0.)
941 {
942 Z0 = 0.;
943 cellData.setTotalConcentration(wCompIdx, 0.);
944 problem().transportModel().totalConcentration(wCompIdx, eIdxGlobal) = 0.;
945 Dune::dgrave << "Regularize totalConcentration(wCompIdx) = "
946 << cellData.totalConcentration(wCompIdx)<< std::endl;
947 }
948 else
949 {
950 Z0 = 1.;
951 cellData.setTotalConcentration(nCompIdx, 0.);
952 problem().transportModel().totalConcentration(nCompIdx,eIdxGlobal) = 0.;
953 Dune::dgrave << "Regularize totalConcentration(eIdxGlobal, nCompIdx) = "
954 << cellData.totalConcentration(nCompIdx)<< std::endl;
955 }
956 }
957
958 PhaseVector pressure;
959 CompositionalFlash<Scalar, FluidSystem> flashSolver;
960
961 // old material law interface is deprecated: Replace this by
962 // const auto& fluidMatrixInteraction = problem().spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
963 // after the release of 3.3, when the deprecated interface is no longer supported
964 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem().spatialParams(), element);
965
966
967 if(getPropValue<TypeTag, Properties::EnableCapillarity>()) // iterate capillary pressure and saturation
968 {
969 unsigned int maxiter = 6;
970 Scalar pc = cellData.capillaryPressure(); // initial guess for pc from last TS
971 // start iteration loop
972 for (unsigned int iter = 0; iter < maxiter; iter++)
973 {
974 switch (pressureType)
975 {
976 case pw:
977 {
978 pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal);
979 pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal) + pc;
980 break;
981 }
982 case pn:
983 {
984 pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal) - pc;
985 pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal);
986 break;
987 }
988 }
989
990 // complete fluid state
991 flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, temperature_);
992
993 // calculate new pc
994 Scalar oldPc = pc;
995 pc = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
996
997 if (fabs(oldPc-pc)<10 && iter != 0)
998 break;
999
1000 if (iter++ == maxiter)
1001 Dune::dinfo << iter << "times iteration of pc was applied at Idx " << eIdxGlobal
1002 << ", pc delta still " << fabs(oldPc-pc) << std::endl;
1003 }
1004 }
1005 else // capillary pressure neglected
1006 {
1007 pressure[wPhaseIdx] = pressure[nPhaseIdx] = asImp_().pressure()[eIdxGlobal];
1008 flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, temperature_);
1009 }
1010
1011 // initialize mobilities
1012 cellData.setMobility(wPhaseIdx, fluidMatrixInteraction.krw(fluidState.saturation(wPhaseIdx))
1013 / cellData.viscosity(wPhaseIdx));
1014 cellData.setMobility(nPhaseIdx, fluidMatrixInteraction.krn(fluidState.saturation(wPhaseIdx))
1015 / cellData.viscosity(nPhaseIdx));
1016
1017 // determine volume mismatch between actual fluid volume and pore volume
1018 Scalar sumConc = (cellData.totalConcentration(wCompIdx)
1019 + cellData.totalConcentration(nCompIdx));
1020 Scalar massw = sumConc * fluidState.phaseMassFraction(wPhaseIdx);
1021 Scalar massn = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
1022
1023 if (Dune::FloatCmp::eq<Scalar>((cellData.density(wPhaseIdx)*cellData.density(nPhaseIdx)), 0))
1024 DUNE_THROW(Dune::MathError, "Sequential2p2c::postProcessUpdate: try to divide by 0 density");
1025 Scalar vol = massw / cellData.density(wPhaseIdx) + massn / cellData.density(nPhaseIdx);
1026 if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30))
1027 {
1028 cellData.volumeError()=(vol - problem().spatialParams().porosity(element));
1029
1030 using std::isnan;
1031 if (isnan(cellData.volumeError()))
1032 {
1033 DUNE_THROW(Dune::MathError, "Sequential2p2c::postProcessUpdate:\n"
1034 << "volErr[" << eIdxGlobal << "] isnan: vol = " << vol
1035 << ", massw = " << massw << ", rho_l = " << cellData.density(wPhaseIdx)
1036 << ", massn = " << massn << ", rho_g = " << cellData.density(nPhaseIdx)
1037 << ", poro = " << problem().spatialParams().porosity(element)
1038 << ", dt = " << problem().timeManager().timeStepSize());
1039 }
1040 }
1041 else
1042 cellData.volumeError()=0.;
1043
1044 return;
1045}
1046
1047}//end namespace Dumux
1048#endif
deprecated.hh
Helpers for deprecation.
math.hh
Define some often used mathematical functions.
vtkmultiwriter.hh
Simplifies writing multi-file VTK datasets.
compositionalflash.hh
Determines the pressures and saturations of all fluid phases given the total mass of all components.
fvpressurecompositional.hh
Base class for compositional pressure equations.
Dumux::harmonicMean
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:68
Dumux::harmonicMeanMatrix
void harmonicMeanMatrix(Dune::FieldMatrix< Scalar, m, n > &K, const Dune::FieldMatrix< Scalar, m, n > &Ki, const Dune::FieldMatrix< Scalar, m, n > &Kj)
Calculate the harmonic mean of a fixed-size matrix.
Definition: math.hh:108
Dumux
Definition: adapt.hh:29
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::pressure
std::string pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::IOName::porosity
std::string porosity() noexcept
I/O name of porosity.
Definition: name.hh:139
Dumux::CompositionalFlash
Flash calculation routines for compositional sequential models.
Definition: compositionalflash.hh:51
Dumux::CompositionalFlash::concentrationFlash2p2c
static void concentrationFlash2p2c(FluidState &fluidState, const Scalar &Z0, const PhaseVector &phasePressure, const Scalar &porosity, const Scalar &temperature)
Definition: compositionalflash.hh:71
Dumux::FVPressure2P2C
The finite volume model for the solution of the compositional pressure equation.
Definition: fvpressure.hh:75
Dumux::FVPressure2P2C::enableVolumeIntegral
bool enableVolumeIntegral
Enables the volume integral of the pressure equation.
Definition: fvpressure.hh:184
Dumux::FVPressure2P2C::ErrorTermLowerBound_
Scalar ErrorTermLowerBound_
Handling of error term: lower bound for error dampening.
Definition: fvpressure.hh:188
Dumux::FVPressure2P2C::getSource
void getSource(EntryType &sourceEntry, const Element &elementI, const CellData &cellDataI, const bool first)
Assembles the source term.
Definition: fvpressure.hh:214
Dumux::FVPressure2P2C::ErrorTermFactor_
Scalar ErrorTermFactor_
Handling of error term: relaxation factor.
Definition: fvpressure.hh:187
Dumux::FVPressure2P2C::getFlux
void getFlux(EntryType &entries, const Intersection &intersection, const CellData &cellDataI, const bool first)
Get flux at an interface between two cells.
Definition: fvpressure.hh:356
Dumux::FVPressure2P2C::getFluxOnBoundary
void getFluxOnBoundary(EntryType &entries, const Intersection &intersection, const CellData &cellDataI, const bool first)
Get flux on Boundary.
Definition: fvpressure.hh:647
Dumux::FVPressure2P2C::problem
const Problem & problem() const
Definition: fvpressure.hh:143
Dumux::FVPressure2P2C::problem_
Problem & problem_
Definition: fvpressure.hh:183
Dumux::FVPressure2P2C::pressureType
static constexpr int pressureType
gives kind of pressure used ( , , )
Definition: fvpressure.hh:191
Dumux::FVPressure2P2C::regulateBoundaryPermeability
bool regulateBoundaryPermeability
Enables regulation of permeability in the direction of a Dirichlet Boundary Condition.
Definition: fvpressure.hh:185
Dumux::FVPressure2P2C::problem
Problem & problem()
Definition: fvpressure.hh:139
Dumux::FVPressure2P2C::minimalBoundaryPermeability
Scalar minimalBoundaryPermeability
Minimal limit for the boundary permeability.
Definition: fvpressure.hh:186
Dumux::FVPressure2P2C::EntryType
Dune::FieldVector< Scalar, 2 > EntryType
Definition: fvpressure.hh:137
Dumux::FVPressure2P2C::getStorage
void getStorage(EntryType &storageEntry, const Element &elementI, const CellData &cellDataI, const bool first)
Assembles the storage term.
Definition: fvpressure.hh:264
Dumux::FVPressure2P2C::ErrorTermUpperBound_
Scalar ErrorTermUpperBound_
Definition: fvpressure.hh:189
Dumux::FVPressure2P2C::FVPressure2P2C
FVPressure2P2C(Problem &problem)
Constructs a FVPressure2P2C object.
Definition: fvpressure.hh:165
Dumux::FVPressure2P2C::updateMaterialLawsInElement
void updateMaterialLawsInElement(const Element &elementI, bool postTimeStep)
Updates secondary variables of one cell.
Definition: fvpressure.hh:908
Dumux::FVPressureCompositional
The finite volume model for the solution of the compositional pressure equation.
Definition: fvpressurecompositional.hh:69
properties.hh
Defines the properties required for the sequential 2p2c models.
Generated by doxygen 1.9.3