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