3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
operatoradaptive.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_MIMETICOPERATOR2PADAPTIVE_HH
25#define DUMUX_MIMETICOPERATOR2PADAPTIVE_HH
26
27#include "croperatoradaptive.hh"
30
31namespace Dumux {
32
43template<class TypeTag>
45{
47
48 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
49 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
50 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
51
52 enum
53 {
54 dim=GridView::dimension,
55 dimWorld=GridView::dimensionworld
56 };
57 using LocalStiffness = typename GET_PROP_TYPE(TypeTag, LocalStiffness);
58
59 using Element = typename GridView::template Codim<0>::Entity;
60
61 using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
62 using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
63 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
64 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
65 using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
66 using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
67 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
68
69 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
70
71 enum
72 {
73 pw = Indices::pressureW,
74 pn = Indices::pressureNw,
75 pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation),
76 wPhaseIdx = Indices::wPhaseIdx,
77 nPhaseIdx = Indices::nPhaseIdx,
78 saturationIdx = Indices::saturationIdx,
79 satEqIdx = Indices::satEqIdx
80 };
81
82 using FieldVector = Dune::FieldVector<Scalar, dimWorld>;
83public:
84
85 MimeticOperatorAssemblerTwoPAdaptive (const GridView& gridView)
86 : ParentType(gridView)
87 {}
88
89 // TODO doc me!
90 template<class Vector>
91 void calculatePressure(LocalStiffness& loc, Vector& u, Problem& problem)
92 {
93 Dune::DynamicVector<Scalar> velocityW(2*dim);
94 Dune::DynamicVector<Scalar> velocityNw(2*dim);
95 Dune::DynamicVector<Scalar> pressTraceW(2*dim);
96 Dune::DynamicVector<Scalar> pressTraceNw(2*dim);
97
98 const auto firstElement = *problem.gridView().template begin<0>();
99 FluidState fluidState;
100 fluidState.setPressure(wPhaseIdx, problem.referencePressure(firstElement));
101 fluidState.setPressure(nPhaseIdx, problem.referencePressure(firstElement));
102 fluidState.setTemperature(problem.temperature(firstElement));
103 fluidState.setSaturation(wPhaseIdx, 1.);
104 fluidState.setSaturation(nPhaseIdx, 0.);
105 Scalar densityDiff = FluidSystem::density(fluidState, nPhaseIdx) - FluidSystem::density(fluidState, wPhaseIdx);
106 Scalar viscosityW = FluidSystem::viscosity(fluidState, wPhaseIdx);
107 Scalar viscosityNw = FluidSystem::viscosity(fluidState, nPhaseIdx);
108
109 // reset velocity
110 for (int i = 0; i < problem.gridView().size(0); i++)
111 {
112 problem.variables().cellData(i).fluxData().resetVelocity();
113 }
114
115 // run over all level elements
116 for (const auto& element : elements(this->gridView_))
117 {
118 int eIdxGlobal = problem.variables().index(element);
119
120 unsigned int numFaces = this->intersectionMapper_.size(eIdxGlobal);
121
122 // get local to global id map and pressure traces
123 velocityW.resize(numFaces);
124 velocityNw.resize(numFaces);
125 pressTraceW.resize(numFaces);
126 pressTraceNw.resize(numFaces);
127
128 CellData& cellData = problem.variables().cellData(eIdxGlobal);
129 FieldVector globalPos = element.geometry().center();
130
131 int intersectionIdx = -1;
132 // get local to global id map and pressure traces
133 for (const auto& intersection : intersections(problem.gridView(), element))
134 {
135 ++intersectionIdx;
136
137 int fIdxGlobal = this->intersectionMapper_.subIndex(element, intersectionIdx);
138
139 Scalar pcPotFace = (problem.bBoxMax() - intersection.geometry().center()) * problem.gravity() * densityDiff;
140
141 switch (pressureType)
142 {
143 case pw:
144 {
145 pressTraceW[intersectionIdx] = u[fIdxGlobal];
146 pressTraceNw[intersectionIdx] = u[fIdxGlobal] + pcPotFace;
147 break;
148 }
149 case pn:
150 {
151 pressTraceNw[intersectionIdx] = u[fIdxGlobal];
152 pressTraceW[intersectionIdx] = u[fIdxGlobal] - pcPotFace;
153
154 break;
155 }
156 }
157 }
158
159 switch (pressureType)
160 {
161 case pw:
162 {
163 Scalar potW = loc.constructPressure(element, pressTraceW);
164 Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
165 Scalar potNw = potW + gravPot;
166
167 cellData.setPotential(wPhaseIdx, potW);
168 cellData.setPotential(nPhaseIdx, potNw);
169
170 gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
171
172 cellData.setPressure(wPhaseIdx, potW - gravPot);
173
174 gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
175
176 cellData.setPressure(nPhaseIdx, potNw - gravPot);
177
178 break;
179 }
180 case pn:
181 {
182 Scalar potNw = loc.constructPressure(element, pressTraceNw);
183 Scalar gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * densityDiff;
184 Scalar potW = potNw - gravPot;
185
186 cellData.setPotential(nPhaseIdx, potNw);
187 cellData.setPotential(wPhaseIdx, potW);
188
189 gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, wPhaseIdx);
190
191 cellData.setPressure(wPhaseIdx, potW - gravPot);
192
193 gravPot = (problem.bBoxMax() - globalPos) * problem.gravity() * FluidSystem::density(fluidState, nPhaseIdx);
194
195 cellData.setPressure(nPhaseIdx, potNw - gravPot);
196
197 break;
198 }
199 }
200
201 // velocity reconstruction: !!! The velocity which is not reconstructed from the primary
202 // pressure variable can be slightly wrong and not conservative!!!!
203 // -> Should not be used for transport!!
204 loc.constructVelocity(element, velocityW, pressTraceW, cellData.potential(wPhaseIdx));
205 loc.constructVelocity(element, velocityNw, pressTraceNw, cellData.potential(nPhaseIdx));
206
207 intersectionIdx = -1;
208 for (const auto& intersection : intersections(problem.gridView(), element))
209 {
210 ++intersectionIdx;
211 int idxInInside = intersection.indexInInside();
212
213 cellData.fluxData().addUpwindPotential(wPhaseIdx, idxInInside, velocityW[intersectionIdx]);
214 cellData.fluxData().addUpwindPotential(nPhaseIdx, idxInInside, velocityNw[intersectionIdx]);
215
216 Scalar mobilityW = 0;
217 Scalar mobilityNw = 0;
218
219 if (intersection.neighbor())
220 {
221 int neighborIdx = problem.variables().index(intersection.outside());
222
223 CellData& cellDataNeighbor = problem.variables().cellData(neighborIdx);
224
225 mobilityW =
226 (velocityW[intersectionIdx] >= 0.) ? cellData.mobility(wPhaseIdx) :
227 cellDataNeighbor.mobility(wPhaseIdx);
228 mobilityNw =
229 (velocityNw[intersectionIdx] >= 0.) ? cellData.mobility(nPhaseIdx) :
230 cellDataNeighbor.mobility(nPhaseIdx);
231
232 if (velocityW[intersectionIdx] >= 0.)
233 {
234 FieldVector velocity(intersection.centerUnitOuterNormal());
235 velocity *= mobilityW/(mobilityW+mobilityNw) * velocityW[intersectionIdx];
236 cellData.fluxData().addVelocity(wPhaseIdx, idxInInside, velocity);
237 cellDataNeighbor.fluxData().addVelocity(wPhaseIdx, intersection.indexInOutside(), velocity);
238 }
239 if (velocityNw[intersectionIdx] >= 0.)
240 {
241 FieldVector velocity(intersection.centerUnitOuterNormal());
242 velocity *= mobilityNw/(mobilityW+mobilityNw) * velocityNw[intersectionIdx];
243 cellData.fluxData().addVelocity(nPhaseIdx, idxInInside, velocity);
244 cellDataNeighbor.fluxData().addVelocity(nPhaseIdx, intersection.indexInOutside(), velocity);
245 }
246
247 cellData.fluxData().setVelocityMarker(idxInInside);
248 }
249 else
250 {
251 BoundaryTypes bctype;
252 problem.boundaryTypes(bctype, intersection);
253 if (bctype.isDirichlet(satEqIdx))
254 {
255 PrimaryVariables boundValues(0.0);
256 problem.dirichlet(boundValues, intersection);
257
258 if (velocityW[intersectionIdx] >= 0.)
259 {
260 mobilityW = cellData.mobility(wPhaseIdx);
261 }
262 else
263 {
264 mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(element),
265 boundValues[saturationIdx]) / viscosityW;
266 }
267
268 if (velocityNw[intersectionIdx] >= 0.)
269 {
270 mobilityNw = cellData.mobility(nPhaseIdx);
271 }
272 else
273 {
274 mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(element),
275 boundValues[saturationIdx]) / viscosityNw;
276 }
277 }
278 else
279 {
280 mobilityW = cellData.mobility(wPhaseIdx);
281 mobilityNw = cellData.mobility(nPhaseIdx);
282 }
283
284 FieldVector velocity(intersection.centerUnitOuterNormal());
285 velocity *= mobilityW/(mobilityW+mobilityNw) * velocityW[intersectionIdx];
286 cellData.fluxData().addVelocity(wPhaseIdx, idxInInside, velocity);
287
288
289 velocity = intersection.centerUnitOuterNormal();
290 velocity *= mobilityNw/(mobilityW+mobilityNw) * velocityNw[intersectionIdx];
291 cellData.fluxData().addVelocity(nPhaseIdx, idxInInside, velocity);
292 cellData.fluxData().setVelocityMarker(idxInInside);
293 }
294 }
295 }
296 }
297};
298} // end namespace Dumux
299#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:281
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Defines a class for Crozieux-Raviart piecewise linear finite element functions.
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
int subIndex(int elemIdx, int fIdx)
Map interface fIdx'th interface of element index to array index.
Definition: intersectionmapper.hh:236
unsigned int size() const
Definition: intersectionmapper.hh:292
Extends CROperatorBase by a generic methods to assemble global stiffness matrix from local stiffness ...
Definition: croperatoradaptive.hh:69
const GridView gridView_
Definition: croperatoradaptive.hh:163
IntersectionMapper intersectionMapper_
Definition: croperatoradaptive.hh:165
Levelwise assembler.
Definition: operatoradaptive.hh:45
void calculatePressure(LocalStiffness &loc, Vector &u, Problem &problem)
Definition: operatoradaptive.hh:91
MimeticOperatorAssemblerTwoPAdaptive(const GridView &gridView)
Definition: operatoradaptive.hh:85
Specifies the properties for immiscible 2p diffusion/pressure models.
Defines the basic properties required for a mimetic method.