3.2-git
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#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
31
32namespace Dumux {
55template<class TypeTag> class FvMpfaO2dVelocity2P
56{
57 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
58 enum
59 {
60 dim = GridView::dimension, dimWorld = GridView::dimensionworld
61 };
62
63 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
64 using Problem = GetPropType<TypeTag, Properties::Problem>;
65
66 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
67
68 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
69 using MaterialLaw = typename SpatialParams::MaterialLaw;
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::BoundaryTypes>;
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 = ReferenceElements::general(element.type()).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, "non-wetting-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 referenceElement = ReferenceElements::general(element.type());
591
592 const LocalPosition& localPos = referenceElement.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 Scalar pcBound = MaterialLaw::pc(
624 problem_.spatialParams().materialLawParams(element), satWBound);
625
626 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
627 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
628
629 pcBound += gravityDiffBound;
630
631 Dune::FieldVector < Scalar, numPhases
632 > lambdaBound(
633 MaterialLaw::krw(
634 problem_.spatialParams().materialLawParams(element),
635 satWBound));
636 lambdaBound[nPhaseIdx] = MaterialLaw::krn(
637 problem_.spatialParams().materialLawParams(element), satWBound);
638 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
639 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
640
641 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
642 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
643 Scalar potentialBoundNw = potentialBoundW;
644
645 //calculate potential gradients
646 switch (pressureType_)
647 {
648 case pw:
649 {
650 potentialBoundNw += pcBound;
651 break;
652 }
653 case pn:
654 {
655 //calculate potential gradients
656 potentialBoundW -= pcBound;
657 break;
658 }
659 }
660
661 Scalar potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
662 Scalar potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
663
664 //store potentials for further calculations (saturation, ...)
665 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
666 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
667
668 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
669 DimVector velocityW(0);
670 DimVector velocityNw(0);
671
672 // calculate capillary pressure gradient
673 DimVector pressGradient = unitDistVec;
674 pressGradient *= (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
675 permeability.mv(pressGradient, velocityW);
676
677 pressGradient = unitDistVec;
678 pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
679 permeability.mv(pressGradient, velocityNw);
680
681 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
682 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
683
684 //velocity is calculated from two vertices of one intersection!
685 velocityW *= 0.5;
686 velocityNw *= 0.5;
687
688 //store velocities
689 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
690 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
691 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
692 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
693 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
694 }
695 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
696 {
697 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
698
699 const auto referenceElement = ReferenceElements::general(element.type());
700
701 const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
702
703 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
704
705 DimVector distVec(globalPosFace - globalPos);
706 Scalar dist = distVec.two_norm();
707 DimVector unitDistVec(distVec);
708 unitDistVec /= dist;
709
710 // get neumann boundary value
711 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
712
713 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
714 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
715
716 DimVector velocityW(unitDistVec);
717 DimVector velocityNw(unitDistVec);
718
719 velocityW *= boundValues[wPhaseIdx] / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
720 velocityNw *= boundValues[nPhaseIdx]
721 / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
722
723 //store potentials for further calculations (saturation, ...)
724 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
725 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
726
727 //store velocities
728 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
729 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
730 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
731 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
732 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
733 }
734 else
735 {
736 DUNE_THROW(Dune::NotImplemented,
737 "No valid boundary condition type defined for pressure equation!");
738 }
739 }
740 }
741}
742
743} // end namespace Dumux
744#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 (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:56
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