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