3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
2dvelocityadaptive.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_FVMPFAL2DVELOCITY2P_ADAPTIVE_HH
25#define DUMUX_FVMPFAL2DVELOCITY2P_ADAPTIVE_HH
26
27#include "2dvelocity.hh"
28
29namespace Dumux {
30
53template<class TypeTag> class FvMpfaL2dVelocity2pAdaptive : public FvMpfaL2dVelocity2p<TypeTag>
54{
55 using ParentType = FvMpfaL2dVelocity2p<TypeTag>;
56 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
57 enum
58 {
59 dim = GridView::dimension, dimWorld = GridView::dimensionworld
60 };
61
62 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
63 using Problem = GetPropType<TypeTag, Properties::Problem>;
64
65 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
66
67 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
68
69 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
70 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
71
72 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
73 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
74 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
75 using CellData = GetPropType<TypeTag, Properties::CellData>;
76
77 using Element = typename GridView::Traits::template Codim<0>::Entity;
78 using Grid = typename GridView::Grid;
79 using IndexSet = typename GridView::IndexSet;
80
81 using Geometry = typename Element::Geometry;
82 using JacobianTransposed = typename Geometry::JacobianTransposed ;
83
84 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
85
86 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
87 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
88 using TransmissibilityCalculator = FvMpfaL2dTransmissibilityCalculator<TypeTag>;
89
90 enum
91 {
92 pw = Indices::pressureW,
93 pn = Indices::pressureNw,
94 pGlobal = Indices::pressureGlobal,
95 sw = Indices::saturationW,
96 sn = Indices::saturationNw,
97 vw = Indices::velocityW,
98 vn = Indices::velocityNw,
99 vt = Indices::velocityTotal
100 };
101 enum
102 {
103 wPhaseIdx = Indices::wPhaseIdx,
104 nPhaseIdx = Indices::nPhaseIdx,
105 pressureIdx = Indices::pressureIdx,
106 saturationIdx = Indices::saturationIdx,
107 pressureEqIdx = Indices::pressureEqIdx,
108 satEqIdx = Indices::satEqIdx,
109 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
110 };
111
112 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
113 using DimVector = Dune::FieldVector<Scalar, dim>;
114
115public:
121 FvMpfaL2dVelocity2pAdaptive(Problem& problem) :
122 ParentType(problem), problem_(problem)
123 {}
124
126 void calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1,
127 CellData& cellData2, CellData& cellData4,
128 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
129
130private:
131 Problem& problem_;
132};
133// end of template
134
147template<class TypeTag>
148void FvMpfaL2dVelocity2pAdaptive<TypeTag>::calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume,
149 CellData& cellData1, CellData& cellData2,
150 CellData& cellData4,
151 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
152{
153 // cell index
154 int globalIdx1 = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
155 int globalIdx2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
156
157 // get pressure values
158 Dune::FieldVector < Scalar, 2 * dim > potW(0);
159 Dune::FieldVector < Scalar, 2 * dim > potNw(0);
160
161 potW[0] = cellData1.potential(wPhaseIdx);
162 potW[1] = cellData2.potential(wPhaseIdx);
163 potW[2] = cellData4.potential(wPhaseIdx);
164
165 potNw[0] = cellData1.potential(nPhaseIdx);
166 potNw[1] = cellData2.potential(nPhaseIdx);
167 potNw[2] = cellData4.potential(nPhaseIdx);
168
169 //get mobilities of the phases
170 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
171 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
172
173 //compute total mobility of cell 1
174 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
175
176 //get mobilities of the phases
177 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
178 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
179
180 //compute total mobility of cell 1
181 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
182
183 //get mobilities of the phases
184 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
185 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
186
187 //compute total mobility of cell 1
188 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
189
190 std::vector < DimVector > lambda(4);
191
192 lambda[0][0] = lambdaTotal1;
193 lambda[0][1] = lambdaTotal1;
194 lambda[1][0] = lambdaTotal2;
195 lambda[1][1] = lambdaTotal2;
196 lambda[3][0] = lambdaTotal4;
197 lambda[3][1] = lambdaTotal4;
198
199 Scalar potentialDiffW12 = 0;
200 Scalar potentialDiffW14 = 0;
201 Scalar potentialDiffW24 = 0;
202
203 Scalar potentialDiffNw12 = 0;
204 Scalar potentialDiffNw14 = 0;
205 Scalar potentialDiffNw24 = 0;
206
207 //flux vector
208 Dune::FieldVector<Scalar, 3> fluxW(0);
209 Dune::FieldVector<Scalar, 3> fluxNw(0);
210
211 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> T(0);
212 DimVector Tu(0);
213 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
214
215 int transmissibilityType = this->transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 3, 3);
216
217 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
218 {
219 u[0] = potW[1];
220 u[1] = potW[2];
221 u[2] = potW[0];
222
223 T.mv(u, Tu);
224
225 fluxW[0] = Tu[1];
226 potentialDiffW12 = Tu[1];
227
228 u[0] = potNw[1];
229 u[1] = potNw[2];
230 u[2] = potNw[0];
231
232 T.mv(u, Tu);
233
234 fluxNw[0] = Tu[1];
235 potentialDiffNw12 = Tu[1];
236 }
237 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
238 {
239 u[0] = potW[0];
240 u[1] = potW[2];
241 u[2] = potW[1];
242
243 T.mv(u, Tu);
244
245 fluxW[0] = Tu[1];
246 potentialDiffW12 = Tu[1];
247
248 u[0] = potNw[0];
249 u[1] = potNw[2];
250 u[2] = potNw[1];
251
252 T.mv(u, Tu);
253
254 fluxNw[0] = Tu[1];
255 potentialDiffNw12 = Tu[1];
256 }
257
258 transmissibilityType = this->transmissibilityCalculator_.calculateLeftHNTransmissibility(T, interactionVolume, lambda, 1, 3, 0);
259
260 if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
261 {
262 u[0] = potW[1];
263 u[1] = potW[0];
264 u[2] = potW[2];
265
266 T.mv(u, Tu);
267
268 fluxW[1] = Tu[1];
269 potentialDiffW24 = Tu[1];
270
271 u[0] = potNw[1];
272 u[1] = potNw[0];
273 u[2] = potNw[2];
274
275 T.mv(u, Tu);
276
277 fluxNw[1] = Tu[1];
278 potentialDiffNw24 = Tu[1];
279 }
280
281 transmissibilityType = this->transmissibilityCalculator_.calculateRightHNTransmissibility(T, interactionVolume, lambda, 3, 0, 1);
282
283 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
284 {
285 u[0] = potW[0];
286 u[1] = potW[1];
287 u[2] = potW[2];
288
289 T.mv(u, Tu);
290
291 fluxW[2] = Tu[1];
292 potentialDiffW14 = -Tu[1];
293
294 u[0] = potNw[0];
295 u[1] = potNw[1];
296 u[2] = potNw[2];
297
298 T.mv(u, Tu);
299
300 fluxNw[2] = Tu[1];
301 potentialDiffNw14 = -Tu[1];
302 }
303
304 //store potentials for further calculations (saturation, ...) -> maybe add new potential to old one!!
305 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW12);
306 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw12);
307 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffW14);
308 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNw14);
309 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffW24);
310 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffNw24);
311 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW12);
312 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw12);
313 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffW14);
314 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNw14);
315 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW24);
316 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw24);
317
318 //compute mobilities of face 1
319 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
320 lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
321 lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
322
323 //compute mobilities of face 4
324 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
325 lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
326 lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
327
328 //compute mobilities of face 2
329 Dune::FieldVector<Scalar, numPhases> lambda24Upw(0.0);
330 lambda24Upw[wPhaseIdx] = (potentialDiffW24 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
331 lambda24Upw[nPhaseIdx] = (potentialDiffNw24 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
332
333 for (int i = 0; i < numPhases; i++)
334 {
335 // evaluate parts of velocity --> always take the normal for which the flux is calculated!
336 DimVector vel12 = interactionVolume.getNormal(0, 0);
337 DimVector vel14 = interactionVolume.getNormal(3, 0);
338 DimVector vel24 = interactionVolume.getNormal(1, 0);
339 DimVector vel21 = interactionVolume.getNormal(0, 0);
340 DimVector vel41 = interactionVolume.getNormal(3, 0);
341 DimVector vel42 = interactionVolume.getNormal(1, 0);
342
343 Dune::FieldVector<Scalar, 3> flux(0);
344 switch (i)
345 {
346 case wPhaseIdx:
347 {
348 flux = fluxW;
349 break;
350 }
351 case nPhaseIdx:
352 {
353 flux = fluxNw;
354 break;
355 }
356 }
357
358 vel12 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0)); //divide by 2 because the flux is related to the half face!
359 vel14 *= flux[2] / (2 * interactionVolume.getFaceArea(3, 0));
360 vel24 *= flux[1] / (2 * interactionVolume.getFaceArea(1, 0));
361 vel21 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0));
362 vel41 *= flux[2] / (4 * interactionVolume.getFaceArea(3, 0));
363 vel42 *= flux[1] / (4 * interactionVolume.getFaceArea(1, 0));
364
365 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
366 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
367 Scalar lambdaT24 = lambda24Upw[wPhaseIdx] + lambda24Upw[nPhaseIdx];
368 Scalar fracFlow12 = (lambdaT12 > ParentType::threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
369 Scalar fracFlow14 = (lambdaT14 > ParentType::threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
370 Scalar fracFlow24 = (lambdaT24 > ParentType::threshold_) ? lambda24Upw[i] / (lambdaT24) : 0.0;
371
372 vel12 *= fracFlow12;
373 vel14 *= fracFlow14;
374 vel24 *= fracFlow24;
375 vel21 *= fracFlow12;
376 vel41 *= fracFlow14;
377 vel42 *= fracFlow24;
378
379 if (innerBoundaryVolumeFaces[globalIdx1][interactionVolume.getIndexOnElement(0, 0)])
380 {
381 vel12 *= 2;
382 vel21 *= 2;
383 }
384 if (innerBoundaryVolumeFaces[globalIdx1][interactionVolume.getIndexOnElement(0, 1)])
385 {
386 vel14 *= 2;
387 vel41 *= 2;
388 }
389 if (innerBoundaryVolumeFaces[globalIdx2][interactionVolume.getIndexOnElement(1, 0)])
390 {
391 vel24 *= 2;
392 vel42 *= 2;
393 }
394
395 //store velocities
396 //set velocity
397 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
398 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel14);
399 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel24);
400 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
401 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel41);
402 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel42);
403 }
404 //set velocity marker
405 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
406 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
407 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
408 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
409 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
410 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
411}
412
413} // end namespace Dumux
414#endif
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::FvMpfaL2dTransmissibilityCalculator
Provides methods for transmissibility calculation in 2-d.
Definition: 2dtransmissibilitycalculator.hh:44
Dumux::FvMpfaL2dVelocity2p
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: lmethod/2dvelocity.hh:58
Dumux::FvMpfaL2dVelocity2pAdaptive
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: 2dvelocityadaptive.hh:54
Dumux::FvMpfaL2dVelocity2pAdaptive::calculateHangingNodeInteractionVolumeVelocity
void calculateHangingNodeInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData4, InnerBoundaryVolumeFaces &innerBoundaryVolumeFaces)
Calculatex velocities for flux faces of a hanging node interaction volume.
Definition: 2dvelocityadaptive.hh:148
Dumux::FvMpfaL2dVelocity2pAdaptive::FvMpfaL2dVelocity2pAdaptive
FvMpfaL2dVelocity2pAdaptive(Problem &problem)
Constructs a FvMpfaL2dVelocity2pAdaptive object.
Definition: 2dvelocityadaptive.hh:121
Dumux::FVMPFALInteractionVolume
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
Dumux::FVMPFALInteractionVolume::getIndexOnElement
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering to numbering of the Dune reference element.
Definition: linteractionvolume.hh:258
Dumux::FVMPFALInteractionVolume::getNormal
DimVector & getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face normal.
Definition: linteractionvolume.hh:364
Dumux::FVMPFALInteractionVolume::getSubVolumeElement
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: linteractionvolume.hh:281
Dumux::FVMPFALInteractionVolume::getFaceArea
Scalar & getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face area.
Definition: linteractionvolume.hh:388
2dvelocity.hh
Velocity calculation using a 2-d MPFA O-method.
Generated by doxygen 1.9.3