3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
cellcentered/pressureadaptive.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_FVPRESSURE2P_ADAPTIVE_HH
25#define DUMUX_FVPRESSURE2P_ADAPTIVE_HH
26
27#include <dune/common/float_cmp.hh>
28
29// dumux environment
33
34namespace Dumux {
35
44template<class TypeTag> class FVPressure2PAdaptive: public FVPressure2P<TypeTag>
45{
47
51
53
56
58 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
60
61 enum
62 {
63 dim = GridView::dimension, dimWorld = GridView::dimensionworld
64 };
65 enum
66 {
67 pw = Indices::pressureW,
68 pn = Indices::pressureNw,
69 pGlobal = Indices::pressureGlobal,
70 sw = Indices::saturationW,
71 sn = Indices::saturationNw
72 };
73 enum
74 {
75 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = getPropValue<TypeTag, Properties::NumPhases>()
76 };
77
78 using Intersection = typename GridView::Intersection;
79
80 using Element = typename GridView::template Codim<0>::Entity;
81
82 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
83 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
84 using FieldMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
85
86protected:
88 using EntryType = typename ParentType::EntryType;
89 enum
90 {
92 };
94
95public:
96 // Function which calculates the flux entry
97 void getFlux(EntryType&, const Intersection&, const CellData&, const bool);
98
105 {
107
108 if (!compressibility_)
109 {
110 const auto element = *problem_.gridView().template begin<0>();
111 FluidState fluidState;
112 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
113 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
114 fluidState.setTemperature(problem_.temperature(element));
115 fluidState.setSaturation(wPhaseIdx, 1.);
116 fluidState.setSaturation(nPhaseIdx, 0.);
117 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
118 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
119 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
120 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
121 }
122
123 velocity_.initialize();
124 }
125
134 void update()
135 {
136 int gridSize = problem_.gridView().size(0);
137 // update RHS vector, matrix
138 if (problem_.gridAdapt().wasAdapted())
139 {
140 this->A_.setSize(gridSize, gridSize); //
141 this->f_.resize(gridSize);
142 this->pressure().resize(gridSize);
143
144
145 for (int i = 0; i < gridSize; i++)
146 {
147 CellData& cellData = problem_.variables().cellData(i);
148
149 switch (pressureType_)
150 {
151 case pw:
152 this->pressure()[i] = cellData.pressure(wPhaseIdx);
153 break;
154 case pn:
155 this->pressure()[i] = cellData.pressure(nPhaseIdx);
156 break;
157 case pGlobal:
158 this->pressure()[i] = cellData.globalPressure();
159 break;
160 }
161 }
162
164 }
165
166
168
170
171 return;
172 }
173
178 {
179 velocity_.calculateVelocity();
180 }
181
186 {
188
190 }
191
197 template<class MultiWriter>
198 void addOutputVtkFields(MultiWriter &writer)
199 {
201
202 velocity_.addOutputVtkFields(writer);
203 }
204
210 FVPressure2PAdaptive(Problem& problem) :
211 ParentType(problem), problem_(problem), velocity_(problem), gravity_(problem.gravity())
212 {
213 if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pGlobal)
214 {
215 DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
216 }
217 if (pressureType_ == pGlobal && compressibility_)
218 {
219 DUNE_THROW(Dune::NotImplemented, "Compressibility not supported for global pressure!");
220 }
221 if (saturationType_ != sw && saturationType_ != sn)
222 {
223 DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
224 }
225
226 density_[wPhaseIdx] = 0.;
227 density_[nPhaseIdx] = 0.;
228 viscosity_[wPhaseIdx] = 0.;
229 viscosity_[nPhaseIdx] = 0.;
230 }
231
232private:
233 Problem& problem_;
235 const GravityVector& gravity_;
236
237 Scalar density_[numPhases];
238 Scalar viscosity_[numPhases];
239
240 static const bool compressibility_ = getPropValue<TypeTag, Properties::EnableCompressibility>();
242 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
244 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
245};
246
255template<class TypeTag>
256void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection& intersection
257 , const CellData& cellData, const bool first)
258{
259 auto elementI = intersection.inside();
260 auto elementJ = intersection.outside();
261
262 if (elementI.level() == elementJ.level() || dim == 3)
263 {
264 ParentType::getFlux(entry, intersection, cellData, first);
265
266 // add the entry only once in case the VisitFacesOnlyOnce option is enabled!!!
267 if (getPropValue<TypeTag, Properties::VisitFacesOnlyOnce>() && elementI.level() < elementJ.level())
268 {
269 entry = 0.;
270 }
271 }
272 // hanging node situation: neighbor has higher level
273 else if (elementJ.level() == elementI.level() + 1)
274 {
275 const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
276
277 // get global coordinates of cell centers
278 const GlobalPosition& globalPosI = elementI.geometry().center();
279 const GlobalPosition& globalPosJ = elementJ.geometry().center();
280
281 int globalIdxI = problem_.variables().index(elementI);
282 int globalIdxJ = problem_.variables().index(elementJ);
283
284 // get mobilities and fractional flow factors
285 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
286 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
287 Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx);
288 Scalar fractionalNwI = cellData.fracFlowFunc(nPhaseIdx);
289 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
290 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
291 Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx);
292 Scalar fractionalNwJ = cellDataJ.fracFlowFunc(nPhaseIdx);
293
294 // get capillary pressure
295 Scalar pcI = cellData.capillaryPressure();
296 Scalar pcJ = cellDataJ.capillaryPressure();
297
298 // get face index
299 int isIndexI = intersection.indexInInside();
300
301 // get face normal
302 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
303
304 // get face area
305 Scalar faceArea = intersection.geometry().volume();
306
307 // Count number of hanging nodes
308 // not really necessary
309 int globalIdxK = 0;
310 auto elementK = intersection.outside();
311 // We are looking for two things:
312 // IsIndexJ, the index of the interface from the neighbor-cell point of view
313 // GlobalIdxK, the index of the third cell
314 // Intersectioniterator around cell I
315 for (const auto& intersectionI : intersections(problem_.gridView(), elementI))
316 {
317 if (intersectionI.neighbor())
318 {
319 auto neighbor2 = intersectionI.outside();
320
321 // make sure we do not choose elemntI as third element
322 // -> faces with hanging node have more than one intersection but only one face index!
323 if (neighbor2 != elementJ && intersectionI.indexInInside() == isIndexI)
324 {
325 globalIdxK = problem_.variables().index(neighbor2);
326 elementK = neighbor2;
327 break;
328 }
329 }
330 }
331
332 const CellData& cellDataK = problem_.variables().cellData(globalIdxK);
333
334 // neighbor cell center in global coordinates
335 const GlobalPosition& globalPosInterface = intersection.geometry().center();
336
337 Dune::FieldVector < Scalar, dimWorld > distVec = globalPosInterface - globalPosI;
338 Scalar lI = distVec * unitOuterNormal;
339 distVec = globalPosJ - globalPosInterface;
340 Scalar lJ = distVec * unitOuterNormal;
341 Scalar l = lI + lJ;
342
343 FieldMatrix permeabilityI(0);
344 FieldMatrix permeabilityJ(0);
345 FieldMatrix permeabilityK(0);
346
347 problem_.spatialParams().meanK(permeabilityI,
348 problem_.spatialParams().intrinsicPermeability(elementI));
349 problem_.spatialParams().meanK(permeabilityJ,
350 problem_.spatialParams().intrinsicPermeability(elementJ));
351 problem_.spatialParams().meanK(permeabilityK,
352 problem_.spatialParams().intrinsicPermeability(elementK));
353
354 // Calculate permeablity component normal to interface
355 Scalar kI, kJ, kK, kMean, ng;
356 Dune::FieldVector < Scalar, dim > permI(0);
357 Dune::FieldVector < Scalar, dim > permJ(0);
358 Dune::FieldVector < Scalar, dim > permK(0);
359
360 permeabilityI.mv(unitOuterNormal, permI);
361 permeabilityJ.mv(unitOuterNormal, permJ);
362 permeabilityK.mv(unitOuterNormal, permK);
363
364 // kI,kJ,kK=(n^T)Kn
365 kI = unitOuterNormal * permI;
366 kJ = unitOuterNormal * permJ;
367 kK = unitOuterNormal * permK;
368
369 // See Diplomarbeit Michael Sinsbeck
370 kMean = kI * kJ * kK * l / (kJ * kK * lI + kI * (kJ + kK) / 2 * lJ);
371
372 ng = gravity_ * unitOuterNormal;
373
374 Scalar fractionalWK = cellDataK.fracFlowFunc(wPhaseIdx);
375 Scalar fractionalNwK = cellDataK.fracFlowFunc(nPhaseIdx);
376
377 Scalar pcK = cellDataK.capillaryPressure();
378 Scalar pcJK = (pcJ + pcK) / 2;
379
380 // Potentials from previous time step are not available.
381 // Instead, calculate mean density, then find potentials,
382 // then upwind density.
383 // pressure from previous time step might also be incorrect.
384
385 Scalar rhoMeanWIJ = density_[wPhaseIdx];
386 Scalar rhoMeanNwIJ = density_[nPhaseIdx];
387 Scalar rhoMeanWIK = density_[wPhaseIdx];
388 Scalar rhoMeanNwIK = density_[nPhaseIdx];
389
390 if (compressibility_)
391 {
392 rhoMeanWIJ = (lJ * cellData.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l;
393 rhoMeanNwIJ = (lJ * cellData.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l;
394 rhoMeanWIK = (lJ * cellData.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l;
395 rhoMeanNwIK = (lJ * cellData.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l;
396 }
397
398 Scalar densityWIJ = density_[wPhaseIdx];
399 Scalar densityNwIJ = density_[nPhaseIdx];
400 Scalar densityWIK = density_[wPhaseIdx];
401 Scalar densityNwIK = density_[nPhaseIdx];
402
403 Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l;
404 Scalar fMeanNwIJ = (lJ * fractionalNwI + lI * fractionalNwJ) / l;
405 Scalar fMeanWIK = (lJ * fractionalWI + lI * fractionalWK) / l;
406 Scalar fMeanNwIK = (lJ * fractionalNwI + lI * fractionalNwK) / l;
407
408 Scalar potentialWIJ = 0;
409 Scalar potentialNwIJ = 0;
410 Scalar potentialWIK = 0;
411 Scalar potentialNwIK = 0;
412
413 // if we are at the very first iteration we can't calculate phase potentials
414 if (!first)
415 {
416 // potentials from previous time step no available.
417 if (compressibility_)
418 {
419 densityWIJ = rhoMeanWIJ;
420 densityNwIJ = rhoMeanNwIJ;
421 densityWIK = rhoMeanWIK;
422 densityNwIK = rhoMeanNwIK;
423 }
424
425 switch (pressureType_)
426 {
427 case pGlobal:
428 {
429 Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
430
431 potentialWIJ = (cellData.globalPressure() - fMeanNwIJ * pcI - (pressJK - fMeanNwIJ * pcJK)) / l
432 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
433 potentialNwIJ = (cellData.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l
434 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
435 potentialWIK = (cellData.globalPressure() - fMeanNwIK * pcI - (pressJK - fMeanNwIK * pcJK)) / l
436 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
437 potentialNwIK = (cellData.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l
438 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
439 break;
440 }
441 default:
442 {
443 potentialWIJ = (cellData.pressure(wPhaseIdx)
444 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
445 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
446 potentialNwIJ = (cellData.pressure(nPhaseIdx)
447 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
448 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
449 potentialWIK = (cellData.pressure(wPhaseIdx)
450 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
451 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
452 potentialNwIK = (cellData.pressure(nPhaseIdx)
453 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
454 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
455 break;
456 }
457 }
458 }
459
460 // do the upwinding of the mobility depending on the phase potentials
461 Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
462 lambdaWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
463 Scalar lambdaNwIJ = (potentialNwIJ > 0) ? lambdaNwI : lambdaNwJ;
464 lambdaNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
465
466 if (compressibility_)
467 {
468 densityWIJ = (potentialWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
469 densityNwIJ = (potentialNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
470 densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
471 densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
472 densityWIK = (potentialWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
473 densityNwIK = (potentialNwIK > 0.) ? cellData.density(nPhaseIdx) : densityNwIK;
474 densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
475 densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
476 }
477
478 // update diagonal entry and right hand side
479 entry[matrix] = (lambdaWIJ + lambdaNwIJ) * kMean / l * faceArea;
480 entry[rhs] = faceArea * lambdaWIJ * kMean * ng
481 * ((densityWIJ) - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2);
482 entry[rhs] += faceArea * lambdaNwIJ * kMean * ng
483 * ((densityNwIJ) - (lJ / l) * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2);
484
485 switch (pressureType_)
486 {
487 case pw:
488 {
489 entry[rhs] += faceArea * lambdaNwIJ * kMean * (pcJK - pcI) / l;
490 break;
491 }
492 case pn:
493 {
494 entry[rhs] -= faceArea * lambdaWIJ * kMean * (pcJK - pcI) / l;
495 break;
496 }
497 }
498
499 // write hanging-node-specific stuff directly into matrix and rhs!
500 this->f_[globalIdxI] -= entry[rhs];
501 this->f_[globalIdxJ] += entry[rhs];
502
503 // set diagonal entry
504 this->A_[globalIdxI][globalIdxI] += entry[matrix];
505 // set off-diagonal
506 this->A_[globalIdxI][globalIdxJ] -= entry[matrix];
507
508 // set entry for cell J
509 this->A_[globalIdxJ][globalIdxI] -= entry[matrix];
510 this->A_[globalIdxJ][globalIdxJ] += entry[matrix];
511
512 // set entry to zero -> matrix already written!
513 entry = 0.;
514
515// std::cout<<"finished hanging node!\n";
516 }
517 else
518 {
519 entry = 0;
520 }
521
522 return;
523}
524
525} // end namespace Dumux
526#endif
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
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
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Finite Volume discretization of a two-phase flow pressure equation of the sequential IMPES model.
Definition: 2p/sequential/diffusion/cellcentered/pressure.hh:115
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 2p/sequential/diffusion/cellcentered/pressure.hh:416
void updateVelocity()
Velocity update.
Definition: 2p/sequential/diffusion/cellcentered/pressure.hh:308
void update()
Pressure update.
Definition: 2p/sequential/diffusion/cellcentered/pressure.hh:264
Finite volume discretization of a two-phase flow pressure equation of the sequential IMPES model.
Definition: cellcentered/pressureadaptive.hh:45
void initialize()
Initializes the adaptive pressure model.
Definition: cellcentered/pressureadaptive.hh:104
FVPressure2PAdaptive(Problem &problem)
Constructs a FVPressure2PAdaptive object.
Definition: cellcentered/pressureadaptive.hh:210
void update()
Pressure update.
Definition: cellcentered/pressureadaptive.hh:134
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: cellcentered/pressureadaptive.hh:198
void getFlux(EntryType &, const Intersection &, const CellData &, const bool)
Function which calculates the flux entry.
Definition: cellcentered/pressureadaptive.hh:256
void updateVelocity()
Velocity update.
Definition: cellcentered/pressureadaptive.hh:185
void calculateVelocity()
Velocity calculation.
Definition: cellcentered/pressureadaptive.hh:177
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
@ rhs
index for the right hand side entry
Definition: sequential/cellcentered/pressure.hh:88
@ matrix
index for the global matrix entry
Definition: sequential/cellcentered/pressure.hh:89
Matrix A_
Global stiffnes matrix (sparse matrix which is build by the initializeMatrix() function)
Definition: sequential/cellcentered/pressure.hh:324
RHSVector f_
Right hand side vector.
Definition: sequential/cellcentered/pressure.hh:325
Dune::FieldVector< Scalar, 2 > EntryType
Definition: sequential/cellcentered/pressure.hh:78
void initializeMatrix()
Initialize the global matrix of the system of equations to solve.
Definition: sequential/cellcentered/pressure.hh:332
Base class for finite volume velocity reconstruction.
Definition: sequential/cellcentered/velocity.hh:48
Defines the properties required for (immiscible) two-phase sequential models.
Finite Volume discretization of a two-phase flow pressure equation.
Finite volume velocity reconstruction.