3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dpressurevelocityadaptive.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2/*****************************************************************************
3 * See the file COPYING for full copying permissions. *
4 * *
5 * This program is free software: you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation, either version 3 of the License, or *
8 * (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU General Public License *
16 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
17 *****************************************************************************/
23#ifndef DUMUX_FVMPFALPRESSUREVELOCITY2P_ADAPTIVE_HH
24#define DUMUX_FVMPFALPRESSUREVELOCITY2P_ADAPTIVE_HH
25
26#include <dune/common/float_cmp.hh>
27#include "3dpressureadaptive.hh"
28#include "3dvelocityadaptive.hh"
29
30namespace Dumux {
31
45template<class TypeTag> class FvMpfaL3dPressureVelocity2pAdaptive: public FvMpfaL3dPressure2pAdaptive<TypeTag>
46{
48 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
49
50 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
51
52 enum
53 {
54 dim = GridView::dimension, dimWorld = GridView::dimensionworld
55 };
56
57 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
58
59 enum
60 {
61 pw = Indices::pressureW,
62 pn = Indices::pressureNw,
63 sw = Indices::saturationW,
64 sn = Indices::saturationNw,
65 wPhaseIdx = Indices::wPhaseIdx,
66 nPhaseIdx = Indices::nPhaseIdx,
67 pressureIdx = Indices::pressureIdx,
68 saturationIdx = Indices::saturationIdx,
69 pressEqIdx = Indices::pressureEqIdx,
70 satEqIdx = Indices::satEqIdx,
71 numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
72 };
73
74 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
75 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
76 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
77 using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
78 using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
79 using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
80 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
81 using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
82 using MaterialLaw = typename SpatialParams::MaterialLaw;
83
84 using InteractionVolume = typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume);
85 using Intersection = typename GridView::Intersection;
86
87 using Element = typename GridView::template Codim<0>::Entity;
88
89 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
90 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
91 using DimVector = Dune::FieldVector<Scalar, dim>;
92
93public:
100 ParentType(problem), problem_(problem), velocity_(problem)
101{
102 density_[wPhaseIdx] = 0.;
103 density_[nPhaseIdx] = 0.;
104 viscosity_[wPhaseIdx] = 0.;
105 viscosity_[nPhaseIdx] = 0.;
106
107 calcVelocityInTransport_ = getParam<bool>("MPFA.CalcVelocityInTransport");
108}
109
110 void calculateVelocity();
111
113 void calculateVelocity(const Intersection&, CellData&);
114 void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
115
118 {
119 this->updateMaterialLaws();
120
121 this->storePressureSolution();
122
125 }
126
132 void initialize(bool solveTwice = true)
133 {
134 const auto element = *problem_.gridView().template begin<0>();
135 FluidState fluidState;
136 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
137 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
138 fluidState.setTemperature(problem_.temperature(element));
139 fluidState.setSaturation(wPhaseIdx, 1.);
140 fluidState.setSaturation(nPhaseIdx, 0.);
141 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
142 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
143 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
144 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
145
147 velocity_.initialize();
148
151
152 return;
153 }
154
160 void update()
161 {
163
166 }
167
174 {
175 return calcVelocityInTransport_;
176 }
177
188 template<class MultiWriter>
189 void addOutputVtkFields(MultiWriter &writer)
190 {
192 velocity_.addOutputVtkFields(writer);
193 }
194
195private:
196 Problem& problem_;
198
199 Scalar density_[numPhases];
200 Scalar viscosity_[numPhases];
201 bool calcVelocityInTransport_;
202
204 static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
206 static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
207};
208// end of template
209
215template<class TypeTag>
217{
218 // run through all vertices
219 for (const auto& vertex : vertices(problem_.gridView()))
220 {
221 int vIdxGlobal = problem_.variables().index(vertex);
222
223 InteractionVolume& interactionVolume = this->interactionVolumes_.interactionVolume(vIdxGlobal);
224
225 // inner interactionvolume
226 if (interactionVolume.isInnerVolume())
227 {
228 // cell index
229 int eIdxGlobal[8];
230 for (int i = 0; i < 8; i++)
231 {
232 eIdxGlobal[i] = problem_.variables().index(interactionVolume.getSubVolumeElement(i));
233 }
234
235 //get the cell Data
236 CellData & cellData1 = problem_.variables().cellData(eIdxGlobal[0]);
237 CellData & cellData2 = problem_.variables().cellData(eIdxGlobal[1]);
238 CellData & cellData3 = problem_.variables().cellData(eIdxGlobal[2]);
239 CellData & cellData4 = problem_.variables().cellData(eIdxGlobal[3]);
240 CellData & cellData5 = problem_.variables().cellData(eIdxGlobal[4]);
241 CellData & cellData6 = problem_.variables().cellData(eIdxGlobal[5]);
242 CellData & cellData7 = problem_.variables().cellData(eIdxGlobal[6]);
243 CellData & cellData8 = problem_.variables().cellData(eIdxGlobal[7]);
244
245 if (!interactionVolume.isHangingNodeVolume())
246 {
247 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume,
248 cellData1, cellData2, cellData3, cellData4,
249 cellData5, cellData6, cellData7, cellData8,
250 this->interactionVolumes_, this->transmissibilityCalculator_);
251 }
252 else
253 {
254 velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume,
255 cellData1, cellData2, cellData3, cellData4,
256 cellData5, cellData6, cellData7, cellData8);
257 }
258 }
259 // at least one face on boundary! (boundary interactionvolume)
260 else
261 {
262 for (int elemIdx = 0; elemIdx < 8; elemIdx++)
263 {
264 if (!interactionVolume.hasSubVolumeElement(elemIdx))
265 {
266 continue;
267 }
268 bool isOutside = false;
269 for (int fIdx = 0; fIdx < dim; fIdx++)
270 {
271 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
272 if (interactionVolume.isOutsideFace(intVolFaceIdx))
273 {
274 isOutside = true;
275 break;
276 }
277 }
278 if (isOutside)
279 {
280 continue;
281 }
282
283 // cell index
284 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(elemIdx));
285 //get the cell Data
286 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
287
288 velocity_.calculateBoundaryInteractionVolumeVelocity(interactionVolume, cellData, elemIdx);
289 }
290 } // end boundaries
291
292 } // end vertex iterator
293
294 return;
295}
296
305template<class TypeTag>
306void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
307{
308 auto elementI = intersection.inside();
309 auto elementJ = intersection.outside();
310
311 int levelI = elementI.level();
312 int levelJ = elementJ.level();
313
314 int eIdxGlobalI = problem_.variables().index(elementI);
315 int eIdxGlobalJ = problem_.variables().index(elementJ);
316
317 CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
318
319 int indexInInside = intersection.indexInInside();
320 int indexInOutside = intersection.indexInOutside();
321
322 Dune::FieldVector<CellData, 8> cellDataTemp;
323
324 if (levelI != levelJ)
325 {
326 DimVector vel(0);
327 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside, vel);
328 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside, vel);
329 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside, 0);
330 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside, 0);
331
332 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside, vel);
333 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside, vel);
334 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside, 0);
335 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside, 0);
336 }
337
338 std::set<int> vIdxGlobal;
339
340 if (levelI >= levelJ)
341 {
342 vIdxGlobal = this->interactionVolumes_.faceVerticeIndices(eIdxGlobalI, indexInInside);
343 }
344 else
345 {
346 vIdxGlobal = this->interactionVolumes_.faceVerticeIndices(eIdxGlobalJ, indexInOutside);
347 }
348
349 std::set<int>::iterator itEnd = vIdxGlobal.end();
350
351 for (std::set<int>::iterator vIdxIt = vIdxGlobal.begin(); vIdxIt != itEnd; ++vIdxIt)
352 {
353 InteractionVolume& interactionVolume = this->interactionVolumes_.interactionVolume(*vIdxIt);
354
355 if (interactionVolume.isInnerVolume())
356 {
357 // cell index
358 std::vector<std::pair<int,int> > localMpfaElemIdx(0);
359
360 int eIdxGlobal[8];
361 for (int i = 0; i < 8; i++)
362 {
363 auto elem = interactionVolume.getSubVolumeElement(i);
364
365 if (interactionVolume.isHangingNodeVolume())
366 {
367 if (elem == elementI)
368 {
369 for (int j = 0; j < 3; j++)
370 {
371 for (int k = 0; k < 8; k++)
372 {
373 if (interactionVolume.getSubVolumeElement(k) == elementJ &&
374 IndexTranslator::getFaceIndexFromSubVolume(i,j)
375 == IndexTranslator::getFaceIndexFromElements(i,k))
376 {
377 std::pair<int,int> localIdx = std::make_pair(i, k);
378 localMpfaElemIdx.push_back(localIdx);
379 }
380 }
381 }
382 }
383 }
384 else
385 {
386 if (localMpfaElemIdx.size() != 1)
387 localMpfaElemIdx.resize(1);
388
389 if (elem == elementI)
390 localMpfaElemIdx[0].first = i;
391 else if (elem == elementJ)
392 localMpfaElemIdx[0].second = i;
393 }
394
395 eIdxGlobal[i] = problem_.variables().index(elem);
396 cellDataTemp[i] = problem_.variables().cellData(eIdxGlobal[i]);
397 }
398
399 int size = localMpfaElemIdx.size();
400
401 // if (size > 1)
402 // std::cout<<"size = "<<size<<"\n";
403
404 for (int i = 0; i < size; i++)
405 {
406 int mpfaFaceIdx = IndexTranslator::getFaceIndexFromElements(localMpfaElemIdx[i].first, localMpfaElemIdx[i].second);
407
408 if (!interactionVolume.isHangingNodeVolume())
409 {
410 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume,
411 cellDataTemp[0], cellDataTemp[1], cellDataTemp[2], cellDataTemp[3],
412 cellDataTemp[4], cellDataTemp[5], cellDataTemp[6], cellDataTemp[7],
413 this->interactionVolumes_, this->transmissibilityCalculator_, mpfaFaceIdx);
414 }
415 else
416 {
417 velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume,
418 cellDataTemp[0], cellDataTemp[1], cellDataTemp[2], cellDataTemp[3],
419 cellDataTemp[4], cellDataTemp[5], cellDataTemp[6], cellDataTemp[7], mpfaFaceIdx);
420 }
421
422 if (levelI >= levelJ)
423 {
424 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
425 cellDataTemp[localMpfaElemIdx[i].first].fluxData().velocity(wPhaseIdx, indexInInside));
426 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
427 cellDataTemp[localMpfaElemIdx[i].first].fluxData().velocity(nPhaseIdx, indexInInside));
428 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
429 cellDataTemp[localMpfaElemIdx[i].first].fluxData().upwindPotential(wPhaseIdx, indexInInside));
430 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
431 cellDataTemp[localMpfaElemIdx[i].first].fluxData().upwindPotential(nPhaseIdx, indexInInside));
432 }
433 if (levelJ >= levelI)
434 {
435 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
436 cellDataTemp[localMpfaElemIdx[i].second].fluxData().velocity(wPhaseIdx, indexInOutside));
437 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
438 cellDataTemp[localMpfaElemIdx[i].second].fluxData().velocity(nPhaseIdx, indexInOutside));
439 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
440 cellDataTemp[localMpfaElemIdx[i].second].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
441 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
442 cellDataTemp[localMpfaElemIdx[i].second].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
443 }
444
445 if (size > 1)
446 {
447 for (int j = 0; j < 8; j++)
448 {
449 cellDataTemp[j] = problem_.variables().cellData(eIdxGlobal[j]);
450 }
451 }
452 }
453 }
454 }
455
456 if (levelI == levelJ)
457 {
458 cellData.fluxData().setVelocityMarker(indexInInside);
459 cellDataJ.fluxData().setVelocityMarker(indexInOutside);
460 }
461 else if (levelI > levelJ)
462 {
463 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside, cellData.fluxData().velocity(wPhaseIdx, indexInInside));
464 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside, cellData.fluxData().velocity(nPhaseIdx, indexInInside));
465 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside, -1*cellData.fluxData().upwindPotential(wPhaseIdx, indexInInside));
466 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside, -1*cellData.fluxData().upwindPotential(nPhaseIdx, indexInInside));
467 }
468 else if (levelJ > levelI)
469 {
470 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside, cellDataJ.fluxData().velocity(wPhaseIdx, indexInOutside));
471 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside, cellDataJ.fluxData().velocity(nPhaseIdx, indexInOutside));
472 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside, -1*cellDataJ.fluxData().upwindPotential(wPhaseIdx, indexInOutside));
473 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside, -1*cellDataJ.fluxData().upwindPotential(nPhaseIdx, indexInOutside));
474 }
475}
476
485template<class TypeTag>
486void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
487{
488 auto element = intersection.inside();
489
490 //get face index
491 int isIndex = intersection.indexInInside();
492
493 //get face normal
494 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
495
496 BoundaryTypes bcType;
497 //get boundary type
498 problem_.boundaryTypes(bcType, intersection);
499 PrimaryVariables boundValues(0.0);
500
501 if (bcType.isDirichlet(pressEqIdx))
502 {
503 problem_.dirichlet(boundValues, intersection);
504
505 // get global coordinates of cell centers
506 const GlobalPosition& globalPosI = element.geometry().center();
507
508 // center of face in global coordinates
509 const GlobalPosition& globalPosJ = intersection.geometry().center();
510
511 // get mobilities and fractional flow factors
512 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
513 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
514
515 // get capillary pressure
516 Scalar pcI = cellData.capillaryPressure();
517
518 // distance vector between barycenters
519 GlobalPosition distVec = globalPosJ - globalPosI;
520
521 // compute distance between cell centers
522 Scalar dist = distVec.two_norm();
523
524 //permeability vector at boundary
525 // compute vectorized permeabilities
526 DimMatrix meanPermeability(0);
527
528 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
529
530 Dune::FieldVector<Scalar, dim> permeability(0);
531 meanPermeability.mv(unitOuterNormal, permeability);
532
533 //determine saturation at the boundary -> if no saturation is known directly at the boundary use the cell saturation
534 Scalar satW = 0;
535 if (bcType.isDirichlet(satEqIdx))
536 {
537 switch (saturationType_)
538 {
539 case sw:
540 {
541 satW = boundValues[saturationIdx];
542 break;
543 }
544 case sn:
545 {
546 satW = 1 - boundValues[saturationIdx];
547 break;
548 }
549 }
550 }
551 else
552 {
553 satW = cellData.saturation(wPhaseIdx);
554 }
555
556 Scalar pressBound = boundValues[pressureIdx];
557 Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
558
559 //determine phase pressures from primary pressure variable
560 Scalar pressWBound = 0;
561 Scalar pressNwBound = 0;
562 if (pressureType_ == pw)
563 {
564 pressWBound = pressBound;
565 pressNwBound = pressBound + pcBound;
566 }
567 else if (pressureType_ == pn)
568 {
569 pressWBound = pressBound - pcBound;
570 pressNwBound = pressBound;
571 }
572
573 Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
574 / viscosity_[wPhaseIdx];
575 Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
576 / viscosity_[nPhaseIdx];
577
578 Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
579 Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
580
581 //calculate potential gradient
582 potentialDiffW = (cellData.pressure(wPhaseIdx) - pressWBound);
583 potentialDiffNw = (cellData.pressure(nPhaseIdx) - pressNwBound);
584
585 potentialDiffW += density_[wPhaseIdx] * (distVec * problem_.gravity());
586 potentialDiffNw += density_[nPhaseIdx] * (distVec * problem_.gravity());
587
588 //store potential gradients for further calculations
589 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, potentialDiffW);
590 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, potentialDiffNw);
591
592 //do the upwinding of the mobility depending on the phase potentials
593 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
594 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
595 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
596 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
597
598
599 Scalar scalarPerm = permeability.two_norm();
600
601 //calculate the gravity term
602 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
603 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
604
605 //calculate unit distVec
606 distVec /= dist;
607 Scalar areaScaling = (unitOuterNormal * distVec);
608 //this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
609 Scalar gravityTermW = (problem_.gravity() * distVec) * density_[wPhaseIdx] * areaScaling;
610 Scalar gravityTermNw = (problem_.gravity() * distVec) * density_[nPhaseIdx] * areaScaling;
611
612 //calculate velocity depending on the pressure used -> use pc = pn - pw
613 switch (pressureType_)
614 {
615 case pw:
616 {
617 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
618 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNw)
619 + 0.5 * (lambdaNwI + lambdaNwBound) * scalarPerm * (pcI - pcBound) / dist;
620 break;
621 }
622 case pn:
623 {
624 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
625 - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
626 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNw);
627 break;
628 }
629 }
630
631 //store velocities
632 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
633 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
634 cellData.fluxData().setVelocityMarker(isIndex);
635
636 } //end dirichlet boundary
637
638 else if (bcType.isNeumann(pressEqIdx))
639 {
640 problem_.neumann(boundValues, intersection);
641
642 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
643 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
644
645 velocityW *= boundValues[wPhaseIdx];
646 velocityNw *= boundValues[nPhaseIdx];
647
648 velocityW /= density_[wPhaseIdx];
649 velocityNw /= density_[nPhaseIdx];
650
651 //store potential gradients for further calculations
652 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
653 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
654
655 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
656 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
657 cellData.fluxData().setVelocityMarker(isIndex);
658 } //end neumann boundary
659 else
660 {
661 DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
662 }
663}
664
665} // end namespace Dumux
666#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:281
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
3-d finite Volume-MPFAL implementation of a two-phase pressure equation on h-adaptive grids
3-d velocity calculation on adaptive grids using a 3-d MPFA L-method.
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag SaturationFormulation
The formulation of the saturation model.
Definition: porousmediumflow/2p/sequential/properties.hh:53
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
Property tag MPFAInteractionVolume
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:96
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
void storePressureSolution()
Globally stores the pressure solution.
Definition: 3dpressure.hh:206
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 3dpressure.hh:385
typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume
Type for storing interaction volume information.
Definition: 3dpressure.hh:155
void updateMaterialLaws()
constitutive functions are initialized and stored in the variables object
Definition: 3dpressure.hh:2481
3-d finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequentia...
Definition: 3dpressureadaptive.hh:78
void update()
updates the pressure field (analog to update function in IMPET)
Definition: 3dpressureadaptive.hh:195
Class for the calculation of 3d velocities from the pressure solution of an IMPES scheme using a grid...
Definition: 3dpressurevelocityadaptive.hh:46
void calculateVelocity()
Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
Definition: 3dpressurevelocityadaptive.hh:216
FvMpfaL3dPressureVelocity2pAdaptive(Problem &problem)
Constructs a FvMpfaL3dPressureVelocity2pAdaptive object.
Definition: 3dpressurevelocityadaptive.hh:99
void updateVelocity()
Function for updating the velocity field if iterations are necessary in the transport solution.
Definition: 3dpressurevelocityadaptive.hh:117
void calculateVelocityOnBoundary(const Intersection &intersection, CellData &cellData)
Calculates the velocity at a boundary.
Definition: 3dpressurevelocityadaptive.hh:486
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: 3dpressurevelocityadaptive.hh:173
void initialize(bool solveTwice=true)
Initializes pressure and velocity.
Definition: 3dpressurevelocityadaptive.hh:132
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: 3dpressurevelocityadaptive.hh:189
void update()
Pressure and velocity update.
Definition: 3dpressurevelocityadaptive.hh:160
Class for calculating 3-d velocities from cell-wise constant pressure values.
Definition: 3dvelocityadaptive.hh:54
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:48
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:212