3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
lmethod/2dvelocity.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_FVMPFAL2DVELOCITY2P_HH
25#define DUMUX_FVMPFAL2DVELOCITY2P_HH
26
27#include <dune/grid/common/gridenums.hh>
28#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
31#include "2dtransmissibilitycalculator.hh"
32
33namespace Dumux {
34
57template<class TypeTag> class FvMpfaL2dVelocity2p
58{
59 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
60 enum
61 {
62 dim = GridView::dimension, dimWorld = GridView::dimensionworld
63 };
64
65 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
66 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
67
68 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
69
70 using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
71 using MaterialLaw = typename SpatialParams::MaterialLaw;
72
73 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
74
75 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
76 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
77
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 CellData = typename GET_PROP_TYPE(TypeTag, CellData);
82
83 using Element = typename GridView::Traits::template Codim<0>::Entity;
84 using Grid = typename GridView::Grid;
85 using IndexSet = typename GridView::IndexSet;
86
87 using Geometry = typename Element::Geometry;
88 using JacobianTransposed = typename Geometry::JacobianTransposed ;
89
90 using GridTypeIndices = typename GET_PROP_TYPE(TypeTag, GridTypeIndices);
91
92 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
93 using TransmissibilityCalculator = FvMpfaL2dTransmissibilityCalculator<TypeTag>;
94 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
95
96 enum
97 {
98 pw = Indices::pressureW,
99 pn = Indices::pressureNw,
100 pGlobal = Indices::pressureGlobal,
101 sw = Indices::saturationW,
102 sn = Indices::saturationNw,
103 vw = Indices::velocityW,
104 vn = Indices::velocityNw,
105 vt = Indices::velocityTotal
106 };
107 enum
108 {
109 wPhaseIdx = Indices::wPhaseIdx,
110 nPhaseIdx = Indices::nPhaseIdx,
111 pressureIdx = Indices::pressureIdx,
112 saturationIdx = Indices::saturationIdx,
113 pressureEqIdx = Indices::pressureEqIdx,
114 satEqIdx = Indices::satEqIdx,
115 numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
116 };
117
118 using LocalPosition = Dune::FieldVector<Scalar, dim>;
119 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
120 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
121 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
122 using DimVector = Dune::FieldVector<Scalar, dim>;
123
124public:
130 FvMpfaL2dVelocity2p(Problem& problem) :
131 problem_(problem), transmissibilityCalculator_(problem), gravity_(problem.gravity())
132 {
133 density_[wPhaseIdx] = 0.;
134 density_[nPhaseIdx] = 0.;
135 viscosity_[wPhaseIdx] = 0.;
136 viscosity_[nPhaseIdx] = 0.;
137
138 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
139 }
140
142 void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
143 CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4,
144 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
145 void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
146 CellData& cellData, int elemIdx);
147
149 void initialize()
150 {
151 const auto element = *problem_.gridView().template begin<0>();
152 FluidState fluidState;
153 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
154 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
155 fluidState.setTemperature(problem_.temperature(element));
156 fluidState.setSaturation(wPhaseIdx, 1.);
157 fluidState.setSaturation(nPhaseIdx, 0.);
158 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
159 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
160 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
161 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
162
163 return;
164 }
165
176 template<class MultiWriter>
177 void addOutputVtkFields(MultiWriter &writer)
178 {
179 if (vtkOutputLevel_ > 0)
180 {
181 Dune::BlockVector < DimVector > &velocityWetting
182 = *(writer.template allocateManagedBuffer<Scalar,dim>(problem_.gridView().size(0)));
183 Dune::BlockVector < DimVector > &velocityNonwetting
184 = *(writer.template allocateManagedBuffer<Scalar,dim>(problem_.gridView().size(0)));
185
186 // compute update vector
187 for (const auto& element : elements(problem_.gridView()))
188 {
189 // cell index
190 int eIdxGlobal = problem_.variables().index(element);
191
192 CellData & cellData = problem_.variables().cellData(eIdxGlobal);
193
194 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
195 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
196
197 // run through all intersections with neighbors and boundary
198 for (const auto& intersection : intersections(problem_.gridView(), element))
199 {
200 int isIndex = intersection.indexInInside();
201
202 fluxW[isIndex] += intersection.geometry().volume()
203 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
204 fluxNw[isIndex] += intersection.geometry().volume()
205 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
206 }
207
208 DimVector refVelocity(0);
209 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
210 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
211
212 const DimVector& localPos = ReferenceElements::general(element.type()).position(0, 0);
213
214 // get the transposed Jacobian of the element mapping
215 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
216
217 // calculate the element velocity by the Piola transformation
218 DimVector elementVelocity(0);
219 jacobianT.umtv(refVelocity, elementVelocity);
220 elementVelocity /= element.geometry().integrationElement(localPos);
221
222 velocityWetting[eIdxGlobal] = elementVelocity;
223
224 refVelocity = 0;
225 refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
226 refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
227
228 // calculate the element velocity by the Piola transformation
229 elementVelocity = 0;
230 jacobianT.umtv(refVelocity, elementVelocity);
231 elementVelocity /= element.geometry().integrationElement(localPos);
232
233 velocityNonwetting[eIdxGlobal] = elementVelocity;
234 }
235
236 writer.attachCellData(velocityWetting, "wetting-velocity", dim);
237 writer.attachCellData(velocityNonwetting, "non-wetting-velocity", dim);
238 }
239
240 return;
241 }
242
243private:
244 Problem& problem_;
245protected:
246 TransmissibilityCalculator transmissibilityCalculator_;
247
248 const GravityVector& gravity_; // vector including the gravity constant
249
250 Scalar density_[numPhases];
251 Scalar viscosity_[numPhases];
252
253 int vtkOutputLevel_;
254
255 static constexpr Scalar threshold_ = 1e-15;
257 static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);
259 static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
261 static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
262};
263// end of template
264
278template<class TypeTag>
279void FvMpfaL2dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
280 CellData& cellData1, CellData& cellData2,
281 CellData& cellData3, CellData& cellData4,
282 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
283{
284 auto element1 = interactionVolume.getSubVolumeElement(0);
285 auto element2 = interactionVolume.getSubVolumeElement(1);
286 auto element3 = interactionVolume.getSubVolumeElement(2);
287 auto element4 = interactionVolume.getSubVolumeElement(3);
288
289 int level1 = element1.level();
290 int level2 = element2.level();
291 int level3 = element3.level();
292 int level4 = element4.level();
293
294 // cell index
295 int eIdxGlobal1 = problem_.variables().index(element1);
296 int eIdxGlobal2 = problem_.variables().index(element2);
297 int eIdxGlobal3 = problem_.variables().index(element3);
298 int eIdxGlobal4 = problem_.variables().index(element4);
299
300 // get pressure values
301 Dune::FieldVector < Scalar, 2 * dim > potW(0);
302 Dune::FieldVector < Scalar, 2 * dim > potNw(0);
303
304 potW[0] = cellData1.potential(wPhaseIdx);
305 potW[1] = cellData2.potential(wPhaseIdx);
306 potW[2] = cellData3.potential(wPhaseIdx);
307 potW[3] = cellData4.potential(wPhaseIdx);
308
309 potNw[0] = cellData1.potential(nPhaseIdx);
310 potNw[1] = cellData2.potential(nPhaseIdx);
311 potNw[2] = cellData3.potential(nPhaseIdx);
312 potNw[3] = cellData4.potential(nPhaseIdx);
313
314 //get mobilities of the phases
315 Dune::FieldVector < Scalar, numPhases > lambda1(cellData1.mobility(wPhaseIdx));
316 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
317
318 //compute total mobility of cell 1
319 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
320
321 //get mobilities of the phases
322 Dune::FieldVector < Scalar, numPhases > lambda2(cellData2.mobility(wPhaseIdx));
323 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
324
325 //compute total mobility of cell 1
326 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
327
328 //get mobilities of the phases
329 Dune::FieldVector < Scalar, numPhases > lambda3(cellData3.mobility(wPhaseIdx));
330 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
331
332 //compute total mobility of cell 1
333 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
334
335 //get mobilities of the phases
336 Dune::FieldVector < Scalar, numPhases > lambda4(cellData4.mobility(wPhaseIdx));
337 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
338
339 //compute total mobility of cell 1
340 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
341
342
343 std::vector < DimVector > lambda(2 * dim);
344
345 lambda[0][0] = lambdaTotal1;
346 lambda[0][1] = lambdaTotal1;
347 lambda[1][0] = lambdaTotal2;
348 lambda[1][1] = lambdaTotal2;
349 lambda[2][0] = lambdaTotal3;
350 lambda[2][1] = lambdaTotal3;
351 lambda[3][0] = lambdaTotal4;
352 lambda[3][1] = lambdaTotal4;
353
354 Scalar potentialDiffW12 = 0;
355 Scalar potentialDiffW14 = 0;
356 Scalar potentialDiffW32 = 0;
357 Scalar potentialDiffW34 = 0;
358
359 Scalar potentialDiffNw12 = 0;
360 Scalar potentialDiffNw14 = 0;
361 Scalar potentialDiffNw32 = 0;
362 Scalar potentialDiffNw34 = 0;
363
364 //flux vector
365 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
366 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
367
368 Dune::FieldMatrix < Scalar, dim, 2 * dim - dim + 1 > T(0);
369 DimVector Tu(0);
370 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
371
372 int lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 2, 3);
373
374 if (lType == TransmissibilityCalculator::rightTriangle)
375 {
376 u[0] = potW[1];
377 u[1] = potW[2];
378 u[2] = potW[0];
379
380 T.mv(u, Tu);
381
382 fluxW[0] = Tu[1];
383 potentialDiffW12 = Tu[1];
384
385 u[0] = potNw[1];
386 u[1] = potNw[2];
387 u[2] = potNw[0];
388
389 T.mv(u, Tu);
390
391 fluxNw[0] = Tu[1];
392 potentialDiffNw12 = Tu[1];
393 }
394 else
395 {
396 u[0] = potW[0];
397 u[1] = potW[3];
398 u[2] = potW[1];
399
400 T.mv(u, Tu);
401
402 fluxW[0] = Tu[1];
403 potentialDiffW12 = Tu[1];
404
405 u[0] = potNw[0];
406 u[1] = potNw[3];
407 u[2] = potNw[1];
408
409 T.mv(u, Tu);
410
411 fluxNw[0] = Tu[1];
412 potentialDiffNw12 = Tu[1];
413 }
414
415 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 1, 2, 3, 0);
416
417 if (lType == TransmissibilityCalculator::rightTriangle)
418 {
419 u[0] = potW[2];
420 u[1] = potW[3];
421 u[2] = potW[1];
422
423 T.mv(u, Tu);
424
425 fluxW[1] = Tu[1];
426 potentialDiffW32 = -Tu[1];
427
428 u[0] = potNw[2];
429 u[1] = potNw[3];
430 u[2] = potNw[1];
431
432 T.mv(u, Tu);
433
434 fluxNw[1] = Tu[1];
435 potentialDiffNw32 = -Tu[1];
436 }
437 else
438 {
439 u[0] = potW[1];
440 u[1] = potW[0];
441 u[2] = potW[2];
442
443 T.mv(u, Tu);
444
445 fluxW[1] = Tu[1];
446 potentialDiffW32 = -Tu[1];
447
448 u[0] = potNw[1];
449 u[1] = potNw[0];
450 u[2] = potNw[2];
451
452 T.mv(u, Tu);
453
454 fluxNw[1] = Tu[1];
455 potentialDiffNw32 = -Tu[1];
456 }
457
458 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 2, 3, 0, 1);
459
460 if (lType == TransmissibilityCalculator::rightTriangle)
461 {
462 u[0] = potW[3];
463 u[1] = potW[0];
464 u[2] = potW[2];
465
466 T.mv(u, Tu);
467
468 fluxW[2] = Tu[1];
469 potentialDiffW34 = Tu[1];
470
471 u[0] = potNw[3];
472 u[1] = potNw[0];
473 u[2] = potNw[2];
474
475 T.mv(u, Tu);
476
477 fluxNw[2] = Tu[1];
478 potentialDiffNw34 = Tu[1];
479 }
480 else
481 {
482 u[0] = potW[2];
483 u[1] = potW[1];
484 u[2] = potW[3];
485
486 T.mv(u, Tu);
487
488 fluxW[2] = Tu[1];
489 potentialDiffW34 = Tu[1];
490
491 u[0] = potNw[2];
492 u[1] = potNw[1];
493 u[2] = potNw[3];
494
495 T.mv(u, Tu);
496
497 fluxNw[2] = Tu[1];
498 potentialDiffNw34 = Tu[1];
499 }
500
501 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 3, 0, 1, 2);
502
503 if (lType == TransmissibilityCalculator::rightTriangle)
504 {
505 u[0] = potW[0];
506 u[1] = potW[1];
507 u[2] = potW[3];
508
509 T.mv(u, Tu);
510
511 fluxW[3] = Tu[1];
512 potentialDiffW14 = -Tu[1];
513
514 u[0] = potNw[0];
515 u[1] = potNw[1];
516 u[2] = potNw[3];
517
518 T.mv(u, Tu);
519
520 fluxNw[3] = Tu[1];
521 potentialDiffNw14 = -Tu[1];
522 }
523 else
524 {
525 u[0] = potW[3];
526 u[1] = potW[2];
527 u[2] = potW[0];
528
529 T.mv(u, Tu);
530
531 fluxW[3] = Tu[1];
532 potentialDiffW14 = -Tu[1];
533
534 u[0] = potNw[3];
535 u[1] = potNw[2];
536 u[2] = potNw[0];
537
538 T.mv(u, Tu);
539
540 fluxNw[3] = Tu[1];
541 potentialDiffNw14 = -Tu[1];
542 }
543
544 //store potentials for further calculations (saturation, ...)
545 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW12);
546 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw12);
547 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffW14);
548 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNw14);
549 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffW32);
550 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffNw32);
551 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW12);
552 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw12);
553 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffW34);
554 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNw34);
555 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffW32);
556 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffNw32);
557 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffW14);
558 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNw14);
559 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW34);
560 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw34);
561
562 //compute mobilities of face 1
563 Dune::FieldVector < Scalar, numPhases > lambda12Upw(0.0);
564 lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
565 lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
566
567 //compute mobilities of face 4
568 Dune::FieldVector < Scalar, numPhases > lambda14Upw(0.0);
569 lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
570 lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
571
572 //compute mobilities of face 2
573 Dune::FieldVector < Scalar, numPhases > lambda32Upw(0.0);
574 lambda32Upw[wPhaseIdx] = (potentialDiffW32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
575 lambda32Upw[nPhaseIdx] = (potentialDiffNw32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
576
577 //compute mobilities of face 3
578 Dune::FieldVector < Scalar, numPhases > lambda34Upw(0.0);
579 lambda34Upw[wPhaseIdx] = (potentialDiffW34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
580 lambda34Upw[nPhaseIdx] = (potentialDiffNw34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
581
582 for (int i = 0; i < numPhases; i++)
583 {
584 // evaluate parts of velocity
585 DimVector vel12 = interactionVolume.getNormal(0, 0);
586 DimVector vel14 = interactionVolume.getNormal(3, 0);
587 DimVector vel23 = interactionVolume.getNormal(1, 0);
588 DimVector vel21 = interactionVolume.getNormal(0, 0);
589 DimVector vel34 = interactionVolume.getNormal(2, 0);
590 DimVector vel32 = interactionVolume.getNormal(1, 0);
591 DimVector vel41 = interactionVolume.getNormal(3, 0);
592 DimVector vel43 = interactionVolume.getNormal(2, 0);
593
594 Dune::FieldVector < Scalar, 2 * dim > flux(0);
595 switch (i)
596 {
597 case wPhaseIdx:
598 {
599 flux = fluxW;
600 break;
601 }
602 case nPhaseIdx:
603 {
604 flux = fluxNw;
605 break;
606 }
607 }
608
609 vel12 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0)); //divide by 2 because the flux is related to the half face!
610 vel14 *= flux[3] / (2 * interactionVolume.getFaceArea(0, 1));
611 vel23 *= flux[1] / (2 * interactionVolume.getFaceArea(1, 0));
612 vel21 *= flux[0] / (2 * interactionVolume.getFaceArea(1, 1));
613 vel34 *= flux[2] / (2 * interactionVolume.getFaceArea(2, 0));
614 vel32 *= flux[1] / (2 * interactionVolume.getFaceArea(2, 1));
615 vel41 *= flux[3] / (2 * interactionVolume.getFaceArea(3, 0));
616 vel43 *= flux[2] / (2 * interactionVolume.getFaceArea(3, 1));
617
618 if (level1 < level2)
619 {
620 vel12 *= 0.5;
621 }
622 else if (level2 < level1)
623 {
624 vel21 *= 0.5;
625 }
626 if (level2 < level3)
627 {
628 vel23 *= 0.5;
629 }
630 else if (level3 < level2)
631 {
632 vel32 *= 0.5;
633 }
634 if (level3 < level4)
635 {
636 vel34 *= 0.5;
637 }
638 else if (level4 < level3)
639 {
640 vel43 *= 0.5;
641 }
642 if (level4 < level1)
643 {
644 vel41 *= 0.5;
645 }
646 else if (level1 < level4)
647 {
648 vel14 *= 0.5;
649 }
650
651 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
652 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
653 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
654 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
655 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
656 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
657 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
658 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
659
660 vel12 *= fracFlow12;
661 vel14 *= fracFlow14;
662 vel23 *= fracFlow32;
663 vel21 *= fracFlow12;
664 vel34 *= fracFlow34;
665 vel32 *= fracFlow32;
666 vel41 *= fracFlow14;
667 vel43 *= fracFlow34;
668
669 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
670 {
671 vel12 *= 2;
672 }
673 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
674 {
675 vel14 *= 2;
676 }
677 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
678 {
679 vel23 *= 2;
680 }
681 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
682 {
683 vel21 *= 2;
684 }
685 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
686 {
687 vel34 *= 2;
688 }
689 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
690 {
691 vel32 *= 2;
692 }
693 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
694 {
695 vel41 *= 2;
696 }
697 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
698 {
699 vel43 *= 2;
700 }
701
702 //store velocities
703 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
704 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel14);
705 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel23);
706 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
707 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 0), vel34);
708 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 1), vel32);
709 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel41);
710 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel43);
711 }
712 //set velocity marker
713 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
714 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
715 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
716 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
717 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 0));
718 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 1));
719 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
720 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
721}
722
732template<class TypeTag>
733void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
734 CellData& cellData, int elemIdx)
735{
736 auto element = interactionVolume.getSubVolumeElement(elemIdx);
737
738 // get global coordinate of cell centers
739 const GlobalPosition& globalPos = element.geometry().center();
740
741 //permeability vector at boundary
742 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
743
744 //get mobilities of the phases
745 Dune::FieldVector < Scalar, numPhases > lambda(cellData.mobility(wPhaseIdx));
746 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
747
748 for (int fIdx = 0; fIdx < dim; fIdx++)
749 {
750 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
751
752 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
753 {
754 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
755 {
756 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
757
758 const auto referenceElement = ReferenceElements::general(element.type());
759
760 const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
761
762 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
763
764 DimVector distVec(globalPosFace - globalPos);
765 Scalar dist = distVec.two_norm();
766 DimVector unitDistVec(distVec);
767 unitDistVec /= dist;
768
769 // get pc and lambda at the boundary
770 Scalar satWBound = cellData.saturation(wPhaseIdx);
771 //check boundary sat at face 1
772 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
773 {
774 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
775 switch (saturationType_)
776 {
777 case sw:
778 {
779 satWBound = satBound;
780 break;
781 }
782 case sn:
783 {
784 satWBound = 1 - satBound;
785 break;
786 }
787 }
788
789 }
790
791 Scalar pcBound = MaterialLaw::pc(
792 problem_.spatialParams().materialLawParams(element), satWBound);
793
794 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
795 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
796
797 pcBound += gravityDiffBound;
798
799 Dune::FieldVector < Scalar, numPhases
800 > lambdaBound(
801 MaterialLaw::krw(
802 problem_.spatialParams().materialLawParams(element),
803 satWBound));
804 lambdaBound[nPhaseIdx] = MaterialLaw::krn(
805 problem_.spatialParams().materialLawParams(element), satWBound);
806 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
807 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
808
809 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
810 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
811 Scalar potentialBoundNw = potentialBoundW;
812
813 //calculate potential gradients
814 switch (pressureType_)
815 {
816 case pw:
817 {
818 potentialBoundNw += pcBound;
819 break;
820 }
821 case pn:
822 {
823 //calculate potential gradients
824 potentialBoundW -= pcBound;
825 break;
826 }
827 }
828
829 Scalar potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
830 Scalar potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
831
832 //store potentials for further calculations (saturation, ...)
833 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
834 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
835
836 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
837 DimVector velocityW(0);
838 DimVector velocityNw(0);
839
840 // calculate capillary pressure gradient
841 DimVector pressGradient = unitDistVec;
842 pressGradient *= (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
843 permeability.mv(pressGradient, velocityW);
844
845 pressGradient = unitDistVec;
846 pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
847 permeability.mv(pressGradient, velocityNw);
848
849 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
850 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
851
852 //velocity is calculated from two vertices of one intersection!
853 velocityW *= 0.5;
854 velocityNw *= 0.5;
855
856 //store velocities
857 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
858 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
859 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
860 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
861 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
862 }
863 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
864 {
865 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
866
867 const auto referenceElement = ReferenceElements::general(element.type());
868
869 const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
870
871 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
872
873 DimVector distVec(globalPosFace - globalPos);
874 Scalar dist = distVec.two_norm();
875 DimVector unitDistVec(distVec);
876 unitDistVec /= dist;
877
878 // get neumann boundary value
879 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
880
881 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
882 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
883
884 DimVector velocityW(unitDistVec);
885 DimVector velocityNw(unitDistVec);
886
887 velocityW *= boundValues[wPhaseIdx] / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
888 velocityNw *= boundValues[nPhaseIdx]
889 / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
890
891 //store potentials for further calculations (saturation, ...)
892 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
893 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
894
895 //store velocities
896 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
897 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
898 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
899 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
900 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
901 }
902 else
903 {
904 DUNE_THROW(Dune::NotImplemented,
905 "No valid boundary condition type defined for pressure equation!");
906 }
907 }
908 }
909}
910
911} // end namespace Dumux
912#endif
GET_PROP_VALUE
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
GET_PROP
#define GET_PROP(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:281
GET_PROP_TYPE
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
2dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation 2-d.
linteractionvolume.hh
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Dumux::intersections
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
Dumux
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Dumux::Properties::SaturationFormulation
Property tag SaturationFormulation
The formulation of the saturation model.
Definition: porousmediumflow/2p/sequential/properties.hh:53
Dumux::Properties::NumPhases
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
Dumux::Properties::VelocityFormulation
Property tag VelocityFormulation
The type of velocity reconstructed for the transport model.
Definition: porousmediumflow/2p/sequential/properties.hh:54
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::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::BoundaryTypes::isNeumann
bool isNeumann(unsigned eqIdx) const
Returns true if an equation is used to specify a Neumann condition.
Definition: common/boundarytypes.hh:269
Dumux::BoundaryTypes::isDirichlet
bool isDirichlet(unsigned eqIdx) const
Returns true if an equation is used to specify a Dirichlet condition.
Definition: common/boundarytypes.hh:236
Dumux::FvMpfaL2dTransmissibilityCalculator
Provides methods for transmissibility calculation in 2-d.
Definition: 2dtransmissibilitycalculator.hh:44
Dumux::FvMpfaL2dVelocity2p
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: lmethod/2dvelocity.hh:58
Dumux::FvMpfaL2dVelocity2p::viscosity_
Scalar viscosity_[numPhases]
Definition: lmethod/2dvelocity.hh:251
Dumux::FvMpfaL2dVelocity2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: lmethod/2dvelocity.hh:177
Dumux::FvMpfaL2dVelocity2p::gravity_
const GravityVector & gravity_
Definition: lmethod/2dvelocity.hh:248
Dumux::FvMpfaL2dVelocity2p::vtkOutputLevel_
int vtkOutputLevel_
Definition: lmethod/2dvelocity.hh:253
Dumux::FvMpfaL2dVelocity2p::calculateBoundaryInteractionVolumeVelocity
void calculateBoundaryInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData, int elemIdx)
Calculates the velocity at a boundary flux faces.
Definition: lmethod/2dvelocity.hh:733
Dumux::FvMpfaL2dVelocity2p::initialize
void initialize()
Initializes the velocity model.
Definition: lmethod/2dvelocity.hh:149
Dumux::FvMpfaL2dVelocity2p::velocityType_
static const int velocityType_
gives kind of velocity used ( , , )
Definition: lmethod/2dvelocity.hh:257
Dumux::FvMpfaL2dVelocity2p::saturationType_
static const int saturationType_
gives kind of saturation used ( , )
Definition: lmethod/2dvelocity.hh:261
Dumux::FvMpfaL2dVelocity2p::transmissibilityCalculator_
TransmissibilityCalculator transmissibilityCalculator_
Definition: lmethod/2dvelocity.hh:246
Dumux::FvMpfaL2dVelocity2p::calculateInnerInteractionVolumeVelocity
void calculateInnerInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData3, CellData &cellData4, InnerBoundaryVolumeFaces &innerBoundaryVolumeFaces)
Calculate velocities for flux faces of an interaction volume.
Definition: lmethod/2dvelocity.hh:279
Dumux::FvMpfaL2dVelocity2p::pressureType_
static const int pressureType_
gives kind of pressure used ( , , )
Definition: lmethod/2dvelocity.hh:259
Dumux::FvMpfaL2dVelocity2p::density_
Scalar density_[numPhases]
Definition: lmethod/2dvelocity.hh:250
Dumux::FvMpfaL2dVelocity2p::threshold_
static constexpr Scalar threshold_
Definition: lmethod/2dvelocity.hh:255
Dumux::FvMpfaL2dVelocity2p::FvMpfaL2dVelocity2p
FvMpfaL2dVelocity2p(Problem &problem)
Constructs a FvMpfaL2dVelocity2p object.
Definition: lmethod/2dvelocity.hh:130
Dumux::FVMPFALInteractionVolume
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
Dumux::FVMPFALInteractionVolume::getIndexOnElement
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering to numbering of the Dune reference element.
Definition: linteractionvolume.hh:258
Dumux::FVMPFALInteractionVolume::getNormal
DimVector & getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face normal.
Definition: linteractionvolume.hh:364
Dumux::FVMPFALInteractionVolume::getSubVolumeElement
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: linteractionvolume.hh:281
Dumux::FVMPFALInteractionVolume::getBoundaryType
BoundaryTypes & getBoundaryType(int subVolumeFaceIdx)
Get boundary condtion types for a flux face.
Definition: linteractionvolume.hh:292
Dumux::FVMPFALInteractionVolume::getFaceArea
Scalar & getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face area.
Definition: linteractionvolume.hh:388
Dumux::FVMPFALInteractionVolume::getNeumannValues
PrimaryVariables & getNeumannValues(int subVolumeFaceIdx)
Get the Neumann boundary condtions for a flux face.
Definition: linteractionvolume.hh:352
Dumux::FVMPFALInteractionVolume::isBoundaryFace
bool isBoundaryFace(int subVolumeFaceIdx)
Returns true if an interaction volume flux face is a boundary face.
Definition: linteractionvolume.hh:330
Dumux::FVMPFALInteractionVolume::getFaceIndexFromSubVolume
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering on element to numbering on interaction volume.
Definition: linteractionvolume.hh:270
Dumux::FVMPFALInteractionVolume::getDirichletValues
PrimaryVariables & getDirichletValues(int subVolumeFaceIdx)
Get the Dirichlet boundary condtions for a flux face.
Definition: linteractionvolume.hh:341
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
Generated by doxygen 1.9.3