diff --git a/dumux/flux/cctpfa/forchheimerslaw.hh b/dumux/flux/cctpfa/forchheimerslaw.hh
index b0b4efe83598f4bb9d79e0b551cdc4c3624bd105..f787f60a5f439fe16ca32d20b1d61cc9d18bbb55 100644
--- a/dumux/flux/cctpfa/forchheimerslaw.hh
+++ b/dumux/flux/cctpfa/forchheimerslaw.hh
@@ -50,7 +50,7 @@ class ForchheimersLawImplementation;
* \tparam GridGeometry the grid geometry
* \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension
*/
-template<class Scalar, class GridGeometry, bool isNetwork>
+template<class Scalar, class GridGeometry, class FluxVariables, bool isNetwork>
class CCTpfaForchheimersLaw;
/*!
@@ -61,8 +61,9 @@ class CCTpfaForchheimersLaw;
template <class TypeTag>
class ForchheimersLawImplementation<TypeTag, DiscretizationMethod::cctpfa>
: public CCTpfaForchheimersLaw<GetPropType<TypeTag, Properties::Scalar>,
- GetPropType<TypeTag, Properties::GridGeometry>,
- (GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension < GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld)>
+ GetPropType<TypeTag, Properties::GridGeometry>,
+ GetPropType<TypeTag, Properties::FluxVariables>,
+ (GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension < GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld)>
{};
/*!
@@ -136,10 +137,10 @@ private:
* \ingroup CCTpfaFlux
* \brief Specialization of the CCTpfaForchheimersLaw grids where dim=dimWorld
*/
-template<class ScalarType, class GridGeometry>
-class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
+template<class ScalarType, class GridGeometry, class FluxVariables>
+class CCTpfaForchheimersLaw<ScalarType, GridGeometry, FluxVariables, /*isNetwork*/ false>
{
- using ThisType = CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>;
+ using ThisType = CCTpfaForchheimersLaw<ScalarType, GridGeometry, FluxVariables, /*isNetwork*/ false>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
@@ -155,6 +156,7 @@ class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
using DimWorldMatrix = Dune::FieldMatrix<ScalarType, dimWorld, dimWorld>;
using DarcysLaw = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>;
+ using UpwindScheme = typename FluxVariables::UpwindScheme;
public:
//! state the scalar type of the law
@@ -354,10 +356,9 @@ private:
{
// 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);
+ Scalar flux = 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;
+ const auto darcyFlux = UpwindScheme::apply(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
// Convert to Darcy velocity
DimWorldVector darcyVelocity = scvf.unitOuterNormal();
@@ -412,8 +413,8 @@ private:
// 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*/);
+ const auto fluxSign = std::copysign(1.0, velocity * scvf.unitOuterNormal());
+ const auto forchheimerUpwindMobility = fluxSign* UpwindScheme::apply(elemVolVars, scvf, upwindTerm, fluxSign, phaseIdx);
// Do not divide by zero. If the mobility is zero, the resulting flux with upwinding will be zero anyway.
if (forchheimerUpwindMobility > 0.0)
@@ -440,8 +441,8 @@ private:
// 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);
+ const auto fluxSign = std::copysign(1.0, oldVelocity * scvf.unitOuterNormal());
+ const Scalar rhoOverMu = fluxSign * UpwindScheme::apply(elemVolVars, scvf, upwindTerm, fluxSign, phaseIdx);
sqrtK.usmv(forchCoeff * rhoOverMu * oldVelocity.two_norm(), oldVelocity, residual);
}
@@ -461,7 +462,6 @@ private:
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:
@@ -473,7 +473,8 @@ private:
// 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);
+ const auto fluxSign = std::copysign(1.0, velocity * scvf.unitOuterNormal());
+ const Scalar rhoOverMu = fluxSign * UpwindScheme::apply(elemVolVars, scvf, upwindTerm, fluxSign, phaseIdx);
if (absV > 1e-20)
{
for (int i = 0; i < dim; i++)
@@ -495,25 +496,6 @@ private:
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>("Flux.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.
*
@@ -554,8 +536,8 @@ private:
* \ingroup CCTpfaFlux
* \brief Specialization of the CCTpfaForchheimersLaw grids where dim<dimWorld
*/
-template<class ScalarType, class GridGeometry>
-class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ true>
+template<class ScalarType, class GridGeometry, class FluxVariables>
+class CCTpfaForchheimersLaw<ScalarType, GridGeometry, FluxVariables, /*isNetwork*/ true>
{
static_assert(AlwaysFalse<ScalarType>::value, "Forchheimer not implemented for network grids");
};
diff --git a/dumux/porousmediumflow/fluxvariables.hh b/dumux/porousmediumflow/fluxvariables.hh
index 1ae4fdbc6c5fcd08d68b926f5cfb2eff5e9f8a1e..3a3478f3138e39948cb1886702ab77332ffc8c2c 100644
--- a/dumux/porousmediumflow/fluxvariables.hh
+++ b/dumux/porousmediumflow/fluxvariables.hh
@@ -40,11 +40,11 @@ namespace Dumux {
* molecular diffusive and heat conduction fluxes.
*
* \param TypeTag The type tag for access to type traits
- * \param UpwindScheme The upwind scheme to be applied to advective fluxes
+ * \param UpScheme The upwind scheme to be applied to advective fluxes
* \note Not all specializations are currently implemented
*/
template<class TypeTag,
- class UpwindScheme = UpwindScheme<GetPropType<TypeTag, Properties::GridGeometry>> >
+ class UpScheme = UpwindScheme<GetPropType<TypeTag, Properties::GridGeometry>> >
class PorousMediumFluxVariables
: public FluxVariablesBase<GetPropType<TypeTag, Properties::Problem>,
typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView,
@@ -61,6 +61,7 @@ class PorousMediumFluxVariables
};
public:
+ using UpwindScheme = UpScheme;
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
using MolecularDiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;