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