3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
omethod/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_FVMPFAO2DVELOCITY2P_HH
25#define DUMUX_FVMPFAO2DVELOCITY2P_HH
26
27#include <dune/grid/common/gridenums.hh>
28
29#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
31#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
32
33#include <dumux/common/deprecated.hh>
34
35namespace Dumux {
58template<class TypeTag> class FvMpfaO2dVelocity2P
59{
60 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
61 enum
62 {
63 dim = GridView::dimension, dimWorld = GridView::dimensionworld
64 };
65
66 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
67 using Problem = GetPropType<TypeTag, Properties::Problem>;
68
69 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
70
71 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
72
73 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
74 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
75
76 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
77 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
78 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
79 using CellData = GetPropType<TypeTag, Properties::CellData>;
80
81 using Element = typename GridView::Traits::template Codim<0>::Entity;
82 using Grid = typename GridView::Grid;
83 using IndexSet = typename GridView::IndexSet;
84
85 using Geometry = typename Element::Geometry;
86 using JacobianTransposed = typename Geometry::JacobianTransposed ;
87
88 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
89
90 using InteractionVolume = FVMPFAOInteractionVolume<TypeTag>;
91 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
92
93 enum
94 {
95 pw = Indices::pressureW,
96 pn = Indices::pressureNw,
97 pGlobal = Indices::pressureGlobal,
98 sw = Indices::saturationW,
99 sn = Indices::saturationNw,
100 vw = Indices::velocityW,
101 vn = Indices::velocityNw,
102 vt = Indices::velocityTotal
103 };
104 enum
105 {
106 wPhaseIdx = Indices::wPhaseIdx,
107 nPhaseIdx = Indices::nPhaseIdx,
108 pressureIdx = Indices::pressureIdx,
109 saturationIdx = Indices::saturationIdx,
110 pressureEqIdx = Indices::pressureEqIdx,
111 satEqIdx = Indices::satEqIdx,
112 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
113 };
114
115 using LocalPosition = Dune::FieldVector<Scalar, dim>;
116 using GlobalPosition = typename Geometry::GlobalCoordinate;
117 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
118 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
119 using DimVector = Dune::FieldVector<Scalar, dim>;
120
121public:
122 /*
123 * \brief Constructs a FvMpfaO2dVelocity2P object
124 *
125 * \param problem A problem class object
126 */
127 FvMpfaO2dVelocity2P(Problem& problem) :
128 problem_(problem), gravity_(problem.gravity())
129 {
130 density_[wPhaseIdx] = 0.;
131 density_[nPhaseIdx] = 0.;
132 viscosity_[wPhaseIdx] = 0.;
133 viscosity_[nPhaseIdx] = 0.;
134
135 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
136 }
137
139 void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1,
140 CellData& cellData2, CellData& cellData3, CellData& cellData4,
141 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
142 void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
143 CellData& cellData, int elemIdx);
144
146 void initialize()
147 {
148 const auto element = *problem_.gridView().template begin<0>();
149 FluidState fluidState;
150 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
151 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
152 fluidState.setTemperature(problem_.temperature(element));
153 fluidState.setSaturation(wPhaseIdx, 1.);
154 fluidState.setSaturation(nPhaseIdx, 0.);
155 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
156 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
157 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
158 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
159
160 return;
161 }
162
173 template<class MultiWriter>
174 void addOutputVtkFields(MultiWriter &writer)
175 {
176 if (vtkOutputLevel_ > 0)
177 {
178 Dune::BlockVector < DimVector > &velocityWetting = *(writer.template allocateManagedBuffer<Scalar,
179 dim>(problem_.gridView().size(0)));
180 Dune::BlockVector < DimVector > &velocityNonwetting = *(writer.template allocateManagedBuffer<Scalar,
181 dim>(problem_.gridView().size(0)));
182
183 // compute update vector
184 for (const auto& element : elements(problem_.gridView()))
185 {
186 // cell index
187 int eIdxGlobal = problem_.variables().index(element);
188
189 CellData & cellData = problem_.variables().cellData(eIdxGlobal);
190
191 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
192 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
193
194 // run through all intersections with neighbors and boundary
195 for (const auto& intersection : intersections(problem_.gridView(), element))
196 {
197 int isIndex = intersection.indexInInside();
198
199 fluxW[isIndex] += intersection.geometry().volume()
200 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
201 fluxNw[isIndex] += intersection.geometry().volume()
202 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
203 }
204
205 DimVector refVelocity(0);
206 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
207 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
208
209 const DimVector& localPos = referenceElement(element).position(0, 0);
210
211 // get the transposed Jacobian of the element mapping
212 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
213
214 // calculate the element velocity by the Piola transformation
215 DimVector elementVelocity(0);
216 jacobianT.umtv(refVelocity, elementVelocity);
217 elementVelocity /= element.geometry().integrationElement(localPos);
218
219 velocityWetting[eIdxGlobal] = elementVelocity;
220
221 refVelocity = 0;
222 refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
223 refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
224
225 // calculate the element velocity by the Piola transformation
226 elementVelocity = 0;
227 jacobianT.umtv(refVelocity, elementVelocity);
228 elementVelocity /= element.geometry().integrationElement(localPos);
229
230 velocityNonwetting[eIdxGlobal] = elementVelocity;
231 }
232
233 writer.attachCellData(velocityWetting, "wetting-velocity", dim);
234 writer.attachCellData(velocityNonwetting, "nonwetting-velocity", dim);
235 }
236
237 return;
238 }
239
240private:
241 Problem& problem_;
242 const GravityVector& gravity_;
243
244 Scalar density_[numPhases];
245 Scalar viscosity_[numPhases];
246
247 int vtkOutputLevel_;
248
249 static constexpr Scalar threshold_ = 1e-15;
251 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
253 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
255 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
256};
257// end of template
258
274template<class TypeTag>
275void FvMpfaO2dVelocity2P<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
276 CellData& cellData1, CellData& cellData2,
277 CellData& cellData3, CellData& cellData4,
278 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
279{
280 // cell index
281 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
282 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
283 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
284 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
285
286 // get pressure values
287 Dune::FieldVector < Scalar, 2 * dim > potW(0);
288 Dune::FieldVector < Scalar, 2 * dim > potNw(0);
289
290 potW[0] = cellData1.potential(wPhaseIdx);
291 potW[1] = cellData2.potential(wPhaseIdx);
292 potW[2] = cellData3.potential(wPhaseIdx);
293 potW[3] = cellData4.potential(wPhaseIdx);
294
295 potNw[0] = cellData1.potential(nPhaseIdx);
296 potNw[1] = cellData2.potential(nPhaseIdx);
297 potNw[2] = cellData3.potential(nPhaseIdx);
298 potNw[3] = cellData4.potential(nPhaseIdx);
299
300 //get mobilities of the phases
301 Dune::FieldVector < Scalar, numPhases > lambda1(cellData1.mobility(wPhaseIdx));
302 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
303
304 //compute total mobility of cell 1
305 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
306
307 //get mobilities of the phases
308 Dune::FieldVector < Scalar, numPhases > lambda2(cellData2.mobility(wPhaseIdx));
309 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
310
311 //compute total mobility of cell 1
312 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
313
314 //get mobilities of the phases
315 Dune::FieldVector < Scalar, numPhases > lambda3(cellData3.mobility(wPhaseIdx));
316 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
317
318 //compute total mobility of cell 1
319 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
320
321 //get mobilities of the phases
322 Dune::FieldVector < Scalar, numPhases > lambda4(cellData4.mobility(wPhaseIdx));
323 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
324
325 //compute total mobility of cell 1
326 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
327
328 Scalar gn12nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 1);
329 Scalar gn12nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 0);
330 Scalar gn14nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 1);
331 Scalar gn14nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 0);
332 Scalar gn12nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 0);
333 Scalar gn12nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 1);
334 Scalar gn23nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 0);
335 Scalar gn23nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 1);
336 Scalar gn43nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 1);
337 Scalar gn43nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 0);
338 Scalar gn23nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 1);
339 Scalar gn23nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 0);
340 Scalar gn43nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 0);
341 Scalar gn43nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 1);
342 Scalar gn14nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 0);
343 Scalar gn14nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 1);
344
345 // compute transmissibility matrix T = CA^{-1}B+F
346 Dune::FieldMatrix < Scalar, 2 * dim, 2 * dim > C(0), F(0), A(0), B(0);
347
348 // evaluate matrix C, F, A, B
349 C[0][0] = -gn12nu12;
350 C[0][3] = -gn12nu14;
351 C[1][0] = gn23nu21;
352 C[1][1] = -gn23nu23;
353 C[2][1] = gn43nu32;
354 C[2][2] = gn43nu34;
355 C[3][2] = -gn14nu43;
356 C[3][3] = gn14nu41;
357
358 F[0][0] = gn12nu12 + gn12nu14;
359 F[1][1] = -gn23nu21 + gn23nu23;
360 F[2][2] = -gn43nu34 - gn43nu32;
361 F[3][3] = gn14nu43 - gn14nu41;
362
363 A[0][0] = gn12nu12 + gn12nu21;
364 A[0][1] = -gn12nu23;
365 A[0][3] = gn12nu14;
366 A[1][0] = -gn23nu21;
367 A[1][1] = gn23nu23 + gn23nu32;
368 A[1][2] = gn23nu34;
369 A[2][1] = -gn43nu32;
370 A[2][2] = -gn43nu34 - gn43nu43;
371 A[2][3] = gn43nu41;
372 A[3][0] = -gn14nu12;
373 A[3][2] = gn14nu43;
374 A[3][3] = -gn14nu41 - gn14nu14;
375
376 // std::cout << A << "\n";
377
378 B[0][0] = gn12nu12 + gn12nu14;
379 B[0][1] = gn12nu21 - gn12nu23;
380 B[1][1] = -gn23nu21 + gn23nu23;
381 B[1][2] = gn23nu34 + gn23nu32;
382 B[2][2] = -gn43nu34 - gn43nu32;
383 B[2][3] = -gn43nu43 + gn43nu41;
384 B[3][0] = -gn14nu12 - gn14nu14;
385 B[3][3] = gn14nu43 - gn14nu41;
386
387 //flux vector
388 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
389 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
390
391 // compute T
392 A.invert();
393 F += C.rightmultiply(B.leftmultiply(A));
394 Dune::FieldMatrix < Scalar, 2 * dim, 2 * dim > T(F);
395
396 T.mv(potW, fluxW);
397 T.mv(potNw, fluxNw);
398
399 Scalar potentialDiffW12 = fluxW[0];
400 Scalar potentialDiffW14 = fluxW[3];
401 Scalar potentialDiffW32 = -fluxW[1];
402 Scalar potentialDiffW34 = -fluxW[2];
403
404 Scalar potentialDiffNw12 = fluxNw[0];
405 Scalar potentialDiffNw14 = fluxNw[3];
406 Scalar potentialDiffNw32 = -fluxNw[1];
407 Scalar potentialDiffNw34 = -fluxNw[2];
408
409 //store potentials for further calculations (saturation, ...)
410 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), fluxW[0]);
411 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), fluxNw[0]);
412 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), fluxW[3]);
413 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), fluxNw[3]);
414 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), fluxW[1]);
415 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), fluxNw[1]);
416 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -fluxW[0]);
417 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -fluxNw[0]);
418 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), -fluxW[2]);
419 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), -fluxNw[2]);
420 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -fluxW[1]);
421 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -fluxNw[1]);
422 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -fluxW[3]);
423 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -fluxNw[3]);
424 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), fluxW[2]);
425 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), fluxNw[2]);
426
427 //compute mobilities of face 1
428 Dune::FieldVector < Scalar, numPhases > lambda12Upw(0.0);
429 lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
430 lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
431
432 //compute mobilities of face 4
433 Dune::FieldVector < Scalar, numPhases > lambda14Upw(0.0);
434 lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
435 lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
436
437 //compute mobilities of face 2
438 Dune::FieldVector < Scalar, numPhases > lambda32Upw(0.0);
439 lambda32Upw[wPhaseIdx] = (potentialDiffW32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
440 lambda32Upw[nPhaseIdx] = (potentialDiffNw32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
441
442 //compute mobilities of face 3
443 Dune::FieldVector < Scalar, numPhases > lambda34Upw(0.0);
444 lambda34Upw[wPhaseIdx] = (potentialDiffW34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
445 lambda34Upw[nPhaseIdx] = (potentialDiffNw34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
446
447 for (int i = 0; i < numPhases; i++)
448 {
449 // evaluate parts of velocity
450 DimVector vel12 = interactionVolume.getNormal(0, 0);
451 DimVector vel14 = interactionVolume.getNormal(0, 1);
452 DimVector vel23 = interactionVolume.getNormal(1, 0);
453 DimVector vel21 = interactionVolume.getNormal(1, 1);
454 DimVector vel34 = interactionVolume.getNormal(2, 0);
455 DimVector vel32 = interactionVolume.getNormal(2, 1);
456 DimVector vel41 = interactionVolume.getNormal(3, 0);
457 DimVector vel43 = interactionVolume.getNormal(3, 1);
458
459 Dune::FieldVector < Scalar, 2 * dim > flux(0);
460 switch (i)
461 {
462 case wPhaseIdx:
463 {
464 flux = fluxW;
465 break;
466 }
467 case nPhaseIdx:
468 {
469 flux = fluxNw;
470 break;
471 }
472 }
473
474 vel12 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0)); //divide by 2 because the flux is related to the half face!
475 vel14 *= flux[3] / (2 * interactionVolume.getFaceArea(0, 1));
476 vel23 *= flux[1] / (2 * interactionVolume.getFaceArea(1, 0));
477 vel21 *= flux[0] / (2 * interactionVolume.getFaceArea(1, 1));
478 vel34 *= flux[2] / (2 * interactionVolume.getFaceArea(2, 0));
479 vel32 *= flux[1] / (2 * interactionVolume.getFaceArea(2, 1));
480 vel41 *= flux[3] / (2 * interactionVolume.getFaceArea(3, 0));
481 vel43 *= flux[2] / (2 * interactionVolume.getFaceArea(3, 1));
482
483 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
484 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
485 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
486 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
487 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
488 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
489 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
490 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
491
492 vel12 *= fracFlow12;
493 vel14 *= fracFlow14;
494 vel23 *= fracFlow32;
495 vel21 *= fracFlow12;
496 vel34 *= fracFlow34;
497 vel32 *= fracFlow32;
498 vel41 *= fracFlow14;
499 vel43 *= fracFlow34;
500
501 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
502 {
503 vel12 *= 2;
504 }
505 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
506 {
507 vel14 *= 2;
508 }
509 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
510 {
511 vel23 *= 2;
512 }
513 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
514 {
515 vel21 *= 2;
516 }
517 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
518 {
519 vel34 *= 2;
520 }
521 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
522 {
523 vel32 *= 2;
524 }
525 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
526 {
527 vel41 *= 2;
528 }
529 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
530 {
531 vel43 *= 2;
532 }
533
534 //store velocities
535 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
536 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel14);
537 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel23);
538 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
539 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 0), vel34);
540 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 1), vel32);
541 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel41);
542 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel43);
543 }
544 //set velocity marker
545 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
546 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
547 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
548 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
549 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 0));
550 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 1));
551 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
552 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
553}
554
564template<class TypeTag>
565void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
566 CellData& cellData, int elemIdx)
567{
568 auto element = interactionVolume.getSubVolumeElement(elemIdx);
569
570 // get global coordinate of cell centers
571 const GlobalPosition& globalPos = element.geometry().center();
572
573 //permeability vector at boundary
574 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
575
576 //get mobilities of the phases
577 Dune::FieldVector < Scalar, numPhases > lambda(cellData.mobility(wPhaseIdx));
578 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
579
580 for (int fIdx = 0; fIdx < dim; fIdx++)
581 {
582 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
583
584 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
585 {
586 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
587 {
588 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
589
590 const auto refElement = referenceElement(element);
591
592 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
593
594 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
595
596 DimVector distVec(globalPosFace - globalPos);
597 Scalar dist = distVec.two_norm();
598 DimVector unitDistVec(distVec);
599 unitDistVec /= dist;
600
601 // get pc and lambda at the boundary
602 Scalar satWBound = cellData.saturation(wPhaseIdx);
603 //check boundary sat at face 1
604 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
605 {
606 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
607 switch (saturationType_)
608 {
609 case sw:
610 {
611 satWBound = satBound;
612 break;
613 }
614 case sn:
615 {
616 satWBound = 1 - satBound;
617 break;
618 }
619 }
620
621 }
622
623 // old material law interface is deprecated: Replace this by
624 // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
625 // after the release of 3.3, when the deprecated interface is no longer supported
626 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
627
628 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
629
630 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
631 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
632
633 pcBound += gravityDiffBound;
634
635 Dune::FieldVector < Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
636 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
637 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
638 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
639
640 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
641 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
642 Scalar potentialBoundNw = potentialBoundW;
643
644 //calculate potential gradients
645 switch (pressureType_)
646 {
647 case pw:
648 {
649 potentialBoundNw += pcBound;
650 break;
651 }
652 case pn:
653 {
654 //calculate potential gradients
655 potentialBoundW -= pcBound;
656 break;
657 }
658 }
659
660 Scalar potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
661 Scalar potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
662
663 //store potentials for further calculations (saturation, ...)
664 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
665 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
666
667 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
668 DimVector velocityW(0);
669 DimVector velocityNw(0);
670
671 // calculate capillary pressure gradient
672 DimVector pressGradient = unitDistVec;
673 pressGradient *= (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
674 permeability.mv(pressGradient, velocityW);
675
676 pressGradient = unitDistVec;
677 pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
678 permeability.mv(pressGradient, velocityNw);
679
680 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
681 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
682
683 //velocity is calculated from two vertices of one intersection!
684 velocityW *= 0.5;
685 velocityNw *= 0.5;
686
687 //store velocities
688 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
689 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
690 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
691 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
692 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
693 }
694 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
695 {
696 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
697
698 const auto refElement = referenceElement(element);
699
700 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
701
702 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
703
704 DimVector distVec(globalPosFace - globalPos);
705 Scalar dist = distVec.two_norm();
706 DimVector unitDistVec(distVec);
707 unitDistVec /= dist;
708
709 // get neumann boundary value
710 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
711
712 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
713 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
714
715 DimVector velocityW(unitDistVec);
716 DimVector velocityNw(unitDistVec);
717
718 velocityW *= boundValues[wPhaseIdx] / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
719 velocityNw *= boundValues[nPhaseIdx]
720 / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
721
722 //store potentials for further calculations (saturation, ...)
723 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
724 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
725
726 //store velocities
727 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
728 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
729 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
730 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
731 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
732 }
733 else
734 {
735 DUNE_THROW(Dune::NotImplemented,
736 "No valid boundary condition type defined for pressure equation!");
737 }
738 }
739 }
740}
741
742} // end namespace Dumux
743#endif
deprecated.hh
Helpers for deprecation.
ointeractionvolume.hh
Class including the information of an interaction volume of a MPFA O-method that does not change with...
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::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:273
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::FvMpfaO2dVelocity2P
Class for calculating velocities from cell-wise constant pressure values.
Definition: omethod/2dvelocity.hh:59
Dumux::FvMpfaO2dVelocity2P::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: omethod/2dvelocity.hh:174
Dumux::FvMpfaO2dVelocity2P::FvMpfaO2dVelocity2P
FvMpfaO2dVelocity2P(Problem &problem)
Definition: omethod/2dvelocity.hh:127
Dumux::FvMpfaO2dVelocity2P::calculateInnerInteractionVolumeVelocity
void calculateInnerInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData3, CellData &cellData4, InnerBoundaryVolumeFaces &innerBoundaryVolumeFaces)
Calculates velocities for all flux faces of an interaction volume.
Definition: omethod/2dvelocity.hh:275
Dumux::FvMpfaO2dVelocity2P::initialize
void initialize()
Initializes the velocity model.
Definition: omethod/2dvelocity.hh:146
Dumux::FvMpfaO2dVelocity2P::calculateBoundaryInteractionVolumeVelocity
void calculateBoundaryInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData, int elemIdx)
Calculates the velocity at a boundary flux faces.
Definition: omethod/2dvelocity.hh:565
Dumux::FVMPFAOInteractionVolume
Class including the information of an interaction volume of a MPFA O-method that does not change with...
Definition: ointeractionvolume.hh:38
Dumux::FVMPFAOInteractionVolume::getBoundaryType
BoundaryTypes & getBoundaryType(int subVolumeFaceIdx)
Get boundary condtion types for a flux face.
Definition: ointeractionvolume.hh:267
Dumux::FVMPFAOInteractionVolume::isBoundaryFace
bool isBoundaryFace(int subVolumeFaceIdx)
Returns true if an interaction volume flux face is a boundary face.
Definition: ointeractionvolume.hh:305
Dumux::FVMPFAOInteractionVolume::getDirichletValues
PrimaryVariables & getDirichletValues(int subVolumeFaceIdx)
Get the Dirichlet boundary condtions for a flux face.
Definition: ointeractionvolume.hh:316
Dumux::FVMPFAOInteractionVolume::getFaceArea
Scalar & getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face area.
Definition: ointeractionvolume.hh:351
Dumux::FVMPFAOInteractionVolume::getNtkrkNu_df
Scalar getNtkrkNu_df(Scalar &relPerm, int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
Definition: ointeractionvolume.hh:434
Dumux::FVMPFAOInteractionVolume::getNeumannValues
PrimaryVariables & getNeumannValues(int subVolumeFaceIdx)
Get the Neumann boundary condtions for a flux face.
Definition: ointeractionvolume.hh:327
Dumux::FVMPFAOInteractionVolume::getNormal
DimVector & getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face normal.
Definition: ointeractionvolume.hh:339
Dumux::FVMPFAOInteractionVolume::getFaceIndexFromSubVolume
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering on element to numbering on interaction volume.
Definition: ointeractionvolume.hh:245
Dumux::FVMPFAOInteractionVolume::getSubVolumeElement
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: ointeractionvolume.hh:256
Dumux::FVMPFAOInteractionVolume::getIndexOnElement
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering to numbering of the Dune reference element.
Definition: ointeractionvolume.hh:233
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
Generated by doxygen 1.9.3