Loading...
Searching...
No Matches
forchheimerfluxmodule.hh
Go to the documentation of this file.
// -*- 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.
*/
#ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH
#define EWOMS_FORCHHEIMER_FLUX_MODULE_HH
#include <opm/common/Exceptions.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/common/darcyfluxmodule.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <stdexcept>
#include <string>
namespace Opm {
template <class TypeTag>
class ForchheimerIntensiveQuantities;
template <class TypeTag>
class ForchheimerExtensiveQuantities;
template <class TypeTag>
class ForchheimerBaseProblem;
template <class TypeTag>
struct ForchheimerFluxModule
{
using FluxIntensiveQuantities = ForchheimerIntensiveQuantities<TypeTag>;
using FluxExtensiveQuantities = ForchheimerExtensiveQuantities<TypeTag>;
using FluxBaseProblem = ForchheimerBaseProblem<TypeTag>;
static void registerParameters()
{}
};
template <class TypeTag>
class ForchheimerBaseProblem
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
template <class Context>
Scalar ergunCoefficient(const Context&,
unsigned,
unsigned) const
{
throw std::logic_error("Not implemented: Problem::ergunCoefficient()");
}
template <class Context>
Evaluation mobilityPassabilityRatio(Context& context,
unsigned spaceIdx,
unsigned timeIdx,
unsigned phaseIdx) const
{
return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx);
}
};
template <class TypeTag>
class ForchheimerIntensiveQuantities
{
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
public:
const Evaluation& ergunCoefficient() const
{ return ergunCoefficient_; }
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_;
std::array<Evaluation, numPhases> mobilityPassabilityRatio_;
};
template <class TypeTag>
class ForchheimerExtensiveQuantities
: public DarcyExtensiveQuantities<TypeTag>
{
using DarcyExtQuants = DarcyExtensiveQuantities<TypeTag>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
using Toolbox = MathToolbox<Evaluation>;
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using DimEvalMatrix = Dune::FieldMatrix<Evaluation, dimWorld, dimWorld>;
public:
const Evaluation& ergunCoefficient() const
{ return ergunCoefficient_; }
Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
{ return mobilityPassabilityRatio_[phaseIdx]; }
protected:
void calculateGradients_(const ElementContext& elemCtx,
unsigned faceIdx,
unsigned timeIdx)
{
DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx);
const auto focusDofIdx = elemCtx.focusDofIndex();
const unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
const 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() +
getValue(intQuantsEx.ergunCoefficient())) / 2;
}
else if (focusDofIdx == j) {
ergunCoefficient_ =
(getValue(intQuantsIn.ergunCoefficient()) +
intQuantsEx.ergunCoefficient()) / 2;
}
else {
ergunCoefficient_ =
(getValue(intQuantsIn.ergunCoefficient()) +
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;
}
const 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] = getValue(up.fluidState().density(phaseIdx));
mobilityPassabilityRatio_[phaseIdx] = 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);
const auto focusDofIdx = elemCtx.focusDofIndex();
const 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_ = 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] = getValue(intQuantsIn.fluidState().density(phaseIdx));
mobilityPassabilityRatio_[phaseIdx] =
getValue(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
}
}
}
void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
const auto focusDofIdx = elemCtx.focusDofIndex();
const auto i = asImp_().interiorIndex();
const 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();
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() +
getValue(intQuantsJ.ergunCoefficient())) / 2;
}
else if (focusDofIdx == j) {
ergunCoefficient_ = (getValue(intQuantsI.ergunCoefficient()) +
intQuantsJ.ergunCoefficient()) / 2;
}
else {
ergunCoefficient_ = (getValue(intQuantsI.ergunCoefficient()) +
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];
}
}
}
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 NumericalProblem("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& mobPassabilityRatio = 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 * mobPassabilityRatio * ergunCoefficient_ * absVel;
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
residual[dimIdx] += sqrtK_[dimIdx] * alpha * velocity[dimIdx];
}
Valgrind::CheckDefined(residual);
}
void gradForchheimerResid_(DimEvalVector& residual,
DimEvalMatrix& gradResid,
unsigned phaseIdx)
{
// TODO (?) use AD for this.
DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
forchheimerResid_(residual, phaseIdx);
constexpr Scalar eps = 1e-11;
DimEvalVector tmp;
for (unsigned i = 0; i < dimWorld; ++i) {
const 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;
}
}
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
std::array<Evaluation, numPhases> mobilityPassabilityRatio_;
// Density of all phases at the integration point
std::array<Evaluation, numPhases> density_;
};
} // namespace Opm
#endif
Opm::DarcyExtensiveQuantities
Provides the Darcy flux module.
Definition darcyfluxmodule.hh:122
Opm::DarcyExtensiveQuantities::calculateBoundaryGradients_
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx, const FluidState &fluidState)
Calculate the gradients at the grid boundary which are required to determine the volumetric fluxes.
Definition darcyfluxmodule.hh:342
Opm::DarcyExtensiveQuantities::calculateGradients_
void calculateGradients_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition darcyfluxmodule.hh:189
Opm::ForchheimerBaseProblem
Provides the defaults for the parameters required by the Forchheimer velocity approach.
Definition forchheimerfluxmodule.hh:85
Opm::ForchheimerBaseProblem::mobilityPassabilityRatio
Evaluation mobilityPassabilityRatio(Context &context, unsigned spaceIdx, unsigned timeIdx, unsigned phaseIdx) const
Returns the ratio between the phase mobility and its passability for a given fluid phase .
Definition forchheimerfluxmodule.hh:119
Opm::ForchheimerBaseProblem::ergunCoefficient
Scalar ergunCoefficient(const Context &, unsigned, unsigned) const
Returns the Ergun coefficient.
Definition forchheimerfluxmodule.hh:100
Opm::ForchheimerExtensiveQuantities
Provides the Forchheimer flux module.
Definition forchheimerfluxmodule.hh:219
Opm::ForchheimerExtensiveQuantities::ergunCoefficient
const Evaluation & ergunCoefficient() const
Return the Ergun coefficent at the face's integration point.
Definition forchheimerfluxmodule.hh:241
Opm::ForchheimerExtensiveQuantities::calculateBoundaryFluxes_
void calculateBoundaryFluxes_(const ElementContext &elemCtx, unsigned bfIdx, unsigned timeIdx)
Calculate the volumetric flux rates of all phases at the domain boundary.
Definition forchheimerfluxmodule.hh:414
Opm::ForchheimerExtensiveQuantities::calculateFluxes_
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Calculate the volumetric fluxes of all phases.
Definition forchheimerfluxmodule.hh:364
Opm::ForchheimerExtensiveQuantities::isDiagonal_
bool isDiagonal_(const DimMatrix &K) const
Check whether all off-diagonal entries of a tensor are zero.
Definition forchheimerfluxmodule.hh:553
Opm::ForchheimerExtensiveQuantities::mobilityPassabilityRatio
Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Return the ratio of the mobility divided by the passability at the face's integration point for a giv...
Definition forchheimerfluxmodule.hh:250
Opm::ForchheimerIntensiveQuantities
Provides the intensive quantities for the Forchheimer module.
Definition forchheimerfluxmodule.hh:134
Opm::ForchheimerIntensiveQuantities::ergunCoefficient
const Evaluation & ergunCoefficient() const
Returns the Ergun coefficient.
Definition forchheimerfluxmodule.hh:148
Opm::ForchheimerIntensiveQuantities::mobilityPassabilityRatio
const Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Returns the passability of a phase.
Definition forchheimerfluxmodule.hh:154
darcyfluxmodule.hh
This file contains the necessary classes to calculate the volumetric fluxes out of a pressure potenti...
fvbaseproperties.hh
Declare the properties used by the infrastructure code of the finite volume discretizations.
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
Opm::ForchheimerFluxModule
Specifies a flux module which uses the Forchheimer relation.
Definition forchheimerfluxmodule.hh:66
Opm::ForchheimerFluxModule::registerParameters
static void registerParameters()
Register all run-time parameters for the flux module.
Definition forchheimerfluxmodule.hh:74