diff --git a/dumux/porousmediumflow/implicit/CMakeLists.txt b/dumux/porousmediumflow/implicit/CMakeLists.txt
index 0540c9217d90efe8131f8e1473fcfb4aaf03a202..e80d7ffae907969199d9f5e8a38a5f650f8fb3e5 100644
--- a/dumux/porousmediumflow/implicit/CMakeLists.txt
+++ b/dumux/porousmediumflow/implicit/CMakeLists.txt
@@ -1,9 +1,7 @@
#install headers
install(FILES
-cpdarcyfluxvariables.hh
-darcyfluxvariables.hh
-forchheimerfluxvariables.hh
-problem.hh
+fluxvariables.hh
+fluxvariablescache.hh
velocityoutput.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/implicit)
diff --git a/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh b/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh
deleted file mode 100644
index dc9aa2c0a3479d757f9ec8091285dedb563df1d4..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh
+++ /dev/null
@@ -1,292 +0,0 @@
-// -*- 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
- * \brief This file contains the data which is required to calculate
- * volume fluxes of fluid phases over a face of a finite volume by means
- * of the Darcy approximation.
- *
- */
-#ifndef DUMUX_CP_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_CP_DARCY_FLUX_VARIABLES_HH
-
-#include <dune/common/float_cmp.hh>
-
-#include <dumux/common/math.hh>
-#include <dumux/common/parameters.hh>
-
-#include <dumux/implicit/properties.hh>
-
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup ImplicitFluxVariables
- * \brief Evaluates the normal component of the Darcy velocity
- * on a (sub)control volume face.
- */
-template <class TypeTag>
-class CpDarcyFluxVariables
-{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldVector<Scalar, dim> DimVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*!
- * \brief Compute / update the flux variables
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- * \todo The fvGeometry should be better initialized, passed and stored as an std::shared_ptr
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- fvGeometryPtr_ = &fvGeometry;
- onBoundary_ = onBoundary;
- faceIdx_ = fIdx;
-
- mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
- calculateVolumeFlux_(problem, element, elemVolVars);
- }
-
-public:
- /*!
- * \brief Return the volumetric flux over a face of a given phase.
- *
- * This is the calculated velocity multiplied by the unit normal
- * and the area of the face. face().normal has already the
- * magnitude of the area.
- *
- * \param phaseIdx index of the phase
- */
- Scalar volumeFlux(const unsigned int phaseIdx) const
- { return volumeFlux_[phaseIdx]; }
-
- /*!
- * \brief Return the velocity of a given phase.
- *
- * This is the full velocity vector on the
- * face (without being multiplied with normal).
- *
- * \param phaseIdx index of the phase
- */
- GlobalPosition velocity(const unsigned int phaseIdx) const
- { return velocity_[phaseIdx] ; }
-
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- unsigned int downstreamIdx(const unsigned phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- unsigned int upstreamIdx(const unsigned phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the SCV (sub-control-volume) face. This may be either
- * a face within the element or a face on the element boundary,
- * depending on the value of onBoundary_.
- */
- const SCVFace &face() const
- {
- if (onBoundary_)
- return fvGeometry_().boundaryFace[faceIdx_];
- else
- return fvGeometry_().subContVolFace[faceIdx_];
- }
-
-protected:
- /*!
- * \brief Actual calculation of the volume flux.
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateVolumeFlux_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the transmissibilities
- const SpatialParams &spatialParams = problem.spatialParams();
-
- const Element& elementI = fvGeometry_().neighbors[face().i];
- FVElementGeometry fvGeometryI;
- fvGeometryI.subContVol[0].global = elementI.geometry().center();
- auto ki = spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0);
- Dune::FieldVector<Scalar, dimWorld> kin;
- ki.mv(face().normal, kin);
- kin /= face().area;
- auto di = face().ipGlobal;
- di -= elementI.geometry().center();
- using std::abs;
- auto ti = abs(di*kin*face().area/(2*di.two_norm2()));
-
- auto tij = 2*ti;
- if (!onBoundary_)
- {
- const Element& elementJ = fvGeometry_().neighbors[face().j];
- FVElementGeometry fvGeometryJ;
- fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
- auto kj = spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0);
- Dune::FieldVector<Scalar, dimWorld> kjn;
- kj.mv(face().normal, kjn);
- kjn /= face().area;
- auto dj = face().ipGlobal;
- dj -= elementJ.geometry().center();
- using std::abs;
- auto tj = abs(dj*kjn*face().area/(2*dj.two_norm2()));
- tij = harmonicMean(ti, tj);
- }
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- const auto& volVarsI = elemVolVars[face().i];
- auto potentialI = volVarsI.pressure(phaseIdx);
-
- const auto& volVarsJ = elemVolVars[face().j];
- auto potentialJ = volVarsJ.pressure(phaseIdx);
-
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
- {
- // calculate the phase density at the integration point. we
- // only do this if the phase is present in both cells
- Scalar SI = volVarsI.fluidState().saturation(phaseIdx);
- Scalar SJ = volVarsJ.fluidState().saturation(phaseIdx);
- Scalar rhoI = volVarsI.fluidState().density(phaseIdx);
- Scalar rhoJ = volVarsJ.fluidState().density(phaseIdx);
- using std::max;
- using std::min;
- Scalar fI = max(0.0, min(SI/1e-5, 0.5));
- Scalar fJ = max(0.0, min(SJ/1e-5, 0.5));
- if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
- // doesn't matter because no wetting phase is present in
- // both cells!
- fI = fJ = 0.5;
- const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
-
- auto globalPosI = elementI.geometry().center();
- potentialI -= density*(problem.gravityAtPos(globalPosI)*globalPosI);
-
- if (onBoundary_)
- {
- potentialJ -= density*(problem.gravityAtPos(face().ipGlobal)*face().ipGlobal);
- }
- else
- {
- const Element& elementJ = fvGeometry_().neighbors[face().j];
- auto globalPosJ = elementJ.geometry().center();
- potentialJ -= density*(problem.gravityAtPos(globalPosJ)*globalPosJ);
- }
- }
-
- auto potentialDiff = potentialI - potentialJ;
-
- // determine the upwind direction
- if (potentialDiff > 0)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
- }
- else
- {
- upstreamIdx_[phaseIdx] = face().j;
- downstreamIdx_[phaseIdx] = face().i;
- }
-
- // obtain the upwind volume variables
- const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
- const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
-
- // set the volume flux
- volumeFlux_[phaseIdx] = tij*potentialDiff;
- volumeFlux_[phaseIdx] *= mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
- + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx);
-
- velocity_[phaseIdx] = face().normal;
- velocity_[phaseIdx] *= volumeFlux_[phaseIdx]/face().area;
- } // over loop all phases
- }
-
- // return const reference to the fvGeometry
- const FVElementGeometry& fvGeometry_() const
- { return *fvGeometryPtr_; }
-
- unsigned int faceIdx_; //!< The index of the sub control volume face
- bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
- unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
- Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
- GlobalPosition velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
- Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
-
-private:
- const FVElementGeometry* fvGeometryPtr_; //!< Information about the geometry of discretization
-};
-
-} // end namespace
-
-#endif // DUMUX_CP_DARCY_FLUX_VARIABLES_HH
diff --git a/dumux/porousmediumflow/implicit/darcyfluxvariables.hh b/dumux/porousmediumflow/implicit/darcyfluxvariables.hh
deleted file mode 100644
index 0257cf95d80354973d6abbf0f065b2811032d1c3..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/implicit/darcyfluxvariables.hh
+++ /dev/null
@@ -1,349 +0,0 @@
-// -*- 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
- * \brief This file contains the data which is required to calculate
- * volume fluxes of fluid phases over a face of a finite volume by means
- * of the Darcy approximation.
- *
- */
-#ifndef DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
-
-#include <dune/common/float_cmp.hh>
-
-#include <dumux/common/math.hh>
-#include <dumux/common/parameters.hh>
-
-#include <dumux/implicit/properties.hh>
-
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup ImplicitFluxVariables
- * \brief Evaluates the normal component of the Darcy velocity
- * on a (sub)control volume face.
- */
-template <class TypeTag>
-class ImplicitDarcyFluxVariables
-{
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*!
- * \brief Compute / update the flux variables
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- * \todo The fvGeometry should be better initialized, passed and stored as an std::shared_ptr
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- fvGeometryPtr_ = &fvGeometry;
- onBoundary_ = onBoundary;
- faceIdx_ = fIdx;
-
- mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
- asImp_().calculateGradients_(problem, element, elemVolVars);
- asImp_().calculateNormalVelocity_(problem, element, elemVolVars);
- }
-
-
- /*!
- * \brief Return the volumetric flux over a face of a given phase.
- *
- * This is the calculated velocity multiplied by the unit normal
- * and the area of the face. face().normal has already the
- * magnitude of the area.
- *
- * \param phaseIdx index of the phase
- */
- Scalar volumeFlux(const unsigned int phaseIdx) const
- { return volumeFlux_[phaseIdx]; }
-
- /*!
- * \brief Return the velocity of a given phase.
- *
- * This is the full velocity vector on the
- * face (without being multiplied with normal).
- *
- * \param phaseIdx index of the phase
- */
- GlobalPosition velocity(const unsigned int phaseIdx) const
- { return velocity_[phaseIdx] ; }
-
- /*!
- * \brief Return intrinsic permeability multiplied with potential
- * gradient multiplied with normal.
- * I.e. everything that does not need upwind decisions.
- *
- * \param phaseIdx index of the phase
- */
- Scalar kGradPNormal(const unsigned int phaseIdx) const
- { return kGradPNormal_[phaseIdx] ; }
-
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- unsigned int downstreamIdx(const unsigned phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- unsigned int upstreamIdx(const unsigned phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the SCV (sub-control-volume) face. This may be either
- * a face within the element or a face on the element boundary,
- * depending on the value of onBoundary_.
- */
- const SCVFace &face() const
- {
- if (onBoundary_)
- return fvGeometry_().boundaryFace[faceIdx_];
- else
- return fvGeometry_().subContVolFace[faceIdx_];
- }
-
-protected:
-
- //! Returns the implementation of the flux variables (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \copydoc asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
-
- /*!
- * \brief Calculation of the potential gradients
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- potentialGrad_[phaseIdx]= 0.0;
-
- for (unsigned int idx = 0;
- idx < face().numFap;
- idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const GlobalPosition &feGrad = face().grad[idx];
-
- // index for the element volume variables
- int volVarsIdx = face().fapIndices[idx];
-
- // the pressure gradient
- GlobalPosition tmp(feGrad);
- tmp *= elemVolVars[volVarsIdx].fluidState().pressure(phaseIdx);
- potentialGrad_[phaseIdx] += tmp;
- }
-
- // correct the pressure gradient by the gravitational acceleration
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
- {
- // average the phase density at the integration point.
- Scalar SI = elemVolVars[face().i].fluidState().saturation(phaseIdx);
- Scalar SJ = elemVolVars[face().j].fluidState().saturation(phaseIdx);
- Scalar rhoI = elemVolVars[face().i].fluidState().density(phaseIdx);
- Scalar rhoJ = elemVolVars[face().j].fluidState().density(phaseIdx);
- // reduce influence if saturation is very small
- using std::max;
- using std::min;
- Scalar fI = max(0.0, min(SI/1e-5, 0.5));
- Scalar fJ = max(0.0, min(SJ/1e-5, 0.5));
- // check whether the phase is not present in both phase
- if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
- fI = fJ = 0.5;
-
- // make gravity acceleration a force
- GlobalPosition f(problem.gravityAtPos(face().ipGlobal));
- f *= (fI*rhoI + fJ*rhoJ)/(fI + fJ); // gravity times averaged density
-
- // calculate the final potential gradient
- potentialGrad_[phaseIdx] -= f;
- } // gravity
- } // loop over all phases
- }
-
- /*!
- * \brief Actual calculation of the normal Darcy velocities.
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateNormalVelocity_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the mean intrinsic permeability
- const SpatialParams &spatialParams = problem.spatialParams();
- DimWorldMatrix K;
- if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvGeometry_(),
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvGeometry_(),
- face().j));
- }
- else
- {
- const Element& elementI = fvGeometry_().neighbors[face().i];
- FVElementGeometry fvGeometryI;
- fvGeometryI.subContVol[0].global = elementI.geometry().center();
-
- const Element& elementJ = fvGeometry_().neighbors[face().j];
- FVElementGeometry fvGeometryJ;
- fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
-
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0),
- spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0));
- }
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // calculate the flux in the normal direction of the
- // current sub control volume face:
- //
- // v = - (K_f grad phi) * n
- // with K_f = rho g / mu K
- //
- // Mind, that the normal has the length of it's area.
- // This means that we are actually calculating
- // Q = - (K grad phi) dot n /|n| * A
-
-
- K.mv(potentialGrad_[phaseIdx], kGradP_[phaseIdx]);
- kGradPNormal_[phaseIdx] = kGradP_[phaseIdx]*face().normal;
-
- // determine the upwind direction
- if (kGradPNormal_[phaseIdx] < 0)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
- }
- else
- {
- upstreamIdx_[phaseIdx] = face().j;
- downstreamIdx_[phaseIdx] = face().i;
- }
-
- // obtain the upwind volume variables
- const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
- const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
-
- // the minus comes from the Darcy relation which states that
- // the flux is from high to low potentials.
- // set the velocity
- velocity_[phaseIdx] = kGradP_[phaseIdx];
- velocity_[phaseIdx] *= - ( mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
- + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx)) ;
-
- // set the volume flux
- volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal;
- }// loop all phases
- }
-
- // set const reference to the fvGeometry
- void setFVGeometryPtr_(const FVElementGeometry& fvGeometry)
- { fvGeometryPtr_ = &fvGeometry; }
-
- // return const reference to the fvGeometry
- const FVElementGeometry& fvGeometry_() const
- { return *fvGeometryPtr_; }
-
- unsigned int faceIdx_; //!< The index of the sub control volume face
- bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
- unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
- Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
- GlobalPosition velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
- Scalar kGradPNormal_[numPhases] ; //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
- GlobalPosition kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
- GlobalPosition potentialGrad_[numPhases] ; //!< Gradient of potential, which drives flow
- Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
-
-private:
- const FVElementGeometry* fvGeometryPtr_; //!< Information about the geometry of discretization
-
-};
-
-} // end namespace
-
-#endif // DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
diff --git a/dumux/porousmediumflow/implicit/forchheimerfluxvariables.hh b/dumux/porousmediumflow/implicit/forchheimerfluxvariables.hh
deleted file mode 100644
index 5fa90c81a374e2e1db3a48caf18a771d50b2b791..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/implicit/forchheimerfluxvariables.hh
+++ /dev/null
@@ -1,410 +0,0 @@
-// -*- 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
- * \brief This file contains the data which is required to calculate
- * volume fluxes of fluid phases over a face of a finite volume by means
- * of the Forchheimer approximation.
- */
-#ifndef DUMUX_IMPLICIT_FORCHHEIMER_FLUX_VARIABLES_HH
-#define DUMUX_IMPLICIT_FORCHHEIMER_FLUX_VARIABLES_HH
-
-#include <dumux/implicit/properties.hh>
-
-#include <dumux/common/parameters.hh>
-#include <dumux/common/math.hh>
-#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(MobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup ImplicitFluxVariables
- * \brief Evaluates the normal component of the Forchheimer velocity
- * on a (sub)control volume face.
- *
- * The commonly used Darcy relation looses it's validity for \f$ Re > 1\f$.
- * If one encounters flow velocities in porous media above this Reynolds number,
- * the Forchheimer relation can be used. Like the Darcy relation, it relates
- * the gradient in potential to velocity.
- * However, this relation is not linear (as in the Darcy case) any more.
- *
- * Therefore, a Newton scheme is implemented in this class in order to calculate a
- * velocity from the current set of variables. This velocity can subsequently be used
- * e.g. by the localresidual.
- *
- * For Reynolds numbers above \f$ 500\f$ the (Standard) forchheimer relation also
- * looses it's validity.
- *
- * The Forchheimer equation looks as follows:
- * \f[ \nabla \left({p_\alpha + \rho_\alpha g z} \right)
- * = - \frac{\mu_\alpha}{k_{r \alpha} K} \underline{v_\alpha}
- * - \frac{c_F}{\eta_{\alpha r} \sqrt{K}} \rho_\alpha
- * |\underline{v_\alpha}| \underline{v_\alpha}
- * \f]
- *
- * For the formulation that is actually used in this class, see the documentation of the
- * function calculating the residual.
- *
- * - This algorithm does not find a solution if the fluid is incompressible and the
- * initial pressure distribution is uniform.
- * - This algorithm assumes the relative passabilities to have the same functional
- * form as the relative permeabilities.
- */
-template <class TypeTag>
-class ImplicitForchheimerFluxVariables
- : public ImplicitDarcyFluxVariables<TypeTag>
-{
- friend class ImplicitDarcyFluxVariables<TypeTag>; // be friends with parent
- typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
-public:
- /*! \copydoc ParentType::update() */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- ParentType::setFVGeometryPtr_(fvGeometry);
- ParentType::onBoundary_ = onBoundary;
- ParentType::faceIdx_ = fIdx;
-
- ParentType::mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
- asImp_().calculateGradients_(problem, element, elemVolVars);
- asImp_().calculateNormalVelocity_(problem, element, elemVolVars);
- }
-
-private:
- //! Returns the implementation of the flux variables (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \copydoc asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
-
-protected:
- /*!
- * \brief Function for calculation of velocities.
- * \note this overloads the method of the darcy flux variables base class
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateNormalVelocity_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // for initial velocity guess calculate base class velocity
- ParentType::calculateNormalVelocity_(problem, element, elemVolVars);
-
- // calculate the mean intrinsic permeability
- const SpatialParams &spatialParams = problem.spatialParams();
- DimWorldMatrix K;
- if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- this->fvGeometry_(),
- this->face().i),
- spatialParams.intrinsicPermeability(element,
- this->fvGeometry_(),
- this->face().j));
- }
- else
- {
- const Element& elementI = this->fvGeometry_().neighbors[this->face().i];
- FVElementGeometry fvGeometryI;
- fvGeometryI.subContVol[0].global = elementI.geometry().center();
-
- const Element& elementJ = this->fvGeometry_().neighbors[this->face().j];
- FVElementGeometry fvGeometryJ;
- fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
-
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0),
- spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0));
- }
-
- // obtain the Forchheimer coefficient from the spatial parameters
- const Scalar forchCoeff = spatialParams.forchCoeff(element,
- this->fvGeometry_(),
- this->face().i);
-
- // Make sure the permeability matrix does not have off-diagonal entries
- assert( isDiagonal_(K) );
-
- DimWorldMatrix sqrtK(0.0);
- using std::sqrt;
- for (int i = 0; i < dim; ++i)
- sqrtK[i][i] = sqrt(K[i][i]);
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // initial guess of velocity: Darcy relation
- // first taken from base class, later overwritten in base class
- GlobalPosition velocity = this->velocity(phaseIdx);
-
- GlobalPosition deltaV; // the change in velocity between Newton iterations
- GlobalPosition residual(10e10); // the residual (function value that is to be minimized)
- GlobalPosition tmp; // temporary variable for numerical differentiation
- DimWorldMatrix gradF; // slope of equation that is to be solved
-
- // search by means of the Newton method for a root of Forchheimer equation
- for (int k = 0; residual.two_norm() > 1e-12 ; ++k) {
-
- if (k >= 30)
- DUNE_THROW(NumericalProblem, "could not determine forchheimer velocity within "
- << k << " iterations");
-
- // calculate current value of Forchheimer equation
- forchheimerResidual_(residual,
- forchCoeff,
- sqrtK,
- K,
- velocity,
- elemVolVars,
- this->potentialGrad_[phaseIdx],
- phaseIdx);
-
- // newton's method: slope (gradF) and current value (residual) of function is known,
- forchheimerDerivative_(gradF,
- forchCoeff,
- sqrtK,
- velocity,
- elemVolVars,
- phaseIdx);
-
- // solve for change in velocity ("x-Axis intercept")
- gradF.solve(deltaV, residual);
- velocity -= deltaV;
- }
-
- this->velocity_[phaseIdx] = velocity ; // store the converged velocity solution
- // velocity dot normal times cross sectional area is volume flux:
- this->volumeFlux_[phaseIdx] = velocity * this->face().normal;
- }// end loop over all phases
- }
-
- /*! \brief Function for calculation of Forchheimer relation.
- *
- * 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[ \underline{v_\alpha} + c_F \sqrt{K} \frac{\rho_\alpha}{\mu_\alpha }
- * |\underline{v_\alpha}| \underline{v_\alpha}
- * + \frac{k_{r \alpha}}{\mu_\alpha} 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$.
- *
- * As long as the correct velocity is not found yet, the above equation is not fulfilled, but
- * the left hand side yields some residual. This function calculates the left hand side of the
- * equation and returns the residual.
- *
- * \param residual The current function value for the given velocity
- * \param forchCoeff The Forchheimer coefficient
- * \param sqrtK A diagonal matrix whose entries are square roots of the permeabilities
- * \param K The permeability matrix
- * \param velocity The current velocity approximation.
- * \param elemVolVars The volume variables of the current element
- * \param potentialGrad The gradient in potential
- * \param phaseIdx The index of the currently considered phase
- */
- void forchheimerResidual_(GlobalPosition & residual,
- const Scalar forchCoeff,
- const DimWorldMatrix & sqrtK,
- const DimWorldMatrix & K,
- const GlobalPosition & velocity,
- const ElementVolumeVariables & elemVolVars,
- const GlobalPosition & potentialGrad,
- const unsigned int phaseIdx) const
- {
- const VolumeVariables upVolVars = elemVolVars[this->upstreamIdx(phaseIdx)];
- const VolumeVariables downVolVars = elemVolVars[this->downstreamIdx(phaseIdx)];
-
- // Obtaining the physical quantities
- Scalar alpha = this->mobilityUpwindWeight_;
- const Scalar mobility = alpha * upVolVars.mobility(phaseIdx)
- + (1. - alpha)* downVolVars.mobility(phaseIdx) ;
-
- const Scalar viscosity = alpha * upVolVars.fluidState().viscosity(phaseIdx)
- + (1. - alpha)* downVolVars.fluidState().viscosity(phaseIdx) ;
-
- const Scalar density = alpha * upVolVars.fluidState().density(phaseIdx)
- + (1. - alpha) * downVolVars.fluidState().density(phaseIdx) ;
-
- // residual = v
- residual = velocity;
-
- // residual += k_r/mu K grad p
- K.usmv(mobility, potentialGrad, residual);
-
- // residual += c_F rho / mu abs(v) sqrt(K) v
- sqrtK.usmv(forchCoeff * density / viscosity * velocity.two_norm(), velocity, residual);
- }
-
-
- /*! \brief Function for calculation of the gradient of Forchheimer relation with respect
- * to velocity.
- *
- * This function already exploits that \f$ \sqrt{K}\f$ is a diagonal matrix. Therefore,
- * we only have to deal with the main entries.
- * The gradient of the Forchheimer relations looks as follows (mind that \f$ \sqrt{K}\f$
- * is a tensor):
- *
- * \f[ f\left(\underline{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} |\underline{v}_\alpha| \sqrt{K}
- * +
- * c_F \frac{\rho_\alpha}{\mu_\alpha}\frac{1}{|\underline{v}_\alpha|} \sqrt{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]
- *
- * \param derivative The gradient of the forchheimer equation
- * \param forchCoeff Forchheimer coefficient
- * \param sqrtK A diagonal matrix whose entries are square roots of the permeabilities.
- * \param velocity The current velocity approximation.
- * \param elemVolVars The volume variables of the current element
- * \param phaseIdx The index of the currently considered phase
- */
- void forchheimerDerivative_(DimWorldMatrix & derivative,
- const Scalar forchCoeff,
- const DimWorldMatrix & sqrtK,
- const GlobalPosition & velocity,
- const ElementVolumeVariables & elemVolVars,
- const unsigned int phaseIdx) const
- {
- const VolumeVariables upVolVars = elemVolVars[this->upstreamIdx(phaseIdx)];
- const VolumeVariables downVolVars = elemVolVars[this->downstreamIdx(phaseIdx)];
-
- // Obtaining the physical quantities
- Scalar alpha = this->mobilityUpwindWeight_;
- const Scalar viscosity = alpha * upVolVars.fluidState().viscosity(phaseIdx)
- + (1. - alpha)* downVolVars.fluidState().viscosity(phaseIdx) ;
-
- const Scalar density = alpha * upVolVars.fluidState().density(phaseIdx)
- + (1. - alpha) * downVolVars.fluidState().density(phaseIdx) ;
-
- // Initialize because for low velocities, we add and do not set the entries.
- derivative = 0.;
-
- 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)
- if(absV > 1e-20)
- for (int i=0; i<dim; i++)
- for (int k=0; k<dim; k++)
- derivative[i][k]= sqrtK[i][i] * velocity[i] * velocity[k] * forchCoeff
- * density / absV / viscosity;
-
- // 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 * density * absV / viscosity ) ) ;
- }
-
- /*!
- * \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.
- *
- */
- bool isDiagonal_(const DimWorldMatrix & K) const
- {
- 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;
- }
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/implicit/problem.hh b/dumux/porousmediumflow/implicit/problem.hh
deleted file mode 100644
index e7696e2617537ab87ba9927a13ae7a2a3b6c6269..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/implicit/problem.hh
+++ /dev/null
@@ -1,193 +0,0 @@
-// -*- 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
- * \brief Base class for all fully implicit porous media problems
- */
-#ifndef DUMUX_IMPLICIT_POROUS_MEDIA_PROBLEM_HH
-#define DUMUX_IMPLICIT_POROUS_MEDIA_PROBLEM_HH
-
-#include <dumux/implicit/problem.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ImplicitBaseProblems
- * \brief Base class for all fully implicit porous media problems
- */
-template<class TypeTag>
-class ImplicitPorousMediaProblem : public ImplicitProblem<TypeTag>
-{
- using ParentType = ImplicitProblem<TypeTag>;
-
- using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- enum {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
-
- using CoordScalar = typename GridView::ctype;
- using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-
-public:
- /*!
- * \brief The constructor
- *
- * \param timeManager The time manager
- * \param gridView The grid view
- * \param verbose Turn verbosity on or off
- */
- ImplicitPorousMediaProblem(TimeManager &timeManager,
- const GridView &gridView,
- const bool verbose = true)
- : ParentType(timeManager, gridView),
- gravity_(0)
- {
- spatialParams_ = std::make_shared<SpatialParams>(asImp_(), gridView);
-
- if (getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.AddVelocity"))
- gravity_[dimWorld-1] = -9.81;
- }
-
- /*!
- * \brief Called by the TimeManager in order to
- * initialize the problem.
- *
- * We initialize the spatial parameters here.
- */
- void init()
- {
- ParentType::init();
- spatialParams_->init();
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief Evaluate the initial phase state at an element (for cc models) or vertex (for box models)
- *
- * \param entity The entity (element or vertex)
- */
- template<class Entity>
- int initialPhasePresence(const Entity& entity)
- {
- static_assert(int(Entity::codimension) == 0 || int(Entity::codimension) == dim, "Entity must be element or vertex");
-
- return asImp_().initialPhasePresenceAtPos(entity.geometry().center());
- }
-
- /*!
- * \brief Evaluate the initial phase state at a given position
- *
- * \param globalPos The global position
- */
- int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
- {
- DUNE_THROW(Dune::InvalidStateException,
- "The problem does not provide a initialPhasePresence() "
- "or initialPhasePresenceAtPos() method.");
- return 0;
- }
-
- /*!
- * \brief Returns the temperature \f$\mathrm{[K]}\f$ at a given global position.
- *
- * This is not specific to the discretization. By default it just
- * calls temperature().
- *
- * \param globalPos The position in global coordinates where the temperature should be specified.
- */
- Scalar temperatureAtPos(const GlobalPosition &globalPos) 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 used.
- */
- Scalar temperature() const
- { DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem"); }
-
- /*!
- * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
- *
- * This is discretization independent interface. By default it
- * just calls gravity().
- */
- const GlobalPosition &gravityAtPos(const GlobalPosition &pos) 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>ProblemEnableGravity</tt>
- * property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds,
- * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$.
- */
- const GlobalPosition &gravity() const
- { return gravity_; }
-
- /*!
- * \brief Returns the spatial parameters object.
- */
- SpatialParams &spatialParams()
- { return *spatialParams_; }
-
- /*!
- * \brief Returns the spatial parameters object.
- */
- const SpatialParams &spatialParams() const
- { return *spatialParams_; }
-
- // \}
-
-protected:
- //! 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); }
-
- GlobalPosition gravity_;
-
- // fluids and material properties
- std::shared_ptr<SpatialParams> spatialParams_;
-};
-
-}
-
-#endif