From 7f50b009b80e632d208f9b0981ce508b714ed25b Mon Sep 17 00:00:00 2001 From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de> Date: Tue, 5 Jul 2016 19:15:14 +0200 Subject: [PATCH] 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. --- dumux/CMakeLists.txt | 1 + dumux/discretization/CMakeLists.txt | 2 + dumux/discretization/box/CMakeLists.txt | 8 + dumux/discretization/box/darcyslaw.hh | 163 ++++++++ dumux/discretization/box/fickslaw.hh | 179 +++++++++ .../box/fluxvariablescachevector.hh | 4 +- .../box/fvelementgeometryvector.hh | 7 +- .../box/stencils.hh | 4 +- dumux/discretization/box/subcontrolvolume.hh | 84 ++++ .../box/subcontrolvolumeface.hh | 61 +++ .../box/volumevariablesvector.hh | 4 +- .../cellcentered/CMakeLists.txt | 7 + .../cellcentered/fluxvariablescachevector.hh | 4 +- .../cellcentered/stencils.hh | 4 +- .../cellcentered/tpfa/CMakeLists.txt | 5 + .../cellcentered/tpfa}/darcyslaw.hh | 123 +----- .../cellcentered/tpfa/fickslaw.hh | 178 ++++++++ .../tpfa/fvelementgeometryvector.hh | 26 +- .../cellcentered/tpfa/subcontrolvolume.hh | 73 ++++ .../cellcentered/tpfa/subcontrolvolumeface.hh | 62 +++ .../cellcentered/volumevariablesvector.hh | 4 +- dumux/discretization/darcyslaw.hh | 47 +++ dumux/discretization/fickslaw.hh | 59 +++ .../fluxvariablesbase.hh | 6 +- .../fvelementgeometry.hh | 9 +- .../subcontrolvolumebase.hh} | 0 .../subcontrolvolumefacebase.hh} | 23 +- .../volumevariables.hh | 6 +- dumux/implicit/box/propertydefaults.hh | 26 +- .../implicit/cellcentered/propertydefaults.hh | 6 +- .../cellcentered/tpfa/propertydefaults.hh | 9 +- dumux/implicit/propertydefaults.hh | 8 +- .../1p/implicit/volumevariables.hh | 2 +- .../2p/implicit/volumevariables.hh | 2 +- dumux/porousmediumflow/CMakeLists.txt | 3 +- .../constitutivelaws/CMakeLists.txt | 5 - .../constitutivelaws/fickslaw.hh | 379 ------------------ .../implicit/fluxvariables.hh | 2 +- 38 files changed, 1009 insertions(+), 586 deletions(-) create mode 100644 dumux/discretization/CMakeLists.txt create mode 100644 dumux/discretization/box/CMakeLists.txt create mode 100644 dumux/discretization/box/darcyslaw.hh create mode 100644 dumux/discretization/box/fickslaw.hh rename dumux/{implicit => discretization}/box/fluxvariablescachevector.hh (98%) rename dumux/{implicit => discretization}/box/fvelementgeometryvector.hh (99%) rename dumux/{implicit => discretization}/box/stencils.hh (98%) create mode 100644 dumux/discretization/box/subcontrolvolume.hh create mode 100644 dumux/discretization/box/subcontrolvolumeface.hh rename dumux/{implicit => discretization}/box/volumevariablesvector.hh (99%) create mode 100644 dumux/discretization/cellcentered/CMakeLists.txt rename dumux/{implicit => discretization}/cellcentered/fluxvariablescachevector.hh (98%) rename dumux/{implicit => discretization}/cellcentered/stencils.hh (98%) create mode 100644 dumux/discretization/cellcentered/tpfa/CMakeLists.txt rename dumux/{porousmediumflow/constitutivelaws => discretization/cellcentered/tpfa}/darcyslaw.hh (60%) create mode 100644 dumux/discretization/cellcentered/tpfa/fickslaw.hh rename dumux/{implicit => discretization}/cellcentered/tpfa/fvelementgeometryvector.hh (93%) create mode 100644 dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh create mode 100644 dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh rename dumux/{implicit => discretization}/cellcentered/volumevariablesvector.hh (99%) create mode 100644 dumux/discretization/darcyslaw.hh create mode 100644 dumux/discretization/fickslaw.hh rename dumux/{implicit => discretization}/fluxvariablesbase.hh (97%) rename dumux/{implicit => discretization}/fvelementgeometry.hh (97%) rename dumux/{implicit/subcontrolvolume.hh => discretization/subcontrolvolumebase.hh} (100%) rename dumux/{implicit/subcontrolvolumeface.hh => discretization/subcontrolvolumefacebase.hh} (87%) rename dumux/{implicit => discretization}/volumevariables.hh (97%) delete mode 100644 dumux/porousmediumflow/constitutivelaws/CMakeLists.txt delete mode 100644 dumux/porousmediumflow/constitutivelaws/fickslaw.hh diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt index 8901582125..a97251962f 100644 --- a/dumux/CMakeLists.txt +++ b/dumux/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory("common") +add_subdirectory("discretization") add_subdirectory("freeflow") add_subdirectory("geomechanics") add_subdirectory("implicit") diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt new file mode 100644 index 0000000000..edb1a46d26 --- /dev/null +++ b/dumux/discretization/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory("box") +add_subdirectory("cellcentered") diff --git a/dumux/discretization/box/CMakeLists.txt b/dumux/discretization/box/CMakeLists.txt new file mode 100644 index 0000000000..b0883c5d5a --- /dev/null +++ b/dumux/discretization/box/CMakeLists.txt @@ -0,0 +1,8 @@ +install(FILES +fluxvariablescachevector.hh +fvelementgeometryvector.hh +stencils.hh +subcontrolvolumeface.hh +subcontrolvolume.hh +volumevariablesvector.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/box) diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh new file mode 100644 index 0000000000..d36eb7f0af --- /dev/null +++ b/dumux/discretization/box/darcyslaw.hh @@ -0,0 +1,163 @@ +// -*- 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 diff --git a/dumux/discretization/box/fickslaw.hh b/dumux/discretization/box/fickslaw.hh new file mode 100644 index 0000000000..8406f61411 --- /dev/null +++ b/dumux/discretization/box/fickslaw.hh @@ -0,0 +1,179 @@ +// -*- 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 + * diffusive mass fluxes due to molecular diffusion with Fick's law. + */ +#ifndef DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH +#define DUMUX_DISCRETIZATION_BOX_FICKS_LAW_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(NumPhases); +NEW_PROP_TAG(FluidSystem); +NEW_PROP_TAG(EffectiveDiffusivityModel); +} + +/*! + * \ingroup CCTpfaFicksLaw + * \brief Specialization of Fick's Law 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 > +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel; + typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; + typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GridView::IndexSet::IndexType IndexType; + typedef typename std::vector<IndexType> Stencil; + + using Element = typename GridView::template Codim<0>::Entity; + + enum { dim = GridView::dimension} ; + enum { dimWorld = GridView::dimensionworld} ; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ; + + typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix; + typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + +public: + + static Scalar flux(const Problem& problem, + const Element& element, + const SubControlVolumeFace& scvFace, + const int phaseIdx, + const int compIdx) + { + // diffusion tensors are always solution dependent + Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx); + + // 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 auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx()); + + // 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)); + + return rho*tij*(xInside - xOutside); + } + + static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { + std::vector<IndexType> stencil; + stencil.clear(); + if (!scvFace.boundary()) + { + stencil.push_back(scvFace.insideScvIdx()); + stencil.push_back(scvFace.outsideScvIdx()); + } + else + stencil.push_back(scvFace.insideScvIdx()); + + return stencil; + } + +private: + + + static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx) + { + Scalar tij; + + 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_(problem, scvFace, 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_(problem, scvFace, outsideD, outsideScv); + + tij = scvFace.area()*(ti * tj)/(ti + tj); + } + else + { + tij = scvFace.area()*ti; + } + + return tij; + } + + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv) + { + 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; + } + + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv) + { + 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; + } +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/box/fluxvariablescachevector.hh b/dumux/discretization/box/fluxvariablescachevector.hh similarity index 98% rename from dumux/implicit/box/fluxvariablescachevector.hh rename to dumux/discretization/box/fluxvariablescachevector.hh index 1208662561..83e69beb5d 100644 --- a/dumux/implicit/box/fluxvariablescachevector.hh +++ b/dumux/discretization/box/fluxvariablescachevector.hh @@ -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> diff --git a/dumux/implicit/box/fvelementgeometryvector.hh b/dumux/discretization/box/fvelementgeometryvector.hh similarity index 99% rename from dumux/implicit/box/fvelementgeometryvector.hh rename to dumux/discretization/box/fvelementgeometryvector.hh index b2da799aca..48c11d229f 100644 --- a/dumux/implicit/box/fvelementgeometryvector.hh +++ b/dumux/discretization/box/fvelementgeometryvector.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> diff --git a/dumux/implicit/box/stencils.hh b/dumux/discretization/box/stencils.hh similarity index 98% rename from dumux/implicit/box/stencils.hh rename to dumux/discretization/box/stencils.hh index 3185f8a55f..8b67b95636 100644 --- a/dumux/implicit/box/stencils.hh +++ b/dumux/discretization/box/stencils.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> diff --git a/dumux/discretization/box/subcontrolvolume.hh b/dumux/discretization/box/subcontrolvolume.hh new file mode 100644 index 0000000000..df886d74d3 --- /dev/null +++ b/dumux/discretization/box/subcontrolvolume.hh @@ -0,0 +1,84 @@ +// -*- 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 diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh new file mode 100644 index 0000000000..b96d5085c5 --- /dev/null +++ b/dumux/discretization/box/subcontrolvolumeface.hh @@ -0,0 +1,61 @@ +// -*- 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, + IndexType scvfIndex, + IndexType indexInElement, + const std::vector<IndexType>& scvIndices, + bool boundary = false) + : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary) {} +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/box/volumevariablesvector.hh b/dumux/discretization/box/volumevariablesvector.hh similarity index 99% rename from dumux/implicit/box/volumevariablesvector.hh rename to dumux/discretization/box/volumevariablesvector.hh index 0b97fdd5c6..4d5e15a686 100644 --- a/dumux/implicit/box/volumevariablesvector.hh +++ b/dumux/discretization/box/volumevariablesvector.hh @@ -20,8 +20,8 @@ * \file * \brief Base class for the volume variables vector */ -#ifndef DUMUX_IMPLICIT_BOX_VOLVARSVECTOR_HH -#define DUMUX_IMPLICIT_BOX_VOLVARSVECTOR_HH +#ifndef DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH +#define DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH #include <dumux/implicit/properties.hh> diff --git a/dumux/discretization/cellcentered/CMakeLists.txt b/dumux/discretization/cellcentered/CMakeLists.txt new file mode 100644 index 0000000000..313cc2adfa --- /dev/null +++ b/dumux/discretization/cellcentered/CMakeLists.txt @@ -0,0 +1,7 @@ +add_subdirectory("tpfa") + +install(FILES +fluxvariablescachevector.hh +stencils.hh +volumevariablesvector.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered) \ No newline at end of file diff --git a/dumux/implicit/cellcentered/fluxvariablescachevector.hh b/dumux/discretization/cellcentered/fluxvariablescachevector.hh similarity index 98% rename from dumux/implicit/cellcentered/fluxvariablescachevector.hh rename to dumux/discretization/cellcentered/fluxvariablescachevector.hh index ffae3dbf96..13026e45e9 100644 --- a/dumux/implicit/cellcentered/fluxvariablescachevector.hh +++ b/dumux/discretization/cellcentered/fluxvariablescachevector.hh @@ -20,8 +20,8 @@ * \file * \brief Base class for the volume variables vector */ -#ifndef DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH -#define DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH +#ifndef DUMUX_DISCRETIZATION_CC_FLUXVARSCACHEVECTOR_HH +#define DUMUX_DISCRETIZATION_CC_FLUXVARSCACHEVECTOR_HH #include <dumux/implicit/properties.hh> diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/discretization/cellcentered/stencils.hh similarity index 98% rename from dumux/implicit/cellcentered/stencils.hh rename to dumux/discretization/cellcentered/stencils.hh index 34e156c31f..b5de5588ef 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/discretization/cellcentered/stencils.hh @@ -20,8 +20,8 @@ * \file * \brief Implements the notion of stencils for cell-centered models */ -#ifndef DUMUX_CC_STENCILS_HH -#define DUMUX_CC_STENCILS_HH +#ifndef DUMUX_DISCRETIZATION_CC_STENCILS_HH +#define DUMUX_DISCRETIZATION_CC_STENCILS_HH #include <set> #include <dumux/implicit/cellcentered/properties.hh> diff --git a/dumux/discretization/cellcentered/tpfa/CMakeLists.txt b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt new file mode 100644 index 0000000000..ff7ebeb2ba --- /dev/null +++ b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt @@ -0,0 +1,5 @@ +install(FILES +fvelementgeometryvector.hh +subcontrolvolumeface.hh +subcontrolvolume.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered/tpfa) \ No newline at end of file diff --git a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh similarity index 60% rename from dumux/porousmediumflow/constitutivelaws/darcyslaw.hh rename to dumux/discretization/cellcentered/tpfa/darcyslaw.hh index f5f6fb07fd..c24d224618 100644 --- a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh +++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh @@ -22,8 +22,8 @@ * 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_POROUSMEDIUMFLOW_DARCYS_LAW_HH -#define DUMUX_POROUSMEDIUMFLOW_DARCYS_LAW_HH +#ifndef DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH +#define DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH #include <memory> @@ -47,15 +47,8 @@ NEW_PROP_TAG(ProblemEnableGravity); /*! * \ingroup DarcysLaw - * \brief Evaluates the normal component of the Darcy velocity - * on a (sub)control volume face. Specializations are provided - * for the different discretization methods. + * \brief Specialization of Darcy's Law for the CCTpfa method. */ -template <class TypeTag, typename DiscretizationMethod = void> -class DarcysLaw -{}; - -// Specialization for the CC-Tpfa method template <class TypeTag> class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type > { @@ -198,116 +191,6 @@ private: } }; -// Specialization 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 diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh new file mode 100644 index 0000000000..cddca62731 --- /dev/null +++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh @@ -0,0 +1,178 @@ +// -*- 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 + * diffusive mass fluxes due to molecular diffusion with Fick's law. + */ +#ifndef DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_HH +#define DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_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(NumPhases); +NEW_PROP_TAG(FluidSystem); +NEW_PROP_TAG(EffectiveDiffusivityModel); +} + +/*! + * \ingroup CCTpfaFicksLaw + * \brief Specialization of Fick's Law for the CCTpfa method. + */ +template <class TypeTag> +class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type > +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel; + typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; + typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GridView::IndexSet::IndexType IndexType; + typedef typename std::vector<IndexType> Stencil; + + using Element = typename GridView::template Codim<0>::Entity; + + enum { dim = GridView::dimension} ; + enum { dimWorld = GridView::dimensionworld} ; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ; + + typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix; + typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + +public: + + static Scalar flux(const Problem& problem, + const Element& element, + const SubControlVolumeFace& scvFace, + const int phaseIdx, + const int compIdx) + { + // diffusion tensors are always solution dependent + Scalar tij = calculateTransmissibility_(problem, element, scvFace, phaseIdx, compIdx); + + // 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 auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx()); + + // 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)); + + return rho*tij*(xInside - xOutside); + } + + static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + { + std::vector<IndexType> stencil; + stencil.clear(); + if (!scvFace.boundary()) + { + stencil.push_back(scvFace.insideScvIdx()); + stencil.push_back(scvFace.outsideScvIdx()); + } + else + stencil.push_back(scvFace.insideScvIdx()); + + return stencil; + } + +private: + + + static Scalar calculateTransmissibility_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx) + { + Scalar tij; + + 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_(problem, element, scvFace, 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_(problem, element, scvFace, outsideD, outsideScv); + + tij = scvFace.area()*(ti * tj)/(ti + tj); + } + else + { + tij = scvFace.area()*ti; + } + + return tij; + } + + static Scalar calculateOmega_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv) + { + GlobalPosition Dnormal; + D.mv(scvFace.unitOuterNormal(), Dnormal); + + auto distanceVector = scvFace.center(); + distanceVector -= scv.center(); + distanceVector /= distanceVector.two_norm2(); + + Scalar omega = Dnormal * distanceVector; + omega *= problem.boxExtrusionFactor(element, scv); + + return omega; + } + + static Scalar calculateOmega_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv) + { + auto distanceVector = scvFace.center(); + distanceVector -= scv.center(); + distanceVector /= distanceVector.two_norm2(); + + Scalar omega = D * (distanceVector * scvFace.unitOuterNormal()); + omega *= problem.boxExtrusionFactor(element, scv); + + return omega; + } +}; +} // end namespace + +#endif diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh similarity index 93% rename from dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh rename to dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh index 1320104471..1a239b4e38 100644 --- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh +++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh @@ -22,12 +22,10 @@ * This builds up the sub control volumes and sub control volume faces * for each element. */ -#ifndef DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH -#define DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH +#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH +#define DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH -#include <dumux/implicit/subcontrolvolume.hh> -#include <dumux/implicit/subcontrolvolumeface.hh> -#include <dumux/implicit/fvelementgeometry.hh> +#include <dumux/implicit/cellcentered/tpfa/properties.hh> #include <dumux/common/elementmap.hh> namespace Dumux @@ -286,21 +284,21 @@ public: IndexType numLocalFaces = element.subEntities(1); std::vector<IndexType> scvfsIndexSet(numLocalFaces); std::vector<IndexType> neighborVolVarIndexSet(numLocalFaces); + IndexType localFaceIdx = 0; for (const auto& intersection : intersections(gridView_, element)) { - IndexType localFaceIdx = intersection.indexInInside(); // inner sub control volume faces if (intersection.neighbor()) { scvfsIndexSet[localFaceIdx] = numScvf_++; auto nIdx = problem.elementMapper().index(intersection.outside()); - neighborVolVarIndexSet[localFaceIdx] = nIdx; + neighborVolVarIndexSet[localFaceIdx++] = nIdx; } // boundary sub control volume faces else if (intersection.boundary()) { scvfsIndexSet[localFaceIdx] = numScvf_++; - neighborVolVarIndexSet[localFaceIdx] = numScvs_ + numBoundaryScvf_++; + neighborVolVarIndexSet[localFaceIdx++] = numScvs_ + numBoundaryScvf_++; } } @@ -361,12 +359,12 @@ private: if (intersection.neighbor() || intersection.boundary()) { localScvfs_.push_back(SubControlVolumeFace(intersection.geometry(), - intersection.geometry().center(), - intersection.centerUnitOuterNormal(), - scvFaceIndices[scvfCounter], - intersection.indexInInside(), - std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}), - intersection.boundary())); + intersection.geometry().center(), + intersection.centerUnitOuterNormal(), + scvFaceIndices[scvfCounter], + intersection.indexInInside(), + std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}), + intersection.boundary())); localScvfIndices_.push_back(scvFaceIndices[scvfCounter]); scvfCounter++; diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh new file mode 100644 index 0000000000..aaba6ee555 --- /dev/null +++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh @@ -0,0 +1,73 @@ +// -*- 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_CC_TPFA_SUBCONTROLVOLUME_HH +#define DUMUX_DISCRETIZATION_CC_TPFA_SUBCONTROLVOLUME_HH + +#include <dune/common/fvector.hh> +#include <dumux/discretization/subcontrolvolumebase.hh> + +namespace Dumux +{ +template<class G, typename I> +class CCTpfaSubControlVolume : 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 cc case + CCTpfaSubControlVolume(Geometry geometry, + IndexType elementIndex) + : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)) {} + + IndexType indexInElement() const + { + return IndexType(0); + } + + //! The global index of this scv + IndexType index() const + { + return this->elementIndex(); + } + + IndexType dofIndex() const + { + return this->elementIndex(); + } + + GlobalPosition dofPosition() const + { + return this->geometry().center(); + } +}; +} // end namespace + +#endif diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh new file mode 100644 index 0000000000..9a3d763c5c --- /dev/null +++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh @@ -0,0 +1,62 @@ +// -*- 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_CC_TPFA_SUBCONTROLVOLUMEFACE_HH +#define DUMUX_DISCRETIZATION_CC_TPFA_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 CCTpfaSubControlVolumeFace : 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: + CCTpfaSubControlVolumeFace(const Geometry& geometry, + const GlobalPosition& ipGlobal, + const GlobalPosition& unitOuterNormal, + IndexType scvfIndex, + IndexType indexInElement, + const std::vector<IndexType>& scvIndices, + bool boundary = false) + : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary) + {} +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/cellcentered/volumevariablesvector.hh b/dumux/discretization/cellcentered/volumevariablesvector.hh similarity index 99% rename from dumux/implicit/cellcentered/volumevariablesvector.hh rename to dumux/discretization/cellcentered/volumevariablesvector.hh index bac9996b57..85534726f1 100644 --- a/dumux/implicit/cellcentered/volumevariablesvector.hh +++ b/dumux/discretization/cellcentered/volumevariablesvector.hh @@ -20,8 +20,8 @@ * \file * \brief Base class for the volume variables vector */ -#ifndef DUMUX_IMPLICIT_CC_VOLVARSVECTOR_HH -#define DUMUX_IMPLICIT_CC_VOLVARSVECTOR_HH +#ifndef DUMUX_DISCRETIZATION_CC_VOLVARSVECTOR_HH +#define DUMUX_DISCRETIZATION_CC_VOLVARSVECTOR_HH #include <dumux/implicit/properties.hh> diff --git a/dumux/discretization/darcyslaw.hh b/dumux/discretization/darcyslaw.hh new file mode 100644 index 0000000000..14e35ee628 --- /dev/null +++ b/dumux/discretization/darcyslaw.hh @@ -0,0 +1,47 @@ +// -*- 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_DARCYS_LAW_HH +#define DUMUX_DISCRETIZATION_DARCYS_LAW_HH + +namespace Dumux +{ + +/*! + * \ingroup DarcysLaw + * \brief Evaluates the normal component of the Darcy velocity + * on a (sub)control volume face. Specializations are provided + * for the different discretization methods. These specializations + * are found in the headers included below. + */ +template <class TypeTag, typename DiscretizationMethod = void> +class DarcysLaw +{}; + +} // end namespace + +#include <dumux/discretization/box/darcyslaw.hh> +#include <dumux/discretization/cellcentered/tpfa/darcyslaw.hh> + +#endif diff --git a/dumux/discretization/fickslaw.hh b/dumux/discretization/fickslaw.hh new file mode 100644 index 0000000000..a1de3ebc2d --- /dev/null +++ b/dumux/discretization/fickslaw.hh @@ -0,0 +1,59 @@ +// -*- 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 + * diffusive mass fluxes due to molecular diffusion with Fick's law. + */ +#ifndef DUMUX_DISCRETIZATION_FICKS_LAW_HH +#define DUMUX_DISCRETIZATION_FICKS_LAW_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(NumPhases); +NEW_PROP_TAG(FluidSystem); +NEW_PROP_TAG(EffectiveDiffusivityModel); +} + +/*! + * \ingroup CCTpfaFicksLaw + * \brief Evaluates the diffusive mass flux according to Fick's law + */ +template <class TypeTag, typename DiscretizationMethod = void> +class FicksLaw +{}; + +} // end namespace + +#include <dumux/discretization/cellcentered/tpfa/fickslaw.hh> +#include <dumux/discretization/box/fickslaw.hh> + +#endif diff --git a/dumux/implicit/fluxvariablesbase.hh b/dumux/discretization/fluxvariablesbase.hh similarity index 97% rename from dumux/implicit/fluxvariablesbase.hh rename to dumux/discretization/fluxvariablesbase.hh index 6fe0a78696..de8061254d 100644 --- a/dumux/implicit/fluxvariablesbase.hh +++ b/dumux/discretization/fluxvariablesbase.hh @@ -20,8 +20,8 @@ * \file * \brief Base class for the flux variables */ -#ifndef DUMUX_IMPLICIT_FLUXVARIABLESBASE_HH -#define DUMUX_IMPLICIT_FLUXVARIABLESBASE_HH +#ifndef DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH +#define DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH #include <dumux/implicit/properties.hh> @@ -29,7 +29,7 @@ namespace Dumux { /*! - * \ingroup ImplicitModel + * \ingroup Discretization * \brief Base class for the flux variables * Actual flux variables inherit from this class */ diff --git a/dumux/implicit/fvelementgeometry.hh b/dumux/discretization/fvelementgeometry.hh similarity index 97% rename from dumux/implicit/fvelementgeometry.hh rename to dumux/discretization/fvelementgeometry.hh index df59ef883b..b2b5cf8e9b 100644 --- a/dumux/implicit/fvelementgeometry.hh +++ b/dumux/discretization/fvelementgeometry.hh @@ -21,15 +21,12 @@ * \brief Class providing iterators over sub control volumes and sub control * volume faces of an element. */ -#ifndef DUMUX_FV_ELEMENTGEOMETRY_HH -#define DUMUX_FV_ELEMENTGEOMETRY_HH +#ifndef DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH +#define DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH #include <dune/common/iteratorrange.hh> #include <dune/common/iteratorfacades.hh> -#include <dumux/implicit/subcontrolvolume.hh> -#include <dumux/implicit/subcontrolvolumeface.hh> - namespace Dumux { namespace Properties @@ -41,7 +38,7 @@ NEW_PROP_TAG(FVElementGeometryVector); } /*! - * \ingroup ImplcititModel + * \ingroup Discretization * \brief An iterator over sub control volumes */ template<class SubControlVolume, class Vector, class FVElementGeometryVector> diff --git a/dumux/implicit/subcontrolvolume.hh b/dumux/discretization/subcontrolvolumebase.hh similarity index 100% rename from dumux/implicit/subcontrolvolume.hh rename to dumux/discretization/subcontrolvolumebase.hh diff --git a/dumux/implicit/subcontrolvolumeface.hh b/dumux/discretization/subcontrolvolumefacebase.hh similarity index 87% rename from dumux/implicit/subcontrolvolumeface.hh rename to dumux/discretization/subcontrolvolumefacebase.hh index f0f964e81d..ab732275d7 100644 --- a/dumux/implicit/subcontrolvolumeface.hh +++ b/dumux/discretization/subcontrolvolumefacebase.hh @@ -20,8 +20,8 @@ * \file * \brief Base class for a sub control volume face */ -#ifndef DUMUX_SUBCONTROLVOLUMEFACE_HH -#define DUMUX_SUBCONTROLVOLUMEFACE_HH +#ifndef DUMUX_DISCRETIZATION_SUBCONTROLVOLUMEFACEBASE_HH +#define DUMUX_DISCRETIZATION_SUBCONTROLVOLUMEFACEBASE_HH #include <utility> #include <dune/common/fvector.hh> @@ -30,12 +30,12 @@ namespace Dumux { /*! - * \ingroup ImplicitModel + * \ingroup Discretization * \brief Base class for a sub control volume face, i.e a part of the boundary * of a sub control volume we computing a flux on. */ template<class Geometry, typename IndexType> -class SubControlVolumeFace +class SubControlVolumeFaceBase { using Scalar = typename Geometry::ctype; static const int dim = Geometry::mydimension; @@ -45,13 +45,13 @@ class SubControlVolumeFace using LocalPosition = Dune::FieldVector<Scalar, dim>; public: - SubControlVolumeFace(const Geometry& geometry, - const GlobalPosition& ipGlobal, - const GlobalPosition& unitOuterNormal, - IndexType scvfIndex, - IndexType indexInElement, - const std::vector<IndexType>& scvIndices, - bool boundary = false) + SubControlVolumeFaceBase(const Geometry& geometry, + const GlobalPosition& ipGlobal, + const GlobalPosition& unitOuterNormal, + IndexType scvfIndex, + IndexType indexInElement, + const std::vector<IndexType>& scvIndices, + bool boundary = false) : geometry_(geometry), ipGlobal_(ipGlobal), ipLocal_(geometry.local(ipGlobal)), @@ -142,4 +142,5 @@ private: } // end namespace + #endif diff --git a/dumux/implicit/volumevariables.hh b/dumux/discretization/volumevariables.hh similarity index 97% rename from dumux/implicit/volumevariables.hh rename to dumux/discretization/volumevariables.hh index fe3b89b594..7acd2069a9 100644 --- a/dumux/implicit/volumevariables.hh +++ b/dumux/discretization/volumevariables.hh @@ -21,10 +21,10 @@ * \brief Base class for the model specific class which provides * access to all volume averaged quantities. */ -#ifndef DUMUX_IMPLICIT_VOLUME_VARIABLES_HH -#define DUMUX_IMPLICIT_VOLUME_VARIABLES_HH +#ifndef DUMUX_DISCRETIZATION_VOLUME_VARIABLES_HH +#define DUMUX_DISCRETIZATION_VOLUME_VARIABLES_HH -#include "properties.hh" +#include <dumux/implicit/properties.hh> #include <dumux/common/valgrind.hh> diff --git a/dumux/implicit/box/propertydefaults.hh b/dumux/implicit/box/propertydefaults.hh index 1eb743f5e8..3681ef036c 100644 --- a/dumux/implicit/box/propertydefaults.hh +++ b/dumux/implicit/box/propertydefaults.hh @@ -27,18 +27,20 @@ #define DUMUX_BOX_PROPERTY_DEFAULTS_HH #include <dumux/implicit/propertydefaults.hh> -#include <dumux/implicit/fvelementgeometry.hh> +#include <dumux/discretization/fvelementgeometry.hh> +#include <dumux/discretization/box/subcontrolvolume.hh> +#include <dumux/discretization/box/subcontrolvolumeface.hh> +#include <dumux/discretization/box/fluxvariablescachevector.hh> +#include <dumux/discretization/box/volumevariablesvector.hh> +#include <dumux/discretization/box/fvelementgeometryvector.hh> +#include <dumux/discretization/box/stencils.hh> #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh> -#include "fluxvariablescachevector.hh" -#include "volumevariablesvector.hh" #include "elementboundarytypes.hh" -#include "fvelementgeometryvector.hh" #include "localresidual.hh" #include "localjacobian.hh" #include "assembler.hh" #include "properties.hh" -#include "stencils.hh" namespace Dumux { @@ -69,7 +71,7 @@ private: using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>; using IndexType = typename GridView::IndexSet::IndexType; public: - using type = SubControlVolume<ScvGeometry, IndexType, /*isBox=*/true>; + using type = BoxSubControlVolume<ScvGeometry, IndexType, /*isBox=*/true>; }; SET_PROP(BoxModel, SubControlVolumeFace) @@ -82,7 +84,7 @@ private: using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>; using IndexType = typename GridView::IndexSet::IndexType; public: - using type = SubControlVolumeFace<ScvfGeometry, IndexType>; + using type = BoxSubControlVolumeFace<ScvfGeometry, IndexType>; }; //! Set the default for the ElementBoundaryTypes @@ -95,22 +97,22 @@ SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper) SET_TYPE_PROP(BoxModel, StencilsVector, BoxStencilsVector<TypeTag>); //! The global current volume variables vector class -SET_TYPE_PROP(BoxModel, CurrentVolumeVariablesVector, Dumux::BoxVolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>); +SET_TYPE_PROP(BoxModel, CurrentVolumeVariablesVector, BoxVolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>); //! The global previous volume variables vector class -SET_TYPE_PROP(BoxModel, PreviousVolumeVariablesVector, Dumux::BoxVolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>); +SET_TYPE_PROP(BoxModel, PreviousVolumeVariablesVector, BoxVolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>); //! The global flux variables cache vector class -SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, Dumux::BoxFluxVariablesCacheVector<TypeTag>); +SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, BoxFluxVariablesCacheVector<TypeTag>); //! Set the BaseLocalResidual to BoxLocalResidual SET_TYPE_PROP(BoxModel, BaseLocalResidual, BoxLocalResidual<TypeTag>); //! Assembler for the global jacobian matrix -SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::BoxAssembler<TypeTag>); +SET_TYPE_PROP(BoxModel, JacobianAssembler, BoxAssembler<TypeTag>); //! The local jacobian operator -SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>); +SET_TYPE_PROP(BoxModel, LocalJacobian, BoxLocalJacobian<TypeTag>); //! indicate that this is a box discretization SET_BOOL_PROP(BoxModel, ImplicitIsBox, true); diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh index 72cfb72f45..ef94323771 100644 --- a/dumux/implicit/cellcentered/propertydefaults.hh +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -28,13 +28,13 @@ #define DUMUX_CC_PROPERTY_DEFAULTS_HH #include <dumux/implicit/propertydefaults.hh> +#include <dumux/discretization/cellcentered/fluxvariablescachevector.hh> +#include <dumux/discretization/cellcentered/volumevariablesvector.hh> +#include <dumux/discretization/cellcentered/stencils.hh> -#include "fluxvariablescachevector.hh" -#include "volumevariablesvector.hh" #include "elementboundarytypes.hh" #include "localresidual.hh" #include "properties.hh" -#include "stencils.hh" #include "localjacobian.hh" #include "assembler.hh" diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index 48af42db4d..cc904e2b56 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -28,9 +28,10 @@ #define DUMUX_CCTPFA_PROPERTY_DEFAULTS_HH #include <dumux/implicit/propertydefaults.hh> -#include <dumux/implicit/fvelementgeometry.hh> #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh> -#include <dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh> +#include <dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh> +#include <dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh> +#include <dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh> #include <dumux/implicit/cellcentered/properties.hh> namespace Dumux { @@ -53,7 +54,7 @@ private: using ScvGeometry = typename Grid::template Codim<0>::Geometry; using IndexType = typename Grid::LeafGridView::IndexSet::IndexType; public: - typedef Dumux::SubControlVolume<ScvGeometry, IndexType, /*isBox=*/false> type; + typedef Dumux::CCTpfaSubControlVolume<ScvGeometry, IndexType> type; }; SET_PROP(CCTpfaModel, SubControlVolumeFace) @@ -63,7 +64,7 @@ private: using ScvfGeometry = typename Grid::template Codim<1>::Geometry; using IndexType = typename Grid::LeafGridView::IndexSet::IndexType; public: - typedef Dumux::SubControlVolumeFace<ScvfGeometry, IndexType> type; + typedef Dumux::CCTpfaSubControlVolumeFace<ScvfGeometry, IndexType> type; }; } // namespace Properties diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 6fa48bf480..2d4219e401 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -38,15 +38,15 @@ #include <dumux/porousmediumflow/implicit/fluxvariables.hh> #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh> -#include <dumux/porousmediumflow/constitutivelaws/darcyslaw.hh> -#include <dumux/porousmediumflow/constitutivelaws/fickslaw.hh> +#include <dumux/discretization/volumevariables.hh> +#include <dumux/discretization/fvelementgeometry.hh> +#include <dumux/discretization/darcyslaw.hh> +#include <dumux/discretization/fickslaw.hh> #include "properties.hh" #include "model.hh" #include "assembler.hh" #include "localjacobian.hh" -#include "volumevariables.hh" -#include "fvelementgeometry.hh" namespace Dumux { diff --git a/dumux/porousmediumflow/1p/implicit/volumevariables.hh b/dumux/porousmediumflow/1p/implicit/volumevariables.hh index f834115929..d888b631cd 100644 --- a/dumux/porousmediumflow/1p/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/1p/implicit/volumevariables.hh @@ -25,7 +25,7 @@ #define DUMUX_1P_VOLUME_VARIABLES_HH #include "properties.hh" -#include <dumux/implicit/volumevariables.hh> +#include <dumux/discretization/volumevariables.hh> #include <dumux/material/fluidstates/immiscible.hh> diff --git a/dumux/porousmediumflow/2p/implicit/volumevariables.hh b/dumux/porousmediumflow/2p/implicit/volumevariables.hh index 1f698eb780..5e001fc651 100644 --- a/dumux/porousmediumflow/2p/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/2p/implicit/volumevariables.hh @@ -27,7 +27,7 @@ #include "properties.hh" -#include <dumux/implicit/volumevariables.hh> +#include <dumux/discretization/volumevariables.hh> #include <dune/common/fvector.hh> diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt index e48629a5b2..de81dec4ff 100644 --- a/dumux/porousmediumflow/CMakeLists.txt +++ b/dumux/porousmediumflow/CMakeLists.txt @@ -11,10 +11,9 @@ add_subdirectory("3p") add_subdirectory("3p3c") add_subdirectory("co2") add_subdirectory("compositional") -add_subdirectory("constitutivelaws") add_subdirectory("immiscible") add_subdirectory("implicit") add_subdirectory("mpnc") add_subdirectory("nonisothermal") add_subdirectory("richards") -add_subdirectory("sequential") \ No newline at end of file +add_subdirectory("sequential") diff --git a/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt b/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt deleted file mode 100644 index e5e8f670be..0000000000 --- a/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt +++ /dev/null @@ -1,5 +0,0 @@ -#install headers -install(FILES -darcyslaw.hh -fickslaw.hh -DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/constitutivelaws) \ No newline at end of file diff --git a/dumux/porousmediumflow/constitutivelaws/fickslaw.hh b/dumux/porousmediumflow/constitutivelaws/fickslaw.hh deleted file mode 100644 index a1829d1898..0000000000 --- a/dumux/porousmediumflow/constitutivelaws/fickslaw.hh +++ /dev/null @@ -1,379 +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 - * diffusive mass fluxes due to molecular diffusion with Fick's law. - */ -#ifndef DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_HH -#define DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_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(NumPhases); -NEW_PROP_TAG(FluidSystem); -NEW_PROP_TAG(EffectiveDiffusivityModel); -} - -/*! - * \ingroup CCTpfaFicksLaw - * \brief Evaluates the diffusive mass flux according to Fick's law - */ -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 > -{ - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GridView::IndexSet::IndexType IndexType; - typedef typename std::vector<IndexType> Stencil; - - using Element = typename GridView::template Codim<0>::Entity; - - enum { dim = GridView::dimension} ; - enum { dimWorld = GridView::dimensionworld} ; - enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ; - - typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix; - typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - -public: - - static Scalar flux(const Problem& problem, - const Element& element, - const SubControlVolumeFace& scvFace, - const int phaseIdx, - const int compIdx) - { - // diffusion tensors are always solution dependent - Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx); - - // 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 auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx()); - - // 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)); - - return rho*tij*(xInside - xOutside); - } - - static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace) - { - std::vector<IndexType> stencil; - stencil.clear(); - if (!scvFace.boundary()) - { - stencil.push_back(scvFace.insideScvIdx()); - stencil.push_back(scvFace.outsideScvIdx()); - } - else - stencil.push_back(scvFace.insideScvIdx()); - - return stencil; - } - -private: - - - static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx) - { - Scalar tij; - - 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_(problem, scvFace, 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_(problem, scvFace, outsideD, outsideScv); - - tij = scvFace.area()*(ti * tj)/(ti + tj); - } - else - { - tij = scvFace.area()*ti; - } - - return tij; - } - - static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv) - { - 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; - } - - static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv) - { - 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; - } -}; - -// 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 diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index 3038322ea0..a528cd9872 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -24,7 +24,7 @@ #define DUMUX_POROUSMEDIUMFLOW_IMPLICIT_FLUXVARIABLES_HH #include <dumux/implicit/properties.hh> -#include <dumux/implicit/fluxvariablesbase.hh> +#include <dumux/discretization/fluxvariablesbase.hh> namespace Dumux { -- GitLab