3.2-git
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 using MaterialLaw = typename SpatialParams::MaterialLaw;
67
68 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
69
70 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
71 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
72
73 using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
74 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
75 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
76 using CellData = GetPropType<TypeTag, Properties::CellData>;
77
78 using Element = typename GridView::Traits::template Codim<0>::Entity;
79 using Grid = typename GridView::Grid;
80 using IndexSet = typename GridView::IndexSet;
81
82 using Geometry = typename Element::Geometry;
83 using JacobianTransposed = typename Geometry::JacobianTransposed ;
84
85 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
86
87 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
88 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
89 using TransmissibilityCalculator = FvMpfaL2dTransmissibilityCalculator<TypeTag>;
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 DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
114 using DimVector = Dune::FieldVector<Scalar, dim>;
115
116public:
122 FvMpfaL2dVelocity2pAdaptive(Problem& problem) :
123 ParentType(problem), problem_(problem)
124 {}
125
127 void calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1,
128 CellData& cellData2, CellData& cellData4,
129 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
130
131private:
132 Problem& problem_;
133};
134// end of template
135
148template<class TypeTag>
149void FvMpfaL2dVelocity2pAdaptive<TypeTag>::calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume,
150 CellData& cellData1, CellData& cellData2,
151 CellData& cellData4,
152 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
153{
154 // cell index
155 int globalIdx1 = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
156 int globalIdx2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
157
158 // get pressure values
159 Dune::FieldVector < Scalar, 2 * dim > potW(0);
160 Dune::FieldVector < Scalar, 2 * dim > potNw(0);
161
162 potW[0] = cellData1.potential(wPhaseIdx);
163 potW[1] = cellData2.potential(wPhaseIdx);
164 potW[2] = cellData4.potential(wPhaseIdx);
165
166 potNw[0] = cellData1.potential(nPhaseIdx);
167 potNw[1] = cellData2.potential(nPhaseIdx);
168 potNw[2] = cellData4.potential(nPhaseIdx);
169
170 //get mobilities of the phases
171 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
172 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
173
174 //compute total mobility of cell 1
175 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
176
177 //get mobilities of the phases
178 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
179 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
180
181 //compute total mobility of cell 1
182 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
183
184 //get mobilities of the phases
185 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
186 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
187
188 //compute total mobility of cell 1
189 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
190
191 std::vector < DimVector > lambda(4);
192
193 lambda[0][0] = lambdaTotal1;
194 lambda[0][1] = lambdaTotal1;
195 lambda[1][0] = lambdaTotal2;
196 lambda[1][1] = lambdaTotal2;
197 lambda[3][0] = lambdaTotal4;
198 lambda[3][1] = lambdaTotal4;
199
200 Scalar potentialDiffW12 = 0;
201 Scalar potentialDiffW14 = 0;
202 Scalar potentialDiffW24 = 0;
203
204 Scalar potentialDiffNw12 = 0;
205 Scalar potentialDiffNw14 = 0;
206 Scalar potentialDiffNw24 = 0;
207
208 //flux vector
209 Dune::FieldVector<Scalar, 3> fluxW(0);
210 Dune::FieldVector<Scalar, 3> fluxNw(0);
211
212 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> T(0);
213 DimVector Tu(0);
214 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
215
216 int transmissibilityType = this->transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 3, 3);
217
218 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
219 {
220 u[0] = potW[1];
221 u[1] = potW[2];
222 u[2] = potW[0];
223
224 T.mv(u, Tu);
225
226 fluxW[0] = Tu[1];
227 potentialDiffW12 = Tu[1];
228
229 u[0] = potNw[1];
230 u[1] = potNw[2];
231 u[2] = potNw[0];
232
233 T.mv(u, Tu);
234
235 fluxNw[0] = Tu[1];
236 potentialDiffNw12 = Tu[1];
237 }
238 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
239 {
240 u[0] = potW[0];
241 u[1] = potW[2];
242 u[2] = potW[1];
243
244 T.mv(u, Tu);
245
246 fluxW[0] = Tu[1];
247 potentialDiffW12 = Tu[1];
248
249 u[0] = potNw[0];
250 u[1] = potNw[2];
251 u[2] = potNw[1];
252
253 T.mv(u, Tu);
254
255 fluxNw[0] = Tu[1];
256 potentialDiffNw12 = Tu[1];
257 }
258
259 transmissibilityType = this->transmissibilityCalculator_.calculateLeftHNTransmissibility(T, interactionVolume, lambda, 1, 3, 0);
260
261 if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
262 {
263 u[0] = potW[1];
264 u[1] = potW[0];
265 u[2] = potW[2];
266
267 T.mv(u, Tu);
268
269 fluxW[1] = Tu[1];
270 potentialDiffW24 = Tu[1];
271
272 u[0] = potNw[1];
273 u[1] = potNw[0];
274 u[2] = potNw[2];
275
276 T.mv(u, Tu);
277
278 fluxNw[1] = Tu[1];
279 potentialDiffNw24 = Tu[1];
280 }
281
282 transmissibilityType = this->transmissibilityCalculator_.calculateRightHNTransmissibility(T, interactionVolume, lambda, 3, 0, 1);
283
284 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
285 {
286 u[0] = potW[0];
287 u[1] = potW[1];
288 u[2] = potW[2];
289
290 T.mv(u, Tu);
291
292 fluxW[2] = Tu[1];
293 potentialDiffW14 = -Tu[1];
294
295 u[0] = potNw[0];
296 u[1] = potNw[1];
297 u[2] = potNw[2];
298
299 T.mv(u, Tu);
300
301 fluxNw[2] = Tu[1];
302 potentialDiffNw14 = -Tu[1];
303 }
304
305 //store potentials for further calculations (saturation, ...) -> maybe add new potential to old one!!
306 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW12);
307 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw12);
308 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffW14);
309 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNw14);
310 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffW24);
311 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffNw24);
312 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW12);
313 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw12);
314 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffW14);
315 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNw14);
316 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW24);
317 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw24);
318
319 //compute mobilities of face 1
320 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
321 lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
322 lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
323
324 //compute mobilities of face 4
325 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
326 lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
327 lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
328
329 //compute mobilities of face 2
330 Dune::FieldVector<Scalar, numPhases> lambda24Upw(0.0);
331 lambda24Upw[wPhaseIdx] = (potentialDiffW24 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
332 lambda24Upw[nPhaseIdx] = (potentialDiffNw24 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
333
334 for (int i = 0; i < numPhases; i++)
335 {
336 // evaluate parts of velocity --> always take the normal for which the flux is calculated!
337 DimVector vel12 = interactionVolume.getNormal(0, 0);
338 DimVector vel14 = interactionVolume.getNormal(3, 0);
339 DimVector vel24 = interactionVolume.getNormal(1, 0);
340 DimVector vel21 = interactionVolume.getNormal(0, 0);
341 DimVector vel41 = interactionVolume.getNormal(3, 0);
342 DimVector vel42 = interactionVolume.getNormal(1, 0);
343
344 Dune::FieldVector<Scalar, 3> flux(0);
345 switch (i)
346 {
347 case wPhaseIdx:
348 {
349 flux = fluxW;
350 break;
351 }
352 case nPhaseIdx:
353 {
354 flux = fluxNw;
355 break;
356 }
357 }
358
359 vel12 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0)); //divide by 2 because the flux is related to the half face!
360 vel14 *= flux[2] / (2 * interactionVolume.getFaceArea(3, 0));
361 vel24 *= flux[1] / (2 * interactionVolume.getFaceArea(1, 0));
362 vel21 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0));
363 vel41 *= flux[2] / (4 * interactionVolume.getFaceArea(3, 0));
364 vel42 *= flux[1] / (4 * interactionVolume.getFaceArea(1, 0));
365
366 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
367 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
368 Scalar lambdaT24 = lambda24Upw[wPhaseIdx] + lambda24Upw[nPhaseIdx];
369 Scalar fracFlow12 = (lambdaT12 > ParentType::threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
370 Scalar fracFlow14 = (lambdaT14 > ParentType::threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
371 Scalar fracFlow24 = (lambdaT24 > ParentType::threshold_) ? lambda24Upw[i] / (lambdaT24) : 0.0;
372
373 vel12 *= fracFlow12;
374 vel14 *= fracFlow14;
375 vel24 *= fracFlow24;
376 vel21 *= fracFlow12;
377 vel41 *= fracFlow14;
378 vel42 *= fracFlow24;
379
380 if (innerBoundaryVolumeFaces[globalIdx1][interactionVolume.getIndexOnElement(0, 0)])
381 {
382 vel12 *= 2;
383 vel21 *= 2;
384 }
385 if (innerBoundaryVolumeFaces[globalIdx1][interactionVolume.getIndexOnElement(0, 1)])
386 {
387 vel14 *= 2;
388 vel41 *= 2;
389 }
390 if (innerBoundaryVolumeFaces[globalIdx2][interactionVolume.getIndexOnElement(1, 0)])
391 {
392 vel24 *= 2;
393 vel42 *= 2;
394 }
395
396 //store velocities
397 //set velocity
398 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
399 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel14);
400 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel24);
401 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
402 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel41);
403 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel42);
404 }
405 //set velocity marker
406 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
407 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
408 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
409 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
410 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
411 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
412}
413
414} // end namespace Dumux
415#endif
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::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:149
Dumux::FvMpfaL2dVelocity2pAdaptive::FvMpfaL2dVelocity2pAdaptive
FvMpfaL2dVelocity2pAdaptive(Problem &problem)
Constructs a FvMpfaL2dVelocity2pAdaptive object.
Definition: 2dvelocityadaptive.hh:122
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