mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
changed: rename ewoms/models/common -> opm/models/common
This commit is contained in:
Arne Morten Kvarving
parent
e01f712294
commit
17e8fa6574
569
opm/models/common/darcyfluxmodule.hh
Normal file
569
opm/models/common/darcyfluxmodule.hh
Normal file
@ -0,0 +1,569 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief This file contains the necessary classes to calculate the
|
||||||
|
* volumetric fluxes out of a pressure potential gradient using the
|
||||||
|
* Darcy relation.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_DARCY_FLUX_MODULE_HH
|
||||||
|
#define EWOMS_DARCY_FLUX_MODULE_HH
|
||||||
|
|
||||||
|
#include "multiphasebaseproperties.hh"
|
||||||
|
#include <opm/models/common/quantitycallbacks.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
#include <opm/material/common/Exceptions.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
#include <dune/common/fmatrix.hh>
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
NEW_PROP_TAG(MaterialLaw);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyIntensiveQuantities;
|
||||||
|
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyExtensiveQuantities;
|
||||||
|
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyBaseProblem;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Specifies a flux module which uses the Darcy relation.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
struct DarcyFluxModule
|
||||||
|
{
|
||||||
|
typedef DarcyIntensiveQuantities<TypeTag> FluxIntensiveQuantities;
|
||||||
|
typedef DarcyExtensiveQuantities<TypeTag> FluxExtensiveQuantities;
|
||||||
|
typedef DarcyBaseProblem<TypeTag> FluxBaseProblem;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the flux module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{ }
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the defaults for the parameters required by the
|
||||||
|
* Darcy velocity approach.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyBaseProblem
|
||||||
|
{ };
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the intensive quantities for the Darcy flux module
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyIntensiveQuantities
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
protected:
|
||||||
|
void update_(const ElementContext& elemCtx OPM_UNUSED,
|
||||||
|
unsigned dofIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{ }
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the Darcy flux module
|
||||||
|
*
|
||||||
|
* The commonly used Darcy relation looses its validity for Reynolds numbers \f$ Re <
|
||||||
|
* 1\f$. If one encounters flow velocities in porous media above this threshold, the
|
||||||
|
* Forchheimer approach can be used.
|
||||||
|
*
|
||||||
|
* The Darcy equation is given by the following relation:
|
||||||
|
*
|
||||||
|
* \f[
|
||||||
|
\vec{v}_\alpha =
|
||||||
|
\left( \nabla p_\alpha - \rho_\alpha \vec{g}\right)
|
||||||
|
\frac{\mu_\alpha}{k_{r,\alpha} K}
|
||||||
|
\f]
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DarcyExtensiveQuantities
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
||||||
|
|
||||||
|
enum { dimWorld = GridView::dimensionworld };
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
|
||||||
|
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
typedef typename FluidSystem::template ParameterCache<Evaluation> ParameterCache;
|
||||||
|
typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
|
||||||
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
||||||
|
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the intrinsic permeability tensor for a given
|
||||||
|
* sub-control volume face.
|
||||||
|
*/
|
||||||
|
const DimMatrix& intrinsicPermability() const
|
||||||
|
{ return K_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the pressure potential gradient of a fluid phase
|
||||||
|
* at the face's integration point [Pa/m]
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase
|
||||||
|
*/
|
||||||
|
const EvalDimVector& potentialGrad(unsigned phaseIdx) const
|
||||||
|
{ return potentialGrad_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the filter velocity of a fluid phase at the
|
||||||
|
* face's integration point [m/s]
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase
|
||||||
|
*/
|
||||||
|
const EvalDimVector& filterVelocity(unsigned phaseIdx) const
|
||||||
|
{ return filterVelocity_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the volume flux of a fluid phase at the face's integration point
|
||||||
|
* \f$[m^3/s / m^2]\f$
|
||||||
|
*
|
||||||
|
* This is the fluid volume of a phase per second and per square meter of face
|
||||||
|
* area.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase
|
||||||
|
*/
|
||||||
|
const Evaluation& volumeFlux(unsigned phaseIdx) const
|
||||||
|
{ return volumeFlux_[phaseIdx]; }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
short upstreamIndex_(unsigned phaseIdx) const
|
||||||
|
{ return upstreamDofIdx_[phaseIdx]; }
|
||||||
|
|
||||||
|
short downstreamIndex_(unsigned phaseIdx) const
|
||||||
|
{ return downstreamDofIdx_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the gradients which are required to determine the volumetric fluxes
|
||||||
|
*
|
||||||
|
* The the upwind directions is also determined by method.
|
||||||
|
*/
|
||||||
|
void calculateGradients_(const ElementContext& elemCtx,
|
||||||
|
unsigned faceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& gradCalc = elemCtx.gradientCalculator();
|
||||||
|
Opm::PressureCallback<TypeTag> pressureCallback(elemCtx);
|
||||||
|
|
||||||
|
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
|
||||||
|
const auto& faceNormal = scvf.normal();
|
||||||
|
|
||||||
|
unsigned i = scvf.interiorIndex();
|
||||||
|
unsigned j = scvf.exteriorIndex();
|
||||||
|
interiorDofIdx_ = static_cast<short>(i);
|
||||||
|
exteriorDofIdx_ = static_cast<short>(j);
|
||||||
|
unsigned focusDofIdx = elemCtx.focusDofIndex();
|
||||||
|
|
||||||
|
// calculate the "raw" pressure gradient
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
pressureCallback.setPhaseIndex(phaseIdx);
|
||||||
|
gradCalc.calculateGradient(potentialGrad_[phaseIdx],
|
||||||
|
elemCtx,
|
||||||
|
faceIdx,
|
||||||
|
pressureCallback);
|
||||||
|
Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// correct the pressure gradients by the gravitational acceleration
|
||||||
|
if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
|
||||||
|
// estimate the gravitational acceleration at a given SCV face
|
||||||
|
// using the arithmetic mean
|
||||||
|
const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
|
||||||
|
const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
|
||||||
|
|
||||||
|
const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
|
||||||
|
const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
|
||||||
|
|
||||||
|
const auto& posIn = elemCtx.pos(i, timeIdx);
|
||||||
|
const auto& posEx = elemCtx.pos(j, timeIdx);
|
||||||
|
const auto& posFace = scvf.integrationPos();
|
||||||
|
|
||||||
|
// the distance between the centers of the control volumes
|
||||||
|
DimVector distVecIn(posIn);
|
||||||
|
DimVector distVecEx(posEx);
|
||||||
|
DimVector distVecTotal(posEx);
|
||||||
|
|
||||||
|
distVecIn -= posFace;
|
||||||
|
distVecEx -= posFace;
|
||||||
|
distVecTotal -= posIn;
|
||||||
|
Scalar absDistTotalSquared = distVecTotal.two_norm2();
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// calculate the hydrostatic pressure at the integration point of the face
|
||||||
|
Evaluation pStatIn;
|
||||||
|
|
||||||
|
if (std::is_same<Scalar, Evaluation>::value ||
|
||||||
|
interiorDofIdx_ == static_cast<int>(focusDofIdx))
|
||||||
|
{
|
||||||
|
const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
|
||||||
|
pStatIn = - rhoIn*(gIn*distVecIn);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
|
||||||
|
pStatIn = - rhoIn*(gIn*distVecIn);
|
||||||
|
}
|
||||||
|
|
||||||
|
// the quantities on the exterior side of the face do not influence the
|
||||||
|
// result for the TPFA scheme, so they can be treated as scalar values.
|
||||||
|
Evaluation pStatEx;
|
||||||
|
|
||||||
|
if (std::is_same<Scalar, Evaluation>::value ||
|
||||||
|
exteriorDofIdx_ == static_cast<int>(focusDofIdx))
|
||||||
|
{
|
||||||
|
const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
|
||||||
|
pStatEx = - rhoEx*(gEx*distVecEx);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
|
||||||
|
pStatEx = - rhoEx*(gEx*distVecEx);
|
||||||
|
}
|
||||||
|
|
||||||
|
// compute the hydrostatic gradient between the two control volumes (this
|
||||||
|
// gradient exhibitis the same direction as the vector between the two
|
||||||
|
// control volume centers and the length (pStaticExterior -
|
||||||
|
// pStaticInterior)/distanceInteriorToExterior
|
||||||
|
Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
|
||||||
|
f *= (pStatEx - pStatIn)/absDistTotalSquared;
|
||||||
|
|
||||||
|
// calculate the final potential gradient
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
|
||||||
|
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
|
||||||
|
if (!Opm::isfinite(potentialGrad_[phaseIdx][dimIdx])) {
|
||||||
|
throw Opm::NumericalIssue("Non-finite potential gradient for phase '"
|
||||||
|
+std::string(FluidSystem::phaseName(phaseIdx))+"'");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
Opm::Valgrind::SetUndefined(K_);
|
||||||
|
elemCtx.problem().intersectionIntrinsicPermeability(K_, elemCtx, faceIdx, timeIdx);
|
||||||
|
Opm::Valgrind::CheckDefined(K_);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
// determine the upstream and downstream DOFs
|
||||||
|
Evaluation tmp = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
|
||||||
|
tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
|
||||||
|
|
||||||
|
if (tmp > 0) {
|
||||||
|
upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
|
||||||
|
downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
|
||||||
|
downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
|
||||||
|
}
|
||||||
|
|
||||||
|
// we only carry the derivatives along if the upstream DOF is the one which
|
||||||
|
// we currently focus on
|
||||||
|
const auto& up = elemCtx.intensiveQuantities(upstreamDofIdx_[phaseIdx], timeIdx);
|
||||||
|
if (upstreamDofIdx_[phaseIdx] == static_cast<int>(focusDofIdx))
|
||||||
|
mobility_[phaseIdx] = up.mobility(phaseIdx);
|
||||||
|
else
|
||||||
|
mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the gradients at the grid boundary which are required to
|
||||||
|
* determine the volumetric fluxes
|
||||||
|
*
|
||||||
|
* The the upwind directions is also determined by method.
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
void calculateBoundaryGradients_(const ElementContext& elemCtx,
|
||||||
|
unsigned boundaryFaceIdx,
|
||||||
|
unsigned timeIdx,
|
||||||
|
const FluidState& fluidState)
|
||||||
|
{
|
||||||
|
const auto& gradCalc = elemCtx.gradientCalculator();
|
||||||
|
Opm::BoundaryPressureCallback<TypeTag, FluidState> pressureCallback(elemCtx, fluidState);
|
||||||
|
|
||||||
|
// calculate the pressure gradient
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
pressureCallback.setPhaseIndex(phaseIdx);
|
||||||
|
gradCalc.calculateBoundaryGradient(potentialGrad_[phaseIdx],
|
||||||
|
elemCtx,
|
||||||
|
boundaryFaceIdx,
|
||||||
|
pressureCallback);
|
||||||
|
Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
|
||||||
|
}
|
||||||
|
|
||||||
|
const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
|
||||||
|
auto i = scvf.interiorIndex();
|
||||||
|
interiorDofIdx_ = static_cast<short>(i);
|
||||||
|
exteriorDofIdx_ = -1;
|
||||||
|
int focusDofIdx = elemCtx.focusDofIndex();
|
||||||
|
|
||||||
|
// calculate the intrinsic permeability
|
||||||
|
const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
|
||||||
|
K_ = intQuantsIn.intrinsicPermeability();
|
||||||
|
|
||||||
|
// correct the pressure gradients by the gravitational acceleration
|
||||||
|
if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
|
||||||
|
// estimate the gravitational acceleration at a given SCV face
|
||||||
|
// using the arithmetic mean
|
||||||
|
const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
|
||||||
|
const auto& posIn = elemCtx.pos(i, timeIdx);
|
||||||
|
const auto& posFace = scvf.integrationPos();
|
||||||
|
|
||||||
|
// the distance between the face center and the center of the control volume
|
||||||
|
DimVector distVecIn(posIn);
|
||||||
|
distVecIn -= posFace;
|
||||||
|
Scalar absDist = distVecIn.two_norm();
|
||||||
|
Scalar gTimesDist = gIn*distVecIn;
|
||||||
|
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// calculate the hydrostatic pressure at the integration point of the face
|
||||||
|
Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
|
||||||
|
Evaluation pStatIn = - gTimesDist*rhoIn;
|
||||||
|
|
||||||
|
Opm::Valgrind::CheckDefined(pStatIn);
|
||||||
|
|
||||||
|
// compute the hydrostatic gradient between the two control volumes (this
|
||||||
|
// gradient exhibitis the same direction as the vector between the two
|
||||||
|
// control volume centers and the length (pStaticExterior -
|
||||||
|
// pStaticInterior)/distanceInteriorToExterior. Note that for the
|
||||||
|
// boundary, 'pStaticExterior' is zero as the boundary pressure is
|
||||||
|
// defined on boundary face's integration point...
|
||||||
|
EvalDimVector f(distVecIn);
|
||||||
|
f *= - pStatIn/absDist;
|
||||||
|
|
||||||
|
// calculate the final potential gradient
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
|
||||||
|
|
||||||
|
Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
|
||||||
|
if (!Opm::isfinite(potentialGrad_[phaseIdx][dimIdx])) {
|
||||||
|
throw Opm::NumericalIssue("Non finite potential gradient for phase '"
|
||||||
|
+std::string(FluidSystem::phaseName(phaseIdx))+"'");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// determine the upstream and downstream DOFs
|
||||||
|
const auto& faceNormal = scvf.normal();
|
||||||
|
|
||||||
|
const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx);
|
||||||
|
|
||||||
|
Scalar kr[numPhases];
|
||||||
|
MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
|
||||||
|
Opm::Valgrind::CheckDefined(kr);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
Evaluation tmp = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
|
||||||
|
tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
|
||||||
|
|
||||||
|
if (tmp > 0) {
|
||||||
|
upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
|
||||||
|
downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
|
||||||
|
downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
|
||||||
|
}
|
||||||
|
|
||||||
|
// take the phase mobility from the DOF in upstream direction
|
||||||
|
if (upstreamDofIdx_[phaseIdx] < 0) {
|
||||||
|
if (interiorDofIdx_ == focusDofIdx)
|
||||||
|
mobility_[phaseIdx] =
|
||||||
|
kr[phaseIdx] / fluidState.viscosity(phaseIdx);
|
||||||
|
else
|
||||||
|
mobility_[phaseIdx] =
|
||||||
|
Toolbox::value(kr[phaseIdx])
|
||||||
|
/ Toolbox::value(fluidState.viscosity(phaseIdx));
|
||||||
|
}
|
||||||
|
else if (upstreamDofIdx_[phaseIdx] != focusDofIdx)
|
||||||
|
mobility_[phaseIdx] = Toolbox::value(intQuantsIn.mobility(phaseIdx));
|
||||||
|
else
|
||||||
|
mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
|
||||||
|
Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the volumetric fluxes of all phases
|
||||||
|
*
|
||||||
|
* The pressure potentials and upwind directions must already be
|
||||||
|
* determined before calling this method!
|
||||||
|
*/
|
||||||
|
void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
|
||||||
|
const DimVector& normal = scvf.normal();
|
||||||
|
Opm::Valgrind::CheckDefined(normal);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
filterVelocity_[phaseIdx] = 0.0;
|
||||||
|
volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
asImp_().calculateFilterVelocity_(phaseIdx);
|
||||||
|
Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
|
||||||
|
|
||||||
|
volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
for (unsigned i = 0; i < normal.size(); ++i)
|
||||||
|
volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the volumetric fluxes at a boundary face of all fluid phases
|
||||||
|
*
|
||||||
|
* The pressure potentials and upwind directions must already be determined before
|
||||||
|
* calling this method!
|
||||||
|
*/
|
||||||
|
void calculateBoundaryFluxes_(const ElementContext& elemCtx,
|
||||||
|
unsigned boundaryFaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
|
||||||
|
const DimVector& normal = scvf.normal();
|
||||||
|
Opm::Valgrind::CheckDefined(normal);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
filterVelocity_[phaseIdx] = 0.0;
|
||||||
|
volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
asImp_().calculateFilterVelocity_(phaseIdx);
|
||||||
|
Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
|
||||||
|
volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
for (unsigned i = 0; i < normal.size(); ++i)
|
||||||
|
volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void calculateFilterVelocity_(unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
#ifndef NDEBUG
|
||||||
|
assert(Opm::isfinite(mobility_[phaseIdx]));
|
||||||
|
for (unsigned i = 0; i < K_.M(); ++ i)
|
||||||
|
for (unsigned j = 0; j < K_.N(); ++ j)
|
||||||
|
assert(std::isfinite(K_[i][j]));
|
||||||
|
#endif
|
||||||
|
|
||||||
|
K_.mv(potentialGrad_[phaseIdx], filterVelocity_[phaseIdx]);
|
||||||
|
filterVelocity_[phaseIdx] *= - mobility_[phaseIdx];
|
||||||
|
|
||||||
|
#ifndef NDEBUG
|
||||||
|
for (unsigned i = 0; i < filterVelocity_[phaseIdx].size(); ++ i)
|
||||||
|
assert(Opm::isfinite(filterVelocity_[phaseIdx][i]));
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
Implementation& asImp_()
|
||||||
|
{ return *static_cast<Implementation*>(this); }
|
||||||
|
|
||||||
|
const Implementation& asImp_() const
|
||||||
|
{ return *static_cast<const Implementation*>(this); }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
// intrinsic permeability tensor and its square root
|
||||||
|
DimMatrix K_;
|
||||||
|
|
||||||
|
// mobilities of all fluid phases [1 / (Pa s)]
|
||||||
|
Evaluation mobility_[numPhases];
|
||||||
|
|
||||||
|
// filter velocities of all phases [m/s]
|
||||||
|
EvalDimVector filterVelocity_[numPhases];
|
||||||
|
|
||||||
|
// the volumetric flux of all fluid phases over the control
|
||||||
|
// volume's face [m^3/s / m^2]
|
||||||
|
Evaluation volumeFlux_[numPhases];
|
||||||
|
|
||||||
|
// pressure potential gradients of all phases [Pa / m]
|
||||||
|
EvalDimVector potentialGrad_[numPhases];
|
||||||
|
|
||||||
|
// upstream, downstream, interior and exterior DOFs
|
||||||
|
short upstreamDofIdx_[numPhases];
|
||||||
|
short downstreamDofIdx_[numPhases];
|
||||||
|
short interiorDofIdx_;
|
||||||
|
short exteriorDofIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
496
opm/models/common/diffusionmodule.hh
Normal file
496
opm/models/common/diffusionmodule.hh
Normal file
@ -0,0 +1,496 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief Classes required for molecular diffusion.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_DIFFUSION_MODULE_HH
|
||||||
|
#define EWOMS_DIFFUSION_MODULE_HH
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
#include <opm/models/common/quantitycallbacks.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
#include <opm/material/common/Exceptions.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
NEW_PROP_TAG(Indices);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Diffusion
|
||||||
|
* \class Opm::DiffusionModule
|
||||||
|
* \brief Provides the auxiliary methods required for consideration of the
|
||||||
|
* diffusion equation.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableDiffusion>
|
||||||
|
class DiffusionModule;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionModule
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionModule<TypeTag, /*enableDiffusion=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the diffusion module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Adds the diffusive mass flux flux to the flux vector over a flux
|
||||||
|
* integration point.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addDiffusiveFlux(RateVector& flux OPM_UNUSED,
|
||||||
|
const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionModule
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionModule<TypeTag, /*enableDiffusion=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||||
|
|
||||||
|
enum { numPhases = FluidSystem::numPhases };
|
||||||
|
enum { numComponents = FluidSystem::numComponents };
|
||||||
|
enum { conti0EqIdx = Indices::conti0EqIdx };
|
||||||
|
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the diffusion module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Adds the mass flux due to molecular diffusion to the flux vector over the
|
||||||
|
* flux integration point.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addDiffusiveFlux(RateVector& flux, const Context& context,
|
||||||
|
unsigned spaceIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
|
||||||
|
|
||||||
|
const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
|
||||||
|
const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
|
||||||
|
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
// arithmetic mean of the phase's molar density
|
||||||
|
Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
|
||||||
|
rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
|
||||||
|
rhoMolar /= 2;
|
||||||
|
|
||||||
|
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||||
|
// mass flux due to molecular diffusion
|
||||||
|
flux[conti0EqIdx + compIdx] +=
|
||||||
|
-rhoMolar
|
||||||
|
* extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
|
||||||
|
* extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Diffusion
|
||||||
|
* \class Opm::DiffusionIntensiveQuantities
|
||||||
|
*
|
||||||
|
* \brief Provides the volumetric quantities required for the
|
||||||
|
* calculation of molecular diffusive fluxes.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableDiffusion>
|
||||||
|
class DiffusionIntensiveQuantities;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionIntensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the tortuousity of the sub-domain of a fluid
|
||||||
|
* phase in the porous medium.
|
||||||
|
*/
|
||||||
|
Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Method tortuosity() does not make sense "
|
||||||
|
"if diffusion is disabled");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the molecular diffusion coefficient for a
|
||||||
|
* component in a phase.
|
||||||
|
*/
|
||||||
|
Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Method diffusionCoefficient() does not "
|
||||||
|
"make sense if diffusion is disabled");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the effective molecular diffusion coefficient of
|
||||||
|
* the porous medium for a component in a phase.
|
||||||
|
*/
|
||||||
|
Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Method effectiveDiffusionCoefficient() "
|
||||||
|
"does not make sense if diffusion is disabled");
|
||||||
|
}
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate diffusive
|
||||||
|
* mass fluxes.
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
void update_(FluidState& fs OPM_UNUSED,
|
||||||
|
typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache OPM_UNUSED,
|
||||||
|
const ElementContext& elemCtx OPM_UNUSED,
|
||||||
|
unsigned dofIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{ }
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionIntensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
|
||||||
|
enum { numPhases = FluidSystem::numPhases };
|
||||||
|
enum { numComponents = FluidSystem::numComponents };
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the molecular diffusion coefficient for a
|
||||||
|
* component in a phase.
|
||||||
|
*/
|
||||||
|
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
|
||||||
|
{ return diffusionCoefficient_[phaseIdx][compIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the tortuousity of the sub-domain of a fluid
|
||||||
|
* phase in the porous medium.
|
||||||
|
*/
|
||||||
|
Evaluation tortuosity(unsigned phaseIdx) const
|
||||||
|
{ return tortuosity_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the effective molecular diffusion coefficient of
|
||||||
|
* the porous medium for a component in a phase.
|
||||||
|
*/
|
||||||
|
Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
|
||||||
|
{ return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate diffusive
|
||||||
|
* mass fluxes.
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
void update_(FluidState& fluidState,
|
||||||
|
typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
|
||||||
|
const ElementContext& elemCtx,
|
||||||
|
unsigned dofIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// TODO: let the problem do this (this is a constitutive
|
||||||
|
// relation of which the model should be free of from the
|
||||||
|
// abstraction POV!)
|
||||||
|
const Evaluation& base =
|
||||||
|
Toolbox::max(0.0001,
|
||||||
|
intQuants.porosity()
|
||||||
|
* intQuants.fluidState().saturation(phaseIdx));
|
||||||
|
tortuosity_[phaseIdx] =
|
||||||
|
1.0 / (intQuants.porosity() * intQuants.porosity())
|
||||||
|
* Toolbox::pow(base, 7.0/3.0);
|
||||||
|
|
||||||
|
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||||
|
diffusionCoefficient_[phaseIdx][compIdx] =
|
||||||
|
FluidSystem::diffusionCoefficient(fluidState,
|
||||||
|
paramCache,
|
||||||
|
phaseIdx,
|
||||||
|
compIdx);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
Evaluation tortuosity_[numPhases];
|
||||||
|
Evaluation diffusionCoefficient_[numPhases][numComponents];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Diffusion
|
||||||
|
* \class Opm::DiffusionExtensiveQuantities
|
||||||
|
*
|
||||||
|
* \brief Provides the quantities required to calculate diffusive mass fluxes.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableDiffusion>
|
||||||
|
class DiffusionExtensiveQuantities;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionExtensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* the diffusive mass fluxes.
|
||||||
|
*/
|
||||||
|
void update_(const ElementContext& elemCtx OPM_UNUSED,
|
||||||
|
unsigned faceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
template <class Context, class FluidState>
|
||||||
|
void updateBoundary_(const Context& context OPM_UNUSED,
|
||||||
|
unsigned bfIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED,
|
||||||
|
const FluidState& fluidState OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief The the gradient of the mole fraction times the face normal.
|
||||||
|
*
|
||||||
|
* \copydoc Doxygen::phaseIdxParam
|
||||||
|
* \copydoc Doxygen::compIdxParam
|
||||||
|
*/
|
||||||
|
const Evaluation& moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED,
|
||||||
|
unsigned compIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("The method moleFractionGradient() does not "
|
||||||
|
"make sense if diffusion is disabled.");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief The effective diffusion coeffcient of a component in a
|
||||||
|
* fluid phase at the face's integration point
|
||||||
|
*
|
||||||
|
* \copydoc Doxygen::phaseIdxParam
|
||||||
|
* \copydoc Doxygen::compIdxParam
|
||||||
|
*/
|
||||||
|
const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED,
|
||||||
|
unsigned compIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("The method effectiveDiffusionCoefficient() "
|
||||||
|
"does not make sense if diffusion is disabled.");
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::DiffusionExtensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
|
||||||
|
enum { dimWorld = GridView::dimensionworld };
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||||
|
|
||||||
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
||||||
|
typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* the diffusive mass fluxes.
|
||||||
|
*/
|
||||||
|
void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& gradCalc = elemCtx.gradientCalculator();
|
||||||
|
Opm::MoleFractionCallback<TypeTag> moleFractionCallback(elemCtx);
|
||||||
|
|
||||||
|
const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
|
||||||
|
const auto& normal = face.normal();
|
||||||
|
const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
|
||||||
|
|
||||||
|
const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
|
||||||
|
const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
moleFractionCallback.setPhaseIndex(phaseIdx);
|
||||||
|
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||||
|
moleFractionCallback.setComponentIndex(compIdx);
|
||||||
|
|
||||||
|
DimEvalVector moleFractionGradient(0.0);
|
||||||
|
gradCalc.calculateGradient(moleFractionGradient,
|
||||||
|
elemCtx,
|
||||||
|
faceIdx,
|
||||||
|
moleFractionCallback);
|
||||||
|
|
||||||
|
moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0;
|
||||||
|
for (unsigned i = 0; i < normal.size(); ++i)
|
||||||
|
moleFractionGradientNormal_[phaseIdx][compIdx] +=
|
||||||
|
normal[i]*moleFractionGradient[i];
|
||||||
|
Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
|
||||||
|
|
||||||
|
// use the arithmetic average for the effective
|
||||||
|
// diffusion coefficients.
|
||||||
|
effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
|
||||||
|
(intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
|
||||||
|
+
|
||||||
|
intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
|
||||||
|
/ 2;
|
||||||
|
|
||||||
|
Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class Context, class FluidState>
|
||||||
|
void updateBoundary_(const Context& context,
|
||||||
|
unsigned bfIdx,
|
||||||
|
unsigned timeIdx,
|
||||||
|
const FluidState& fluidState)
|
||||||
|
{
|
||||||
|
const auto& stencil = context.stencil(timeIdx);
|
||||||
|
const auto& face = stencil.boundaryFace(bfIdx);
|
||||||
|
|
||||||
|
const auto& elemCtx = context.elementContext();
|
||||||
|
unsigned insideScvIdx = face.interiorIndex();
|
||||||
|
const auto& insideScv = stencil.subControlVolume(insideScvIdx);
|
||||||
|
|
||||||
|
const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
|
||||||
|
const auto& fluidStateInside = intQuantsInside.fluidState();
|
||||||
|
|
||||||
|
// distance between the center of the SCV and center of the boundary face
|
||||||
|
DimVector distVec = face.integrationPos();
|
||||||
|
distVec -= context.element().geometry().global(insideScv.localGeometry().center());
|
||||||
|
|
||||||
|
Scalar dist = distVec * face.normal();
|
||||||
|
|
||||||
|
// if the following assertation triggers, the center of the
|
||||||
|
// center of the interior SCV was not inside the element!
|
||||||
|
assert(dist > 0);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||||
|
// calculate mole fraction gradient using two-point
|
||||||
|
// gradients
|
||||||
|
moleFractionGradientNormal_[phaseIdx][compIdx] =
|
||||||
|
(fluidState.moleFraction(phaseIdx, compIdx)
|
||||||
|
-
|
||||||
|
fluidStateInside.moleFraction(phaseIdx, compIdx))
|
||||||
|
/ dist;
|
||||||
|
Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
|
||||||
|
|
||||||
|
// use effective diffusion coefficients of the interior finite
|
||||||
|
// volume.
|
||||||
|
effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
|
||||||
|
intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
|
||||||
|
|
||||||
|
Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief The the gradient of the mole fraction times the face normal.
|
||||||
|
*
|
||||||
|
* \copydoc Doxygen::phaseIdxParam
|
||||||
|
* \copydoc Doxygen::compIdxParam
|
||||||
|
*/
|
||||||
|
const Evaluation& moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
|
||||||
|
{ return moleFractionGradientNormal_[phaseIdx][compIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief The effective diffusion coeffcient of a component in a
|
||||||
|
* fluid phase at the face's integration point
|
||||||
|
*
|
||||||
|
* \copydoc Doxygen::phaseIdxParam
|
||||||
|
* \copydoc Doxygen::compIdxParam
|
||||||
|
*/
|
||||||
|
const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
|
||||||
|
{ return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
|
||||||
|
|
||||||
|
private:
|
||||||
|
Evaluation moleFractionGradientNormal_[numPhases][numComponents];
|
||||||
|
Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents];
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
859
opm/models/common/energymodule.hh
Normal file
859
opm/models/common/energymodule.hh
Normal file
@ -0,0 +1,859 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief Contains the classes required to consider energy as a
|
||||||
|
* conservation quantity in a multi-phase module.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_ENERGY_MODULE_HH
|
||||||
|
#define EWOMS_ENERGY_MODULE_HH
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
#include <opm/models/common/quantitycallbacks.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
#include <opm/material/common/Exceptions.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
|
||||||
|
#include <string>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
NEW_PROP_TAG(Indices);
|
||||||
|
NEW_PROP_TAG(EnableEnergy);
|
||||||
|
NEW_PROP_TAG(ThermalConductionLaw);
|
||||||
|
NEW_PROP_TAG(ThermalConductionLawParams);
|
||||||
|
NEW_PROP_TAG(SolidEnergyLaw);
|
||||||
|
NEW_PROP_TAG(SolidEnergyLawParams);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
/*!
|
||||||
|
* \ingroup Energy
|
||||||
|
* \brief Provides the auxiliary methods required for consideration of
|
||||||
|
* the energy equation.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableEnergy>
|
||||||
|
class EnergyModule;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyModule
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyModule<TypeTag, /*enableEnergy=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
|
||||||
|
|
||||||
|
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||||
|
|
||||||
|
typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the energy module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the name of a primary variable or an empty
|
||||||
|
* string if the specified primary variable index does not belong to
|
||||||
|
* the energy module.
|
||||||
|
*/
|
||||||
|
static std::string primaryVarName(unsigned pvIdx OPM_UNUSED)
|
||||||
|
{ return ""; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the name of an equation or an empty
|
||||||
|
* string if the specified equation index does not belong to
|
||||||
|
* the energy module.
|
||||||
|
*/
|
||||||
|
static std::string eqName(unsigned eqIdx OPM_UNUSED)
|
||||||
|
{ return ""; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the relative weight of a primary variable for
|
||||||
|
* calculating relative errors.
|
||||||
|
*/
|
||||||
|
static Scalar primaryVarWeight(const Model& model OPM_UNUSED,
|
||||||
|
unsigned globalDofIdx OPM_UNUSED,
|
||||||
|
unsigned pvIdx OPM_UNUSED)
|
||||||
|
{ return -1; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the relative weight of a equation of the residual.
|
||||||
|
*/
|
||||||
|
static Scalar eqWeight(const Model& model OPM_UNUSED,
|
||||||
|
unsigned globalDofIdx OPM_UNUSED,
|
||||||
|
unsigned eqIdx OPM_UNUSED)
|
||||||
|
{ return -1; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Given a fluid state, set the temperature in the primary variables
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
static void setPriVarTemperatures(PrimaryVariables& priVars OPM_UNUSED,
|
||||||
|
const FluidState& fs OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Given a fluid state, set the enthalpy rate which emerges
|
||||||
|
* from a volumetric rate.
|
||||||
|
*/
|
||||||
|
template <class RateVector, class FluidState>
|
||||||
|
static void setEnthalpyRate(RateVector& rateVec OPM_UNUSED,
|
||||||
|
const FluidState& fluidState OPM_UNUSED,
|
||||||
|
unsigned phaseIdx OPM_UNUSED,
|
||||||
|
const Evaluation& volume OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the rate of the enthalpy flux to a rate vector.
|
||||||
|
*/
|
||||||
|
static void setEnthalpyRate(RateVector& rateVec OPM_UNUSED,
|
||||||
|
const Evaluation& rate OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the rate of the enthalpy flux to a rate vector.
|
||||||
|
*/
|
||||||
|
static void addToEnthalpyRate(RateVector& rateVec OPM_UNUSED,
|
||||||
|
const Evaluation& rate OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the rate of the conductive energy flux to a rate vector.
|
||||||
|
*/
|
||||||
|
static Scalar thermalConductionRate(const ExtensiveQuantities& extQuants OPM_UNUSED)
|
||||||
|
{ return 0.0; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for a fluid phase to an equation
|
||||||
|
* vector
|
||||||
|
*/
|
||||||
|
template <class LhsEval>
|
||||||
|
static void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage OPM_UNUSED,
|
||||||
|
const IntensiveQuantities& intQuants OPM_UNUSED,
|
||||||
|
unsigned phaseIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for a fluid phase to an equation
|
||||||
|
* vector
|
||||||
|
*/
|
||||||
|
template <class LhsEval, class Scv>
|
||||||
|
static void addFracturePhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage OPM_UNUSED,
|
||||||
|
const IntensiveQuantities& intQuants OPM_UNUSED,
|
||||||
|
const Scv& scv OPM_UNUSED,
|
||||||
|
unsigned phaseIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for the fracture part a fluid phase to an
|
||||||
|
* equation vector
|
||||||
|
*/
|
||||||
|
template <class LhsEval>
|
||||||
|
static void addSolidEnergyStorage(Dune::FieldVector<LhsEval, numEq>& storage OPM_UNUSED,
|
||||||
|
const IntensiveQuantities& intQuants OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Evaluates the advective energy fluxes over a face of a
|
||||||
|
* subcontrol volume and adds the result in the flux vector.
|
||||||
|
*
|
||||||
|
* This method is called by compute flux (base class)
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addAdvectiveFlux(RateVector& flux OPM_UNUSED,
|
||||||
|
const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Evaluates the advective energy fluxes over a fracture
|
||||||
|
* which should be attributed to a face of a subcontrol
|
||||||
|
* volume and adds the result in the flux vector.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void handleFractureFlux(RateVector& flux OPM_UNUSED,
|
||||||
|
const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Adds the diffusive energy flux to the flux vector over the face of a
|
||||||
|
* sub-control volume.
|
||||||
|
*
|
||||||
|
* This method is called by compute flux (base class)
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addDiffusiveFlux(RateVector& flux OPM_UNUSED,
|
||||||
|
const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyModule
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyModule<TypeTag, /*enableEnergy=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
|
||||||
|
|
||||||
|
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||||
|
enum { numPhases = FluidSystem::numPhases };
|
||||||
|
enum { energyEqIdx = Indices::energyEqIdx };
|
||||||
|
enum { temperatureIdx = Indices::temperatureIdx };
|
||||||
|
|
||||||
|
typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the energy module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the name of a primary variable or an empty
|
||||||
|
* string if the specified primary variable index does not belong to
|
||||||
|
* the energy module.
|
||||||
|
*/
|
||||||
|
static std::string primaryVarName(unsigned pvIdx)
|
||||||
|
{
|
||||||
|
if (pvIdx == temperatureIdx)
|
||||||
|
return "temperature";
|
||||||
|
return "";
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the name of an equation or an empty
|
||||||
|
* string if the specified equation index does not belong to
|
||||||
|
* the energy module.
|
||||||
|
*/
|
||||||
|
static std::string eqName(unsigned eqIdx)
|
||||||
|
{
|
||||||
|
if (eqIdx == energyEqIdx)
|
||||||
|
return "energy";
|
||||||
|
return "";
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the relative weight of a primary variable for
|
||||||
|
* calculating relative errors.
|
||||||
|
*/
|
||||||
|
static Scalar primaryVarWeight(const Model& model, unsigned globalDofIdx, unsigned pvIdx)
|
||||||
|
{
|
||||||
|
if (pvIdx != temperatureIdx)
|
||||||
|
return -1;
|
||||||
|
|
||||||
|
// make the weight of the temperature primary variable inversly proportional to its value
|
||||||
|
return std::max(1.0/1000, 1.0/model.solution(/*timeIdx=*/0)[globalDofIdx][temperatureIdx]);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the relative weight of a equation.
|
||||||
|
*/
|
||||||
|
static Scalar eqWeight(const Model& model OPM_UNUSED,
|
||||||
|
unsigned globalDofIdx OPM_UNUSED,
|
||||||
|
unsigned eqIdx)
|
||||||
|
{
|
||||||
|
if (eqIdx != energyEqIdx)
|
||||||
|
return -1;
|
||||||
|
|
||||||
|
// approximate change of internal energy of 1kg of liquid water for a temperature
|
||||||
|
// change of 30K
|
||||||
|
return 1.0 / (4.184e3 * 30.0);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the rate of energy flux of a rate vector.
|
||||||
|
*/
|
||||||
|
static void setEnthalpyRate(RateVector& rateVec, const Evaluation& rate)
|
||||||
|
{ rateVec[energyEqIdx] = rate; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the rate of the enthalpy flux to a rate vector.
|
||||||
|
*/
|
||||||
|
static void addToEnthalpyRate(RateVector& rateVec, const Evaluation& rate)
|
||||||
|
{ rateVec[energyEqIdx] += rate; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the conductive energy flux for a given flux integration point.
|
||||||
|
*/
|
||||||
|
static Evaluation thermalConductionRate(const ExtensiveQuantities& extQuants)
|
||||||
|
{ return -extQuants.temperatureGradNormal() * extQuants.thermalConductivity(); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Given a fluid state, set the enthalpy rate which emerges
|
||||||
|
* from a volumetric rate.
|
||||||
|
*/
|
||||||
|
template <class RateVector, class FluidState>
|
||||||
|
static void setEnthalpyRate(RateVector& rateVec,
|
||||||
|
const FluidState& fluidState,
|
||||||
|
unsigned phaseIdx,
|
||||||
|
const Evaluation& volume)
|
||||||
|
{
|
||||||
|
rateVec[energyEqIdx] =
|
||||||
|
volume
|
||||||
|
* fluidState.density(phaseIdx)
|
||||||
|
* fluidState.enthalpy(phaseIdx);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Given a fluid state, set the temperature in the primary variables
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
static void setPriVarTemperatures(PrimaryVariables& priVars,
|
||||||
|
const FluidState& fs)
|
||||||
|
{
|
||||||
|
priVars[temperatureIdx] = Toolbox::value(fs.temperature(/*phaseIdx=*/0));
|
||||||
|
#ifndef NDEBUG
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
assert(std::abs(Toolbox::value(fs.temperature(/*phaseIdx=*/0))
|
||||||
|
- Toolbox::value(fs.temperature(phaseIdx))) < 1e-30);
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for a fluid phase to an equation
|
||||||
|
* vector
|
||||||
|
*/
|
||||||
|
template <class LhsEval>
|
||||||
|
static void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
|
||||||
|
const IntensiveQuantities& intQuants,
|
||||||
|
unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
const auto& fs = intQuants.fluidState();
|
||||||
|
storage[energyEqIdx] +=
|
||||||
|
Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(fs.internalEnergy(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(intQuants.porosity());
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for a fluid phase to an equation
|
||||||
|
* vector
|
||||||
|
*/
|
||||||
|
template <class Scv, class LhsEval>
|
||||||
|
static void addFracturePhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
|
||||||
|
const IntensiveQuantities& intQuants,
|
||||||
|
const Scv& scv,
|
||||||
|
unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
const auto& fs = intQuants.fractureFluidState();
|
||||||
|
storage[energyEqIdx] +=
|
||||||
|
Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(fs.internalEnergy(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
|
||||||
|
* Toolbox::template decay<LhsEval>(intQuants.fracturePorosity())
|
||||||
|
* Toolbox::template decay<LhsEval>(intQuants.fractureVolume())/scv.volume();
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Add the energy storage term for a fluid phase to an equation
|
||||||
|
* vector
|
||||||
|
*/
|
||||||
|
template <class LhsEval>
|
||||||
|
static void addSolidEnergyStorage(Dune::FieldVector<LhsEval, numEq>& storage,
|
||||||
|
const IntensiveQuantities& intQuants)
|
||||||
|
{ storage[energyEqIdx] += Opm::decay<LhsEval>(intQuants.solidInternalEnergy()); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Evaluates the advective energy fluxes for a flux integration point and adds
|
||||||
|
* the result in the flux vector.
|
||||||
|
*
|
||||||
|
* This method is called by compute flux (base class)
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addAdvectiveFlux(RateVector& flux,
|
||||||
|
const Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
|
||||||
|
|
||||||
|
// advective energy flux in all phases
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!context.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// intensive quantities of the upstream and the downstream DOFs
|
||||||
|
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
|
||||||
|
const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx);
|
||||||
|
|
||||||
|
flux[energyEqIdx] +=
|
||||||
|
extQuants.volumeFlux(phaseIdx)
|
||||||
|
* up.fluidState().enthalpy(phaseIdx)
|
||||||
|
* up.fluidState().density(phaseIdx);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Evaluates the advective energy fluxes over a fracture which should be
|
||||||
|
* attributed to a face of a subcontrol volume and adds the result in the flux
|
||||||
|
* vector.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void handleFractureFlux(RateVector& flux,
|
||||||
|
const Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& scvf = context.stencil(timeIdx).interiorFace(spaceIdx);
|
||||||
|
const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
|
||||||
|
|
||||||
|
// reduce the energy flux in the matrix by the half the width occupied by the
|
||||||
|
// fracture
|
||||||
|
flux[energyEqIdx] *=
|
||||||
|
1 - extQuants.fractureWidth()/(2*scvf.area());
|
||||||
|
|
||||||
|
// advective energy flux in all phases
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!context.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// intensive quantities of the upstream and the downstream DOFs
|
||||||
|
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
|
||||||
|
const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx);
|
||||||
|
|
||||||
|
flux[energyEqIdx] +=
|
||||||
|
extQuants.fractureVolumeFlux(phaseIdx)
|
||||||
|
* up.fluidState().enthalpy(phaseIdx)
|
||||||
|
* up.fluidState().density(phaseIdx);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Adds the diffusive energy flux to the flux vector over the face of a
|
||||||
|
* sub-control volume.
|
||||||
|
*
|
||||||
|
* This method is called by compute flux (base class)
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
static void addDiffusiveFlux(RateVector& flux,
|
||||||
|
const Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
|
||||||
|
|
||||||
|
// diffusive energy flux
|
||||||
|
flux[energyEqIdx] +=
|
||||||
|
- extQuants.temperatureGradNormal()
|
||||||
|
* extQuants.thermalConductivity();
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Energy
|
||||||
|
* \class Opm::EnergyIndices
|
||||||
|
*
|
||||||
|
* \brief Provides the indices required for the energy equation.
|
||||||
|
*/
|
||||||
|
template <unsigned PVOffset, bool enableEnergy>
|
||||||
|
struct EnergyIndices;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyIndices
|
||||||
|
*/
|
||||||
|
template <unsigned PVOffset>
|
||||||
|
struct EnergyIndices<PVOffset, /*enableEnergy=*/false>
|
||||||
|
{
|
||||||
|
//! The index of the primary variable representing temperature
|
||||||
|
enum { temperatureIdx = -1000 };
|
||||||
|
|
||||||
|
//! The index of the equation representing the conservation of energy
|
||||||
|
enum { energyEqIdx = -1000 };
|
||||||
|
|
||||||
|
protected:
|
||||||
|
enum { numEq_ = 0 };
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyIndices
|
||||||
|
*/
|
||||||
|
template <unsigned PVOffset>
|
||||||
|
struct EnergyIndices<PVOffset, /*enableEnergy=*/true>
|
||||||
|
{
|
||||||
|
//! The index of the primary variable representing temperature
|
||||||
|
enum { temperatureIdx = PVOffset };
|
||||||
|
|
||||||
|
//! The index of the equation representing the conservation of energy
|
||||||
|
enum { energyEqIdx = PVOffset };
|
||||||
|
|
||||||
|
protected:
|
||||||
|
enum { numEq_ = 1 };
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Energy
|
||||||
|
* \class Opm::EnergyIntensiveQuantities
|
||||||
|
*
|
||||||
|
* \brief Provides the volumetric quantities required for the energy equation.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableEnergy>
|
||||||
|
class EnergyIntensiveQuantities;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyIntensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyIntensiveQuantities<TypeTag, /*enableEnergy=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the volumetric internal energy \f$\mathrm{[J/(m^3]}\f$ of the
|
||||||
|
* solid matrix in the sub-control volume.
|
||||||
|
*/
|
||||||
|
Evaluation solidInternalEnergy() const
|
||||||
|
{
|
||||||
|
throw std::logic_error("solidInternalEnergy() does not make sense for isothermal models");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the total thermal conductivity \f$\mathrm{[W/m^2 / (K/m)]}\f$ of
|
||||||
|
* the solid matrix in the sub-control volume.
|
||||||
|
*/
|
||||||
|
Evaluation thermalConductivity() const
|
||||||
|
{
|
||||||
|
throw std::logic_error("thermalConductivity() does not make sense for isothermal models");
|
||||||
|
}
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the temperatures of the fluids of a fluid state.
|
||||||
|
*/
|
||||||
|
template <class FluidState, class Context>
|
||||||
|
static void updateTemperatures_(FluidState& fluidState,
|
||||||
|
const Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
Scalar T = context.problem().temperature(context, spaceIdx, timeIdx);
|
||||||
|
fluidState.setTemperature(Toolbox::createConstant(T));
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* energy fluxes.
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
void update_(FluidState& fs OPM_UNUSED,
|
||||||
|
typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache OPM_UNUSED,
|
||||||
|
const ElementContext& elemCtx OPM_UNUSED,
|
||||||
|
unsigned dofIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{ }
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyIntensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyIntensiveQuantities<TypeTag, /*enableEnergy=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ThermalConductionLaw) ThermalConductionLaw;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, SolidEnergyLaw) SolidEnergyLaw;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||||
|
|
||||||
|
enum { numPhases = FluidSystem::numPhases };
|
||||||
|
enum { energyEqIdx = Indices::energyEqIdx };
|
||||||
|
enum { temperatureIdx = Indices::temperatureIdx };
|
||||||
|
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the temperatures of the fluids of a fluid state.
|
||||||
|
*/
|
||||||
|
template <class FluidState, class Context>
|
||||||
|
static void updateTemperatures_(FluidState& fluidState,
|
||||||
|
const Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& priVars = context.primaryVars(spaceIdx, timeIdx);
|
||||||
|
Evaluation val;
|
||||||
|
if (std::is_same<Evaluation, Scalar>::value) // finite differences
|
||||||
|
val = Toolbox::createConstant(priVars[temperatureIdx]);
|
||||||
|
else {
|
||||||
|
// automatic differentiation
|
||||||
|
if (timeIdx == 0)
|
||||||
|
val = Toolbox::createVariable(priVars[temperatureIdx], temperatureIdx);
|
||||||
|
else
|
||||||
|
val = Toolbox::createConstant(priVars[temperatureIdx]);
|
||||||
|
}
|
||||||
|
fluidState.setTemperature(val);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* energy fluxes.
|
||||||
|
*/
|
||||||
|
template <class FluidState>
|
||||||
|
void update_(FluidState& fs,
|
||||||
|
typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
|
||||||
|
const ElementContext& elemCtx,
|
||||||
|
unsigned dofIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
// set the specific enthalpies of the fluids
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
fs.setEnthalpy(phaseIdx,
|
||||||
|
FluidSystem::enthalpy(fs, paramCache, phaseIdx));
|
||||||
|
}
|
||||||
|
|
||||||
|
// compute and set the volumetric internal energy of the solid phase
|
||||||
|
const auto& problem = elemCtx.problem();
|
||||||
|
const auto& solidEnergyParams = problem.solidEnergyLawParams(elemCtx, dofIdx, timeIdx);
|
||||||
|
const auto& thermalCondParams = problem.thermalConductionLawParams(elemCtx, dofIdx, timeIdx);
|
||||||
|
|
||||||
|
solidInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyParams, fs);
|
||||||
|
thermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalCondParams, fs);
|
||||||
|
|
||||||
|
Opm::Valgrind::CheckDefined(solidInternalEnergy_);
|
||||||
|
Opm::Valgrind::CheckDefined(thermalConductivity_);
|
||||||
|
}
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the volumetric internal energy \f$\mathrm{[J/m^3]}\f$ of the
|
||||||
|
* solid matrix in the sub-control volume.
|
||||||
|
*/
|
||||||
|
const Evaluation& solidInternalEnergy() const
|
||||||
|
{ return solidInternalEnergy_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the total conductivity capacity \f$\mathrm{[W/m^2 / (K/m)]}\f$ of
|
||||||
|
* the solid matrix in the sub-control volume.
|
||||||
|
*/
|
||||||
|
const Evaluation& thermalConductivity() const
|
||||||
|
{ return thermalConductivity_; }
|
||||||
|
|
||||||
|
private:
|
||||||
|
Evaluation solidInternalEnergy_;
|
||||||
|
Evaluation thermalConductivity_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Energy
|
||||||
|
* \class Opm::EnergyExtensiveQuantities
|
||||||
|
*
|
||||||
|
* \brief Provides the quantities required to calculate energy fluxes.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, bool enableEnergy>
|
||||||
|
class EnergyExtensiveQuantities;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyExtensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyExtensiveQuantities<TypeTag, /*enableEnergy=*/false>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* energy fluxes.
|
||||||
|
*/
|
||||||
|
void update_(const ElementContext& elemCtx OPM_UNUSED,
|
||||||
|
unsigned faceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
template <class Context, class FluidState>
|
||||||
|
void updateBoundary_(const Context& context OPM_UNUSED,
|
||||||
|
unsigned bfIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED,
|
||||||
|
const FluidState& fs OPM_UNUSED)
|
||||||
|
{}
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief The temperature gradient times the face normal [K m^2 / m]
|
||||||
|
*/
|
||||||
|
Scalar temperatureGradNormal() const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Calling temperatureGradNormal() does not make sense "
|
||||||
|
"for isothermal models");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief The total thermal conductivity at the face \f$\mathrm{[W/m^2 / (K/m)]}\f$
|
||||||
|
*/
|
||||||
|
Scalar thermalConductivity() const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Calling thermalConductivity() does not make sense for "
|
||||||
|
"isothermal models");
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \copydoc Opm::EnergyExtensiveQuantities
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class EnergyExtensiveQuantities<TypeTag, /*enableEnergy=*/true>
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
|
||||||
|
enum { dimWorld = GridView::dimensionworld };
|
||||||
|
typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
|
||||||
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Update the quantities required to calculate
|
||||||
|
* energy fluxes.
|
||||||
|
*/
|
||||||
|
void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& gradCalc = elemCtx.gradientCalculator();
|
||||||
|
Opm::TemperatureCallback<TypeTag> temperatureCallback(elemCtx);
|
||||||
|
|
||||||
|
EvalDimVector temperatureGrad;
|
||||||
|
gradCalc.calculateGradient(temperatureGrad,
|
||||||
|
elemCtx,
|
||||||
|
faceIdx,
|
||||||
|
temperatureCallback);
|
||||||
|
|
||||||
|
// scalar product of temperature gradient and scvf normal
|
||||||
|
const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(faceIdx);
|
||||||
|
|
||||||
|
temperatureGradNormal_ = 0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
temperatureGradNormal_ += (face.normal()[dimIdx]*temperatureGrad[dimIdx]);
|
||||||
|
|
||||||
|
const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
|
||||||
|
const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
|
||||||
|
const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
|
||||||
|
|
||||||
|
// arithmetic mean
|
||||||
|
thermalConductivity_ =
|
||||||
|
0.5 * (intQuantsInside.thermalConductivity() + intQuantsOutside.thermalConductivity());
|
||||||
|
Opm::Valgrind::CheckDefined(thermalConductivity_);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class Context, class FluidState>
|
||||||
|
void updateBoundary_(const Context& context, unsigned bfIdx, unsigned timeIdx, const FluidState& fs)
|
||||||
|
{
|
||||||
|
const auto& stencil = context.stencil(timeIdx);
|
||||||
|
const auto& face = stencil.boundaryFace(bfIdx);
|
||||||
|
|
||||||
|
const auto& elemCtx = context.elementContext();
|
||||||
|
unsigned insideScvIdx = face.interiorIndex();
|
||||||
|
const auto& insideScv = elemCtx.stencil(timeIdx).subControlVolume(insideScvIdx);
|
||||||
|
|
||||||
|
const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
|
||||||
|
const auto& fsInside = intQuantsInside.fluidState();
|
||||||
|
|
||||||
|
// distance between the center of the SCV and center of the boundary face
|
||||||
|
DimVector distVec = face.integrationPos();
|
||||||
|
distVec -= insideScv.geometry().center();
|
||||||
|
|
||||||
|
Scalar tmp = 0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
tmp += distVec[dimIdx] * face.normal()[dimIdx];
|
||||||
|
Scalar dist = tmp;
|
||||||
|
|
||||||
|
// if the following assertation triggers, the center of the
|
||||||
|
// center of the interior SCV was not inside the element!
|
||||||
|
assert(dist > 0);
|
||||||
|
|
||||||
|
// calculate the temperature gradient using two-point gradient
|
||||||
|
// appoximation
|
||||||
|
temperatureGradNormal_ =
|
||||||
|
(fs.temperature(/*phaseIdx=*/0) - fsInside.temperature(/*phaseIdx=*/0)) / dist;
|
||||||
|
|
||||||
|
// take the value for thermal conductivity from the interior finite volume
|
||||||
|
thermalConductivity_ = intQuantsInside.thermalConductivity();
|
||||||
|
}
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief The temperature gradient times the face normal [K m^2 / m]
|
||||||
|
*/
|
||||||
|
const Evaluation& temperatureGradNormal() const
|
||||||
|
{ return temperatureGradNormal_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief The total thermal conductivity at the face \f$\mathrm{[W/m^2 /
|
||||||
|
* (K/m)]}\f$
|
||||||
|
*/
|
||||||
|
const Evaluation& thermalConductivity() const
|
||||||
|
{ return thermalConductivity_; }
|
||||||
|
|
||||||
|
private:
|
||||||
|
Evaluation temperatureGradNormal_;
|
||||||
|
Evaluation thermalConductivity_;
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
35
opm/models/common/flux.hh
Normal file
35
opm/models/common/flux.hh
Normal file
@ -0,0 +1,35 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief This file contains the necessary classes to calculate the
|
||||||
|
* velocity out of a pressure potential gradient.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_VELOCITY_MODULES_HH
|
||||||
|
#define EWOMS_VELOCITY_MODULES_HH
|
||||||
|
|
||||||
|
#include <opm/models/common/darcyfluxmodule.hh>
|
||||||
|
#include <opm/models/common/forchheimerfluxmodule.hh>
|
||||||
|
|
||||||
|
#endif
|
||||||
583
opm/models/common/forchheimerfluxmodule.hh
Normal file
583
opm/models/common/forchheimerfluxmodule.hh
Normal file
@ -0,0 +1,583 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief This file contains the necessary classes to calculate the
|
||||||
|
* volumetric fluxes out of a pressure potential gradient using the
|
||||||
|
* Forchhheimer approach.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH
|
||||||
|
#define EWOMS_FORCHHEIMER_FLUX_MODULE_HH
|
||||||
|
|
||||||
|
#include "darcyfluxmodule.hh"
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
#include <opm/material/common/Exceptions.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
#include <dune/common/fmatrix.hh>
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
NEW_PROP_TAG(MaterialLaw);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerIntensiveQuantities;
|
||||||
|
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerExtensiveQuantities;
|
||||||
|
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerBaseProblem;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Specifies a flux module which uses the Forchheimer relation.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
struct ForchheimerFluxModule
|
||||||
|
{
|
||||||
|
typedef ForchheimerIntensiveQuantities<TypeTag> FluxIntensiveQuantities;
|
||||||
|
typedef ForchheimerExtensiveQuantities<TypeTag> FluxExtensiveQuantities;
|
||||||
|
typedef ForchheimerBaseProblem<TypeTag> FluxBaseProblem;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the flux module.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the defaults for the parameters required by the
|
||||||
|
* Forchheimer velocity approach.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerBaseProblem
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the Ergun coefficient.
|
||||||
|
*
|
||||||
|
* The Ergun coefficient is a measure how much the velocity is
|
||||||
|
* reduced by turbolence. It is a quantity that does not depend on
|
||||||
|
* the fluid phase but only on the porous medium in question. A
|
||||||
|
* value of 0 means that the velocity is not influenced by
|
||||||
|
* turbolence.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Scalar ergunCoefficient(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::ergunCoefficient()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the ratio between the phase mobility
|
||||||
|
* \f$k_{r,\alpha}\f$ and its passability
|
||||||
|
* \f$\eta_{r,\alpha}\f$ for a given fluid phase
|
||||||
|
* \f$\alpha\f$.
|
||||||
|
*
|
||||||
|
* The passability coefficient specifies the influence of the
|
||||||
|
* other fluid phases on the turbolent behaviour of a given fluid
|
||||||
|
* phase. By default it is equal to the relative permeability. The
|
||||||
|
* mobility to passability ratio is the inverse of phase' the viscosity.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Evaluation mobilityPassabilityRatio(Context& context,
|
||||||
|
unsigned spaceIdx,
|
||||||
|
unsigned timeIdx,
|
||||||
|
unsigned phaseIdx) const
|
||||||
|
{
|
||||||
|
return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the intensive quantities for the Forchheimer module
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerIntensiveQuantities
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Returns the Ergun coefficient.
|
||||||
|
*
|
||||||
|
* The Ergun coefficient is a measure how much the velocity is
|
||||||
|
* reduced by turbolence. A value of 0 means that it is not
|
||||||
|
* influenced.
|
||||||
|
*/
|
||||||
|
const Evaluation& ergunCoefficient() const
|
||||||
|
{ return ergunCoefficient_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the passability of a phase.
|
||||||
|
*/
|
||||||
|
const Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
|
||||||
|
{ return mobilityPassabilityRatio_[phaseIdx]; }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& problem = elemCtx.problem();
|
||||||
|
ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||||
|
mobilityPassabilityRatio_[phaseIdx] =
|
||||||
|
problem.mobilityPassabilityRatio(elemCtx,
|
||||||
|
dofIdx,
|
||||||
|
timeIdx,
|
||||||
|
phaseIdx);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
Evaluation ergunCoefficient_;
|
||||||
|
Evaluation mobilityPassabilityRatio_[numPhases];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup FluxModules
|
||||||
|
* \brief Provides the Forchheimer flux module
|
||||||
|
*
|
||||||
|
* The commonly used Darcy relation looses its validity for Reynolds numbers \f$ Re <
|
||||||
|
* 1\f$. If one encounters flow velocities in porous media above this threshold, the
|
||||||
|
* Forchheimer approach can be used. Like the Darcy approach, it is a relation of with
|
||||||
|
* the fluid velocity in terms of the gradient of pressure potential. However, this
|
||||||
|
* relation is not linear (as in the Darcy case) any more.
|
||||||
|
*
|
||||||
|
* Therefore, the Newton scheme is used to solve the Forchheimer equation. This velocity
|
||||||
|
* is then used like the Darcy velocity e.g. by the local residual.
|
||||||
|
*
|
||||||
|
* Note that for Reynolds numbers above \f$\approx 500\f$ the standard Forchheimer
|
||||||
|
* relation also looses it's validity.
|
||||||
|
*
|
||||||
|
* The Forchheimer equation is given by the following relation:
|
||||||
|
*
|
||||||
|
* \f[
|
||||||
|
\nabla p_\alpha - \rho_\alpha \vec{g} =
|
||||||
|
- \frac{\mu_\alpha}{k_{r,\alpha}} K^{-1}\vec{v}_\alpha
|
||||||
|
- \frac{\rho_\alpha C_E}{\eta_{r,\alpha}} \sqrt{K}^{-1}
|
||||||
|
\left| \vec{v}_\alpha \right| \vec{v}_\alpha
|
||||||
|
\f]
|
||||||
|
*
|
||||||
|
* Where \f$C_E\f$ is the modified Ergun parameter and \f$\eta_{r,\alpha}\f$ is the
|
||||||
|
* passability which is given by a closure relation (usually it is assumed to be
|
||||||
|
* identical to the relative permeability). To avoid numerical problems, the relation
|
||||||
|
* implemented by this class multiplies both sides with \f$\frac{k_{r_alpha}}{mu} K\f$,
|
||||||
|
* so we get
|
||||||
|
*
|
||||||
|
* \f[
|
||||||
|
\frac{k_{r_alpha}}{mu} K \left( \nabla p_\alpha - \rho_\alpha \vec{g}\right) =
|
||||||
|
- \vec{v}_\alpha
|
||||||
|
- \frac{\rho_\alpha C_E}{\eta_{r,\alpha}} \frac{k_{r_alpha}}{mu} \sqrt{K}
|
||||||
|
\left| \vec{v}_\alpha \right| \vec{v}_\alpha
|
||||||
|
\f]
|
||||||
|
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class ForchheimerExtensiveQuantities
|
||||||
|
: public DarcyExtensiveQuantities<TypeTag>
|
||||||
|
{
|
||||||
|
typedef DarcyExtensiveQuantities<TypeTag> DarcyExtQuants;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
|
||||||
|
|
||||||
|
enum { dimWorld = GridView::dimensionworld };
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
||||||
|
typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
|
||||||
|
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
||||||
|
typedef Dune::FieldMatrix<Evaluation, dimWorld, dimWorld> DimEvalMatrix;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Return the Ergun coefficent at the face's integration point.
|
||||||
|
*/
|
||||||
|
const Evaluation& ergunCoefficient() const
|
||||||
|
{ return ergunCoefficient_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the ratio of the mobility divided by the passability at the face's
|
||||||
|
* integration point for a given fluid phase.
|
||||||
|
*
|
||||||
|
* Usually, that's the inverse of the viscosity.
|
||||||
|
*/
|
||||||
|
Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
|
||||||
|
{ return mobilityPassabilityRatio_[phaseIdx]; }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
void calculateGradients_(const ElementContext& elemCtx,
|
||||||
|
unsigned faceIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx);
|
||||||
|
|
||||||
|
auto focusDofIdx = elemCtx.focusDofIndex();
|
||||||
|
unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
|
||||||
|
unsigned j = static_cast<unsigned>(this->exteriorDofIdx_);
|
||||||
|
const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
|
||||||
|
const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
|
||||||
|
|
||||||
|
// calculate the square root of the intrinsic permeability
|
||||||
|
assert(isDiagonal_(this->K_));
|
||||||
|
sqrtK_ = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
|
||||||
|
|
||||||
|
// obtain the Ergun coefficient. Lacking better ideas, we use its the arithmetic mean.
|
||||||
|
if (focusDofIdx == i) {
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(intQuantsIn.ergunCoefficient() +
|
||||||
|
Opm::getValue(intQuantsEx.ergunCoefficient()))/2;
|
||||||
|
}
|
||||||
|
else if (focusDofIdx == j)
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(Opm::getValue(intQuantsIn.ergunCoefficient()) +
|
||||||
|
intQuantsEx.ergunCoefficient())/2;
|
||||||
|
else
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(Opm::getValue(intQuantsIn.ergunCoefficient()) +
|
||||||
|
Opm::getValue(intQuantsEx.ergunCoefficient()))/2;
|
||||||
|
|
||||||
|
// obtain the mobility to passability ratio for each phase.
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
unsigned upIdx = static_cast<unsigned>(this->upstreamIndex_(phaseIdx));
|
||||||
|
const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
|
||||||
|
|
||||||
|
if (focusDofIdx == upIdx) {
|
||||||
|
density_[phaseIdx] =
|
||||||
|
up.fluidState().density(phaseIdx);
|
||||||
|
mobilityPassabilityRatio_[phaseIdx] =
|
||||||
|
up.mobilityPassabilityRatio(phaseIdx);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
density_[phaseIdx] =
|
||||||
|
Opm::getValue(up.fluidState().density(phaseIdx));
|
||||||
|
mobilityPassabilityRatio_[phaseIdx] =
|
||||||
|
Opm::getValue(up.mobilityPassabilityRatio(phaseIdx));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class FluidState>
|
||||||
|
void calculateBoundaryGradients_(const ElementContext& elemCtx,
|
||||||
|
unsigned boundaryFaceIdx,
|
||||||
|
unsigned timeIdx,
|
||||||
|
const FluidState& fluidState)
|
||||||
|
{
|
||||||
|
DarcyExtQuants::calculateBoundaryGradients_(elemCtx,
|
||||||
|
boundaryFaceIdx,
|
||||||
|
timeIdx,
|
||||||
|
fluidState);
|
||||||
|
|
||||||
|
auto focusDofIdx = elemCtx.focusDofIndex();
|
||||||
|
unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
|
||||||
|
const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
|
||||||
|
|
||||||
|
// obtain the Ergun coefficient. Because we are on the boundary here, we will
|
||||||
|
// take the Ergun coefficient of the interior
|
||||||
|
if (focusDofIdx == i)
|
||||||
|
ergunCoefficient_ = intQuantsIn.ergunCoefficient();
|
||||||
|
else
|
||||||
|
ergunCoefficient_ = Opm::getValue(intQuantsIn.ergunCoefficient());
|
||||||
|
|
||||||
|
// calculate the square root of the intrinsic permeability
|
||||||
|
assert(isDiagonal_(this->K_));
|
||||||
|
sqrtK_ = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
|
||||||
|
|
||||||
|
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
if (focusDofIdx == i) {
|
||||||
|
density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
|
||||||
|
mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
density_[phaseIdx] =
|
||||||
|
Opm::getValue(intQuantsIn.fluidState().density(phaseIdx));
|
||||||
|
mobilityPassabilityRatio_[phaseIdx] =
|
||||||
|
Opm::getValue(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the volumetric fluxes of all phases
|
||||||
|
*
|
||||||
|
* The pressure potentials and upwind directions must already be
|
||||||
|
* determined before calling this method!
|
||||||
|
*/
|
||||||
|
void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
auto focusDofIdx = elemCtx.focusDofIndex();
|
||||||
|
auto i = asImp_().interiorIndex();
|
||||||
|
auto j = asImp_().exteriorIndex();
|
||||||
|
const auto& intQuantsI = elemCtx.intensiveQuantities(i, timeIdx);
|
||||||
|
const auto& intQuantsJ = elemCtx.intensiveQuantities(j, timeIdx);
|
||||||
|
|
||||||
|
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
|
||||||
|
const auto& normal = scvf.normal();
|
||||||
|
Opm::Valgrind::CheckDefined(normal);
|
||||||
|
|
||||||
|
// obtain the Ergun coefficient from the intensive quantity object. Until a
|
||||||
|
// better method comes along, we use arithmetic averaging.
|
||||||
|
if (focusDofIdx == i)
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(intQuantsI.ergunCoefficient() +
|
||||||
|
Opm::getValue(intQuantsJ.ergunCoefficient())) / 2;
|
||||||
|
else if (focusDofIdx == j)
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(Opm::getValue(intQuantsI.ergunCoefficient()) +
|
||||||
|
intQuantsJ.ergunCoefficient()) / 2;
|
||||||
|
else
|
||||||
|
ergunCoefficient_ =
|
||||||
|
(Opm::getValue(intQuantsI.ergunCoefficient()) +
|
||||||
|
Opm::getValue(intQuantsJ.ergunCoefficient())) / 2;
|
||||||
|
|
||||||
|
///////////////
|
||||||
|
// calculate the weights of the upstream and the downstream control volumes
|
||||||
|
///////////////
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
this->filterVelocity_[phaseIdx] = 0.0;
|
||||||
|
this->volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
calculateForchheimerFlux_(phaseIdx);
|
||||||
|
|
||||||
|
this->volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
|
||||||
|
this->volumeFlux_[phaseIdx] +=
|
||||||
|
this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Calculate the volumetric flux rates of all phases at the domain boundary
|
||||||
|
*/
|
||||||
|
void calculateBoundaryFluxes_(const ElementContext& elemCtx,
|
||||||
|
unsigned bfIdx,
|
||||||
|
unsigned timeIdx)
|
||||||
|
{
|
||||||
|
const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
|
||||||
|
const auto& normal = boundaryFace.normal();
|
||||||
|
|
||||||
|
///////////////
|
||||||
|
// calculate the weights of the upstream and the downstream degrees of freedom
|
||||||
|
///////////////
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
this->filterVelocity_[phaseIdx] = 0.0;
|
||||||
|
this->volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
calculateForchheimerFlux_(phaseIdx);
|
||||||
|
|
||||||
|
this->volumeFlux_[phaseIdx] = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
this->volumeFlux_[phaseIdx] +=
|
||||||
|
this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void calculateForchheimerFlux_(unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
// initial guess: filter velocity is zero
|
||||||
|
DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
|
||||||
|
velocity = 0.0;
|
||||||
|
|
||||||
|
// the change of velocity between two consecutive Newton iterations
|
||||||
|
DimEvalVector deltaV(1e5);
|
||||||
|
// the function value that is to be minimized of the equation that is to be
|
||||||
|
// fulfilled
|
||||||
|
DimEvalVector residual;
|
||||||
|
// derivative of equation that is to be solved
|
||||||
|
DimEvalMatrix gradResid;
|
||||||
|
|
||||||
|
// search by means of the Newton method for a root of Forchheimer equation
|
||||||
|
unsigned newtonIter = 0;
|
||||||
|
while (deltaV.one_norm() > 1e-11) {
|
||||||
|
if (newtonIter >= 50)
|
||||||
|
throw Opm::NumericalIssue("Could not determine Forchheimer velocity within "
|
||||||
|
+std::to_string(newtonIter)+" iterations");
|
||||||
|
++newtonIter;
|
||||||
|
|
||||||
|
// calculate the residual and its Jacobian matrix
|
||||||
|
gradForchheimerResid_(residual, gradResid, phaseIdx);
|
||||||
|
|
||||||
|
// newton method
|
||||||
|
gradResid.solve(deltaV, residual);
|
||||||
|
velocity -= deltaV;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void forchheimerResid_(DimEvalVector& residual, unsigned phaseIdx) const
|
||||||
|
{
|
||||||
|
const DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
|
||||||
|
|
||||||
|
// Obtaining the upstreamed quantities
|
||||||
|
const auto& mobility = this->mobility_[phaseIdx];
|
||||||
|
const auto& density = density_[phaseIdx];
|
||||||
|
const auto& mobilityPassabilityRatio = mobilityPassabilityRatio_[phaseIdx];
|
||||||
|
|
||||||
|
// optain the quantites for the integration point
|
||||||
|
const auto& pGrad = this->potentialGrad_[phaseIdx];
|
||||||
|
|
||||||
|
// residual = v_\alpha
|
||||||
|
residual = velocity;
|
||||||
|
|
||||||
|
// residual += mobility_\alpha K(\grad p_\alpha - \rho_\alpha g)
|
||||||
|
// -> this->K_.usmv(mobility, pGrad, residual);
|
||||||
|
assert(isDiagonal_(this->K_));
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
|
||||||
|
residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx];
|
||||||
|
|
||||||
|
// Forchheimer turbulence correction:
|
||||||
|
//
|
||||||
|
// residual +=
|
||||||
|
// \rho_\alpha
|
||||||
|
// * mobility_\alpha
|
||||||
|
// * C_E / \eta_{r,\alpha}
|
||||||
|
// * abs(v_\alpha) * sqrt(K)*v_\alpha
|
||||||
|
//
|
||||||
|
// -> sqrtK_.usmv(density*mobilityPassabilityRatio*ergunCoefficient_*velocity.two_norm(),
|
||||||
|
// velocity,
|
||||||
|
// residual);
|
||||||
|
Evaluation absVel = 0.0;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
absVel += velocity[dimIdx]*velocity[dimIdx];
|
||||||
|
// the derivatives of the square root of 0 are undefined, so we must guard
|
||||||
|
// against this case
|
||||||
|
if (absVel <= 0.0)
|
||||||
|
absVel = 0.0;
|
||||||
|
else
|
||||||
|
absVel = Toolbox::sqrt(absVel);
|
||||||
|
const auto& alpha = density*mobilityPassabilityRatio*ergunCoefficient_*absVel;
|
||||||
|
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
|
||||||
|
residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx];
|
||||||
|
Opm::Valgrind::CheckDefined(residual);
|
||||||
|
}
|
||||||
|
|
||||||
|
void gradForchheimerResid_(DimEvalVector& residual,
|
||||||
|
DimEvalMatrix& gradResid,
|
||||||
|
unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
// TODO (?) use AD for this.
|
||||||
|
DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
|
||||||
|
forchheimerResid_(residual, phaseIdx);
|
||||||
|
|
||||||
|
Scalar eps = 1e-11;
|
||||||
|
DimEvalVector tmp;
|
||||||
|
for (unsigned i = 0; i < dimWorld; ++i) {
|
||||||
|
Scalar coordEps = std::max(eps, Toolbox::scalarValue(velocity[i]) * (1 + eps));
|
||||||
|
velocity[i] += coordEps;
|
||||||
|
forchheimerResid_(tmp, phaseIdx);
|
||||||
|
tmp -= residual;
|
||||||
|
tmp /= coordEps;
|
||||||
|
gradResid[i] = tmp;
|
||||||
|
velocity[i] -= coordEps;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Check whether all off-diagonal entries of a tensor are zero.
|
||||||
|
*
|
||||||
|
* \param K the tensor that is to be checked.
|
||||||
|
* \return True iff all off-diagonals are zero.
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
bool isDiagonal_(const DimMatrix& K) const
|
||||||
|
{
|
||||||
|
for (unsigned i = 0; i < dimWorld; i++) {
|
||||||
|
for (unsigned j = 0; j < dimWorld; j++) {
|
||||||
|
if (i == j)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
if (std::abs(K[i][j]) > 1e-25)
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
Implementation& asImp_()
|
||||||
|
{ return *static_cast<Implementation *>(this); }
|
||||||
|
|
||||||
|
const Implementation& asImp_() const
|
||||||
|
{ return *static_cast<const Implementation *>(this); }
|
||||||
|
|
||||||
|
protected:
|
||||||
|
// intrinsic permeability tensor and its square root
|
||||||
|
DimVector sqrtK_;
|
||||||
|
|
||||||
|
// Ergun coefficient of all phases at the integration point
|
||||||
|
Evaluation ergunCoefficient_;
|
||||||
|
|
||||||
|
// Passability of all phases at the integration point
|
||||||
|
Evaluation mobilityPassabilityRatio_[numPhases];
|
||||||
|
|
||||||
|
// Density of all phases at the integration point
|
||||||
|
Evaluation density_[numPhases];
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
177
opm/models/common/multiphasebaseextensivequantities.hh
Normal file
177
opm/models/common/multiphasebaseextensivequantities.hh
Normal file
@ -0,0 +1,177 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \copydoc Opm::MultiPhaseBaseExtensiveQuantities
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_MULTI_PHASE_BASE_EXTENSIVE_QUANTITIES_HH
|
||||||
|
#define EWOMS_MULTI_PHASE_BASE_EXTENSIVE_QUANTITIES_HH
|
||||||
|
|
||||||
|
#include "multiphasebaseproperties.hh"
|
||||||
|
|
||||||
|
#include <opm/models/common/quantitycallbacks.hh>
|
||||||
|
#include <ewoms/disc/common/fvbaseextensivequantities.hh>
|
||||||
|
#include <opm/models/utils/parametersystem.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief This class calculates the pressure potential gradients and
|
||||||
|
* the filter velocities for multi-phase flow in porous media
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class MultiPhaseBaseExtensiveQuantities
|
||||||
|
: public GET_PROP_TYPE(TypeTag, DiscExtensiveQuantities)
|
||||||
|
, public GET_PROP_TYPE(TypeTag, FluxModule)::FluxExtensiveQuantities
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, DiscExtensiveQuantities) ParentType;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
|
||||||
|
typedef typename FluxModule::FluxExtensiveQuantities FluxExtensiveQuantities;
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the extensive quantities.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{
|
||||||
|
FluxModule::registerParameters();
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Update the extensive quantities for a given sub-control-volume-face.
|
||||||
|
*
|
||||||
|
* \param elemCtx Reference to the current element context.
|
||||||
|
* \param scvfIdx The local index of the sub-control-volume face for
|
||||||
|
* which the extensive quantities should be calculated.
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
|
||||||
|
{
|
||||||
|
ParentType::update(elemCtx, scvfIdx, timeIdx);
|
||||||
|
|
||||||
|
// compute the pressure potential gradients
|
||||||
|
FluxExtensiveQuantities::calculateGradients_(elemCtx, scvfIdx, timeIdx);
|
||||||
|
|
||||||
|
// Check whether the pressure potential is in the same direction as the face
|
||||||
|
// normal or in the opposite one
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
|
||||||
|
Opm::Valgrind::SetUndefined(upstreamScvIdx_[phaseIdx]);
|
||||||
|
Opm::Valgrind::SetUndefined(downstreamScvIdx_[phaseIdx]);
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
upstreamScvIdx_[phaseIdx] = FluxExtensiveQuantities::upstreamIndex_(phaseIdx);
|
||||||
|
downstreamScvIdx_[phaseIdx] = FluxExtensiveQuantities::downstreamIndex_(phaseIdx);
|
||||||
|
}
|
||||||
|
|
||||||
|
FluxExtensiveQuantities::calculateFluxes_(elemCtx, scvfIdx, timeIdx);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Update the extensive quantities for a given boundary face.
|
||||||
|
*
|
||||||
|
* \param context Reference to the current execution context.
|
||||||
|
* \param bfIdx The local index of the boundary face for which
|
||||||
|
* the extensive quantities should be calculated.
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
* \param fluidState The FluidState on the domain boundary.
|
||||||
|
* \param paramCache The FluidSystem's parameter cache.
|
||||||
|
*/
|
||||||
|
template <class Context, class FluidState>
|
||||||
|
void updateBoundary(const Context& context,
|
||||||
|
unsigned bfIdx,
|
||||||
|
unsigned timeIdx,
|
||||||
|
const FluidState& fluidState)
|
||||||
|
{
|
||||||
|
ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState);
|
||||||
|
|
||||||
|
FluxExtensiveQuantities::calculateBoundaryGradients_(context.elementContext(),
|
||||||
|
bfIdx,
|
||||||
|
timeIdx,
|
||||||
|
fluidState);
|
||||||
|
FluxExtensiveQuantities::calculateBoundaryFluxes_(context.elementContext(),
|
||||||
|
bfIdx,
|
||||||
|
timeIdx);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the local index of the upstream control volume for a given phase as
|
||||||
|
* a function of the normal flux.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase for which the upstream
|
||||||
|
* direction is requested.
|
||||||
|
*/
|
||||||
|
short upstreamIndex(unsigned phaseIdx) const
|
||||||
|
{ return upstreamScvIdx_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the local index of the downstream control volume
|
||||||
|
* for a given phase as a function of the normal flux.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase for which the downstream
|
||||||
|
* direction is requested.
|
||||||
|
*/
|
||||||
|
short downstreamIndex(unsigned phaseIdx) const
|
||||||
|
{ return downstreamScvIdx_[phaseIdx]; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the weight of the upstream control volume
|
||||||
|
* for a given phase as a function of the normal flux.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase
|
||||||
|
*/
|
||||||
|
Scalar upstreamWeight(unsigned phaseIdx OPM_UNUSED) const
|
||||||
|
{ return 1.0; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the weight of the downstream control volume
|
||||||
|
* for a given phase as a function of the normal flux.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase
|
||||||
|
*/
|
||||||
|
Scalar downstreamWeight(unsigned phaseIdx) const
|
||||||
|
{ return 1.0 - upstreamWeight(phaseIdx); }
|
||||||
|
|
||||||
|
private:
|
||||||
|
short upstreamScvIdx_[numPhases];
|
||||||
|
short downstreamScvIdx_[numPhases];
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
253
opm/models/common/multiphasebasemodel.hh
Normal file
253
opm/models/common/multiphasebasemodel.hh
Normal file
@ -0,0 +1,253 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \copydoc Opm::MultiPhaseBaseModel
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_MULTI_PHASE_BASE_MODEL_HH
|
||||||
|
#define EWOMS_MULTI_PHASE_BASE_MODEL_HH
|
||||||
|
|
||||||
|
#include <opm/material/densead/Math.hpp>
|
||||||
|
|
||||||
|
#include "multiphasebaseproperties.hh"
|
||||||
|
#include "multiphasebaseproblem.hh"
|
||||||
|
#include "multiphasebaseextensivequantities.hh"
|
||||||
|
|
||||||
|
#include <opm/models/common/flux.hh>
|
||||||
|
#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
|
||||||
|
|
||||||
|
#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
|
||||||
|
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
|
||||||
|
#include <opm/material/thermal/NullThermalConductionLaw.hpp>
|
||||||
|
#include <opm/material/thermal/NullSolidEnergyLaw.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
template <class TypeTag>
|
||||||
|
class MultiPhaseBaseModel;
|
||||||
|
}
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
//! The generic type tag for problems using the immiscible multi-phase model
|
||||||
|
NEW_TYPE_TAG(MultiPhaseBaseModel, INHERITS_FROM(VtkMultiPhase, VtkTemperature));
|
||||||
|
|
||||||
|
//! Specify the splices of the MultiPhaseBaseModel type tag
|
||||||
|
SET_SPLICES(MultiPhaseBaseModel, SpatialDiscretizationSplice);
|
||||||
|
|
||||||
|
//! Set the default spatial discretization
|
||||||
|
//!
|
||||||
|
//! We use a vertex centered finite volume method by default
|
||||||
|
SET_TAG_PROP(MultiPhaseBaseModel, SpatialDiscretizationSplice, VcfvDiscretization);
|
||||||
|
|
||||||
|
//! set the number of equations to the number of phases
|
||||||
|
SET_INT_PROP(MultiPhaseBaseModel, NumEq, GET_PROP_TYPE(TypeTag, Indices)::numEq);
|
||||||
|
//! The number of phases is determined by the fluid system
|
||||||
|
SET_INT_PROP(MultiPhaseBaseModel, NumPhases, GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases);
|
||||||
|
//! Number of chemical species in the system
|
||||||
|
SET_INT_PROP(MultiPhaseBaseModel, NumComponents, GET_PROP_TYPE(TypeTag, FluidSystem)::numComponents);
|
||||||
|
|
||||||
|
//! The type of the base base class for actual problems
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel, BaseProblem, Opm::MultiPhaseBaseProblem<TypeTag>);
|
||||||
|
|
||||||
|
//! By default, use the Darcy relation to determine the phase velocity
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel, FluxModule, Opm::DarcyFluxModule<TypeTag>);
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the material law to the null law by default.
|
||||||
|
*/
|
||||||
|
SET_PROP(MultiPhaseBaseModel, MaterialLaw)
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef Opm::NullMaterialTraits<Scalar, FluidSystem::numPhases> Traits;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef Opm::NullMaterial<Traits> type;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the property for the material parameters by extracting
|
||||||
|
* it from the material law.
|
||||||
|
*/
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel,
|
||||||
|
MaterialLawParams,
|
||||||
|
typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
|
||||||
|
|
||||||
|
//! set the energy storage law for the solid to the one which assumes zero heat capacity
|
||||||
|
//! by default
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel,
|
||||||
|
SolidEnergyLaw,
|
||||||
|
Opm::NullSolidEnergyLaw<typename GET_PROP_TYPE(TypeTag, Scalar)>);
|
||||||
|
|
||||||
|
//! extract the type of the parameter objects for the solid energy storage law from the
|
||||||
|
//! law itself
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel,
|
||||||
|
SolidEnergyLawParams,
|
||||||
|
typename GET_PROP_TYPE(TypeTag, SolidEnergyLaw)::Params);
|
||||||
|
|
||||||
|
//! set the thermal conduction law to a dummy one by default
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel,
|
||||||
|
ThermalConductionLaw,
|
||||||
|
Opm::NullThermalConductionLaw<typename GET_PROP_TYPE(TypeTag, Scalar)>);
|
||||||
|
|
||||||
|
//! extract the type of the parameter objects for the thermal conduction law from the law
|
||||||
|
//! itself
|
||||||
|
SET_TYPE_PROP(MultiPhaseBaseModel,
|
||||||
|
ThermalConductionLawParams,
|
||||||
|
typename GET_PROP_TYPE(TypeTag, ThermalConductionLaw)::Params);
|
||||||
|
|
||||||
|
//! disable gravity by default
|
||||||
|
SET_BOOL_PROP(MultiPhaseBaseModel, EnableGravity, false);
|
||||||
|
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup MultiPhaseBaseModel
|
||||||
|
* \brief A base class for fully-implicit multi-phase porous-media flow models
|
||||||
|
* which assume multiple fluid phases.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class MultiPhaseBaseModel : public GET_PROP_TYPE(TypeTag, Discretization)
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Discretization) ParentType;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
|
||||||
|
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
|
||||||
|
typedef typename GridView::template Codim<0>::Entity Element;
|
||||||
|
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
enum { numComponents = FluidSystem::numComponents };
|
||||||
|
|
||||||
|
public:
|
||||||
|
MultiPhaseBaseModel(Simulator& simulator)
|
||||||
|
: ParentType(simulator)
|
||||||
|
{ }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the immiscible model.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{
|
||||||
|
ParentType::registerParameters();
|
||||||
|
|
||||||
|
// register runtime parameters of the VTK output modules
|
||||||
|
Opm::VtkMultiPhaseModule<TypeTag>::registerParameters();
|
||||||
|
Opm::VtkTemperatureModule<TypeTag>::registerParameters();
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns true iff a fluid phase is used by the model.
|
||||||
|
*
|
||||||
|
* \param phaseIdx The index of the fluid phase in question
|
||||||
|
*/
|
||||||
|
bool phaseIsConsidered(unsigned phaseIdx OPM_UNUSED) const
|
||||||
|
{ return true; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Compute the total storage inside one phase of all
|
||||||
|
* conservation quantities.
|
||||||
|
*
|
||||||
|
* \copydetails Doxygen::storageParam
|
||||||
|
* \copydetails Doxygen::phaseIdxParam
|
||||||
|
*/
|
||||||
|
void globalPhaseStorage(EqVector& storage, unsigned phaseIdx)
|
||||||
|
{
|
||||||
|
assert(0 <= phaseIdx && phaseIdx < numPhases);
|
||||||
|
|
||||||
|
storage = 0;
|
||||||
|
|
||||||
|
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(this->gridView());
|
||||||
|
std::mutex mutex;
|
||||||
|
#ifdef _OPENMP
|
||||||
|
#pragma omp parallel
|
||||||
|
#endif
|
||||||
|
{
|
||||||
|
// Attention: the variables below are thread specific and thus cannot be
|
||||||
|
// moved in front of the #pragma!
|
||||||
|
unsigned threadId = ThreadManager::threadId();
|
||||||
|
ElementContext elemCtx(this->simulator_);
|
||||||
|
ElementIterator elemIt = threadedElemIt.beginParallel();
|
||||||
|
EqVector tmp;
|
||||||
|
|
||||||
|
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
|
||||||
|
const Element& elem = *elemIt;
|
||||||
|
if (elem.partitionType() != Dune::InteriorEntity)
|
||||||
|
continue; // ignore ghost and overlap elements
|
||||||
|
|
||||||
|
elemCtx.updateStencil(elem);
|
||||||
|
elemCtx.updateIntensiveQuantities(/*timeIdx=*/0);
|
||||||
|
|
||||||
|
const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
|
||||||
|
|
||||||
|
for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
|
||||||
|
const auto& scv = stencil.subControlVolume(dofIdx);
|
||||||
|
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
|
||||||
|
|
||||||
|
tmp = 0;
|
||||||
|
this->localResidual(threadId).addPhaseStorage(tmp,
|
||||||
|
elemCtx,
|
||||||
|
dofIdx,
|
||||||
|
/*timeIdx=*/0,
|
||||||
|
phaseIdx);
|
||||||
|
tmp *= scv.volume()*intQuants.extrusionFactor();
|
||||||
|
|
||||||
|
mutex.lock();
|
||||||
|
storage += tmp;
|
||||||
|
mutex.unlock();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
storage = this->gridView_.comm().sum(storage);
|
||||||
|
}
|
||||||
|
|
||||||
|
void registerOutputModules_()
|
||||||
|
{
|
||||||
|
ParentType::registerOutputModules_();
|
||||||
|
|
||||||
|
// add the VTK output modules which make sense for all multi-phase models
|
||||||
|
this->addOutputModule(new Opm::VtkMultiPhaseModule<TypeTag>(this->simulator_));
|
||||||
|
this->addOutputModule(new Opm::VtkTemperatureModule<TypeTag>(this->simulator_));
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const Implementation& asImp_() const
|
||||||
|
{ return *static_cast<const Implementation *>(this); }
|
||||||
|
};
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
409
opm/models/common/multiphasebaseproblem.hh
Normal file
409
opm/models/common/multiphasebaseproblem.hh
Normal file
@ -0,0 +1,409 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \copydoc Opm::MultiPhaseBaseProblem
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_MULTI_PHASE_BASE_PROBLEM_HH
|
||||||
|
#define EWOMS_MULTI_PHASE_BASE_PROBLEM_HH
|
||||||
|
|
||||||
|
#include "multiphasebaseproperties.hh"
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproblem.hh>
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
|
||||||
|
#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
|
||||||
|
#include <opm/material/common/Means.hpp>
|
||||||
|
#include <opm/material/common/Unused.hpp>
|
||||||
|
#include <opm/material/common/Exceptions.hpp>
|
||||||
|
|
||||||
|
#include <dune/common/fvector.hh>
|
||||||
|
#include <dune/common/fmatrix.hh>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
NEW_PROP_TAG(SolidEnergyLawParams);
|
||||||
|
NEW_PROP_TAG(ThermalConductionLawParams);
|
||||||
|
NEW_PROP_TAG(EnableGravity);
|
||||||
|
NEW_PROP_TAG(FluxModule);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief The base class for the problems of ECFV discretizations which deal
|
||||||
|
* with a multi-phase flow through a porous medium.
|
||||||
|
*/
|
||||||
|
template<class TypeTag>
|
||||||
|
class MultiPhaseBaseProblem
|
||||||
|
: public FvBaseProblem<TypeTag>
|
||||||
|
, public GET_PROP_TYPE(TypeTag, FluxModule)::FluxBaseProblem
|
||||||
|
{
|
||||||
|
//! \cond SKIP_THIS
|
||||||
|
typedef Opm::FvBaseProblem<TypeTag> ParentType;
|
||||||
|
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, SolidEnergyLawParams) SolidEnergyLawParams;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ThermalConductionLawParams) ThermalConductionLawParams;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params MaterialLawParams;
|
||||||
|
|
||||||
|
enum { dimWorld = GridView::dimensionworld };
|
||||||
|
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||||
|
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
|
||||||
|
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
||||||
|
//! \endcond
|
||||||
|
|
||||||
|
public:
|
||||||
|
/*!
|
||||||
|
* \copydoc Problem::FvBaseProblem(Simulator& )
|
||||||
|
*/
|
||||||
|
MultiPhaseBaseProblem(Simulator& simulator)
|
||||||
|
: ParentType(simulator)
|
||||||
|
{ init_(); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Register all run-time parameters for the problem and the model.
|
||||||
|
*/
|
||||||
|
static void registerParameters()
|
||||||
|
{
|
||||||
|
ParentType::registerParameters();
|
||||||
|
|
||||||
|
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableGravity,
|
||||||
|
"Use the gravity correction for the pressure gradients.");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the intrinsic permeability of an intersection.
|
||||||
|
*
|
||||||
|
* This method is specific to the finite volume discretizations. If left unspecified,
|
||||||
|
* it calls the intrinsicPermeability() method for the intersection's interior and
|
||||||
|
* exterior finite volumes and averages them harmonically. Note that if this function
|
||||||
|
* is defined, the intrinsicPermeability() method does not need to be defined by the
|
||||||
|
* problem (if a finite-volume discretization is used).
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
void intersectionIntrinsicPermeability(DimMatrix& result,
|
||||||
|
const Context& context,
|
||||||
|
unsigned intersectionIdx,
|
||||||
|
unsigned timeIdx) const
|
||||||
|
{
|
||||||
|
const auto& scvf = context.stencil(timeIdx).interiorFace(intersectionIdx);
|
||||||
|
|
||||||
|
const DimMatrix& K1 = asImp_().intrinsicPermeability(context, scvf.interiorIndex(), timeIdx);
|
||||||
|
const DimMatrix& K2 = asImp_().intrinsicPermeability(context, scvf.exteriorIndex(), timeIdx);
|
||||||
|
|
||||||
|
// entry-wise harmonic mean. this is almost certainly wrong if
|
||||||
|
// you have off-main diagonal entries in your permeabilities!
|
||||||
|
for (unsigned i = 0; i < dimWorld; ++i)
|
||||||
|
for (unsigned j = 0; j < dimWorld; ++j)
|
||||||
|
result[i][j] = Opm::harmonicMean(K1[i][j], K2[i][j]);
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \name Problem parameters
|
||||||
|
*/
|
||||||
|
// \{
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ at a given position
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::intrinsicPermeability()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the porosity [] of the porous medium for a given
|
||||||
|
* control volume.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Scalar porosity(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::porosity()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the parameter object for the energy storage law of the solid in a
|
||||||
|
* sub-control volume.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
const SolidEnergyLawParams&
|
||||||
|
solidEnergyParams(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::solidEnergyParams()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the parameter object for the thermal conductivity law in a
|
||||||
|
* sub-control volume.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
const ThermalConductionLawParams&
|
||||||
|
thermalConductionParams(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::thermalConductionParams()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Define the tortuosity.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Scalar tortuosity(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::tortuosity()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Define the dispersivity.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Scalar dispersivity(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
throw std::logic_error("Not implemented: Problem::dispersivity()");
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the material law parameters \f$\mathrm{[K]}\f$ within a control volume.
|
||||||
|
*
|
||||||
|
* If you get a compiler error at this method, you set the
|
||||||
|
* MaterialLaw property to something different than
|
||||||
|
* Opm::NullMaterialLaw. In this case, you have to overload the
|
||||||
|
* matererialLaw() method in the derived class!
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
const MaterialLawParams &
|
||||||
|
materialLawParams(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{
|
||||||
|
static MaterialLawParams dummy;
|
||||||
|
return dummy;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
Scalar temperature(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{ return asImp_().temperature(); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
|
||||||
|
*
|
||||||
|
* This is not specific to the discretization. By default it just
|
||||||
|
* throws an exception so it must be overloaded by the problem if
|
||||||
|
* no energy equation is to be used.
|
||||||
|
*/
|
||||||
|
Scalar temperature() const
|
||||||
|
{ throw std::logic_error("Not implemented:temperature() method not implemented by the actual problem"); }
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
|
||||||
|
*
|
||||||
|
* \param context Reference to the object which represents the
|
||||||
|
* current execution context.
|
||||||
|
* \param spaceIdx The local index of spatial entity defined by the context
|
||||||
|
* \param timeIdx The index used by the time discretization.
|
||||||
|
*/
|
||||||
|
template <class Context>
|
||||||
|
const DimVector& gravity(const Context& context OPM_UNUSED,
|
||||||
|
unsigned spaceIdx OPM_UNUSED,
|
||||||
|
unsigned timeIdx OPM_UNUSED) const
|
||||||
|
{ return asImp_().gravity(); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
|
||||||
|
*
|
||||||
|
* This method is used for problems where the gravitational
|
||||||
|
* acceleration does not depend on the spatial position. The
|
||||||
|
* default behaviour is that if the <tt>EnableGravity</tt>
|
||||||
|
* property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds,
|
||||||
|
* else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$.
|
||||||
|
*/
|
||||||
|
const DimVector& gravity() const
|
||||||
|
{ return gravity_; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Mark grid cells for refinement or coarsening
|
||||||
|
*
|
||||||
|
* \return The number of elements marked for refinement or coarsening.
|
||||||
|
*/
|
||||||
|
unsigned markForGridAdaptation()
|
||||||
|
{
|
||||||
|
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||||
|
|
||||||
|
unsigned numMarked = 0;
|
||||||
|
ElementContext elemCtx( this->simulator() );
|
||||||
|
auto gridView = this->simulator().vanguard().gridView();
|
||||||
|
auto& grid = this->simulator().vanguard().grid();
|
||||||
|
auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
|
||||||
|
auto elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
|
||||||
|
for (; elemIt != elemEndIt; ++elemIt)
|
||||||
|
{
|
||||||
|
const auto& element = *elemIt ;
|
||||||
|
elemCtx.updateAll( element );
|
||||||
|
|
||||||
|
// HACK: this should better be part of an AdaptionCriterion class
|
||||||
|
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||||
|
Scalar minSat = 1e100 ;
|
||||||
|
Scalar maxSat = -1e100;
|
||||||
|
size_t nDofs = elemCtx.numDof(/*timeIdx=*/0);
|
||||||
|
for (unsigned dofIdx = 0; dofIdx < nDofs; ++dofIdx)
|
||||||
|
{
|
||||||
|
const auto& intQuant = elemCtx.intensiveQuantities( dofIdx, /*timeIdx=*/0 );
|
||||||
|
minSat = std::min(minSat,
|
||||||
|
Toolbox::value(intQuant.fluidState().saturation(phaseIdx)));
|
||||||
|
maxSat = std::max(maxSat,
|
||||||
|
Toolbox::value(intQuant.fluidState().saturation(phaseIdx)));
|
||||||
|
}
|
||||||
|
|
||||||
|
const Scalar indicator =
|
||||||
|
(maxSat - minSat)/(std::max<Scalar>(0.01, maxSat+minSat)/2);
|
||||||
|
if( indicator > 0.2 && element.level() < 2 ) {
|
||||||
|
grid.mark( 1, element );
|
||||||
|
++ numMarked;
|
||||||
|
}
|
||||||
|
else if ( indicator < 0.025 ) {
|
||||||
|
grid.mark( -1, element );
|
||||||
|
++ numMarked;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
grid.mark( 0, element );
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// get global sum so that every proc is on the same page
|
||||||
|
numMarked = this->simulator().vanguard().grid().comm().sum( numMarked );
|
||||||
|
|
||||||
|
return numMarked;
|
||||||
|
}
|
||||||
|
|
||||||
|
// \}
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/*!
|
||||||
|
* \brief Converts a Scalar value to an isotropic Tensor
|
||||||
|
*
|
||||||
|
* This is convenient e.g. for specifying intrinsic permebilities:
|
||||||
|
* \code{.cpp}
|
||||||
|
* auto permTensor = this->toDimMatrix_(1e-12);
|
||||||
|
* \endcode
|
||||||
|
*
|
||||||
|
* \param val The scalar value which should be expressed as a tensor
|
||||||
|
*/
|
||||||
|
DimMatrix toDimMatrix_(Scalar val) const
|
||||||
|
{
|
||||||
|
DimMatrix ret(0.0);
|
||||||
|
for (unsigned i = 0; i < DimMatrix::rows; ++i)
|
||||||
|
ret[i][i] = val;
|
||||||
|
return ret;
|
||||||
|
}
|
||||||
|
|
||||||
|
DimVector gravity_;
|
||||||
|
|
||||||
|
private:
|
||||||
|
//! Returns the implementation of the problem (i.e. static polymorphism)
|
||||||
|
Implementation& asImp_()
|
||||||
|
{ return *static_cast<Implementation *>(this); }
|
||||||
|
//! \copydoc asImp_()
|
||||||
|
const Implementation& asImp_() const
|
||||||
|
{ return *static_cast<const Implementation *>(this); }
|
||||||
|
|
||||||
|
void init_()
|
||||||
|
{
|
||||||
|
gravity_ = 0.0;
|
||||||
|
if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity))
|
||||||
|
gravity_[dimWorld-1] = -9.81;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
69
opm/models/common/multiphasebaseproperties.hh
Normal file
69
opm/models/common/multiphasebaseproperties.hh
Normal file
@ -0,0 +1,69 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
* \ingroup MultiPhaseBaseModel
|
||||||
|
*
|
||||||
|
* \brief Defines the common properties required by the porous medium
|
||||||
|
* multi-phase models.
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH
|
||||||
|
#define EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
#include <ewoms/io/vtkmultiphasemodule.hh>
|
||||||
|
#include <ewoms/io/vtktemperaturemodule.hh>
|
||||||
|
|
||||||
|
BEGIN_PROPERTIES
|
||||||
|
|
||||||
|
//! The splice to be used for the spatial discretization
|
||||||
|
NEW_PROP_TAG(SpatialDiscretizationSplice);
|
||||||
|
//! Number of fluid phases in the system
|
||||||
|
NEW_PROP_TAG(NumPhases);
|
||||||
|
//! Number of chemical species in the system
|
||||||
|
NEW_PROP_TAG(NumComponents);
|
||||||
|
//! Enumerations used by the model
|
||||||
|
NEW_PROP_TAG(Indices);
|
||||||
|
//! The material law which ought to be used (extracted from the spatial parameters)
|
||||||
|
NEW_PROP_TAG(MaterialLaw);
|
||||||
|
//! The context material law (extracted from the spatial parameters)
|
||||||
|
NEW_PROP_TAG(MaterialLawParams);
|
||||||
|
//! The material law for the energy stored in the solid matrix
|
||||||
|
NEW_PROP_TAG(SolidEnergyLaw);
|
||||||
|
//! The parameters of the material law for energy storage of the solid
|
||||||
|
NEW_PROP_TAG(SolidEnergyLawParams);
|
||||||
|
//! The material law for thermal conduction
|
||||||
|
NEW_PROP_TAG(ThermalConductionLaw);
|
||||||
|
//! The parameters of the material law for thermal conduction
|
||||||
|
NEW_PROP_TAG(ThermalConductionLawParams);
|
||||||
|
//!The fluid systems including the information about the phases
|
||||||
|
NEW_PROP_TAG(FluidSystem);
|
||||||
|
//! Specifies the relation used for velocity
|
||||||
|
NEW_PROP_TAG(FluxModule);
|
||||||
|
|
||||||
|
//! Returns whether gravity is considered in the problem
|
||||||
|
NEW_PROP_TAG(EnableGravity);
|
||||||
|
|
||||||
|
END_PROPERTIES
|
||||||
|
|
||||||
|
#endif
|
||||||
482
opm/models/common/quantitycallbacks.hh
Normal file
482
opm/models/common/quantitycallbacks.hh
Normal file
@ -0,0 +1,482 @@
|
|||||||
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||||
|
// vi: set et ts=4 sw=4 sts=4:
|
||||||
|
/*
|
||||||
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
|
OPM is free software: you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation, either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
OPM is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License
|
||||||
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
Consult the COPYING file in the top-level source directory of this
|
||||||
|
module for the precise wording of the license and the list of
|
||||||
|
copyright holders.
|
||||||
|
*/
|
||||||
|
/*!
|
||||||
|
* \file
|
||||||
|
*
|
||||||
|
* \brief This method contains all callback classes for quantities
|
||||||
|
* that are required by some extensive quantities
|
||||||
|
*/
|
||||||
|
#ifndef EWOMS_QUANTITY_CALLBACKS_HH
|
||||||
|
#define EWOMS_QUANTITY_CALLBACKS_HH
|
||||||
|
|
||||||
|
#include <ewoms/disc/common/fvbaseproperties.hh>
|
||||||
|
|
||||||
|
#include <opm/material/common/MathToolbox.hpp>
|
||||||
|
#include <opm/material/common/Valgrind.hpp>
|
||||||
|
|
||||||
|
#include <type_traits>
|
||||||
|
#include <utility>
|
||||||
|
|
||||||
|
namespace Opm {
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for temperature.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class TemperatureCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
typedef decltype(std::declval<IQFluidState>().temperature(0)) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
TemperatureCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the temperature given the index of a degree of freedom within an
|
||||||
|
* element context.
|
||||||
|
*
|
||||||
|
* In this context, we assume that thermal equilibrium applies, i.e. that the
|
||||||
|
* temperature of all phases is equal.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{ return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().temperature(/*phaseIdx=*/0); }
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for a phase pressure.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class PressureCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
typedef decltype(std::declval<IQFluidState>().pressure(0)) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
PressureCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{ Opm::Valgrind::SetUndefined(phaseIdx_); }
|
||||||
|
|
||||||
|
PressureCallback(const ElementContext& elemCtx, unsigned phaseIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which the pressure
|
||||||
|
* should be returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the pressure of the specified phase given the index of a degree of
|
||||||
|
* freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().pressure(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for a phase pressure.
|
||||||
|
*/
|
||||||
|
template <class TypeTag, class FluidState>
|
||||||
|
class BoundaryPressureCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQRawFluidState;
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<IQRawFluidState>::type>::type IQFluidState;
|
||||||
|
typedef typename IQFluidState::Scalar IQScalar;
|
||||||
|
typedef Opm::MathToolbox<IQScalar> Toolbox;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef IQScalar ResultType;
|
||||||
|
|
||||||
|
BoundaryPressureCallback(const ElementContext& elemCtx, const FluidState& boundaryFs)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, boundaryFs_(boundaryFs)
|
||||||
|
{ Opm::Valgrind::SetUndefined(phaseIdx_); }
|
||||||
|
|
||||||
|
BoundaryPressureCallback(const ElementContext& elemCtx,
|
||||||
|
const FluidState& boundaryFs,
|
||||||
|
unsigned phaseIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, boundaryFs_(boundaryFs)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which the pressure
|
||||||
|
* should be returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the pressure of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().pressure(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
IQScalar boundaryValue() const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return boundaryFs_.pressure(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
const FluidState& boundaryFs_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for the density of a phase.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class DensityCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
typedef decltype(std::declval<IQFluidState>().density(0)) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
DensityCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{ Opm::Valgrind::SetUndefined(phaseIdx_); }
|
||||||
|
|
||||||
|
DensityCallback(const ElementContext& elemCtx, unsigned phaseIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which the density
|
||||||
|
* should be returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the density of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().density(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for the molar density of a phase.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class MolarDensityCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef decltype(std::declval<IQFluidState>().molarDensity(0)) ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
MolarDensityCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{ Opm::Valgrind::SetUndefined(phaseIdx_); }
|
||||||
|
|
||||||
|
MolarDensityCallback(const ElementContext& elemCtx, unsigned phaseIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which the molar
|
||||||
|
* density should be returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the molar density of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().molarDensity(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for the viscosity of a phase.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class ViscosityCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
typedef decltype(std::declval<IQFluidState>().viscosity(0)) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
ViscosityCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{ Opm::Valgrind::SetUndefined(phaseIdx_); }
|
||||||
|
|
||||||
|
ViscosityCallback(const ElementContext& elemCtx, unsigned phaseIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which the viscosity
|
||||||
|
* should be returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the viscosity of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().viscosity(phaseIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for the velocity of a phase at the center of a DOF.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class VelocityCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||||
|
|
||||||
|
typedef decltype(IntensiveQuantities().velocityCenter()) ResultRawType;
|
||||||
|
|
||||||
|
enum { dim = GridView::dimensionworld };
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename ResultType::field_type ResultFieldType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultFieldType>::ValueType ResultFieldValueType;
|
||||||
|
|
||||||
|
VelocityCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the velocity of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{ return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).velocityCenter(); }
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for the velocity of a phase at the center of a DOF.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class VelocityComponentCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(IntensiveQuantities().velocityCenter()[0]) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
VelocityComponentCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{ Opm::Valgrind::SetUndefined(dimIdx_); }
|
||||||
|
|
||||||
|
VelocityComponentCallback(const ElementContext& elemCtx, unsigned dimIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, dimIdx_(dimIdx)
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the component of the velocity
|
||||||
|
* which should be returned.
|
||||||
|
*/
|
||||||
|
void setDimIndex(unsigned dimIdx)
|
||||||
|
{ dimIdx_ = dimIdx; }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the velocity of a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(dimIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).velocityCenter()[dimIdx_];
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned dimIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \ingroup Discretization
|
||||||
|
*
|
||||||
|
* \brief Callback class for a mole fraction of a component in a phase.
|
||||||
|
*/
|
||||||
|
template <class TypeTag>
|
||||||
|
class MoleFractionCallback
|
||||||
|
{
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||||
|
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||||
|
|
||||||
|
typedef decltype(std::declval<IntensiveQuantities>().fluidState()) IQFluidState;
|
||||||
|
typedef decltype(std::declval<IQFluidState>().moleFraction(0, 0)) ResultRawType;
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef typename std::remove_const<typename std::remove_reference<ResultRawType>::type>::type ResultType;
|
||||||
|
typedef typename Opm::MathToolbox<ResultType>::ValueType ResultValueType;
|
||||||
|
|
||||||
|
MoleFractionCallback(const ElementContext& elemCtx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
{
|
||||||
|
Opm::Valgrind::SetUndefined(phaseIdx_);
|
||||||
|
Opm::Valgrind::SetUndefined(compIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
MoleFractionCallback(const ElementContext& elemCtx, unsigned phaseIdx, unsigned compIdx)
|
||||||
|
: elemCtx_(elemCtx)
|
||||||
|
, phaseIdx_(static_cast<unsigned short>(phaseIdx))
|
||||||
|
, compIdx_(static_cast<unsigned short>(compIdx))
|
||||||
|
{}
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the fluid phase for which a mole fraction should be
|
||||||
|
* returned.
|
||||||
|
*/
|
||||||
|
void setPhaseIndex(unsigned phaseIdx)
|
||||||
|
{ phaseIdx_ = static_cast<unsigned short>(phaseIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Set the index of the component for which the mole fraction should be
|
||||||
|
* returned.
|
||||||
|
*/
|
||||||
|
void setComponentIndex(unsigned compIdx)
|
||||||
|
{ compIdx_ = static_cast<unsigned short>(compIdx); }
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* \brief Return the mole fraction of a component in a phase given the index of a
|
||||||
|
* degree of freedom within an element context.
|
||||||
|
*/
|
||||||
|
ResultType operator()(unsigned dofIdx) const
|
||||||
|
{
|
||||||
|
Opm::Valgrind::CheckDefined(phaseIdx_);
|
||||||
|
Opm::Valgrind::CheckDefined(compIdx_);
|
||||||
|
return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().moleFraction(phaseIdx_, compIdx_);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
const ElementContext& elemCtx_;
|
||||||
|
unsigned short phaseIdx_;
|
||||||
|
unsigned short compIdx_;
|
||||||
|
};
|
||||||
|
|
||||||
|
} // namespace Opm
|
||||||
|
|
||||||
|
#endif
|
||||||
Loading…
Reference in New Issue
Block a user