Commit 7f50b009 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

introduce discretization folder

A new folder named discretization is introduced on the /dumux level.
In this folder there are all the geometry classes as well as vol vars, flux var cache
and physical laws (darcy, fick etc), since they could be possibly reused in decoupled
models and are not restricted to implicit models. NOTE: The fluxvariables and the flux
variables cache are still in the folder /porousmediumflow/implicit/. It is to be discussed
whether or not they should be moved into the discretization folder as well.
parent fd58b779
add_subdirectory("common")
add_subdirectory("discretization")
add_subdirectory("freeflow")
add_subdirectory("geomechanics")
add_subdirectory("implicit")
......
add_subdirectory("box")
add_subdirectory("cellcentered")
install(FILES
fluxvariablescachevector.hh
fvelementgeometryvector.hh
stencils.hh
subcontrolvolumeface.hh
subcontrolvolume.hh
volumevariablesvector.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/box)
// -*- 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 and mass fluxes of fluid phases over a face of a finite volume by means
* of the Darcy approximation. Specializations are provided for the different discretization methods.
*/
#ifndef DUMUX_DISCRETIZATION_BOX_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_BOX_DARCYS_LAW_HH
#include <memory>
#include <dune/common/float_cmp.hh>
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
namespace Dumux
{
namespace Properties
{
// forward declaration of properties
NEW_PROP_TAG(ProblemEnableGravity);
}
/*!
* \ingroup DarcysLaw
* \brief Specialization of Darcy's Law for the box method.
*/
template <class TypeTag>
class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IndexSet::IndexType IndexType;
typedef typename GridView::ctype CoordScalar;
typedef std::vector<IndexType> Stencil;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> FeCache;
typedef typename FeCache::FiniteElementType::Traits::LocalBasisType FeLocalBasis;
typedef typename FluxVariablesCache::FaceData FaceData;
public:
static Scalar flux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvFace,
const IndexType phaseIdx)
{
// get the precalculated local jacobian and shape values at the integration point
const auto& faceData = problem.model().fluxVarsCache()[scvFace].faceData();
const auto& fvGeometry = problem.model().fvGeometries(element);
const auto& insideScv = problem.model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
const auto extrusionFactor = problem.model().curVolVars(insideScv).extrusionFactor();
const auto K = problem.spatialParams().intrinsicPermeability(insideScv);
// evaluate gradP - rho*g at integration point
DimVector gradP(0.0);
for (const auto& scv : fvGeometry.scvs())
{
// the global shape function gradient
DimVector gradI;
faceData.jacInvT.mv(faceData.localJacobian[scv.indexInElement()][0], gradI);
gradI *= problem.model().curVolVars(scv).pressure(phaseIdx);
gradP += gradI;
}
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
{
// gravitational acceleration
DimVector g(problem.gravityAtPos(scvFace.center()));
// interpolate the density at the IP
const auto& insideVolVars = problem.model().curVolVars(insideScv);
const auto& outsideScv = problem.model().fvGeometries().subControlVolume(scvFace.outsideScvIdx());
const auto& outsideVolVars = problem.model().curVolVars(outsideScv);
Scalar rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
// turn gravity into a force
g *= rho;
// subtract from pressure gradient
gradP -= g;
}
// apply the permeability and return the flux
auto KGradP = applyPermeability(K, gradP);
return -1.0*(KGradP*scvFace.unitOuterNormal())*scvFace.area()*extrusionFactor;
}
static DimVector applyPermeability(const DimWorldMatrix& K, const DimVector& gradI)
{
DimVector result(0.0);
K.mv(gradI, result);
return result;
}
static DimVector applyPermeability(const Scalar k, const DimVector& gradI)
{
DimVector result(gradI);
result *= k;
return result;
}
// This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method.
static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
{
return Stencil(0);
}
static FaceData calculateFaceData(const Problem& problem, const Element& element, const typename Element::Geometry& geometry, const FeLocalBasis& localBasis, const SubControlVolumeFace& scvFace)
{
FaceData faceData;
// evaluate shape functions and gradients at the integration point
const auto ipLocal = geometry.local(scvFace.center());
faceData.jacInvT = geometry.jacobianInverseTransposed(ipLocal);
localBasis.evaluateJacobian(ipLocal, faceData.localJacobian);
localBasis.evaluateFunction(ipLocal, faceData.shapeValues);
return std::move(faceData);
}
};
} // end namespace
#endif
......@@ -21,8 +21,8 @@
* \brief This file contains the data which is required to calculate
* diffusive mass fluxes due to molecular diffusion with Fick's law.
*/
#ifndef DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_HH
#define DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_HH
#ifndef DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
#define DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
#include <dune/common/float_cmp.hh>
......@@ -45,15 +45,10 @@ NEW_PROP_TAG(EffectiveDiffusivityModel);
/*!
* \ingroup CCTpfaFicksLaw
* \brief Evaluates the diffusive mass flux according to Fick's law
* \brief Specialization of Fick's Law for the box method.
*/
template <class TypeTag, typename DiscretizationMethod = void>
class FicksLaw
{};
// Specialization for the CC-Tpfa Method
template <class TypeTag>
class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type >
class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
......@@ -179,201 +174,6 @@ private:
}
};
// Specialization for the Box Method
template <class TypeTag>
class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
void update(const Problem &problem,
const Element& element,
const SubControlVolumeFace &scvFace,
int phaseIdx, int compIdx)
{
DUNE_THROW(Dune::NotImplemented, "Fick's law for the Box method is not yet implemented!");
problemPtr_ = &problem;
scvFacePtr_ = &scvFace;
phaseIdx_ = phaseIdx;
compIdx_ = compIdx;
updateStencil_();
// TODO for non solution dependent diffusion tensors...
}
void update(const Problem &problem, const SubControlVolumeFace &scvFace,
int phaseIdx, int compIdx,
VolumeVariables* boundaryVolVars)
{
boundaryVolVars_ = boundaryVolVars;
update(problem, scvFace, phaseIdx, compIdx);
}
void beginFluxComputation(bool boundaryVolVarsUpdated = false)
{
// diffusion tensors are always solution dependent
updateTransmissibilities_();
// Get the inside volume variables
const auto insideScvIdx = scvFace_().insideScvIdx();
const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
const auto* insideVolVars = &problem_().model().curVolVars(insideScv);
// and the outside volume variables
const VolumeVariables* outsideVolVars;
if (!scvFace_().boundary())
outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx());
else
{
outsideVolVars = boundaryVolVars_;
if (!boundaryVolVarsUpdated)
{
// update the boudary volvars for Dirichlet boundaries
const auto element = problem_().model().fvGeometries().element(insideScv);
const auto dirichletPriVars = problem_().dirichlet(element, scvFace_());
boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv);
}
}
// compute the diffusive flux
const auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_);
const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_);
const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_));
rhoDGradXNormal_ = rho*tij_*(xInside - xOutside);
}
/*!
* \brief A function to calculate the mass flux over a sub control volume face
*
* \param phaseIdx The index of the phase of which the flux is to be calculated
* \param compIdx The index of the transported component
*/
Scalar flux() const
{
return rhoDGradXNormal_;
}
std::set<IndexType> stencil() const
{
return stencil_;
}
protected:
void updateTransmissibilities_()
{
const auto insideScvIdx = scvFace_().insideScvIdx();
const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
const auto& insideVolVars = problem_().model().curVolVars(insideScvIdx);
auto insideD = insideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx_), insideD);
Scalar ti = calculateOmega_(insideD, insideScv);
if (!scvFace_().boundary())
{
const auto outsideScvIdx = scvFace_().outsideScvIdx();
const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
const auto& outsideVolVars = problem_().model().curVolVars(outsideScvIdx);
auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx_), outsideD);
Scalar tj = -1.0*calculateOmega_(outsideD, outsideScv);
tij_ = scvFace_().area()*(ti * tj)/(ti + tj);
}
else
{
tij_ = scvFace_().area()*ti;
}
}
void updateStencil_()
{
// fill the stencil
if (!scvFace_().boundary())
stencil_= {scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()};
else
// fill the stencil
stencil_ = {scvFace_().insideScvIdx()};
}
Scalar calculateOmega_(const DimWorldMatrix &D, const SubControlVolume &scv) const
{
GlobalPosition Dnormal;
D.mv(scvFace_().unitOuterNormal(), Dnormal);
auto distanceVector = scvFace_().center();
distanceVector -= scv.center();
distanceVector /= distanceVector.two_norm2();
Scalar omega = Dnormal * distanceVector;
omega *= problem_().model().curVolVars(scv).extrusionFactor();
return omega;
}
Scalar calculateOmega_(Scalar D, const SubControlVolume &scv) const
{
auto distanceVector = scvFace_().center();
distanceVector -= scv.center();
distanceVector /= distanceVector.two_norm2();
Scalar omega = D * (distanceVector * scvFace_().unitOuterNormal());
omega *= problem_().model().curVolVars(scv).extrusionFactor();
return omega;
}
const Problem &problem_() const
{
return *problemPtr_;
}
const SubControlVolumeFace& scvFace_() const
{
return *scvFacePtr_;
}
const Problem *problemPtr_;
const SubControlVolumeFace *scvFacePtr_; //!< Pointer to the sub control volume face for which the flux variables are created
std::set<IndexType> stencil_; //!< Indices of the cells of which the pressure is needed for the flux calculation
//! Boundary volume variables (they only get updated on Dirichlet boundaries)
VolumeVariables* boundaryVolVars_;
IndexType phaseIdx_;
IndexType compIdx_;
//! Precomputed values
Scalar tij_; //!< transmissibility for the flux calculation
Scalar rhoDGradXNormal_; //! rho*D(grad(x))*n
};
} // end namespace
#endif
......@@ -20,8 +20,8 @@
* \file
* \brief Base class for the volume variables vector
*/
#ifndef DUMUX_IMPLICIT_BOX_FLUXVARSCACHEVECTOR_HH
#define DUMUX_IMPLICIT_BOX_FLUXVARSCACHEVECTOR_HH
#ifndef DUMUX_DISCRETIZATION_BOX_FLUXVARSCACHEVECTOR_HH
#define DUMUX_DISCRETIZATION_BOX_FLUXVARSCACHEVECTOR_HH
#include <dumux/implicit/properties.hh>
......
......@@ -22,17 +22,14 @@
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
#ifndef DUMUX_IMPLICIT_BOX_FV_GEOMETRY_VECTOR_HH
#define DUMUX_IMPLICIT_BOX_FV_GEOMETRY_VECTOR_HH
#ifndef DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
#define DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dumux/implicit/box/properties.hh>
#include <dumux/implicit/subcontrolvolume.hh>
#include <dumux/implicit/subcontrolvolumeface.hh>
#include <dumux/implicit/fvelementgeometry.hh>
#include <dumux/common/elementmap.hh>
#include <dumux/common/math.hh>
......
......@@ -20,8 +20,8 @@
* \file
* \brief Implements the notion of stencils for vertex-centered models
*/
#ifndef DUMUX_BOX_STENCILS_HH
#define DUMUX_BOX_STENCILS_HH
#ifndef DUMUX_DISCRETIZATION_BOX_STENCILS_HH
#define DUMUX_DISCRETIZATION_BOX_STENCILS_HH
#include <set>
#include <dumux/implicit/box/properties.hh>
......
// -*- 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 a sub control volume
*/
#ifndef DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUME_HH
#define DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUME_HH
#include <dune/common/fvector.hh>
#include <dumux/discretization/subcontrolvolumebase.hh>
namespace Dumux
{
template<class G, typename I>
class BoxSubControlVolume : public SubControlVolumeBase<G, I>
{
public:
// exported types
using Geometry = G;
using IndexType = I;
private:
using Scalar = typename Geometry::ctype;
enum { dimworld = Geometry::coorddimension };
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
// the contructor in the box case
BoxSubControlVolume(Geometry geometry,
IndexType scvIdx,
IndexType elementIndex,
IndexType indexInElement,
IndexType dofIndex)
: SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)),
scvIdx_(scvIdx), indexInElement_(std::move(indexInElement)),
dofIndex_(std::move(dofIndex)) {}
IndexType indexInElement() const
{
return indexInElement_;
}
//! The global index of this scv
IndexType index() const
{
return scvIdx_;
}
IndexType dofIndex() const
{
return dofIndex_;
}
GlobalPosition dofPosition() const
{
return this->geometry().corner(indexInElement_);
}
private:
IndexType scvIdx_;
IndexType indexInElement_;
IndexType dofIndex_;
};
} // end namespace
#endif
// -*- 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 a sub control volume face
*/
#ifndef DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUMEFACE_HH
#define DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUMEFACE_HH
#include <utility>
#include <dune/common/fvector.hh>
#include <dumux/discretization/subcontrolvolumefacebase.hh>
namespace Dumux
{
/*!
* \ingroup Discretization
* \brief Class for a sub control volume face in the box method, i.e a part of the boundary
* of a sub control volume we compute fluxes on. We simply use the base class here.
*/
template<class Geometry, typename IndexType>
class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexType>
{
using Scalar = typename Geometry::ctype;
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
using LocalPosition = Dune::FieldVector<Scalar, dim>;
public:
BoxSubControlVolumeFace(const Geometry& geometry,
const GlobalPosition& ipGlobal,
const GlobalPosition& unitOuterNormal,