mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-27 09:40:59 -06:00
578 lines
21 KiB
C++
578 lines
21 KiB
C++
// -*- 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 <opm/common/Exceptions.hpp>
|
|
|
|
#include <opm/models/discretization/common/fvbaseproperties.hh>
|
|
|
|
#include <opm/material/common/Valgrind.hpp>
|
|
|
|
#include <dune/common/fvector.hh>
|
|
#include <dune/common/fmatrix.hh>
|
|
|
|
#include <cmath>
|
|
|
|
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
|
|
{
|
|
using FluxIntensiveQuantities = ForchheimerIntensiveQuantities<TypeTag>;
|
|
using FluxExtensiveQuantities = ForchheimerExtensiveQuantities<TypeTag>;
|
|
using FluxBaseProblem = ForchheimerBaseProblem<TypeTag>;
|
|
|
|
/*!
|
|
* \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
|
|
{
|
|
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
|
|
using Evaluation = GetPropType<TypeTag, Properties::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&,
|
|
unsigned,
|
|
unsigned) 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
|
|
{
|
|
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
|
|
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
|
|
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
|
|
|
|
enum { numPhases = getPropValue<TypeTag, Properties::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>
|
|
{
|
|
using DarcyExtQuants = DarcyExtensiveQuantities<TypeTag>;
|
|
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
|
|
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
|
|
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:
|
|
/*!
|
|
* \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() +
|
|
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;
|
|
|
|
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);
|
|
|
|
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_ = 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));
|
|
}
|
|
}
|
|
}
|
|
|
|
/*!
|
|
* \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();
|
|
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];
|
|
}
|
|
}
|
|
|
|
/*!
|
|
* \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 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& 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];
|
|
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
|