transfluxmodule.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 This file is part of the Open Porous Media project (OPM).
5
6 OPM 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 2 of the License, or
9 (at your option) any later version.
10
11 OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
31#ifndef EWOMS_TRANS_FLUX_MODULE_HH
32#define EWOMS_TRANS_FLUX_MODULE_HH
33
36#include <opm/material/common/Valgrind.hpp>
37#include <dune/common/fvector.hh>
38#include <dune/common/fmatrix.hh>
39
40namespace Opm {
41
42template <class TypeTag>
43class TransIntensiveQuantities;
44
45template <class TypeTag>
46class TransExtensiveQuantities;
47
48template <class TypeTag>
49class TransBaseProblem;
50
54template <class TypeTag>
56{
63 static void registerParameters()
64 { }
65};
66
71template <class TypeTag>
73{ };
74
78template <class TypeTag>
80{
82protected:
83 void update_(const ElementContext&, unsigned, unsigned)
84 { }
85};
86
90template <class TypeTag>
92{
94
102
103 enum { dimWorld = GridView::dimensionworld };
104 enum { numPhases = FluidSystem::numPhases };
105
106 typedef MathToolbox<Evaluation> Toolbox;
107 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
108 typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
109 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
110
111public:
115 const DimMatrix& intrinsicPermeability() const
116 {
117 throw std::logic_error("The ECL transmissibility module does not provide an explicit intrinsic permeability");
118 }
119
126 const EvalDimVector& potentialGrad(unsigned) const
127 {
128 throw std::logic_error("The ECL transmissibility module does not provide explicit potential gradients");
129 }
130
137 const Evaluation& pressureDifference(unsigned phaseIdx) const
138 { return pressureDifference_[phaseIdx]; }
139
146 const EvalDimVector& filterVelocity(unsigned) const
147 {
148 throw std::logic_error("The ECL transmissibility module does not provide explicit filter velocities");
149 }
150
160 const Evaluation& volumeFlux(unsigned phaseIdx) const
161 { return volumeFlux_[phaseIdx]; }
162
163protected:
171 unsigned upstreamIndex_(unsigned phaseIdx) const
172 {
173 assert(phaseIdx < numPhases);
174
175 return upIdx_[phaseIdx];
176 }
177
185 unsigned downstreamIndex_(unsigned phaseIdx) const
186 {
187 assert(phaseIdx < numPhases);
188
189 return dnIdx_[phaseIdx];
190 }
191
192 void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
193 { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
194
195 void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
196 { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
197
201 void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
202 {
203 Valgrind::SetUndefined(*this);
204
205 // only valied for element center finite volume discretization
206 static const bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
207
208 static_assert(isEcfv);
209
210 const auto& stencil = elemCtx.stencil(timeIdx);
211 const auto& scvf = stencil.interiorFace(scvfIdx);
212
213 interiorDofIdx_ = scvf.interiorIndex();
214 exteriorDofIdx_ = scvf.exteriorIndex();
215 assert(interiorDofIdx_ != exteriorDofIdx_);
216
217 unsigned I = stencil.globalSpaceIndex(interiorDofIdx_);
218 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_);
219
220 Scalar trans = transmissibility_(elemCtx, scvfIdx, timeIdx);
221
222 // estimate the gravity correction: for performance reasons we use a simplified
223 // approach for this flux module that assumes that gravity is constant and always
224 // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
225 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
226
227 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
228 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx_, timeIdx);
229
230 Scalar zIn = dofCenterDepth_(elemCtx, interiorDofIdx_, timeIdx);
231 Scalar zEx = dofCenterDepth_(elemCtx, exteriorDofIdx_, timeIdx);
232 // the distances from the DOF's depths. (i.e., the additional depth of the
233 // exterior DOF)
234 Scalar distZ = zIn - zEx;
235
236 for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
237 if (!FluidSystem::phaseIsActive(phaseIdx))
238 continue;
239
240 // check shortcut: if the mobility of the phase is zero in the interior as
241 // well as the exterior DOF, we can skip looking at the phase.
242 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
243 intQuantsEx.mobility(phaseIdx) <= 0.0)
244 {
245 upIdx_[phaseIdx] = interiorDofIdx_;
246 dnIdx_[phaseIdx] = exteriorDofIdx_;
247 pressureDifference_[phaseIdx] = 0.0;
248 volumeFlux_[phaseIdx] = 0.0;
249 continue;
250 }
251
252 // do the gravity correction: compute the hydrostatic pressure for the
253 // external at the depth of the internal one
254 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
255 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
256 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
257
258 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
259 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
260
261 pressureExterior += rhoAvg*(distZ*g);
262
263 pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
264
265 // decide the upstream index for the phase. for this we make sure that the
266 // degree of freedom which is regarded upstream if both pressures are equal
267 // is always the same: if the pressure is equal, the DOF with the lower
268 // global index is regarded to be the upstream one.
269 if (pressureDifference_[phaseIdx] > 0.0) {
270 upIdx_[phaseIdx] = exteriorDofIdx_;
271 dnIdx_[phaseIdx] = interiorDofIdx_;
272 }
273 else if (pressureDifference_[phaseIdx] < 0.0) {
274 upIdx_[phaseIdx] = interiorDofIdx_;
275 dnIdx_[phaseIdx] = exteriorDofIdx_;
276 }
277 else {
278 // if the pressure difference is zero, we chose the DOF which has the
279 // larger volume associated to it as upstream DOF
280 Scalar Vin = elemCtx.dofVolume(interiorDofIdx_, /*timeIdx=*/0);
281 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx_, /*timeIdx=*/0);
282 if (Vin > Vex) {
283 upIdx_[phaseIdx] = interiorDofIdx_;
284 dnIdx_[phaseIdx] = exteriorDofIdx_;
285 }
286 else if (Vin < Vex) {
287 upIdx_[phaseIdx] = exteriorDofIdx_;
288 dnIdx_[phaseIdx] = interiorDofIdx_;
289 }
290 else {
291 assert(Vin == Vex);
292 // if the volumes are also equal, we pick the DOF which exhibits the
293 // smaller global index
294 if (I < J) {
295 upIdx_[phaseIdx] = interiorDofIdx_;
296 dnIdx_[phaseIdx] = exteriorDofIdx_;
297 }
298 else {
299 upIdx_[phaseIdx] = exteriorDofIdx_;
300 dnIdx_[phaseIdx] = interiorDofIdx_;
301 }
302 }
303 }
304
305 // this is slightly hacky because in the automatic differentiation case, it
306 // only works for the element centered finite volume method. for ebos this
307 // does not matter, though.
308 unsigned upstreamIdx = upstreamIndex_(phaseIdx);
309 const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
310
311 if (upstreamIdx == interiorDofIdx_)
312 volumeFlux_[phaseIdx] =
313 pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans);
314 else
315 volumeFlux_[phaseIdx] =
316 pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*(-trans));
317 }
318 }
319
323 template <class FluidState>
324 void calculateBoundaryGradients_(const ElementContext& elemCtx,
325 unsigned scvfIdx,
326 unsigned timeIdx,
327 const FluidState& exFluidState)
328 {
329 const auto& stencil = elemCtx.stencil(timeIdx);
330 const auto& scvf = stencil.boundaryFace(scvfIdx);
331
332 interiorDofIdx_ = scvf.interiorIndex();
333
334 Scalar trans = transmissibilityBoundary_(elemCtx, scvfIdx, timeIdx);
335
336 // estimate the gravity correction: for performance reasons we use a simplified
337 // approach for this flux module that assumes that gravity is constant and always
338 // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
339 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
340
341 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
342
343 // this is quite hacky because the dune grid interface does not provide a
344 // cellCenterDepth() method (so we ask the problem to provide it). The "good"
345 // solution would be to take the Z coordinate of the element centroids, but since
346 // ECL seems to like to be inconsistent on that front, it needs to be done like
347 // here...
348 Scalar zIn = dofCenterDepth_(elemCtx, interiorDofIdx_, timeIdx);
349 Scalar zEx = scvf.integrationPos()[dimWorld - 1];
350
351 // the distances from the DOF's depths. (i.e., the additional depth of the
352 // exterior DOF)
353 Scalar distZ = zIn - zEx;
354
355 for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
356 if (!FluidSystem::phaseIsActive(phaseIdx))
357 continue;
358
359 // do the gravity correction: compute the hydrostatic pressure for the
360 // integration position
361 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
362 const auto& rhoEx = exFluidState.density(phaseIdx);
363 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
364
365 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
366 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
367 pressureExterior += rhoAvg*(distZ*g);
368
369 pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
370
371 // decide the upstream index for the phase. for this we make sure that the
372 // degree of freedom which is regarded upstream if both pressures are equal
373 // is always the same: if the pressure is equal, the DOF with the lower
374 // global index is regarded to be the upstream one.
375 if (pressureDifference_[phaseIdx] > 0.0) {
376 upIdx_[phaseIdx] = -1;
377 dnIdx_[phaseIdx] = interiorDofIdx_;
378 }
379 else {
380 upIdx_[phaseIdx] = interiorDofIdx_;
381 dnIdx_[phaseIdx] = -1;
382 }
383
384 short upstreamIdx = upstreamIndex_(phaseIdx);
385 if (upstreamIdx == interiorDofIdx_) {
386
387 // this is slightly hacky because in the automatic differentiation case, it
388 // only works for the element centered finite volume method. for ebos this
389 // does not matter, though.
390 const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
391
392 volumeFlux_[phaseIdx] =
393 pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans);
394
395 }
396 else {
397 // compute the phase mobility using the material law parameters of the
398 // interior element. \todo {this could probably be done more efficiently}
399 const auto& matParams =
400 elemCtx.problem().materialLawParams(elemCtx,
401 interiorDofIdx_,
402 /*timeIdx=*/0);
403 std::array<typename FluidState::Scalar,numPhases> kr;
404 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
405
406 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
407 volumeFlux_[phaseIdx] =
408 pressureDifference_[phaseIdx]*mob*(-trans);
409
410 }
411 }
412 }
413
417 void calculateFluxes_(const ElementContext&, unsigned, unsigned)
418 { }
419
420 void calculateBoundaryFluxes_(const ElementContext&, unsigned, unsigned)
421 {}
422
423private:
424
425 Scalar transmissibility_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) const
426 {
427 const auto& stencil = elemCtx.stencil(timeIdx);
428 const auto& face = stencil.interiorFace(scvfIdx);
429 const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).globalPos();
430 const auto& exteriorPos = stencil.subControlVolume(face.exteriorIndex()).globalPos();
431 auto distVec0 = face.integrationPos() - interiorPos;
432 auto distVec1 = face.integrationPos() - exteriorPos;
433 Scalar ndotDistIn = std::abs(face.normal() * distVec0);
434 Scalar ndotDistExt = std::abs(face.normal() * distVec1);
435
436 Scalar distSquaredIn = distVec0 * distVec0;
437 Scalar distSquaredExt = distVec1 * distVec1;
438 const auto& K0mat = elemCtx.problem().intrinsicPermeability(elemCtx, face.interiorIndex(), timeIdx);
439 const auto& K1mat = elemCtx.problem().intrinsicPermeability(elemCtx, face.exteriorIndex(), timeIdx);
440 // the permeability per definition aligns with the grid
441 // we only support diagonal permeability tensor
442 // and can therefore neglect off-diagonal values
443 int idx = 0;
444 Scalar val = 0.0;
445 for (unsigned i = 0; i < dimWorld; ++ i){
446 if (std::abs(face.normal()[i]) > val) {
447 val = std::abs(face.normal()[i]);
448 idx = i;
449 }
450 }
451 const Scalar& K0 = K0mat[idx][idx];
452 const Scalar& K1 = K1mat[idx][idx];
453 const Scalar T0 = K0 * ndotDistIn / distSquaredIn;
454 const Scalar T1 = K1 * ndotDistExt / distSquaredExt;
455 return T0 * T1 / (T0 + T1);
456 }
457 Scalar transmissibilityBoundary_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) const
458 {
459 const auto& stencil = elemCtx.stencil(timeIdx);
460 const auto& face = stencil.interiorFace(scvfIdx);
461 const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).globalPos();
462 auto distVec0 = face.integrationPos() - interiorPos;
463 Scalar ndotDistIn = face.normal() * distVec0;
464 Scalar distSquaredIn = distVec0 * distVec0;
465 const auto& K0mat = elemCtx.problem().intrinsicPermeability(elemCtx, face.interiorIndex(), timeIdx);
466 // the permeability per definition aligns with the grid
467 // we only support diagonal permeability tensor
468 // and can therefore neglect off-diagonal values
469 int idx = 0;
470 Scalar val = 0.0;
471 for (unsigned i = 0; i < dimWorld; ++ i){
472 if (std::abs(face.normal()[i]) > val) {
473 val = std::abs(face.normal()[i]);
474 idx = i;
475 }
476 }
477 const Scalar& K0 = K0mat[idx][idx];
478 const Scalar T0 = K0 * ndotDistIn / distSquaredIn;
479 return T0;
480 }
481
482 template <class Context>
483 Scalar dofCenterDepth_(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
484 {
485 const auto& pos = context.pos(spaceIdx, timeIdx);
486 return pos[dimWorld-1];
487 }
488
489 Implementation& asImp_()
490 { return *static_cast<Implementation*>(this); }
491
492 const Implementation& asImp_() const
493 { return *static_cast<const Implementation*>(this); }
494
495 // the volumetric flux of all phases [m^3/s]
496 Evaluation volumeFlux_[numPhases];
497
498 // the difference in effective pressure between the exterior and the interior degree
499 // of freedom [Pa]
500 Evaluation pressureDifference_[numPhases];
501
502 // the local indices of the interior and exterior degrees of freedom
503 unsigned short interiorDofIdx_;
504 unsigned short exteriorDofIdx_;
505 short upIdx_[numPhases];
506 short dnIdx_[numPhases];
507};
508
509} // namespace Opm
510
511#endif
Provides the defaults for the parameters required by the transmissibility based volume flux calculati...
Definition: transfluxmodule.hh:73
Provides the transmissibility based flux module.
Definition: transfluxmodule.hh:92
const Evaluation & volumeFlux(unsigned phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: transfluxmodule.hh:160
unsigned downstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in downstream direction.
Definition: transfluxmodule.hh:185
const EvalDimVector & filterVelocity(unsigned) const
Return the filter velocity of a fluid phase at the face's integration point [m/s].
Definition: transfluxmodule.hh:146
void calculateGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Update the required gradients for interior faces.
Definition: transfluxmodule.hh:201
void updateSolvent(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: transfluxmodule.hh:192
const Evaluation & pressureDifference(unsigned phaseIdx) const
Return the gravity corrected pressure difference between the interior and the exterior of a face.
Definition: transfluxmodule.hh:137
void calculateFluxes_(const ElementContext &, unsigned, unsigned)
Update the volumetric fluxes for all fluid phases on the interior faces of the context.
Definition: transfluxmodule.hh:417
void calculateBoundaryFluxes_(const ElementContext &, unsigned, unsigned)
Definition: transfluxmodule.hh:420
const EvalDimVector & potentialGrad(unsigned) const
Return the pressure potential gradient of a fluid phase at the face's integration point [Pa/m].
Definition: transfluxmodule.hh:126
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx, const FluidState &exFluidState)
Update the required gradients for boundary faces.
Definition: transfluxmodule.hh:324
const DimMatrix & intrinsicPermeability() const
Return the intrinsic permeability tensor at a face [m^2].
Definition: transfluxmodule.hh:115
unsigned upstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in upstream direction.
Definition: transfluxmodule.hh:171
void updatePolymer(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: transfluxmodule.hh:195
Provides the intensive quantities for the transmissibility based flux module.
Definition: transfluxmodule.hh:80
void update_(const ElementContext &, unsigned, unsigned)
Definition: transfluxmodule.hh:83
Defines the common properties required by the porous medium multi-phase models.
Definition: blackoilboundaryratevector.hh:37
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:235
Specifies a flux module which uses transmissibilities.
Definition: transfluxmodule.hh:56
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: transfluxmodule.hh:63