diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index 109ac32bc948b404b1cab200119e4c4dae3ca99c..d6ed9a01e590cd2ef1a882328545962c2520195b 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -40,7 +40,7 @@ class DarcysLawImplementation;
/*!
* \ingroup CCTpfaDiscretization
* \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation
- * \note Darcy's law is speialized for network and surface grids (i.e. if grid dim < dimWorld)
+ * \note Darcy's law is specialized for network and surface grids (i.e. if grid dim < dimWorld)
* \tparam Scalar the scalar type for scalar physical quantities
* \tparam FVGridGeometry the grid geometry
* \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension
@@ -51,7 +51,7 @@ class CCTpfaDarcysLaw;
/*!
* \ingroup CCTpfaDiscretization
* \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation
- * \note Darcy's law is speialized for network and surface grids (i.e. if grid dim < dimWorld)
+ * \note Darcy's law is specialized for network and surface grids (i.e. if grid dim < dimWorld)
*/
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethod::cctpfa>
diff --git a/dumux/discretization/cellcentered/tpfa/forchheimerslaw.hh b/dumux/discretization/cellcentered/tpfa/forchheimerslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8ec28c739e25678391ae91f23293c80129c62890
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/forchheimerslaw.hh
@@ -0,0 +1,558 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCTpfaDiscretization
+ * \brief Forchheimers's law for cell-centered finite volume schemes with two-point flux approximation
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_FORCHHEIMERS_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_FORCHHEIMERS_LAW_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/typetraits/typetraits.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/tpfa/darcyslaw.hh>
+
+namespace Dumux {
+
+// forward declarations
+template<class TypeTag, DiscretizationMethod discMethod>
+class ForchheimersLawImplementation;
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Forchheimer's law for cell-centered finite volume schemes with two-point flux approximation
+ * \note Forchheimer's law is specialized for network and surface grids (i.e. if grid dim < dimWorld)
+ * \tparam Scalar the scalar type for scalar physical quantities
+ * \tparam FVGridGeometry the grid geometry
+ * \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension
+ */
+template<class Scalar, class FVGridGeometry, bool isNetwork>
+class CCTpfaForchheimersLaw;
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Forchheimer's law for cell-centered finite volume schemes with two-point flux approximation
+ * \note Forchheimer's law is specialized for network and surface grids (i.e. if grid dim < dimWorld)
+ */
+template <class TypeTag>
+class ForchheimersLawImplementation<TypeTag, DiscretizationMethod::cctpfa>
+: public CCTpfaForchheimersLaw<typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ (GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView::dimension < GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView::dimensionworld)>
+{};
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Class that fills the cache corresponding to tpfa Forchheimer's Law
+ */
+template<class FVGridGeometry>
+class TpfaForchheimersLawCacheFiller
+{
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+
+public:
+ //! Function to fill a TpfaForchheimersLawCache of a given scvf
+ //! This interface has to be met by any advection-related cache filler class
+ //! TODO: Probably get cache type out of the filler
+ template<class FluxVariablesCache, class Problem, class ElementVolumeVariables, class FluxVariablesCacheFiller>
+ static void fill(FluxVariablesCache& scvfFluxVarsCache,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCacheFiller& fluxVarsCacheFiller)
+ {
+ scvfFluxVarsCache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+};
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief The cache corresponding to tpfa Forchheimer's Law
+ */
+template<class AdvectionType, class FVGridGeometry>
+class TpfaForchheimersLawCache
+{
+ using Scalar = typename AdvectionType::Scalar;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+ static constexpr int dimWorld = FVGridGeometry::GridView::dimensionworld;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+public:
+ using Filler = TpfaForchheimersLawCacheFiller<FVGridGeometry>;
+
+ template<class Problem, class ElementVolumeVariables>
+ void updateAdvection(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
+ harmonicMeanSqrtK_ = AdvectionType::calculateHarmonicMeanSqrtPermeability(problem, elemVolVars, scvf);
+ }
+
+ const Scalar& advectionTij() const
+ { return tij_; }
+
+ const DimWorldMatrix& harmonicMeanSqrtPermeability() const
+ { return harmonicMeanSqrtK_; }
+
+private:
+ Scalar tij_;
+ DimWorldMatrix harmonicMeanSqrtK_;
+};
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Specialization of the CCTpfaForchheimersLaw grids where dim=dimWorld
+ */
+template<class ScalarType, class FVGridGeometry>
+class CCTpfaForchheimersLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
+{
+ using ThisType = CCTpfaForchheimersLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using DimWorldVector = Dune::FieldVector<ScalarType, dimWorld>;
+ using DimWorldMatrix = Dune::FieldMatrix<ScalarType, dimWorld, dimWorld>;
+
+ using DarcysLaw = CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>;
+
+ public:
+ //! state the scalar type of the law
+ using Scalar = ScalarType;
+
+ //! state the discretization method this implementation belongs to
+ static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
+
+ //! state the type for the corresponding cache
+ using Cache = TpfaForchheimersLawCache<ThisType, FVGridGeometry>;
+
+ /*! \brief Compute the advective flux using the Forchheimer equation
+ *
+ * see e.g. Nield & Bejan: Convection in Porous Media \cite nield2006
+ *
+ * The relative passability \f$ \eta_r\f$ is the "Forchheimer-equivalent" to the relative
+ * permeability \f$ k_r\f$.
+ * We use the same function as for \f$ k_r \f$ (VG, BC, linear) other authors use a simple
+ * power law e.g.: \f$\eta_{rw} = S_w^3\f$
+ *
+ * Some rearrangements have been made to the original Forchheimer relation:
+ *
+ * \f[ \mathbf{v_\alpha} + c_F \sqrt{\mathbf{K}} \frac{\rho_\alpha}{\mu_\alpha }
+ * |\mathbf{v_\alpha}| \mathbf{v_\alpha}
+ * + \frac{k_{r \alpha}}{\mu_\alpha} \mathbf{K} \nabla \left(p_\alpha
+ * + \rho_\alpha g z \right)= 0
+ * \f]
+ *
+ * This already includes the assumption \f$ k_r(S_w) = \eta_r(S_w)\f$:
+ * - \f$\eta_{rw} = S_w^x\f$ looks very similar to e.g. Van Genuchten relative permeabilities
+ * - Fichot, et al. (2006), Nuclear Engineering and Design, state that several authors claim
+ * that \f$ k_r(S_w), \eta_r(S_w)\f$ can be chosen equal
+ * - It leads to the equation not degenerating for the case of \f$S_w=1\f$, because I do not
+ * need to multiply with two different functions, and therefore there are terms not being
+ * zero.
+ * - If this assumption is not to be made: Some regularization needs to be introduced ensuring
+ * that not all terms become zero for \f$S_w=1\f$.
+ *
+ * This non-linear equations is solved for \f$\mathbf{v_\alpha}\f$ using Newton's method
+ * and an analytical derivative w.r.t. \f$\mathbf{v_\alpha}\f$.
+ *
+ * The gradient of the Forchheimer relations looks as follows (mind that \f$\sqrt{\mathbf{K}}\f$
+ * is a tensor):
+ *
+ * \f[ f\left(\mathbf{v_\alpha}\right) =
+ * \left(
+ * \begin{array}{ccc}
+ * 1 & 0 &0 \\
+ * 0 & 1 &0 \\
+ * 0 & 0 &1 \\
+ * \end{array}
+ * \right)
+ * +
+ * c_F \frac{\rho_\alpha}{\mu_\alpha} |\mathbf{v}_\alpha| \sqrt{\mathbf{K}}
+ * +
+ * c_F \frac{\rho_\alpha}{\mu_\alpha}\frac{1}{|\mathbf{v}_\alpha|} \sqrt{\mathbf{K}}
+ * \left(
+ * \begin{array}{ccc}
+ * v_x^2 & v_xv_y & v_xv_z \\
+ * v_yv_x & v_{y}^2 & v_yv_z \\
+ * v_zv_x & v_zv_y &v_{z}^2 \\
+ * \end{array}
+ * \right)
+ * \f]
+ *
+ * \note We restrict the use of Forchheimer's law to diagonal permeability tensors so far. This might be changed to
+ * general tensors using eigenvalue decomposition to get \f$\sqrt{\mathbf{K}}\f$
+ */
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ int phaseIdx,
+ const ElementFluxVarsCache& elemFluxVarsCache)
+ {
+ const auto velocity = velocity_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
+ Scalar flux = velocity * scvf.unitOuterNormal();
+ flux *= scvf.area();
+ return flux;
+ }
+
+ //! The flux variables cache has to be bound to an element prior to flux calculations
+ //! During the binding, the transmissibility will be computed and stored using the method below.
+ template<class Problem, class ElementVolumeVariables>
+ static Scalar calculateTransmissibility(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ return DarcysLaw::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+
+ /*! \brief Returns the harmonic mean of \f$\sqrt{K_0}\f$ and \f$\sqrt{K_1}\f$.
+ *
+ * This is a specialization for scalar-valued permeabilities which returns a tensor with identical diagonal entries.
+ */
+ template<class Problem, class ElementVolumeVariables,
+ bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value,
+ std::enable_if_t<scalarPerm, int> = 0>
+ static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using std::sqrt;
+ Scalar harmonicMeanSqrtK(0.0);
+
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideVolVars = elemVolVars[insideScvIdx];
+ const Scalar Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal());
+ const Scalar sqrtKi = sqrt(Ki);
+
+ if (!scvf.boundary())
+ {
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+ const Scalar Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal());
+ const Scalar sqrtKj = sqrt(Kj);
+ harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal());
+ }
+ else
+ harmonicMeanSqrtK = sqrtKi;
+
+ // Convert to tensor
+ DimWorldMatrix result(0.0);
+ for (int i = 0; i < dimWorld; ++i)
+ result[i][i] = harmonicMeanSqrtK;
+
+ return result;
+ }
+
+ /*! \brief Returns the harmonic mean of \f$\sqrt{\mathbf{K_0}}\f$ and \f$\sqrt{\mathbf{K_1}}\f$.
+ *
+ * This is a specialization for tensor-valued permeabilities.
+ */
+ template<class Problem, class ElementVolumeVariables,
+ bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value,
+ std::enable_if_t<!scalarPerm, int> = 0>
+ static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using std::sqrt;
+ DimWorldMatrix sqrtKi(0.0);
+ DimWorldMatrix sqrtKj(0.0);
+ DimWorldMatrix harmonicMeanSqrtK(0.0);
+
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideVolVars = elemVolVars[insideScvIdx];
+ const auto& Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal());
+ // Make sure the permeability matrix does not have off-diagonal entries
+ // TODO implement method using eigenvalues and eigenvectors for general tensors
+ if (!isDiagonal_(Ki))
+ DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported.");
+
+ for (int i = 0; i < dim; ++i)
+ sqrtKi[i][i] = sqrt(Ki[i][i]);
+
+ if (!scvf.boundary())
+ {
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+ const auto& Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal());
+ // Make sure the permeability matrix does not have off-diagonal entries
+ if (!isDiagonal_(Kj))
+ DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported.");
+
+ for (int i = 0; i < dim; ++i)
+ sqrtKj[i][i] = sqrt(Kj[i][i]);
+
+ harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal());
+ }
+ else
+ harmonicMeanSqrtK = sqrtKi;
+
+ return harmonicMeanSqrtK;
+ }
+
+private:
+
+ //! Computes the face velocity based on the Forchheimer equation
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
+ static DimWorldVector velocity_(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ int phaseIdx,
+ const ElementFluxVarsCache& elemFluxVarsCache)
+ {
+ // Get the volume flux based on Darcy's law. The value returned by this method needs to be multiplied with the
+ // mobility (upwinding).
+ Scalar darcyFlux = DarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
+ auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.mobility(phaseIdx); };
+ const Scalar darcyUpwindMobility = doUpwind_(scvf, elemVolVars, upwindTerm, !std::signbit(darcyFlux) /*insideIsUpstream*/);
+ darcyFlux *= darcyUpwindMobility;
+
+ // Convert to Darcy velocity
+ DimWorldVector darcyVelocity = scvf.unitOuterNormal();
+ darcyVelocity *= darcyFlux / scvf.area();
+
+ // Get the harmonic mean of the square root of permeability
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& sqrtK = fluxVarsCache.harmonicMeanSqrtPermeability();
+
+ // Obtain the Forchheimer coefficient from the spatial parameters
+ const Scalar forchCoeff = problem.spatialParams().forchCoeff(scvf);
+
+ // Initial guess of velocity: Darcy relation
+ DimWorldVector velocity = darcyVelocity;
+
+ DimWorldVector deltaV(0.0); // the change in velocity between Newton iterations
+ DimWorldVector residual(10e10); // the residual (function value that is to be minimized)
+ DimWorldMatrix gradF(0.0); // slope of equation that is to be solved
+
+ // Search by means of the Newton method for a root of Forchheimer equation
+ static const Scalar epsilon = getParamFromGroup<Scalar>(problem.paramGroup(), "Forchheimer.NewtonTolerance", 1e-12);
+ static const std::size_t maxNumIter = getParamFromGroup<std::size_t>(problem.paramGroup(), "Forchheimer.MaxIterations", 30);
+ for (int k = 0; residual.two_norm() > epsilon ; ++k)
+ {
+ if (k >= maxNumIter)
+ DUNE_THROW(NumericalProblem, "could not determine forchheimer velocity within "
+ << k << " iterations");
+
+ // calculate current value of Forchheimer equation
+ forchheimerResidual_(residual,
+ forchCoeff,
+ sqrtK,
+ darcyVelocity,
+ velocity,
+ elemVolVars,
+ scvf,
+ phaseIdx);
+
+ // newton's method: slope (gradF) and current value (residual) of function is known,
+ forchheimerDerivative_(gradF,
+ forchCoeff,
+ sqrtK,
+ velocity,
+ elemVolVars,
+ scvf,
+ phaseIdx);
+
+ // solve for change in velocity ("x-Axis intercept")
+ gradF.solve(deltaV, residual);
+ velocity -= deltaV;
+ }
+
+ // The fluxvariables expect a value on which an upwinding of the mobility is performed.
+ // We therefore divide by the mobility here.
+ const Scalar forchheimerUpwindMobility = doUpwind_(scvf, elemVolVars, upwindTerm,
+ !std::signbit(velocity * scvf.unitOuterNormal()) /*insideIsUpstream*/);
+
+ // Do not divide by zero. If the mobility is zero, the resulting flux with upwinding will be zero anyway.
+ if (forchheimerUpwindMobility > 0.0)
+ velocity /= forchheimerUpwindMobility;
+
+ return velocity;
+ }
+
+ //! calculate the residual of the Forchheimer equation
+ template<class ElementVolumeVariables>
+ static void forchheimerResidual_(DimWorldVector& residual,
+ const Scalar forchCoeff,
+ const DimWorldMatrix& sqrtK,
+ const DimWorldVector& darcyVelocity,
+ const DimWorldVector& oldVelocity,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx)
+ {
+ residual = oldVelocity;
+
+ // residual += k_r/mu K grad p
+ residual -= darcyVelocity;
+
+ // residual += c_F rho / mu abs(v) sqrt(K) v
+ auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); };
+ bool insideIsUpstream = !std::signbit(oldVelocity * scvf.unitOuterNormal());
+ const Scalar rhoOverMu = doUpwind_(scvf, elemVolVars, upwindTerm, insideIsUpstream);
+ sqrtK.usmv(forchCoeff * rhoOverMu * oldVelocity.two_norm(), oldVelocity, residual);
+ }
+
+ //! The analytical derivative of the the Forcheimer equation's residual
+ template<class ElementVolumeVariables>
+ static void forchheimerDerivative_(DimWorldMatrix& derivative,
+ const Scalar forchCoeff,
+ const DimWorldMatrix& sqrtK,
+ const DimWorldVector& velocity,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx)
+ {
+
+
+ // Initialize because for low velocities, we add and do not set the entries.
+ derivative = 0.0;
+
+ auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); };
+ bool insideIsUpstream = !std::signbit(velocity * scvf.unitOuterNormal());
+
+ const Scalar absV = velocity.two_norm() ;
+ // This part of the derivative is only calculated if the velocity is sufficiently small:
+ // do not divide by zero!
+ // This should not be very bad because the derivative is only intended to give an
+ // approximation of the gradient of the
+ // function that goes into the Newton scheme.
+ // This is important in the case of a one-phase region in two-phase flow. The non-existing
+ // phase has a velocity of zero (kr=0).
+ // f = sqrtK * vTimesV (this is a matrix)
+ // f *= forchCoeff density / |velocity| / viscosity (multiply all entries with scalars)
+ const Scalar rhoOverMu = doUpwind_(scvf, elemVolVars, upwindTerm, insideIsUpstream);
+ if (absV > 1e-20)
+ {
+ for (int i = 0; i < dim; i++)
+ {
+ for (int k = 0; k <= i; k++)
+ {
+ derivative[i][k] = sqrtK[i][i] * velocity[i] * velocity[k] * forchCoeff
+ / absV * rhoOverMu;
+
+ if (k < i)
+ derivative[k][i] = derivative[i][k];
+ }
+ }
+ }
+
+ // add on the main diagonal:
+ // 1 + sqrtK_i forchCoeff density |v| / viscosity
+ for (int i = 0; i < dim; i++)
+ derivative[i][i] += (1.0 + (sqrtK[i][i]*forchCoeff * absV * rhoOverMu));
+ }
+
+ template<class ElementVolumeVariables, class UpwindTermFunction>
+ static Scalar doUpwind_(const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ const UpwindTermFunction& upwindTerm,
+ const bool insideIsUpstream)
+ {
+ static const Scalar upwindWeight = getParam<Scalar>("Implicit.UpwindWeight");
+
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+ if (!insideIsUpstream) // if sign of flux is negative
+ return (upwindWeight*upwindTerm(outsideVolVars)
+ + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
+ else
+ return (upwindWeight*upwindTerm(insideVolVars)
+ + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
+ }
+
+ /*!
+ * \brief Check whether all off-diagonal entries of a tensor are zero.
+ *
+ * \param K the tensor that is to be checked.
+ *
+ * \return True if all off-diagonals are zero.
+ *
+ */
+ static bool isDiagonal_(const DimWorldMatrix & K)
+ {
+ using std::abs;
+ for (int i = 0; i < dim; i++) {
+ for (int k = 0; k < dim; k++) {
+ if ((i != k) && (abs(K[i][k]) >= 1e-25)) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return volVars.permeability(); }
+
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
+};
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Specialization of the CCTpfaForchheimersLaw grids where dim<dimWorld
+ */
+template<class ScalarType, class FVGridGeometry>
+class CCTpfaForchheimersLaw<ScalarType, FVGridGeometry, /*isNetwork*/ true>
+{
+ static_assert(AlwaysFalse<ScalarType>::value, "Forchheimer not implemented for network grids");
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/forchheimerslaw.hh b/dumux/discretization/forchheimerslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8d4de49360c13277f7de52fe08740873b014a466
--- /dev/null
+++ b/dumux/discretization/forchheimerslaw.hh
@@ -0,0 +1,55 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Discretization
+ * \brief Forchheimer's law specialized for different discretization schemes
+ * This file contains the data which is required to calculate
+ * volume and mass fluxes of fluid phases over a face of a finite volume by means
+ * of the Forchheimer approximation. Specializations are provided for the different discretization methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_FORCHHEIMERS_LAW_HH
+#define DUMUX_DISCRETIZATION_FORCHHEIMERS_LAW_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+// forward declaration
+template <class TypeTag, DiscretizationMethod discMethod>
+class ForchheimersLawImplementation
+{
+ static_assert(discMethod == DiscretizationMethod::cctpfa, "Forchheimer only implemented for cctpfa!");
+};
+
+/*!
+ * \ingroup Discretization
+ * \brief Evaluates the normal component of the Forchheimer velocity on a (sub)control volume face.
+ * \note Specializations are provided for the different discretization methods.
+ * These specializations are found in the headers included below.
+ */
+template <class TypeTag>
+using ForchheimersLaw = ForchheimersLawImplementation<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
+
+} // end namespace Dumux
+
+#include <dumux/discretization/cellcentered/tpfa/forchheimerslaw.hh>
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/1ptestproblem.hh b/test/porousmediumflow/1p/implicit/1ptestproblem.hh
index fdad5160ec2cb9894b09bd0f682e3db11a791c13..dc893c2572a3be07652fc97701f59deb9ee67e85 100644
--- a/test/porousmediumflow/1p/implicit/1ptestproblem.hh
+++ b/test/porousmediumflow/1p/implicit/1ptestproblem.hh
@@ -37,6 +37,10 @@
#include "1ptestspatialparams.hh"
+#if FORCHHEIMER
+#include <dumux/discretization/forchheimerslaw.hh>
+#endif
+
namespace Dumux
{
template <class TypeTag>
@@ -71,6 +75,10 @@ SET_PROP(OnePTestTypeTag, SpatialParams)
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using type = OnePTestSpatialParams<FVGridGeometry, Scalar>;
};
+
+#ifdef FORCHHEIMER
+SET_TYPE_PROP(OnePTestTypeTag, AdvectionType, ForchheimersLaw<TypeTag>);
+#endif
}
/*!
diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index 268cfce8876ce6ec5725970d02e2cb398d5b67a4..fcf04f8f5583f31905c9a1de7a5df05b754cc966 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -37,6 +37,16 @@ dune_add_test(NAME test_1pbox
${CMAKE_CURRENT_BINARY_DIR}/1ptestbox-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1pbox test_1pfv.input -Problem.Name 1ptestbox")
+dune_add_test(NAME test_1pforchheimercctpfa
+ SOURCES test_1pfv.cc
+ COMPILE_DEFINITIONS TYPETAG=OnePTestCCTpfaTypeTag FORCHHEIMER=1
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1pforchheimercctpfa.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1ptestforchheimertpfa-00001.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1pforchheimercctpfa"
+ --zeroThreshold {"velocity_liq \(m/s\)":1e-12})
+
# a gstat test (becaue it's a random permeability field we can't test against a reference solution)
dune_add_test(NAME test_1pwithgstat
SOURCES test_1pfv.cc
diff --git a/test/porousmediumflow/1p/implicit/test_1pforchheimercctpfa.input b/test/porousmediumflow/1p/implicit/test_1pforchheimercctpfa.input
new file mode 100644
index 0000000000000000000000000000000000000000..4a13c84901c1e1ce6d3deae077c80be973d598a8
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_1pforchheimercctpfa.input
@@ -0,0 +1,23 @@
+[TimeLoop]
+DtInitial = 1 # [s]
+TEnd = 1 # [s]
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 2 10
+
+[Problem]
+Name = 1ptestforchheimertpfa # name passed to the output routines
+EnableGravity = false
+
+[SpatialParams]
+LensLowerLeft = 0.2 0.2
+LensUpperRight = 0.8 0.8
+
+Permeability = 1e-7 # [m^2]
+PermeabilityLens = 1e-7 # [m^2]
+ForchCoeff = 0.55
+
+[Vtk]
+AddVelocity = true
diff --git a/test/references/1pforchheimercctpfa.vtu b/test/references/1pforchheimercctpfa.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2137998a675f8842813bc64a88a1983ef05742ca
--- /dev/null
+++ b/test/references/1pforchheimercctpfa.vtu
@@ -0,0 +1,56 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="20" NumberOfPoints="33">
+ <CellData Scalars="p" Vectors="velocity_liq (m/s)">
+ <DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
+ 195000 195000 185000 185000 175000 175000 165000 165000 155000 155000 145000 145000
+ 135000 135000 125000 125000 115000 115000 105000 105000
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
+ 0 0.236926 0 -0 0.236926 0 0 0.236926 0 -0 0.236926 0
+ 0 0.236926 0 -0 0.236926 0 0 0.236926 0 -0 0.236926 0
+ 5.82077e-15 0.236926 0 5.82077e-15 0.236926 0 0 0.236926 0 -0 0.236926 0
+ 0 0.236926 0 -0 0.236926 0 0 0.236926 0 -0 0.236926 0
+ 0 0.236926 0 -0 0.236926 0 0 0.236926 0 -0 0.236926 0
+ </DataArray>
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.5 0 0 0 0.1 0 0.5 0.1 0
+ 1 0 0 1 0.1 0 0 0.2 0 0.5 0.2 0
+ 1 0.2 0 0 0.3 0 0.5 0.3 0 1 0.3 0
+ 0 0.4 0 0.5 0.4 0 1 0.4 0 0 0.5 0
+ 0.5 0.5 0 1 0.5 0 0 0.6 0 0.5 0.6 0
+ 1 0.6 0 0 0.7 0 0.5 0.7 0 1 0.7 0
+ 0 0.8 0 0.5 0.8 0 1 0.8 0 0 0.9 0
+ 0.5 0.9 0 1 0.9 0 0 1 0 0.5 1 0
+ 1 1 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 2 3 7 6
+ 3 5 8 7 6 7 10 9 7 8 11 10
+ 9 10 13 12 10 11 14 13 12 13 16 15
+ 13 14 17 16 15 16 19 18 16 17 20 19
+ 18 19 22 21 19 20 23 22 21 22 25 24
+ 22 23 26 25 24 25 28 27 25 26 29 28
+ 27 28 31 30 28 29 32 31
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>