From 0a7bbb63d0eda97fdb6440de41b18b0f35fb7b5e Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Fri, 26 Jun 2020 09:57:44 +0200 Subject: [PATCH 01/36] [common] Add class for the decomposition of vectors into basis --- dumux/common/vectordecomposition.hh | 203 ++++++++++++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 dumux/common/vectordecomposition.hh diff --git a/dumux/common/vectordecomposition.hh b/dumux/common/vectordecomposition.hh new file mode 100644 index 0000000000..0254c55210 --- /dev/null +++ b/dumux/common/vectordecomposition.hh @@ -0,0 +1,203 @@ +// -*- 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 3 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 . * + *****************************************************************************/ + +#ifndef DUMUX_COMMON_VECTORDECOMPOSITION_HH +#define DUMUX_COMMON_VECTORDECOMPOSITION_HH + +#include + +#include + +#include "math.hh" + +namespace Dumux { + +class VectorDecomposition +{ + +public: + template + static std::tuple, std::vector, bool> calculateVectorDecomposition(const DV x, const std::vector& v) + { + std::size_t numVectors = v.size(); + static constexpr auto dim = DV::dimension; + if(numVectors == 0) + DUNE_THROW(Dune::InvalidStateException, "Can't perform decomposition without any given vectors!"); + + if(isAligned(x, v[0])) + return {std::vector({0}), calculateCoefficients(x, v[0]), true}; + + using Scalar = typename DV::value_type; + std::vector coefficients(dim, std::numeric_limits::max()); + std::vector indices; + bool foundValidSolution = false; + + if constexpr (dim == 3) + { + for (std::size_t i = 0; i < numVectors - 2; ++i) + { + for (std::size_t j = i+1; j < numVectors - 1; ++j) + { + for (std::size_t k = j+1; k < numVectors; ++k) + { + auto [coeff, valid] = calculateCoefficients(x, v[i], v[j], v[k]); + using std::max_element; using std::move; + if(valid && (*max_element(coeff.begin(), coeff.end()) < *max_element(coefficients.begin(), coefficients.end()))) + { + coefficients = move(coeff); + indices = {i, j, k}; + foundValidSolution = true; + } + } + } + } + } + else if constexpr (dim == 2) + { + for (std::size_t i = 0; i < numVectors - 1; ++i) + { + for (std::size_t j = i+1; j < numVectors; ++j) + { + auto [coeff, valid] = calculateCoefficients(x, v[i], v[j]); + using std::max_element; using std::move; + if(valid && (*max_element(coeff.begin(), coeff.end()) < *max_element(coefficients.begin(), coefficients.end()))) + { + coefficients = move(coeff); + indices = {i, j}; + foundValidSolution = true; + } + } + } + } + + return {indices, coefficients, foundValidSolution}; + } + +private: + //Cramer's rule for the solution of 2D system of equation + template + static std::pair, bool> calculateCoefficients(DV x, DV v1, DV v2) + { + std::vector coeff; + const auto v1_norm = v1.two_norm(); + const auto v2_norm = v2.two_norm(); + + v1 /= v1_norm; + v2 /= v2_norm; + const auto xNorm = x.two_norm(); + + x /= xNorm; + + const auto A_f = v1[0]*v2[1] - v1[1]*v2[0]; + const auto A_f_1 = x[0]*v2[1] - x[1]*v2[0]; + const auto A_f_2 = v1[0]*x[1] - v1[1]*x[0]; + + if(std::abs(A_f) < 1.0e-8) + return {coeff, false}; + + coeff.resize(2); + coeff[0] = A_f_1 / A_f ; + coeff[1] = A_f_2 / A_f ; + + for(auto& c : coeff) + if(std::abs(c) < 1.0e-12) + c = 0.0; + + coeff[0] *= xNorm/v1_norm; + coeff[1] *= xNorm/v2_norm; + + if(coeff[0] >= -1.0e-30 && coeff[1] >= -1.0e-30) + return {coeff, true}; + + return {coeff, false}; + } + + template + static std::pair, bool> calculateCoefficients(DV x, DV v1, DV v2, DV v3) + { + std::vector coeff; + const auto v1_norm = v1.two_norm(); + const auto v2_norm = v2.two_norm(); + const auto v3_norm = v3.two_norm(); + const auto xNorm = x.two_norm(); + + v1 /= v1_norm; + v2 /= v2_norm; + v3 /= v3_norm; + + x /= xNorm; + + const auto A_f = v1*crossProduct(v2,v3); + const auto A_f_1 = x * crossProduct(v2,v3); + const auto A_f_2 = v1 * crossProduct(x,v3); + const auto A_f_3 = v1 * crossProduct(v2,x); + + if(std::abs(A_f) < 1.0e-8) + return {coeff, false}; + + coeff.resize(3); + coeff[0] = A_f_1 / A_f ; + coeff[1] = A_f_2 / A_f ; + coeff[2] = A_f_3 / A_f ; + + for(auto& c : coeff) + if(std::abs(c) < 1.0e-12) + c = 0.0; + + coeff[0] *= xNorm/v1_norm; + coeff[1] *= xNorm/v2_norm; + coeff[2] *= xNorm/v3_norm; + + if(coeff[0] >=-1.0e-30 && coeff[1] >=-1.0e-30 && coeff[2] >=-1.0e-30) + return {coeff, true}; + + return {coeff, false}; + } + + template + static bool isAligned(DV x, DV v1) + { + const auto v1_norm = v1.two_norm(); + + v1 /= v1_norm; + const auto xNorm = x.two_norm(); + + x /= xNorm; + + if(std::abs(v1*x - 1.0) < 1.0e-10) + return true; + + return false; + } + + template + static std::vector calculateCoefficients(const DV& x, const DV& v1) + { + std::vector coeff; + coeff.resize(1); + coeff[0] = std::abs(x*v1); + coeff[0] /= v1.two_norm2(); + + return coeff; + } +}; + +} // end namespace Dumux + +#endif -- GitLab From f5a43b541e089398a11244ea0b7307d77e88ee6b Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:07:22 +0200 Subject: [PATCH 02/36] [disc][wmpfa] Add draft of gridgeometry and elementgeometry --- .../cellcentered/wmpfa/fvelementgeometry.hh | 175 +++++++++ .../cellcentered/wmpfa/fvgridgeometry.hh | 331 ++++++++++++++++++ 2 files changed, 506 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh create mode 100644 dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh diff --git a/dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh new file mode 100644 index 0000000000..4b6a2e9067 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/fvelementgeometry.hh @@ -0,0 +1,175 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models + * This builds up the sub control volumes and sub control volume faces + * for each element in the local scope we are restricting to, e.g. stencil or element. + */ +#ifndef DUMUX_DISCRETIZATION_CCWMPFA_FV_ELEMENT_GEOMETRY_HH +#define DUMUX_DISCRETIZATION_CCWMPFA_FV_ELEMENT_GEOMETRY_HH + +#include +#include +#include + +#include +#include +#include +#include + +namespace Dumux { + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models + * This builds up the sub control volumes and sub control volume faces + * for each element in the local scope we are restricting to, e.g. stencil or element. + * \tparam GG the finite volume grid geometry type + * \tparam enableGridGeometryCache if the grid geometry is cached or not + * \note This class is specialized for versions with and without caching the fv geometries on the grid view + */ +template +class CCWMpfaFVElementGeometry; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models + * Specialization for grid caching enabled + * \note The finite volume geometries are stored in the corresponding GridGeometry + */ +template +class CCWMpfaFVElementGeometry +{ + using ThisType = CCWMpfaFVElementGeometry; + using GridView = typename GG::GridView; + using GridIndexType = typename IndexTraits::GridIndex; + +public: + //! export type of the element + using Element = typename GridView::template Codim<0>::Entity; + //! export type of subcontrol volume + using SubControlVolume = typename GG::SubControlVolume; + //! export type of subcontrol volume face + using SubControlVolumeFace = typename GG::SubControlVolumeFace; + //! export type of finite volume grid geometry + using GridGeometry = GG; + + //! the maximum number of scvs per element + static constexpr std::size_t maxNumElementScvs = 1; + //! the maximum number of scvfs per element (use cubes for maximum) + static constexpr std::size_t maxNumElementScvfs = 2*GridView::dimension; + + //! Constructor + CCWMpfaFVElementGeometry(const GridGeometry& gridGeometry) + : gridGeometryPtr_(&gridGeometry) {} + + //! Get an elment sub control volume with a global scv index + //! We separate element and neighbor scvs to speed up mapping + const SubControlVolume& scv(GridIndexType scvIdx) const + { + return gridGeometry().scv(scvIdx); + } + + //! Get an element sub control volume face with a global scvf index + //! We separate element and neighbor scvfs to speed up mapping + const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const + { + return gridGeometry().scvf(scvfIdx); + } + + //! Get the scvf on the same face but from the other side + //! Note that e.g. the normals might be different in the case of surface grids + const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const + { + return gridGeometry().flipScvf(scvfIdx, outsideScvIdx); + } + + //! iterator range for sub control volumes. Iterates over + //! all scvs of the bound element (not including neighbor scvs) + //! This is a free function found by means of ADL + //! To iterate over all sub control volumes of this FVElementGeometry use + //! for (auto&& scv : scvs(fvGeometry)) + friend inline Dune::IteratorRange< ScvIterator, ThisType> > + scvs(const CCWMpfaFVElementGeometry& fvGeometry) + { + using ScvIterator = Dumux::ScvIterator, ThisType>; + return Dune::IteratorRange(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry), + ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry)); + } + + //! iterator range for sub control volumes faces. Iterates over + //! all scvfs of the bound element (not including neighbor scvfs) + //! This is a free function found by means of ADL + //! To iterate over all sub control volume faces of this FVElementGeometry use + //! for (auto&& scvf : scvfs(fvGeometry)) + friend inline Dune::IteratorRange< ScvfIterator, ThisType> > + scvfs(const CCWMpfaFVElementGeometry& fvGeometry) + { + const auto& g = fvGeometry.gridGeometry(); + const auto scvIdx = fvGeometry.scvIndices_[0]; + using ScvfIterator = Dumux::ScvfIterator, ThisType>; + return Dune::IteratorRange(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry), + ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry)); + } + + //! number of sub control volumes in this fv element geometry + std::size_t numScv() const + { + return scvIndices_.size(); + } + + //! number of sub control volumes in this fv element geometry + std::size_t numScvf() const + { + return gridGeometry().scvfIndicesOfScv(scvIndices_[0]).size(); + } + + //! Binding of an element, called by the local jacobian to prepare element assembly + void bind(const Element& element) + { + this->bindElement(element); + } + + //! Bind only element-local + void bindElement(const Element& element) + { + elementPtr_ = &element; + scvIndices_[0] = gridGeometry().elementMapper().index(*elementPtr_); + } + + //! The global finite volume geometry we are a restriction of + const GridGeometry& gridGeometry() const + { return *gridGeometryPtr_; } + + //! Returns whether one of the geometry's scvfs lies on a boundary + bool hasBoundaryScvf() const + { return gridGeometry().hasBoundaryScvf(scvIndices_[0]); } + +private: + + const Element* elementPtr_; + std::array scvIndices_; + const GridGeometry* gridGeometryPtr_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh new file mode 100644 index 0000000000..bde9384594 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh @@ -0,0 +1,331 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered TPFA models on a grid view + * This builds up the sub control volumes and sub control volume faces + * for each element of the grid partition. + */ +#ifndef DUMUX_DISCRETIZATION_CCWMpfa_FV_GRID_GEOMETRY_HH +#define DUMUX_DISCRETIZATION_CCWMpfa_FV_GRID_GEOMETRY_HH + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace Dumux { + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The default traits for the tpfa finite volume grid geometry + * Defines the scv and scvf types and the mapper types + * \tparam the grid view type + */ +template> +struct CCWMpfaDefaultGridGeometryTraits +: public MapperTraits +{ + using SubControlVolume = CCSubControlVolume; + using SubControlVolumeFace = CCWMpfaSubControlVolumeFace; + + template + using ConnectivityMap = CCSimpleConnectivityMap; + + template + using LocalView = CCWMpfaFVElementGeometry; + + //! State the maximum admissible number of neighbors per scvf + //! Per default, we allow for 8 branches on network/surface grids, where + //! conformity is assumed. For normal grids, we allow a maximum of one + //! hanging node per scvf. Use different traits if you need more. + static constexpr int maxNumScvfNeighbors = int(GridView::dimension) > +class CCWMpfaFVGridGeometry; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered weighted MPFA models on a grid view + * This builds up the sub control volumes and sub control volume faces + * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster + */ +template +class CCWMpfaFVGridGeometry +: public BaseGridGeometry +{ + using ThisType = CCWMpfaFVGridGeometry; + using ParentType = BaseGridGeometry; + using ConnectivityMap = typename Traits::template ConnectivityMap; + using GridIndexType = typename IndexTraits::GridIndex; + using LocalIndexType = typename IndexTraits::LocalIndex; + using Element = typename GV::template Codim<0>::Entity; + + static const int dim = GV::dimension; + static const int dimWorld = GV::dimensionworld; + +public: + //! export the type of the fv element geometry (the local view type) + using LocalView = typename Traits::template LocalView; + //! export the type of sub control volume + using SubControlVolume = typename Traits::SubControlVolume; + //! export the type of sub control volume + using SubControlVolumeFace = typename Traits::SubControlVolumeFace; + //! export dof mapper type + using DofMapper = typename Traits::ElementMapper; + + //! export the discretization method this geometry belongs to + static constexpr DiscretizationMethod discMethod = DiscretizationMethod::ccwmpfa; + + //! The maximum admissible stencil size (used for static memory allocation during assembly) + static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1; + + //! export the grid view type + using GridView = GV; + + //! Constructor + CCWMpfaFVGridGeometry(const GridView& gridView) + : ParentType(gridView) + { + // Check if the overlap size is what we expect + if (!CheckOverlapSize::isValid(gridView)) + DUNE_THROW(Dune::InvalidStateException, "The ccwmpfa discretization method needs at least an overlap of 1 for parallel computations. " + << " Set the parameter \"Grid.Overlap\" in the input file."); + } + + //! the element mapper is the dofMapper + //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa... + const DofMapper& dofMapper() const + { return this->elementMapper(); } + + //! The total number of sub control volumes + std::size_t numScv() const + { + return scvs_.size(); + } + + //! The total number of sub control volume faces + std::size_t numScvf() const + { + return scvfs_.size(); + } + + //! The total number of boundary sub control volume faces + std::size_t numBoundaryScvf() const + { + return numBoundaryScvf_; + } + + //! The total number of degrees of freedom + std::size_t numDofs() const + { return this->gridView().size(0); } + + //! update all fvElementGeometries (do this again after grid adaption) + void update() + { + ParentType::update(); + + // clear containers (necessary after grid refinement) + scvs_.clear(); + scvfs_.clear(); + scvfIndicesOfScv_.clear(); + flipScvfIndices_.clear(); + + // determine size of containers + std::size_t numScvs = numDofs(); + std::size_t numScvf = 0; + for (const auto& element : elements(this->gridView())) + numScvf += element.subEntities(1); + + // reserve memory + scvs_.resize(numScvs); + scvfs_.reserve(numScvf); + scvfIndicesOfScv_.resize(numScvs); + hasBoundaryScvf_.resize(numScvs, false); + + // Build the scvs and scv faces + GridIndexType scvfIdx = 0; + numBoundaryScvf_ = 0; + for (const auto& element : elements(this->gridView())) + { + const auto eIdx = this->elementMapper().index(element); + scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx); + + // the element-wise index sets for finite volume geometry + std::vector scvfsIndexSet; + scvfsIndexSet.reserve(element.subEntities(1)); + + // for network grids there might be multiple intersection with the same geometryInInside + // we indentify those by the indexInInside for now (assumes conforming grids at branching facets) + using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage; + + LocalIndexType localScvfIdx = 0; + for (const auto& intersection : intersections(this->gridView(), element)) + { + // inner sub control volume faces (includes periodic boundaries) + if (intersection.neighbor()) + { + // update the grid geometry if we have periodic boundaries + if (intersection.boundary()) + this->setPeriodic(); + + if (dim == dimWorld) + { + const auto nIdx = this->elementMapper().index(intersection.outside()); + scvfs_.emplace_back(intersection, + intersection.geometry(), + scvfIdx, + localScvfIdx, + ScvfGridIndexStorage({eIdx, nIdx}), + false); + scvfsIndexSet.push_back(scvfIdx++); + } + // this is for network grids + // (will be optimized away of dim == dimWorld) + else + DUNE_THROW(Dune::NotImplemented, "Weighted MPFA schemes not implemented for network grids"); + } + // boundary sub control volume faces + else if (intersection.boundary()) + { + scvfs_.emplace_back(intersection, + intersection.geometry(), + scvfIdx, + localScvfIdx, + ScvfGridIndexStorage({eIdx, static_cast(this->gridView().size(0) + numBoundaryScvf_++)}), + true); + scvfsIndexSet.push_back(scvfIdx++); + + hasBoundaryScvf_[eIdx] = true; + } + ++localScvfIdx; + } + + // Save the scvf indices belonging to this scv to build up fv element geometries fast + scvfIndicesOfScv_[eIdx] = scvfsIndexSet; + } + + // Make the flip index set for network, surface, and periodic grids + if (this->isPeriodic()) + { + flipScvfIndices_.resize(scvfs_.size()); + for (auto&& scvf : scvfs_) + { + if (scvf.boundary()) + continue; + + flipScvfIndices_[scvf.index()].resize(scvf.numOutsideScvs()); + const auto insideScvIdx = scvf.insideScvIdx(); + // check which outside scvf has the insideScvIdx index in its outsideScvIndices + for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) + flipScvfIndices_[scvf.index()][i] = findFlippedScvfIndex_(insideScvIdx, scvf.outsideScvIdx(i)); + } + } + + // build the connectivity map for an effecient assembly + connectivityMap_.update(*this); + } + + //! Get a sub control volume with a global scv index + const SubControlVolume& scv(GridIndexType scvIdx) const + { + return scvs_[scvIdx]; + } + + //! Get a sub control volume face with a global scvf index + const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const + { + return scvfs_[scvfIdx]; + } + + //! Get the scvf on the same face but from the other side + //! Note that e.g. the normals might be different in the case of surface grids + const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const + { + return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; + } + + //! Get the sub control volume face indices of an scv by global index + const std::vector& scvfIndicesOfScv(GridIndexType scvIdx) const + { + return scvfIndicesOfScv_[scvIdx]; + } + + /*! + * \brief Returns the connectivity map of which dofs have derivatives with respect + * to a given dof. + */ + const ConnectivityMap &connectivityMap() const + { return connectivityMap_; } + + //! Returns whether one of the geometry's scvfs lies on a boundary + bool hasBoundaryScvf(GridIndexType eIdx) const + { return hasBoundaryScvf_[eIdx]; } + +private: + // find the scvf that has insideScvIdx in its outsideScvIdx list and outsideScvIdx as its insideScvIdx + GridIndexType findFlippedScvfIndex_(GridIndexType insideScvIdx, GridIndexType outsideScvIdx) + { + // go over all potential scvfs of the outside scv + for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx]) + { + const auto& outsideScvf = this->scvf(outsideScvfIndex); + for (unsigned int j = 0; j < outsideScvf.numOutsideScvs(); ++j) + if (outsideScvf.outsideScvIdx(j) == insideScvIdx) + return outsideScvf.index(); + } + + DUNE_THROW(Dune::InvalidStateException, "No flipped version of this scvf found!"); + } + + //! connectivity map for efficient assembly + ConnectivityMap connectivityMap_; + + //! containers storing the global data + std::vector scvs_; + std::vector scvfs_; + std::vector> scvfIndicesOfScv_; + std::size_t numBoundaryScvf_; + std::vector hasBoundaryScvf_; + + //! needed for embedded surface and network grids (dim < dimWorld) + std::vector> flipScvfIndices_; +}; + +} // end namespace Dumux + +#endif -- GitLab From 898afb8f59e195533ecf0bf8980a4580dfdf1e7a Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:07:59 +0200 Subject: [PATCH 03/36] [disc][wmpfa] Add draft of elementvolumevariables --- .../wmpfa/elementvolumevariables.hh | 243 ++++++++++++++++++ 1 file changed, 243 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh diff --git a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh new file mode 100644 index 0000000000..4618ee503f --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh @@ -0,0 +1,243 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered tpfa models + */ +#ifndef DUMUX_DISCRETIZATION_CCWMPFA_ELEMENT_VOLUMEVARIABLES_HH +#define DUMUX_DISCRETIZATION_CCWMPFA_ELEMENT_VOLUMEVARIABLES_HH + +#include +#include +#include + +#include + +namespace Dumux { + +namespace CCWMpfa { + + /*! + * \ingroup CCMpfaDiscretization + * \brief Computes how many boundary vol vars come into play for flux calculations + * on an element (for a given element finite volume geometry). This number here + * is probably always higher than the actually needed number of volume variables. + * However, we want to make sure it is high enough so that enough memory is reserved + * in the element volume variables below. + * \todo TODO What about non-symmetric schemes? Is there a better way for estimating this? + * + * \param fvGeometry the element finite volume geometry + */ + template + std::size_t maxNumBoundaryVolVars(const FVElementGeometry& fvGeometry) + { + const auto& gridGeometry = fvGeometry.gridGeometry(); + + std::size_t numBoundaryVolVars = fvGeometry.numScvf(); + for (const auto& scvf : scvfs(fvGeometry)) + { + if(!scvf.boundary()) + { + const auto outsideScvIdx = scvf.outsideScvIdx(); + const auto outsideElement = gridGeometry.element(outsideScvIdx); + FVElementGeometry fvGeometryJ; + fvGeometryJ.bindElement(outsideElement); + if(fvGeometryJ.hasBoundaryScvf()) + numBoundaryVolVars += fvGeometryJ.numScvf(); + } + } + + return numBoundaryVolVars; + } + + /*! + * \ingroup CCMpfaDiscretization + * \brief Adds the boundary volume variables found within the stencil to the + * provided containers and stores the indices associated with them. + * + * \param volVars The container where the volume variables are stored + * \param volVarIndices The container where the volume variable indices are stored + * \param problem The problem containing the Dirichlet boundary conditions + * \param element The element to which the finite volume geometry was bound + * \param fvGeometry The element finite volume geometry + */ + template + std::pair, std::vector> + boundaryVolVars(const Problem& problem, + const typename FVElemGeom::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElemGeom& fvGeometry) + { + std::vector volVars; volVars.reserve(fvGeometry.numScvf()); + std::vector volVarIndices; volVarIndices.reserve(fvGeometry.numScvf()); + const auto& gridGeometry = fvGeometry.gridGeometry(); + + // treat the BCs inside the element + if (fvGeometry.hasBoundaryScvf()) + { + const auto boundElemIdx = gridGeometry.elementMapper().index(element); + const auto& scvI = fvGeometry.scv(boundElemIdx); + + for (const auto& scvf : scvfs(fvGeometry)) + { + if (!scvf.boundary()) + continue; + + // Only proceed on dirichlet boundaries. On Neumann + // boundaries the "outside" vol vars cannot be properly defined. + if (problem.boundaryTypes(element, scvf).hasOnlyDirichlet()) + { + VolumeVariables dirichletVolVars; + dirichletVolVars.update(elementSolution(problem.dirichlet(element, scvf)), + problem, + element, + scvI); + + volVars.emplace_back(std::move(dirichletVolVars)); + volVarIndices.push_back(scvf.outsideScvIdx()); + } + } + } + return std::make_pair(volVars, volVarIndices); + } +} // end namespace CCMpfa + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered tpfa models + * \note The class is specilized for versions with and without caching + * \tparam GVV the grid volume variables type + * \tparam cachingEnabled if the cache is enabled + */ +template +class CCWMpfaElementVolumeVariables +{}; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered tpfa models with caching + * \note the volume variables are stored for the whole grid view in the corresponding GridVolumeVariables class + */ +template +class CCWMpfaElementVolumeVariables +{ +public: + //! export type of the grid volume variables + using GridVolumeVariables = GVV; + + //! export type of the volume variables + using VolumeVariables = typename GridVolumeVariables::VolumeVariables; + + //! Constructor + CCWMpfaElementVolumeVariables(const GridVolumeVariables& gridVolVars) + : gridVolVarsPtr_(&gridVolVars) + , numScv_(gridVolVars.problem().gridGeometry().numScv()) + {} + + //! operator for the access with an scv + template::value, int> = 0> + const VolumeVariables& operator [](const SubControlVolume& scv) const + { + if (scv.dofIndex() < numScv_) + return gridVolVars().volVars(scv.dofIndex()); + else + return boundaryVolumeVariables_[getLocalIdx_(scv.dofIndex())]; + } + + //! operator for the access with an index + const VolumeVariables& operator [](const std::size_t scvIdx) const + { + if (scvIdx < numScv_) + return gridVolVars().volVars(scvIdx); + else + return boundaryVolumeVariables_[getLocalIdx_(scvIdx)]; + } + + //! precompute all boundary volume variables in a stencil of an element, the remaining ones are cached + template + void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const SolutionVector& sol) + { + if (!fvGeometry.hasBoundaryScvf()) + return; + + clear_(); + + auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), element, fvGeometry); + std::move(volVars.begin(), volVars.end(), std::back_inserter(boundaryVolumeVariables_)); + std::move(indices.begin(), indices.end(), std::back_inserter(boundaryVolVarIndices_)); + + // add boundary volVars of neighbors + for (const auto& scvf : scvfs(fvGeometry)) + { + if (scvf.boundary()) + continue; + + const auto outsideScvIdx = scvf.outsideScvIdx(); + const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx); + auto fvGeometryJ = localView(fvGeometry.gridGeometry()); + fvGeometryJ.bind(outsideElement); + + if (!fvGeometryJ.hasBoundaryScvf()) + continue; + + auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), outsideElement, fvGeometryJ); + std::move(volVars.begin(), volVars.end(), std::back_inserter(boundaryVolumeVariables_)); + std::move(indices.begin(), indices.end(), std::back_inserter(boundaryVolVarIndices_)); + } + } + + //! precompute the volume variables of an element - do nothing: volVars are cached + template + void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const SolutionVector& sol) + {} + + //! The global volume variables object we are a restriction of + const GridVolumeVariables& gridVolVars() const + { return *gridVolVarsPtr_; } + +private: + //! Clear all local storage + void clear_() + { + boundaryVolVarIndices_.clear(); + boundaryVolumeVariables_.clear(); + } + + const GridVolumeVariables* gridVolVarsPtr_; + + //! map a global scv index to the local storage index + int getLocalIdx_(const int volVarIdx) const + { + auto it = std::find(boundaryVolVarIndices_.begin(), boundaryVolVarIndices_.end(), volVarIdx); + assert(it != boundaryVolVarIndices_.end() && "Could not find the current volume variables for volVarIdx!"); + return std::distance(boundaryVolVarIndices_.begin(), it); + } + + std::vector boundaryVolVarIndices_; + std::vector boundaryVolumeVariables_; + const std::size_t numScv_; +}; + +} // end namespace Dumux + +#endif -- GitLab From e3926b39160e8e2680a2f4fb40336efff86ba324 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:08:50 +0200 Subject: [PATCH 04/36] [disc][wmpfa] Interpolation operator class Currently only the harmonic averaging point interpolator is implemented --- .../wmpfa/interpolationoperator.hh | 204 ++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh diff --git a/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh new file mode 100644 index 0000000000..183b11cc66 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh @@ -0,0 +1,204 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCMpfaDiscretization + * \brief Base class for interaction volumes of mpfa methods. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERPOLATIONOPERATOR_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_INTERPOLATIONOPERATOR_HH + +#include + +#include +#include + +#include "dumux/common/math.hh" + +namespace Dumux { + + +//! Empty interpolator class +class EmptyInterpolator {}; + +template< class T, bool Enable> +class HapInterpolatorBase +{ + using GridView = typename T::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using Position = typename T::GlobalPosition; + + using Scalar = typename Position::value_type; + + struct LocalInterpolationData + { + Position point; + std::pair weights; + }; + +public: + //! When global caching is enabled, precompute transmissibilities for all scv faces + void clear() + { + interpolationData_.clear(); + isUpdated_ = false; + } + + template + void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const TF& tensor) + { + clear(); + interpolationData_.resize(fvGeometry.numScvf()); + for (auto&& scvf : scvfs(fvGeometry)) + { + if(!scvf.boundary()) + { + const auto insideScvIdx = scvf.insideScvIdx(); + const auto outsideScvIdx = scvf.outsideScvIdx(); + const auto& insideVolVars = elemVolVars[insideScvIdx]; + const auto& outsideVolVars = elemVolVars[outsideScvIdx]; + + const auto tauK = vtmv(scvf.unitOuterNormal(), tensor(insideVolVars), scvf.unitOuterNormal()); + const auto tauL = vtmv(scvf.unitOuterNormal(), tensor(outsideVolVars), scvf.unitOuterNormal()); + + using std::abs; + const auto centerK = fvGeometry.scv(insideScvIdx).center(); + const auto centerL = fvGeometry.scv(outsideScvIdx).center(); + const auto distK = abs((scvf.center() - fvGeometry.scv(insideScvIdx).center()) * scvf.unitOuterNormal()); + const auto distL = abs((scvf.center() - fvGeometry.scv(outsideScvIdx).center()) * scvf.unitOuterNormal()); + + const auto omegaK = distL*tauK / (distL*tauK + distK*tauL); + const auto omegaL = distK*tauL / (distL*tauK + distK*tauL); + + auto point = (omegaK * centerK) + (omegaL * centerL) + + (distL*distK / (distL*tauK + distK*tauL)) + * mv(tensor(insideVolVars) - tensor(outsideVolVars), scvf.unitOuterNormal()); + + interpolationData_[scvf.localIndex()] = (LocalInterpolationData({point, std::make_pair(omegaK, omegaL)})); + } + else + { + interpolationData_[scvf.localIndex()] = (LocalInterpolationData({scvf.center(), std::make_pair(1.0, 0.0)})); + } + } + isUpdated_ = true; + } + + bool isUpdated() + { + return isUpdated_; + } + + template + const std::vector getDistanceVectors(const FVElementGeometry& fvGeometry) const + { + std::vector distances; + distances.resize(fvGeometry.numScvf()); + for (auto&& scvf : scvfs(fvGeometry)) + distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].point - fvGeometry.scv(scvf.insideScvIdx()).center()); + + return distances; + } + + template + std::vector getDistanceVectors(const FVElementGeometry& fvGeometry) + { + std::vector distances; + distances.resize(fvGeometry.numScvf()); + for (auto&& scvf : scvfs(fvGeometry)) + distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].point - fvGeometry.scv(scvf.insideScvIdx()).center()); + + return distances; + } + + +private: + std::vector interpolationData_; + bool isUpdated_; +}; + +template +class HapInterpolatorBase : public EmptyInterpolator {}; + +template +class HapInterpolationOperator +{ + +public: + //! state the traits type publicly + using Traits = T; + using Position = typename T::GlobalPosition; + using Scalar = typename Position::value_type; + + static constexpr bool advectionEnabled = PT::enableAdvection; + static constexpr bool diffusionEnabled = PT::enableMolecularDiffusion; + static constexpr bool heatConductionEnabled = PT::enableHeatConduction; + + //! export the underlying process-specific interpolator types + using AdvectionInterpolator = HapInterpolatorBase; + using DiffusionInterpolator = HapInterpolatorBase; + using HeatConductionInterpolator = HapInterpolatorBase; + + //! return references to the interpolator containing data related to advection + const AdvectionInterpolator& advectionInterpolator() const { return advectionInterpolator_; } + AdvectionInterpolator& advectionInterpolator() { return advectionInterpolator_; } + + //! return references to the interpolator containing data related to diffusion + const DiffusionInterpolator& diffusionInterpolator() const { return diffusionInterpolator_; } + DiffusionInterpolator& diffusionInterpolator() { return diffusionInterpolator_; } + + //! return references to the interpolator containing data related to heat conduction + const HeatConductionInterpolator& heatConductionInterpolator() const { return heatConductionInterpolator_; } + HeatConductionInterpolator& heatConductionInterpolator() { return heatConductionInterpolator_; } + + template + void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) + { + // (maybe) solve system subject to intrinsic permeability + if constexpr (advectionEnabled) + { + auto getTensor = [] (const auto& volVars) { return volVars.permeability(); }; + advectionInterpolator_.bind(element, fvGeometry, elemVolVars, getTensor); + } + + // (maybe) solve system subject to diffusion tensors + if constexpr (diffusionEnabled) + { + } + + // (maybe) solve system subject to thermal conductivity + if constexpr (heatConductionEnabled) + { + } + } + +private: + AdvectionInterpolator advectionInterpolator_; + DiffusionInterpolator diffusionInterpolator_; + HeatConductionInterpolator heatConductionInterpolator_; +}; + +} // end namespace Dumux + +#endif -- GitLab From a84df56595e3c0c9bc2e55ef400c71529646cc8a Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:09:59 +0200 Subject: [PATCH 05/36] [disc][wmpfa] Face data handle Data handle for the conormal decomposition and for the calculation of one-sided fluxes --- .../cellcentered/wmpfa/facedatahandle.hh | 197 ++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/facedatahandle.hh diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh new file mode 100644 index 0000000000..aefdc2be1c --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -0,0 +1,197 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief Data handle class for face data of wmpfa methods. + */ +#ifndef DUMUX_DISCRETIZATION_CC_WMPFA_FACEDATAHANDLE_HH +#define DUMUX_DISCRETIZATION_CC_WMPFA_FACEDATAHANDLE_HH + +#include +#include + +#include +#include "dumux/common/vectordecomposition.hh" + +namespace Dumux { + +namespace WMpfaHelper{ + +template +void eraseZeros(Coeff& c, Stencil& s) +{ + assert(c.size() == s.size()); + for(std::size_t i=0; i +class AdvectionDataHandle +{ + using Scalar = typename IntOp::Scalar; + using OneSidedStencil = std::vector; + using Stencils = std::pair; + using OneSidedCoefficients = std::vector; + using Coefficients = std::pair; + + using Interpolator = typename IntOp::AdvectionInterpolator; + + static constexpr int numPhases = PhysicsTraits::numPhases; + +public: + + template + void decompose(const Interpolator& intOp, const TF& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + { + boundaryFace_ = scvf.boundary(); + const auto coNormal = mv(tensor(scvf.insideScvIdx()), scvf.unitOuterNormal()); + auto&& [indices, coeff, found] = VectorDecomposition::calculateVectorDecomposition(coNormal, intOp.getDistanceVectors(fvGeometry)); + + if(!found) + DUNE_THROW(Dune::InvalidStateException, "CoNormal decomposition not found"); + else + { + WMpfaHelper::eraseZeros(coeff, indices); + + if(scvf.insideScvIdx() < scvf.outsideScvIdx()) + { + coefficients_[0].first = std::move(coeff); + stencils_[0].first = std::move(indices); + } + else + { + coefficients_[0].second = std::move(coeff); + stencils_[0].second = std::move(indices); + } + } + } + + bool valid() + { + return (!boundaryFace_ && coefficients_[0].first.size() > 0 && coefficients_[0].second.size() > 0) + || (boundaryFace_ && coefficients_[0].first.size() > 0); + } + +private: + std::array coefficients_; + std::array stencils_; + bool boundaryFace_; +}; + +//! Data handle for quantities related to diffusion +template +class DiffusionDataHandle +{ + static constexpr int numPhases = PhysicsTraits::numPhases; + static constexpr int numComponents = PhysicsTraits::numComponents; + + using Interpolator = typename IntOp::DiffusionInterpolator; + +public: + + template + void decompose(const Interpolator& intOp, const TF& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + {} +private: + using Scalar = typename IntOp::Scalar; + using OneSidedCoefficients = std::vector; + using Coefficients = std::pair; + std::array coNormalDecomposition_; +}; + +//! Data handle for quantities related to heat conduction +template +class HeatConductionDataHandle +{ + using Interpolator = typename IntOp::HeatConductionInterpolator; + +public: + + template + void decompose(const Interpolator& intOp, const TF& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + {} + +private: + using Scalar = typename IntOp::Scalar; + using OneSidedStencil = std::vector; + using Stencils = std::pair; + using OneSidedCoefficients = std::vector; + using Coefficients = std::pair; + std::array coefficients_; + std::array stencils; +}; + + +//! Process-dependent data handles when related process is disabled +template +class AdvectionDataHandle : public EmptyDataHandle {}; +template +class DiffusionDataHandle : public EmptyDataHandle {}; +template +class HeatConductionDataHandle : public EmptyDataHandle {}; + +} + +template +class FaceDataHandle +{ + +public: + //! export the underlying process-specific handle types + using AdvectionHandle = WMpfaDataHandle::AdvectionDataHandle; + using DiffusionHandle = WMpfaDataHandle::DiffusionDataHandle; + using HeatConductionHandle = WMpfaDataHandle::HeatConductionDataHandle; + + //! return references to the handle containing data related to advection + const AdvectionHandle& advectionHandle() const { return advectionHandle_; } + AdvectionHandle& advectionHandle() { return advectionHandle_; } + + //! return references to the handle containing data related to diffusion + const DiffusionHandle& diffusionHandle() const { return diffusionHandle_; } + DiffusionHandle& diffusionHandle() { return diffusionHandle_; } + + //! return references to the handle containing data related to heat conduction + const HeatConductionHandle& heatConductionHandle() const { return heatConductionHandle_; } + HeatConductionHandle& heatConductionHandle() { return heatConductionHandle_; } + +private: + AdvectionHandle advectionHandle_; + DiffusionHandle diffusionHandle_; + HeatConductionHandle heatConductionHandle_; +}; + +} // end namespace Dumux + +#endif -- GitLab From 277f88e5674bdb4e7c64af09e025107934018e50 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:13:52 +0200 Subject: [PATCH 06/36] [disc][wmpfa] Add draft for caches --- .../wmpfa/elementfluxvariablescache.hh | 116 +++++++++ .../wmpfa/gridfluxvariablescache.hh | 237 ++++++++++++++++++ 2 files changed, 353 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/elementfluxvariablescache.hh create mode 100644 dumux/discretization/cellcentered/wmpfa/gridfluxvariablescache.hh diff --git a/dumux/discretization/cellcentered/wmpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/wmpfa/elementfluxvariablescache.hh new file mode 100644 index 0000000000..31402a606d --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/elementfluxvariablescache.hh @@ -0,0 +1,116 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The flux variables caches for an element + */ +#ifndef DUMUX_DISCRETIZATION_CCWMPFA_ELEMENT_FLUXVARSCACHE_HH +#define DUMUX_DISCRETIZATION_CCWMPFA_ELEMENT_FLUXVARSCACHE_HH + +#include +#include +#include +#include + +namespace Dumux { + +template +struct DataStorage +{ + std::vector interpolationOperators; + std::vector facesDataHandles; +}; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The flux variables caches for an element + * \note The class is specialized for a version with and without caching + * If grid caching is enabled the flux caches are stored for the whole gridview in the corresponding + * GridFluxVariablesCache which is memory intensive but faster. For caching disabled the + * flux caches are locally computed for each element whenever needed. + */ +template +class CCWMpfaElementFluxVariablesCache; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The flux variables caches for an element with caching enabled + */ +template +class CCWMpfaElementFluxVariablesCache +{ + //! the type of the flux variables cache filler + using FluxVariablesCacheFiller = typename GFVC::Traits::FluxVariablesCacheFiller; + +public: + //! export the type of the grid flux variables cache + using GridFluxVariablesCache = GFVC; + + //! export the type of the flux variables cache + using FluxVariablesCache = typename GFVC::FluxVariablesCache; + + //! export the data handle types used + using InterpolationOperator = typename GFVC::InterpolationOperator; + using FaceDataHandle = typename GFVC::FaceDataHandle; + + CCWMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global) + : gridFluxVarsCachePtr_(&global) {} + + //! Specialization for the global caching being enabled - do nothing here + template + void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) {} + + //! Specialization for the global caching being enabled - do nothing here + template + void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) {} + + //! Specialization for the global caching being enabled - do nothing here + template + void bindScvf(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const typename FVElementGeometry::SubControlVolumeFace& scvf) {} + + //! Specialization for the global caching being enabled - do nothing here + template + void update(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) {} + + //! access operators in the case of caching + template + const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const + { return gridFluxVarsCache()[scvf]; } + + //! The global object we are a restriction of + const GridFluxVariablesCache& gridFluxVarsCache() const + { return *gridFluxVarsCachePtr_; } + +private: + const GridFluxVariablesCache* gridFluxVarsCachePtr_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/cellcentered/wmpfa/gridfluxvariablescache.hh b/dumux/discretization/cellcentered/wmpfa/gridfluxvariablescache.hh new file mode 100644 index 0000000000..c5ee7af6e9 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/gridfluxvariablescache.hh @@ -0,0 +1,237 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief Flux variable caches on a gridview + */ +#ifndef DUMUX_DISCRETIZATION_CCWMPFA_GRID_FLUXVARSCACHE_HH +#define DUMUX_DISCRETIZATION_CCWMPFA_GRID_FLUXVARSCACHE_HH + +// make the local view function available whenever we use this class +#include +#include "elementfluxvariablescache.hh" + +namespace Dumux { + +template +struct DataHandlePhysicsTraits +{ + static constexpr bool enableAdvection = ModelTraits::enableAdvection(); + static constexpr bool enableMolecularDiffusion = ModelTraits::enableMolecularDiffusion(); + static constexpr bool enableHeatConduction = ModelTraits::enableEnergyBalance(); + + static constexpr int numPhases = ModelTraits::numFluidPhases(); + static constexpr int numComponents = ModelTraits::numFluidComponents(); +}; + +template +struct CCWMpfaDefaultGridFluxVariablesCacheTraits +{ + using Problem = P; + using FluxVariablesCache = FVC; + using FluxVariablesCacheFiller = FVCF; + + using InterpolationOperator = IOP; + using FaceDataHandle = FDH; + + template + using LocalView = CCWMpfaElementFluxVariablesCache; +}; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Flux variable caches on a gridview + * \note The class is specialized for a version with and without grid caching + */ +template +class CCWMpfaGridFluxVariablesCache; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Flux variable caches on a gridview with grid caching enabled + * \note The flux caches of the gridview are stored which is memory intensive but faster + */ +template +class CCWMpfaGridFluxVariablesCache +{ + using Problem = typename TheTraits::Problem; + using ThisType = CCWMpfaGridFluxVariablesCache; + + //! the flux variable cache filler type + using FluxVariablesCacheFiller = typename TheTraits::FluxVariablesCacheFiller; +public: + //! export the Traits + using Traits = TheTraits; + + //! export the data handle types used + using InterpolationOperator = typename Traits::InterpolationOperator; + using FaceDataHandle = typename Traits::FaceDataHandle; + + //! export the flux variable cache type + using FluxVariablesCache = typename Traits::FluxVariablesCache; + + //! make it possible to query if caching is enabled + static constexpr bool cachingEnabled = true; + + //! export the type of the local view + using LocalView = typename Traits::template LocalView; + + //! The constructor + CCWMpfaGridFluxVariablesCache(const Problem& problem) + : problemPtr_(&problem) + {} + + //! When global caching is enabled, precompute transmissibilities for all scv faces + template + void update(const GridGeometry& gridGeometry, + const GridVolumeVariables& gridVolVars, + const SolutionVector& sol, + bool forceUpdate = false) + { + // Update only if the filler puts solution-dependent + // stuff into the caches or if update is enforced + if (FluxVariablesCacheFiller::isSolDependent || forceUpdate) + { + // clear previous data if forced update is desired + if (forceUpdate) + { + clear_(); + + dataStorage_.interpolationOperators.resize(gridGeometry.gridView().size(0)); + dataStorage_.facesDataHandles.resize(gridGeometry.numScvf()); + + // reserve memory estimate for caches, interaction volumes and corresponding data + fluxVarsCache_.resize(gridGeometry.numScvf()); + } + + // instantiate helper class to fill the caches + FluxVariablesCacheFiller filler(problem()); + + for (const auto& element : elements(gridGeometry.gridView())) + { + auto fvGeometry = localView(gridGeometry); + fvGeometry.bind(element); + + auto elemVolVars = localView(gridVolVars); + elemVolVars.bind(element, fvGeometry, sol); + + auto& intOperator = dataStorage_.interpolationOperators[gridGeometry.elementMapper().index(element)]; + intOperator.bind(element, fvGeometry, elemVolVars); + + // Prepare all caches of the scvfs inside the corresponding interaction volume. Skip + // those ivs that are touching a boundary, we only store the data on interior ivs here. + for (const auto& scvf : scvfs(fvGeometry)) + //if (!isEmbeddedInBoundaryIV_(scvf, gridGeometry) && !fluxVarsCache_[scvf.index()].isUpdated()) + filler.fill(*this, fluxVarsCache_[scvf.index()], dataStorage_, element, fvGeometry, elemVolVars, scvf, forceUpdate); + } + } + } + + template + void updateElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) + { + // // Update only if the filler puts + // // solution-dependent stuff into the caches + // if (FluxVariablesCacheFiller::isSolDependent) + // { + // const auto& gridGeometry = fvGeometry.gridGeometry(); + // const auto& assemblyMapI = gridGeometry.connectivityMap()[gridGeometry.elementMapper().index(element)]; + + // // helper class to fill flux variables caches + // FluxVariablesCacheFiller filler(problem()); + + // // first, set all the caches to "outdated" + // for (const auto& scvf : scvfs(fvGeometry)) + // fluxVarsCache_[scvf.index()].setUpdateStatus(false); + // for (const auto& dataJ : assemblyMapI) + // for (const auto scvfIdx : dataJ.scvfsJ) + // fluxVarsCache_[scvfIdx].setUpdateStatus(false); + + // // go through the caches maybe update them + // for (const auto& scvf : scvfs(fvGeometry)) + // { + // auto& scvfCache = fluxVarsCache_[scvf.index()]; + // if (!isEmbeddedInBoundaryIV_(scvf, gridGeometry) && !scvfCache.isUpdated()) + // filler.fill(*this, scvfCache, ivDataStorage_, element, fvGeometry, elemVolVars, scvf); + // } + + // for (const auto& dataJ : assemblyMapI) + // { + // const auto elementJ = gridGeometry.element(dataJ.globalJ); + // for (const auto scvfIdx : dataJ.scvfsJ) + // { + // auto& scvfCache = fluxVarsCache_[scvfIdx]; + // const auto& scvf = fvGeometry.scvf(scvfIdx); + // if (!isEmbeddedInBoundaryIV_(scvf, gridGeometry) && !scvfCache.isUpdated()) + // filler.fill(*this, scvfCache, ivDataStorage_, elementJ, fvGeometry, elemVolVars, scvf); + // } + // } + // } + } + + //! access operators in the case of caching + template + const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const + { return fluxVarsCache_[scvf.index()]; } + + //! access operators in the case of caching + template + FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) + { return fluxVarsCache_[scvf.index()]; } + + template + const InterpolationOperator& interpolationOperator(const SubControlVolume& scv) const + { return dataStorage_.interpolationOperators[scv.elementIndex()]; } + + template + const FaceDataHandle& faceDataHandle(const SubControlVolumeFace& scvf) const + { return dataStorage_.facesDataHandles[scvf.index()]; } + + const Problem& problem() const + { return *problemPtr_; } + +private: + + //! clear all containers + void clear_() + { + fluxVarsCache_.clear(); + dataStorage_.interpolationOperators.clear(); + dataStorage_.facesDataHandles.clear(); + } + + const Problem* problemPtr_; + std::vector fluxVarsCache_; + + // stored interaction volumes and handles + using DS = DataStorage; + DS dataStorage_; +}; + +} // end namespace Dumux + +#endif -- GitLab From e1218e2ce7e157273867ef2420b60c02b55ce215 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:14:13 +0200 Subject: [PATCH 07/36] [disc][wmpfa] Add draft for gridvolumevariables --- .../cellcentered/wmpfa/gridvolumevariables.hh | 64 +++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/gridvolumevariables.hh diff --git a/dumux/discretization/cellcentered/wmpfa/gridvolumevariables.hh b/dumux/discretization/cellcentered/wmpfa/gridvolumevariables.hh new file mode 100644 index 0000000000..74105036f6 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/gridvolumevariables.hh @@ -0,0 +1,64 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The grid volume variables class for cell centered tpfa models + */ +#ifndef DUMUX_DISCRETIZATION_CC_WMPFA_GRID_VOLUMEVARIABLES_HH +#define DUMUX_DISCRETIZATION_CC_WMPFA_GRID_VOLUMEVARIABLES_HH + +#include +#include + +namespace Dumux { + +template +struct CCWMpfaDefaultGridVolumeVariablesTraits +{ + using Problem = P; + using VolumeVariables = VV; + + template + using LocalView = CCWMpfaElementVolumeVariables; +}; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Base class for the grid volume variables + * \note This class has a cached version and a non-cached version + * \tparam Problem the type of problem we are solving + * \tparam VolumeVariables the type of volume variables we are using for the model + * \tparam Traits the traits class injecting the problem, volVar and elemVolVars type + * \tparam cachingEnabled if the cache is enabled + */ +template > +class CCWMpfaGridVolumeVariables : public CCGridVolumeVariables +{ +public: + using ParentType = CCGridVolumeVariables; + using ParentType::ParentType; +}; + +} // end namespace Dumux + +#endif -- GitLab From 501f5b3ea3cf823a26cb59a5fe0b83952234b779 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:14:37 +0200 Subject: [PATCH 08/36] [disc][wmpfa] First draft of Subcontrolvolume faces --- .../wmpfa/subcontrolvolumeface.hh | 236 ++++++++++++++++++ 1 file changed, 236 insertions(+) create mode 100644 dumux/discretization/cellcentered/wmpfa/subcontrolvolumeface.hh diff --git a/dumux/discretization/cellcentered/wmpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/wmpfa/subcontrolvolumeface.hh new file mode 100644 index 0000000000..225cf24761 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/subcontrolvolumeface.hh @@ -0,0 +1,236 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The sub control volume face + */ +#ifndef DUMUX_DISCRETIZATION_CC_WEIGHTED_MPFA_SUBCONTROLVOLUMEFACE_HH +#define DUMUX_DISCRETIZATION_CC_WEIGHTED_MPFA_SUBCONTROLVOLUMEFACE_HH + +#include +#include + +#include +#include +#include + +#include +#include +#include + +namespace Dumux { + +/*! + * \ingroup CCWMpfaDiscretization + * \brief Default traits class to be used for the sub-control volume faces + * for the cell-centered finite volume scheme using weighted MPFA + * \tparam GV the type of the grid view + */ +template +struct CCWMpfaDefaultScvfGeometryTraits +{ + using Grid = typename GridView::Grid; + + static constexpr int dim = Grid::dimension; + static constexpr int dimWorld = Grid::dimensionworld; + + using Scalar = typename Grid::ctype; + using GridIndexType = typename IndexTraits::GridIndex; + using LocalIndexType = typename IndexTraits::LocalIndex; + using GridIndexStorage = typename std::vector; + + // we use geometry traits that use static corner vectors to and a fixed geometry type + template + struct ScvfMLGTraits : public Dune::MultiLinearGeometryTraits + { + // we use static vectors to store the corners as we know + // the number of corners in advance (2^(dim-1) corners (1<<(dim-1)) + template< int mydim, int cdim > + struct CornerStorage + { + using Type = Dune::ReservedVector< Dune::FieldVector< ct, cdim >, (1<<(dim-1)) >; + }; + }; + + using Geometry = Dune::MultiLinearGeometry >; + using CornerStorage = typename ScvfMLGTraits::template CornerStorage::Type; + using GlobalPosition = typename CornerStorage::value_type; + using BoundaryFlag = Dumux::BoundaryFlag; +}; + +/*! + * \ingroup CCWMpfaDiscretization + * \brief The sub control volume face + * \tparam GV the type of the grid view + * \tparam T the scvf geometry traits + */ +template > +class CCWMpfaSubControlVolumeFace +: public SubControlVolumeFaceBase, T> +{ + using ThisType = CCWMpfaSubControlVolumeFace; + using ParentType = SubControlVolumeFaceBase; + using GridIndexType = typename T::GridIndexType; + using LocalIndexType = typename T::LocalIndexType; + using Scalar = typename T::Scalar; + using CornerStorage = typename T::CornerStorage; + using GridIndexStorage = typename T::GridIndexStorage; + using Geometry = typename T::Geometry; + using BoundaryFlag = typename T::BoundaryFlag; + +public: + //! export the type used for global coordinates + using GlobalPosition = typename T::GlobalPosition; + //! state the traits public and thus export all types + using Traits = T; + + // the default constructor + CCWMpfaSubControlVolumeFace() = default; + + /*! + * \brief Constructor with intersection + * + * \param is The intersection + * \param isGeometry The geometry of the intersection + * \param scvfIndex The global index of this scv face + * \param scvIndices The inside/outside scv indices connected to this face + * \param isBoundary Bool to specify whether or not the scvf is on an interior or the domain boundary + */ + template + CCWMpfaSubControlVolumeFace(const Intersection& is, + const typename Intersection::Geometry& isGeometry, + GridIndexType scvfIndex, + LocalIndexType localScvfIndex, + const GridIndexStorage& scvIndices, + bool isBoundary) + : ParentType() + , geomType_(isGeometry.type()) + , area_(isGeometry.volume()) + , center_(isGeometry.center()) + , unitOuterNormal_(is.centerUnitOuterNormal()) + , scvfIndex_(scvfIndex) + , localScvfIndex_(localScvfIndex) + , scvIndices_(scvIndices) + , boundary_(isBoundary) + , boundaryFlag_{is} + { + corners_.resize(isGeometry.corners()); + for (int i = 0; i < isGeometry.corners(); ++i) + corners_[i] = isGeometry.corner(i); + } + + //! The center of the sub control volume face + const GlobalPosition& center() const + { + return center_; + } + + //! The integration point for flux evaluations in global coordinates + const GlobalPosition& ipGlobal() const + { + // Return center for now + return center_; + } + + //! The area of the sub control volume face + Scalar area() const + { + return area_; + } + + //! returns bolean if the sub control volume face is on the boundary + bool boundary() const + { + return boundary_; + } + + //! The unit outer normal of the sub control volume face + const GlobalPosition& unitOuterNormal() const + { + return unitOuterNormal_; + } + + //! index of the inside sub control volume for spatial param evaluation + GridIndexType insideScvIdx() const + { + return scvIndices_[0]; + } + + //! index of the outside sub control volume for spatial param evaluation + // This results in undefined behaviour if boundary is true + GridIndexType outsideScvIdx(int i = 0) const + { + return scvIndices_[i+1]; + } + + //! The number of outside scvs connection via this scv face + std::size_t numOutsideScvs() const + { + return scvIndices_.size()-1; + } + + //! The global index of this sub control volume face + GridIndexType index() const + { + return scvfIndex_; + } + + //! The local index of this sub control volume face + LocalIndexType localIndex() const + { + return localScvfIndex_; + } + + //! return the i-th corner of this sub control volume face + const GlobalPosition& corner(int i) const + { + assert(i < corners_.size() && "provided index exceeds the number of corners"); + return corners_[i]; + } + + //! The geometry of the sub control volume face + Geometry geometry() const + { + return Geometry(geomType_, corners_); + } + + //! Return the boundary flag + typename BoundaryFlag::value_type boundaryFlag() const + { + return boundaryFlag_.get(); + } + +private: + Dune::GeometryType geomType_; + CornerStorage corners_; + Scalar area_; + GlobalPosition center_; + GlobalPosition unitOuterNormal_; + GridIndexType scvfIndex_; + LocalIndexType localScvfIndex_; + GridIndexStorage scvIndices_; + bool boundary_; + BoundaryFlag boundaryFlag_; +}; + +} // end namespace Dumux + +#endif -- GitLab From 71a5c7cc56229390f227122115bf7634b3a32123 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 17:15:36 +0200 Subject: [PATCH 09/36] [disc][wmpfa] Add WMpfa model class --- dumux/discretization/ccwmpfa.hh | 139 ++++++++++++++++++ .../cellcentered/wmpfa/methods.hh | 40 +++++ 2 files changed, 179 insertions(+) create mode 100644 dumux/discretization/ccwmpfa.hh create mode 100644 dumux/discretization/cellcentered/wmpfa/methods.hh diff --git a/dumux/discretization/ccwmpfa.hh b/dumux/discretization/ccwmpfa.hh new file mode 100644 index 0000000000..86630c8e97 --- /dev/null +++ b/dumux/discretization/ccwmpfa.hh @@ -0,0 +1,139 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Discretization + * \brief Properties for all models using cell-centered finite volume scheme with WMPFA + * \note Inherit from these properties to use a cell-centered finite volume scheme with WMPFA + */ + +#ifndef DUMUX_DISCRETIZATION_CC_WMPFA_HH +#define DUMUX_DISCRETIZATION_CC_WMPFA_HH + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Dumux { +namespace Properties { + +//! Type tag for the cell-centered wmpfa scheme. +// Create new type tags +namespace TTag { +struct CCWMpfaModel { using InheritsFrom = std::tuple; }; +} // end namespace TTag + +//! Set the default for the grid geometry +template +struct GridGeometry +{ +private: + static constexpr bool enableCache = getPropValue(); + using GridView = typename GetPropType::LeafGridView; +public: + using type = CCWMpfaFVGridGeometry; +}; + +//! The grid volume variables vector class +template +struct GridVolumeVariables +{ +private: + static constexpr bool enableCache = getPropValue(); + using Problem = GetPropType; + using VolumeVariables = GetPropType; +public: + using type = CCWMpfaGridVolumeVariables; +}; + +//! The grid flux variables cache vector class +template +struct GridFluxVariablesCache +{ +private: + static constexpr bool enableCache = getPropValue(); + using Problem = GetPropType; + using FluxVariablesCache = GetPropType; + using FluxVariablesCacheFiller = GetPropType; + + using PhysicsTraits = DataHandlePhysicsTraits>; + + struct InterpolationTraits + { + using GridGeometry = GetPropType; + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + }; + + using InterpolationOperator = HapInterpolationOperator; + using DataHandle = FaceDataHandle; + + using Traits = CCWMpfaDefaultGridFluxVariablesCacheTraits; + +public: + using type = CCWMpfaGridFluxVariablesCache; +}; + +//! Set the default for the ElementBoundaryTypes +template +struct ElementBoundaryTypes { using type = CCElementBoundaryTypes; }; + +//! Set the BaseLocalResidual to CCLocalResidual +template +struct BaseLocalResidual { using type = CCLocalResidual; }; +} // namespace Properties + +namespace Detail { + +template +struct ProblemTraits +{ +private: + using GG = std::decay_t().gridGeometry())>; + using Element = typename GG::GridView::template Codim<0>::Entity; + using SubControlVolumeFace = typename GG::SubControlVolumeFace; +public: + using GridGeometry = GG; + // BoundaryTypes is whatever the problem returns from boundaryTypes(element, scvf) + using BoundaryTypes = std::decay_t().boundaryTypes(std::declval(), std::declval()))>; +}; + +} // end namespace Detail + +} // namespace Dumux + +#endif diff --git a/dumux/discretization/cellcentered/wmpfa/methods.hh b/dumux/discretization/cellcentered/wmpfa/methods.hh new file mode 100644 index 0000000000..29218bdca6 --- /dev/null +++ b/dumux/discretization/cellcentered/wmpfa/methods.hh @@ -0,0 +1,40 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaDiscretization + * \brief The available weighted mpfa schemes in Dumux + */ +#ifndef DUMUX_DISCRETIZATION_CC_WMPFA_METHODS_HH +#define DUMUX_DISCRETIZATION_CC_WMPFA_METHODS_HH + +namespace Dumux { + + /*! + * \brief The available weighted mpfa schemes in Dumux + * \ingroup CCWMpfaDiscretization + */ + enum class WMpfaMethods : unsigned int + { + avgMethod, nltpfa, nlmpfa + }; + +} // end namespace Dumux + +#endif -- GitLab From db384761fa50a56f14d2353884b3ab4167045ede Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 18:12:33 +0200 Subject: [PATCH 10/36] [assembly][wmpfa] Add WMpfa discretization to assembly --- dumux/assembly/jacobianpattern.hh | 3 ++- dumux/assembly/partialreassembler.hh | 13 +++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/dumux/assembly/jacobianpattern.hh b/dumux/assembly/jacobianpattern.hh index d3be1809c7..66f2f0ca39 100644 --- a/dumux/assembly/jacobianpattern.hh +++ b/dumux/assembly/jacobianpattern.hh @@ -85,7 +85,8 @@ Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry) */ template = 0> + || (GridGeometry::discMethod == DiscretizationMethod::ccmpfa) + || (GridGeometry::discMethod == DiscretizationMethod::ccwmpfa)), int> = 0> Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry) { const auto numDofs = gridGeometry.numDofs(); diff --git a/dumux/assembly/partialreassembler.hh b/dumux/assembly/partialreassembler.hh index 52defa5d3c..94e73ee60b 100644 --- a/dumux/assembly/partialreassembler.hh +++ b/dumux/assembly/partialreassembler.hh @@ -407,6 +407,19 @@ public: using ParentType::ParentType; }; +/*! + * \ingroup Assembly + * \brief The partial reassembler engine specialized for the cellcentered WMPFA method + */ +template +class PartialReassemblerEngine +: public PartialReassemblerEngine +{ + using ParentType = PartialReassemblerEngine; +public: + using ParentType::ParentType; +}; + //! helper struct to determine whether the an engine class has vertex colors struct hasVertexColor { -- GitLab From d0c4853eefabfdc96b1365865a3a3a97e4f58c95 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 18:13:02 +0200 Subject: [PATCH 11/36] [disc][fluxstencil] Draft implementation of ccwmpfa fluxstencil Currently all neighbors of neihbors are added, which is not really efficient. --- dumux/discretization/fluxstencil.hh | 56 +++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh index f671e143b3..28ac968c81 100644 --- a/dumux/discretization/fluxstencil.hh +++ b/dumux/discretization/fluxstencil.hh @@ -123,6 +123,62 @@ public: } }; +/* + * \ingroup Discretization + * \brief Flux stencil specialization for the cell-centered wmpfa scheme + * \tparam FVElementGeometry The local view on the finite volume grid geometry + */ +template +class FluxStencil +{ + using GridGeometry = typename FVElementGeometry::GridGeometry; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using GridIndexType = typename IndexTraits::GridIndex; + +public: + //! We don't know yet how many faces couple to a neighboring element + //! Each cell I couples to a cell J always only via one face + using ScvfStencilIForJ = std::vector; + + //! The flux stencil type + using Stencil = typename SubControlVolumeFace::Traits::GridIndexStorage; + + //! Returns the indices of the elements required for flux calculation on an scvf. + static Stencil stencil(const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + const auto& gridGeometry = fvGeometry.gridGeometry(); + + //ToDo correct stencil! At the moment simply all neighors of neighbors are added + Stencil stencil({scvf.insideScvIdx()}); + for (const auto& scvf : scvfs(fvGeometry)) + { + if (scvf.boundary()) + continue; + + stencil.push_back(scvf.outsideScvIdx()); + + const auto outsideScvIdx = scvf.outsideScvIdx(); + const auto outsideElement = gridGeometry.element(outsideScvIdx); + auto fvGeometryJ = localView(gridGeometry); + fvGeometryJ.bind(outsideElement); + + for (const auto& scvfJ : scvfs(fvGeometryJ)) + { + if (scvfJ.boundary() || scvf.insideScvIdx() == scvfJ.outsideScvIdx()) + continue; + + stencil.push_back(scvfJ.outsideScvIdx()); + } + } + + return stencil; + } +}; + } // end namespace Dumux #endif -- GitLab From 920c9e82d8837359d1fc58e77da971b03c0fc113 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 18:15:12 +0200 Subject: [PATCH 12/36] [disc][wmpfa] First draft of fluxvariablescache --- dumux/porousmediumflow/fluxvariablescache.hh | 12 ++ .../fluxvariablescachefiller.hh | 187 ++++++++++++++++++ 2 files changed, 199 insertions(+) diff --git a/dumux/porousmediumflow/fluxvariablescache.hh b/dumux/porousmediumflow/fluxvariablescache.hh index 892d251b17..77390de5e8 100644 --- a/dumux/porousmediumflow/fluxvariablescache.hh +++ b/dumux/porousmediumflow/fluxvariablescache.hh @@ -144,6 +144,18 @@ private: unsigned int idxInOutsideFaces_; // index of scvf among outside scvfs of iv-local "positive" face (only surface grids) }; +// specialization for the cell centered wmpfa method +template +class PorousMediumFluxVariablesCacheImplementation +: public AdvectionCacheChooser::enableAdvection()> +, public DiffusionCacheChooser::enableMolecularDiffusion()> +, public EnergyCacheChooser::enableEnergyBalance()> +{ +public: + //! export type used for scalar values + using Scalar = GetPropType; +}; + } // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index 7091268cd5..a085f6349e 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -702,6 +702,193 @@ private: const SecondaryLocalFaceData* secondaryLocalFaceData_; }; +//! Specialization of the flux variables cache filler for the cell centered mpfa method +template +class PorousMediumFluxVariablesCacheFillerImplementation +{ + using ModelTraits = GetPropType; + using Problem = GetPropType; + using GridView = typename GetPropType::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + using GridGeometry = GetPropType; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using ElementVolumeVariables = typename GetPropType::LocalView; + using ElementFluxVariablesCache = typename GetPropType::LocalView; + + using InterpolationOperator = typename ElementFluxVariablesCache::InterpolationOperator; + using FaceDataHandle = typename ElementFluxVariablesCache::FaceDataHandle; + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + + static constexpr bool advectionEnabled = ModelTraits::enableAdvection(); + static constexpr bool diffusionEnabled = ModelTraits::enableMolecularDiffusion(); + static constexpr bool heatConductionEnabled = ModelTraits::enableEnergyBalance(); + + static constexpr bool advectionIsSolDependent = getPropValue(); + static constexpr bool diffusionIsSolDependent = getPropValue(); + static constexpr bool heatConductionIsSolDependent = getPropValue(); + +public: + //! This cache filler is always solution-dependent, as it updates the + //! vectors of cell unknowns with which the transmissibilities have to be + //! multiplied in order to obtain the fluxes. + static constexpr bool isSolDependent = true; + + //! The constructor. Sets problem pointer. + PorousMediumFluxVariablesCacheFillerImplementation(const Problem& problem) + : problemPtr_(&problem) {} + + /*! + * \brief function to fill the flux variables caches + * + * \param fluxVarsCacheStorage Class that holds the scvf flux vars caches + * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf + * \param DataStorage Class that stores the interpolation operator & handles + * \param element The element + * \param fvGeometry The finite volume geometry + * \param elemVolVars The element volume variables (primary/secondary variables) + * \param scvf The corresponding sub-control volume face + * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones) + */ + template + void fill(FluxVarsCacheStorage& fluxVarsCacheStorage, + FluxVariablesCache& scvfFluxVarsCache, + DataStorage& dataStorage, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + bool forceUpdateAll = false) + { + // Set pointers + const auto& gridGeometry = fvGeometry.gridGeometry(); + + const auto& intOp = fluxVarsCacheStorage.interpolationOperator(fvGeometry.scv(scvf.insideScvIdx())); + faceHandlePtr_ = &dataStorage.facesDataHandles[scvf.index()]; + + // fill the physics-related quantities of the caches + if (forceUpdateAll) + { + if constexpr (advectionEnabled) + { + // lambda to obtain the permeability tensor + auto getTensor = [&elemVolVars] (typename GridView::IndexSet::IndexType index) { return elemVolVars[index].permeability(); }; + prepareHandle_(faceHandlePtr_->advectionHandle(), intOp.advectionInterpolator(), getTensor, fvGeometry, scvf); + fillAdvection_(scvfFluxVarsCache, faceHandlePtr_->advectionHandle(), intOp.advectionInterpolator(), element, fvGeometry, elemVolVars, scvf); + } + // if constexpr (diffusionEnabled) + // { + // fillDiffusion_(scvfFluxVarsCache, intOp, faceDataHandle(), element, fvGeometry, elemVolVars, scvf); + // } + // if constexpr (heatConductionEnabled) + // { + // fillHeatConduction_(scvfFluxVarsCache, intOp, faceDataHandle(), element, fvGeometry, elemVolVars, scvf); + // } + } + else + { + if constexpr (advectionEnabled && advectionIsSolDependent) + { + // lambda to obtain the permeability tensor + auto getTensor = [&elemVolVars] (typename GridView::IndexSet::IndexType index) { return elemVolVars[index].permeability(); }; + prepareHandle_(faceHandlePtr_->advectionHandle(), intOp.advectionInterpolator(), getTensor, fvGeometry, scvf); + fillAdvection_(scvfFluxVarsCache, faceHandlePtr_->advectionHandle(), intOp.advectionInterpolator(), element, fvGeometry, elemVolVars, scvf); + } + // if constexpr (diffusionEnabled && diffusionIsSolDependent) + // { + // fillDiffusion_(scvfFluxVarsCache, intOp, faceDataHandle().diffusionHandle(), element, fvGeometry, elemVolVars, scvf); + // } + // if constexpr (heatConductionEnabled && heatConductionIsSolDependent) + // { + // fillHeatConduction_(scvfFluxVarsCache, intOp, faceDataHandle(), element, fvGeometry, elemVolVars, scvf); + // } + } + } + + const FaceDataHandle& faceDataHandle() const + { + return *faceHandlePtr_; + } + +private: + template< class Handle, class IntOP, class TensorFunction > + void prepareHandle_(Handle& handle, const IntOP& intOp, const TensorFunction& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + { + handle.decompose(intOp, tensor, fvGeometry, scvf); + } + + //! method to fill the advective quantities + template + void fillAdvection_(FluxVariablesCache& scvfFluxVarsCache, + AdvectionHandle& handle, + const AdvectionInterpolator& intOp, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) + { + // using AdvectionType = GetPropType; + // using AdvectionFiller = typename AdvectionType::Cache::Filler; + + // // forward to the filler for the advective quantities + // AdvectionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this); + } + + // //! method to fill the diffusive quantities + // template + // void fillDiffusion_(FluxVariablesCache& scvfFluxVarsCache, + // typename FaceDataHandle::DiffusionHandle& handle, + // const typename InterpolationOperator::DiffusionInterpolator& intOp, + // const Element& element, + // const FVElementGeometry& fvGeometry, + // const ElementVolumeVariables& elemVolVars, + // const SubControlVolumeFace& scvf) + // { + // // using DiffusionType = GetPropType; + // // using DiffusionFiller = typename DiffusionType::Cache::Filler; + // // using FluidSystem = GetPropType; + + // // static constexpr int numPhases = ModelTraits::numFluidPhases(); + // // static constexpr int numComponents = ModelTraits::numFluidComponents(); + + // // // forward to the filler of the diffusive quantities + // // if constexpr (FluidSystem::isTracerFluidSystem()) + // // for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + // // for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx) + // // DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this); + // // else + // // for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + // // for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx) + // // if (compIdx != FluidSystem::getMainComponent(phaseIdx)) + // // DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this); + // } + + // //! method to fill the quantities related to heat conduction + // template + // void fillHeatConduction_(FluxVariablesCache& scvfFluxVarsCache, + // typename FaceDataHandle::HeatConductionHandle& handle, + // const typename InterpolationOperator::HeatConductionInterpolator& intOp, + // const Element& element, + // const FVElementGeometry& fvGeometry, + // const ElementVolumeVariables& elemVolVars, + // const SubControlVolumeFace& scvf) + // { + // // using HeatConductionType = GetPropType; + // // using HeatConductionFiller = typename HeatConductionType::Cache::Filler; + + // // // forward to the filler of the diffusive quantities + // // HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this); + // } + + const Problem& problem() const { return *problemPtr_; } + + const Problem* problemPtr_; + FaceDataHandle* faceHandlePtr_; +}; + } // end namespace Dumux #endif -- GitLab From d3c9070a9b382e0a3d34046212d2e93979ba4f76 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 1 Jul 2020 18:15:47 +0200 Subject: [PATCH 13/36] [disc][wmpfa] Add ccwmpfa as possible discretization method --- dumux/discretization/cellcentered/elementsolution.hh | 12 ++++++++---- dumux/discretization/method.hh | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/dumux/discretization/cellcentered/elementsolution.hh b/dumux/discretization/cellcentered/elementsolution.hh index 5aabd0934e..8ffaafdf94 100644 --- a/dumux/discretization/cellcentered/elementsolution.hh +++ b/dumux/discretization/cellcentered/elementsolution.hh @@ -112,7 +112,8 @@ private: template auto elementSolution(const Element& element, const SolutionVector& sol, const GridGeometry& gg) -> std::enable_if_t()[0])>> > @@ -128,7 +129,8 @@ auto elementSolution(const Element& element, const SolutionVector& sol, const Gr template auto elementSolution(const Element& element, const ElementVolumeVariables& elemVolVars, const FVElementGeometry& gg) -> std::enable_if_t> { using PrimaryVariables = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables; @@ -143,7 +145,8 @@ auto elementSolution(const Element& element, const ElementVolumeVariables& elemV template auto elementSolution(PrimaryVariables&& priVars) -> std::enable_if_t> { return CCElementSolution(std::move(priVars)); @@ -157,7 +160,8 @@ auto elementSolution(PrimaryVariables&& priVars) template auto elementSolution(const PrimaryVariables& priVars) -> std::enable_if_t> { return CCElementSolution(priVars); diff --git a/dumux/discretization/method.hh b/dumux/discretization/method.hh index 15478a56bd..92221391ab 100644 --- a/dumux/discretization/method.hh +++ b/dumux/discretization/method.hh @@ -35,7 +35,7 @@ namespace Dumux { */ enum class DiscretizationMethod { - none, box, cctpfa, ccmpfa, staggered, fem + none, box, cctpfa, ccmpfa, ccwmpfa, staggered, fem }; } // end namespace Dumux -- GitLab From 2a6f4b3503a411b262cc3f6bc36190b1e5d75819 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 2 Jul 2020 16:43:31 +0200 Subject: [PATCH 14/36] [disc][wmpfa] More appropriate face data handle --- .../cellcentered/wmpfa/facedatahandle.hh | 72 ++++++++++++++----- dumux/flux/darcyslaw.hh | 1 + dumux/flux/upwindscheme.hh | 5 ++ 3 files changed, 60 insertions(+), 18 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index aefdc2be1c..92f1c0c1f5 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -61,19 +61,32 @@ template class AdvectionDataHandle { using Scalar = typename IntOp::Scalar; - using OneSidedStencil = std::vector; - using Stencils = std::pair; - using OneSidedCoefficients = std::vector; - using Coefficients = std::pair; + //ToDo get the index types + using SubControlVolume = typename IntOp::GridGeometry::SubControlVolume; + using GridIndexType = typename SubControlVolume::Traits::GridIndexType; + using LocalIndexType = typename SubControlVolume::Traits::LocalIndexType; + using Position = typename IntOp::Position; + + struct Entry{ + Scalar coefficient; + GridIndexType index; + Position position; + }; + using Entries = std::vector; using Interpolator = typename IntOp::AdvectionInterpolator; static constexpr int numPhases = PhysicsTraits::numPhases; public: + void clear() + { + visited_ = 0; + entries_.clear(); + } - template - void decompose(const Interpolator& intOp, const TF& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + template + void decompose(const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf) { boundaryFace_ = scvf.boundary(); const auto coNormal = mv(tensor(scvf.insideScvIdx()), scvf.unitOuterNormal()); @@ -84,29 +97,52 @@ public: else { WMpfaHelper::eraseZeros(coeff, indices); + updateEntries(intOp, indices, coeff, fvGeometry, scvf); + } - if(scvf.insideScvIdx() < scvf.outsideScvIdx()) - { - coefficients_[0].first = std::move(coeff); - stencils_[0].first = std::move(indices); - } - else + ++visited_; + } + + template + void updateEntries(const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf) + { + const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); + auto& entries = entries_[visited_]; + entries.push_back({0.0, scv.dofIndex(), scv.center()}); + for(std::size_t i=0; i 0 && coefficients_[0].second.size() > 0) - || (boundaryFace_ && coefficients_[0].first.size() > 0); + return ((!boundaryFace_ && entries_[0].size() > 0 && entries_[1].size() > 0) + || (boundaryFace_ && entries_[0].size() > 0)) && visited_ == 2; } private: - std::array coefficients_; - std::array stencils_; + int visited_; + std::array scvfIndices_; + std::array scvIndices_; + std::array entries_; bool boundaryFace_; }; diff --git a/dumux/flux/darcyslaw.hh b/dumux/flux/darcyslaw.hh index de930e5604..4966fe4f16 100644 --- a/dumux/flux/darcyslaw.hh +++ b/dumux/flux/darcyslaw.hh @@ -51,5 +51,6 @@ using DarcysLaw = DarcysLawImplementation #include #include +#include #endif diff --git a/dumux/flux/upwindscheme.hh b/dumux/flux/upwindscheme.hh index 3fdf09d943..866382303e 100644 --- a/dumux/flux/upwindscheme.hh +++ b/dumux/flux/upwindscheme.hh @@ -190,6 +190,11 @@ template class UpwindSchemeImpl : public UpwindSchemeImpl {}; +//! Upwind scheme for cell-centered wmpfa schemes +template +class UpwindSchemeImpl +: public UpwindSchemeImpl {}; + } // end namespace Dumux #endif -- GitLab From 93b4bb34d930de8e2bb21e4d4fa08c93862c285a Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 2 Jul 2020 16:44:18 +0200 Subject: [PATCH 15/36] [disc][wmpfa] More general interpolators --- .../wmpfa/interpolationoperator.hh | 63 ++++++++++++++++--- 1 file changed, 56 insertions(+), 7 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh index 183b11cc66..30757172b7 100644 --- a/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh +++ b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh @@ -46,10 +46,45 @@ class HapInterpolatorBase using Scalar = typename Position::value_type; + struct Entry + { + Entry() + {entry = {-1,0.0};} + + Entry(std::size_t idx, Scalar weight) + { + entry = {idx, weight}; + } + + Scalar weight() const + {return entry.second;} + + std::size_t dofIndex() const + {return entry.first;} + + private: + std::pair entry; + }; + struct LocalInterpolationData { - Position point; - std::pair weights; + LocalInterpolationData() = default; + + LocalInterpolationData(Position p, Entry e1, Entry e2) + { + position_ = p; + entries_ = {e1,e2}; + } + + const std::array& entries() const + { return entries_;} + + Position position() const + {return position_;} + + private: + Position position_ = {}; + std::array entries_ = {Entry(), Entry()}; }; public: @@ -70,6 +105,8 @@ public: interpolationData_.resize(fvGeometry.numScvf()); for (auto&& scvf : scvfs(fvGeometry)) { + const auto insideScvIdx = scvf.insideScvIdx(); + const auto outsideScvIdx = scvf.outsideScvIdx(); if(!scvf.boundary()) { const auto insideScvIdx = scvf.insideScvIdx(); @@ -89,15 +126,15 @@ public: const auto omegaK = distL*tauK / (distL*tauK + distK*tauL); const auto omegaL = distK*tauL / (distL*tauK + distK*tauL); - auto point = (omegaK * centerK) + (omegaL * centerL) + auto position = (omegaK * centerK) + (omegaL * centerL) + (distL*distK / (distL*tauK + distK*tauL)) * mv(tensor(insideVolVars) - tensor(outsideVolVars), scvf.unitOuterNormal()); - interpolationData_[scvf.localIndex()] = (LocalInterpolationData({point, std::make_pair(omegaK, omegaL)})); + interpolationData_[scvf.localIndex()] = LocalInterpolationData(position, Entry(insideScvIdx, omegaK), Entry(outsideScvIdx, omegaL) ); } else { - interpolationData_[scvf.localIndex()] = (LocalInterpolationData({scvf.center(), std::make_pair(1.0, 0.0)})); + interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.center(), Entry(insideScvIdx, 0.0), Entry(outsideScvIdx, 1.0) ); } } isUpdated_ = true; @@ -114,7 +151,7 @@ public: std::vector distances; distances.resize(fvGeometry.numScvf()); for (auto&& scvf : scvfs(fvGeometry)) - distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].point - fvGeometry.scv(scvf.insideScvIdx()).center()); + distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].position() - fvGeometry.scv(scvf.insideScvIdx()).center()); return distances; } @@ -125,11 +162,21 @@ public: std::vector distances; distances.resize(fvGeometry.numScvf()); for (auto&& scvf : scvfs(fvGeometry)) - distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].point - fvGeometry.scv(scvf.insideScvIdx()).center()); + distances[scvf.localIndex()]= (interpolationData_[scvf.localIndex()].position() - fvGeometry.scv(scvf.insideScvIdx()).center()); return distances; } + const LocalInterpolationData getInterpolationData(std::size_t localIdx) const + { + return interpolationData_[localIdx]; + } + + LocalInterpolationData getInterpolationData(std::size_t localIdx) + { + return interpolationData_[localIdx]; + } + private: std::vector interpolationData_; @@ -146,6 +193,8 @@ class HapInterpolationOperator public: //! state the traits type publicly using Traits = T; + using GridGeometry = typename T::GridGeometry; + using GridView = typename T::GridView; using Position = typename T::GlobalPosition; using Scalar = typename Position::value_type; -- GitLab From 08d303cd3f03f0e43dcbbd96a24fe40360caefbd Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 9 Jul 2020 14:21:31 +0200 Subject: [PATCH 16/36] [ccwmpfa] Minor fixes --- .../cellcentered/wmpfa/facedatahandle.hh | 20 ++++++++++++++++--- .../cellcentered/wmpfa/fvgridgeometry.hh | 2 +- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index 92f1c0c1f5..32e5370f85 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -82,7 +82,14 @@ public: void clear() { visited_ = 0; - entries_.clear(); + entries_[0].clear(); + entries_[1].clear(); + } + + void prepare() + { + if(visited_ == 2 || boundaryFace_) + clear(); } template @@ -116,7 +123,7 @@ public: { if(e.dofIndex() != scvf.insideScvIdx()) { - auto c = coeff[scvf.localIndex()]*e.weight(); + auto c = coeff[i]*e.weight(); entries[0].coefficient += c; //ToDo pos is not needed for each entry!!! entries.push_back({c, e.dofIndex(), intData.position()}); @@ -125,6 +132,13 @@ public: } scvfIndices_[visited_] = scvf.index(); scvIndices_[visited_] = scvf.insideScvIdx(); + + if(boundaryFace_) + { + ++visited_; + scvfIndices_[visited_] = scvf.index(); + scvIndices_[visited_] = scvf.outsideScvIdx(); + } } const auto& subFluxData(GridIndexType idx) const @@ -143,7 +157,7 @@ private: std::array scvfIndices_; std::array scvIndices_; std::array entries_; - bool boundaryFace_; + bool boundaryFace_ = {false}; }; //! Data handle for quantities related to diffusion diff --git a/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh index bde9384594..dd079ec776 100644 --- a/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh +++ b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh @@ -112,7 +112,7 @@ public: static constexpr DiscretizationMethod discMethod = DiscretizationMethod::ccwmpfa; //! The maximum admissible stencil size (used for static memory allocation during assembly) - static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1; + static constexpr int maxElementStencilSize = LocalView::maxNumElementScvfs*LocalView::maxNumElementScvfs*Traits::maxNumScvfNeighbors + 1; //! export the grid view type using GridView = GV; -- GitLab From 54e9f2e4bcd40baeeee0afa3194438cc90178b87 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 9 Jul 2020 14:22:24 +0200 Subject: [PATCH 17/36] [ccwmpfa] Update stencils and related boundary volVars --- .../wmpfa/elementvolumevariables.hh | 45 ++++++++++++------- dumux/discretization/fluxstencil.hh | 9 ++-- .../fluxvariablescachefiller.hh | 16 +++---- 3 files changed, 44 insertions(+), 26 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh index 4618ee503f..67f54c2ceb 100644 --- a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh @@ -112,6 +112,17 @@ namespace CCWMpfa { volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(scvf.outsideScvIdx()); } + else + { + VolumeVariables dirichletVolVars; + dirichletVolVars.update(elementSolution(problem.dirichlet(element, scvf)), + problem, + element, + scvI); + + volVars.emplace_back(std::move(dirichletVolVars)); + volVarIndices.push_back(scvf.outsideScvIdx()); + } } } return std::make_pair(volVars, volVarIndices); @@ -175,30 +186,30 @@ public: const FVElementGeometry& fvGeometry, const SolutionVector& sol) { - if (!fvGeometry.hasBoundaryScvf()) - return; - clear_(); - auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), element, fvGeometry); - std::move(volVars.begin(), volVars.end(), std::back_inserter(boundaryVolumeVariables_)); - std::move(indices.begin(), indices.end(), std::back_inserter(boundaryVolVarIndices_)); + if (fvGeometry.hasBoundaryScvf()) + { + auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), element, fvGeometry); + std::move(volVars.begin(), volVars.end(), std::back_inserter(boundaryVolumeVariables_)); + std::move(indices.begin(), indices.end(), std::back_inserter(boundaryVolVarIndices_)); + } + + const auto& gridGeometry = fvGeometry.gridGeometry(); + const auto globalI = gridGeometry.elementMapper().index(element); + const auto& assemblyMapI = gridGeometry.connectivityMap()[globalI]; // add boundary volVars of neighbors - for (const auto& scvf : scvfs(fvGeometry)) + for (const auto& dataJ : assemblyMapI) { - if (scvf.boundary()) - continue; - - const auto outsideScvIdx = scvf.outsideScvIdx(); - const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx); - auto fvGeometryJ = localView(fvGeometry.gridGeometry()); - fvGeometryJ.bind(outsideElement); + const auto& elementJ = gridGeometry.element(dataJ.globalJ); + auto fvGeometryJ = localView(gridGeometry); + fvGeometryJ.bind(elementJ); if (!fvGeometryJ.hasBoundaryScvf()) continue; - auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), outsideElement, fvGeometryJ); + auto&& [volVars, indices] = CCWMpfa::boundaryVolVars(gridVolVars().problem(), elementJ, fvGeometryJ); std::move(volVars.begin(), volVars.end(), std::back_inserter(boundaryVolumeVariables_)); std::move(indices.begin(), indices.end(), std::back_inserter(boundaryVolVarIndices_)); } @@ -229,6 +240,10 @@ private: int getLocalIdx_(const int volVarIdx) const { auto it = std::find(boundaryVolVarIndices_.begin(), boundaryVolVarIndices_.end(), volVarIdx); + if(it == boundaryVolVarIndices_.end()) + { + int shit = 1; + } assert(it != boundaryVolVarIndices_.end() && "Could not find the current volume variables for volVarIdx!"); return std::distance(boundaryVolVarIndices_.begin(), it); } diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh index 28ac968c81..80fde0ccc5 100644 --- a/dumux/discretization/fluxstencil.hh +++ b/dumux/discretization/fluxstencil.hh @@ -147,12 +147,12 @@ public: //! Returns the indices of the elements required for flux calculation on an scvf. static Stencil stencil(const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) { const auto& gridGeometry = fvGeometry.gridGeometry(); - //ToDo correct stencil! At the moment simply all neighors of neighbors are added + //ToDo correct stencil! At the moment simply all neighors of the neighbor are added Stencil stencil({scvf.insideScvIdx()}); for (const auto& scvf : scvfs(fvGeometry)) { @@ -160,7 +160,10 @@ public: continue; stencil.push_back(scvf.outsideScvIdx()); + } + if (!scvf.boundary()) + { const auto outsideScvIdx = scvf.outsideScvIdx(); const auto outsideElement = gridGeometry.element(outsideScvIdx); auto fvGeometryJ = localView(gridGeometry); diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index a085f6349e..a21773189e 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -732,10 +732,9 @@ class PorousMediumFluxVariablesCacheFillerImplementation(); public: - //! This cache filler is always solution-dependent, as it updates the - //! vectors of cell unknowns with which the transmissibilities have to be - //! multiplied in order to obtain the fluxes. - static constexpr bool isSolDependent = true; + static constexpr bool isSolDependent = (advectionEnabled && advectionIsSolDependent) || + (diffusionEnabled && diffusionIsSolDependent) || + (heatConductionEnabled && heatConductionIsSolDependent); //! The constructor. Sets problem pointer. PorousMediumFluxVariablesCacheFillerImplementation(const Problem& problem) @@ -817,6 +816,7 @@ private: template< class Handle, class IntOP, class TensorFunction > void prepareHandle_(Handle& handle, const IntOP& intOp, const TensorFunction& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) { + handle.prepare(); handle.decompose(intOp, tensor, fvGeometry, scvf); } @@ -830,11 +830,11 @@ private: const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) { - // using AdvectionType = GetPropType; - // using AdvectionFiller = typename AdvectionType::Cache::Filler; + using AdvectionType = GetPropType; + using AdvectionFiller = typename AdvectionType::Cache::Filler; - // // forward to the filler for the advective quantities - // AdvectionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this); + // forward to the filler for the advective quantities + AdvectionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this); } // //! method to fill the diffusive quantities -- GitLab From 1b03b307b5731a65fcbd5d50e6bc8283324afa9e Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 09:54:16 +0200 Subject: [PATCH 18/36] [wmpfa] Only store one-sided information in facedatahandle --- .../wmpfa/elementvolumevariables.hh | 4 -- .../cellcentered/wmpfa/facedatahandle.hh | 44 +++++++------------ 2 files changed, 15 insertions(+), 33 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh index 67f54c2ceb..7ee0572a7d 100644 --- a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh @@ -240,10 +240,6 @@ private: int getLocalIdx_(const int volVarIdx) const { auto it = std::find(boundaryVolVarIndices_.begin(), boundaryVolVarIndices_.end(), volVarIdx); - if(it == boundaryVolVarIndices_.end()) - { - int shit = 1; - } assert(it != boundaryVolVarIndices_.end() && "Could not find the current volume variables for volVarIdx!"); return std::distance(boundaryVolVarIndices_.begin(), it); } diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index 32e5370f85..62e9ff6cf2 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -72,24 +72,22 @@ class AdvectionDataHandle GridIndexType index; Position position; }; + +public: using Entries = std::vector; using Interpolator = typename IntOp::AdvectionInterpolator; static constexpr int numPhases = PhysicsTraits::numPhases; -public: void clear() { - visited_ = 0; - entries_[0].clear(); - entries_[1].clear(); + entries_.clear(); } void prepare() { - if(visited_ == 2 || boundaryFace_) - clear(); + clear(); } template @@ -106,16 +104,13 @@ public: WMpfaHelper::eraseZeros(coeff, indices); updateEntries(intOp, indices, coeff, fvGeometry, scvf); } - - ++visited_; } template void updateEntries(const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf) { const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); - auto& entries = entries_[visited_]; - entries.push_back({0.0, scv.dofIndex(), scv.center()}); + entries_.push_back({0.0, scv.dofIndex(), scv.center()}); for(std::size_t i=0; i 0 && entries_[1].size() > 0) - || (boundaryFace_ && entries_[0].size() > 0)) && visited_ == 2; + return entries_.size() > 0; } private: - int visited_; - std::array scvfIndices_; - std::array scvIndices_; - std::array entries_; + GridIndexType scvfIndices_; + GridIndexType scvIndices_; + Entries entries_; bool boundaryFace_ = {false}; }; -- GitLab From 7523aa3751b3eab61568168b13f2df25c2439663 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 09:55:58 +0200 Subject: [PATCH 19/36] [wmpfa] Always make the flip index set --- .../cellcentered/wmpfa/fvgridgeometry.hh | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh index dd079ec776..c8aa116d97 100644 --- a/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh +++ b/dumux/discretization/cellcentered/wmpfa/fvgridgeometry.hh @@ -239,21 +239,18 @@ public: scvfIndicesOfScv_[eIdx] = scvfsIndexSet; } - // Make the flip index set for network, surface, and periodic grids - if (this->isPeriodic()) + // Make the flip index set + flipScvfIndices_.resize(scvfs_.size()); + for (auto&& scvf : scvfs_) { - flipScvfIndices_.resize(scvfs_.size()); - for (auto&& scvf : scvfs_) - { - if (scvf.boundary()) - continue; - - flipScvfIndices_[scvf.index()].resize(scvf.numOutsideScvs()); - const auto insideScvIdx = scvf.insideScvIdx(); - // check which outside scvf has the insideScvIdx index in its outsideScvIndices - for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) - flipScvfIndices_[scvf.index()][i] = findFlippedScvfIndex_(insideScvIdx, scvf.outsideScvIdx(i)); - } + if (scvf.boundary()) + continue; + + flipScvfIndices_[scvf.index()].resize(scvf.numOutsideScvs()); + const auto insideScvIdx = scvf.insideScvIdx(); + // check which outside scvf has the insideScvIdx index in its outsideScvIndices + for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) + flipScvfIndices_[scvf.index()][i] = findFlippedScvfIndex_(insideScvIdx, scvf.outsideScvIdx(i)); } // build the connectivity map for an effecient assembly -- GitLab From ca9b6d64ece14470e87cd2f867ab275c9643e039 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 10:03:05 +0200 Subject: [PATCH 20/36] [wmpfa] Use correct index type for interpolator --- .../cellcentered/wmpfa/interpolationoperator.hh | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh index 30757172b7..2dc88703d4 100644 --- a/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh +++ b/dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh @@ -40,10 +40,13 @@ class EmptyInterpolator {}; template< class T, bool Enable> class HapInterpolatorBase { + using GridGeometry = typename T::GridGeometry; using GridView = typename T::GridView; using Element = typename GridView::template Codim<0>::Entity; using Position = typename T::GlobalPosition; + using GridIndexType = typename GridGeometry::SubControlVolume::Traits::GridIndexType; + using Scalar = typename Position::value_type; struct Entry @@ -51,7 +54,7 @@ class HapInterpolatorBase Entry() {entry = {-1,0.0};} - Entry(std::size_t idx, Scalar weight) + Entry(GridIndexType idx, Scalar weight) { entry = {idx, weight}; } @@ -59,11 +62,11 @@ class HapInterpolatorBase Scalar weight() const {return entry.second;} - std::size_t dofIndex() const + GridIndexType dofIndex() const {return entry.first;} private: - std::pair entry; + std::pair entry; }; struct LocalInterpolationData @@ -167,12 +170,12 @@ public: return distances; } - const LocalInterpolationData getInterpolationData(std::size_t localIdx) const + const LocalInterpolationData getInterpolationData(GridIndexType localIdx) const { return interpolationData_[localIdx]; } - LocalInterpolationData getInterpolationData(std::size_t localIdx) + LocalInterpolationData getInterpolationData(GridIndexType localIdx) { return interpolationData_[localIdx]; } -- GitLab From 1a09cb4b08047637bf4f6655249cd32d86a112f8 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 10:29:06 +0200 Subject: [PATCH 21/36] [wmpfa] Unused variables --- dumux/porousmediumflow/fluxvariablescachefiller.hh | 3 --- 1 file changed, 3 deletions(-) diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index a21773189e..adef8995ea 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -762,9 +762,6 @@ public: const SubControlVolumeFace& scvf, bool forceUpdateAll = false) { - // Set pointers - const auto& gridGeometry = fvGeometry.gridGeometry(); - const auto& intOp = fluxVarsCacheStorage.interpolationOperator(fvGeometry.scv(scvf.insideScvIdx())); faceHandlePtr_ = &dataStorage.facesDataHandles[scvf.index()]; -- GitLab From a7aab6e854ca59a7d1e1d95b2677daebf98ec68d Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 10:30:41 +0200 Subject: [PATCH 22/36] [flux][wmpfa] Implementation of darcyslaw --- dumux/flux/ccwmpfa/CMakeLists.txt | 3 + dumux/flux/ccwmpfa/darcyslaw.hh | 237 ++++++++++++++++++++++++++++++ 2 files changed, 240 insertions(+) create mode 100644 dumux/flux/ccwmpfa/CMakeLists.txt create mode 100644 dumux/flux/ccwmpfa/darcyslaw.hh diff --git a/dumux/flux/ccwmpfa/CMakeLists.txt b/dumux/flux/ccwmpfa/CMakeLists.txt new file mode 100644 index 0000000000..2860457238 --- /dev/null +++ b/dumux/flux/ccwmpfa/CMakeLists.txt @@ -0,0 +1,3 @@ +install(FILES +darcyslaw.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/flux/ccwmpfa) diff --git a/dumux/flux/ccwmpfa/darcyslaw.hh b/dumux/flux/ccwmpfa/darcyslaw.hh new file mode 100644 index 0000000000..66fec35463 --- /dev/null +++ b/dumux/flux/ccwmpfa/darcyslaw.hh @@ -0,0 +1,237 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCWMpfaFlux + * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation + */ +#ifndef DUMUX_DISCRETIZATION_CC_WMPFA_DARCYS_LAW_HH +#define DUMUX_DISCRETIZATION_CC_WMPFA_DARCYS_LAW_HH + +#include +#include +#include + +#include +#include + +namespace Dumux { + +//! forward declaration of the method-specific implementation +template +class DarcysLawImplementation; + +/*! + * \ingroup CCWMpfaFlux + * \brief Class that fills the cache corresponding to wmpfa Darcy's Law + */ +template +class WMpfaDarcysLawCacheFiller +{ + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + +public: + //! Function to fill a WMpfaDarcysLawCache of a given scvf + //! This interface has to be met by any advection-related cache filler class + //! TODO: Probably get cache type out of the filler + template + static void fill(FluxVariablesCache& scvfFluxVarsCache, + const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const FluxVariablesCacheFiller& fluxVarsCacheFiller) + { + scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.faceDataHandle()); + } +}; + +/*! + * \ingroup CCWMpfaFlux + * \brief The cache corresponding to tpfa Darcy's Law + */ +template +class WMpfaDarcysLawCache +{ + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + + using DataHandle = typename ElementFluxVariablesCache::FaceDataHandle; + using AdvectionDataHandle = typename DataHandle::AdvectionHandle; + +public: + using Filler = WMpfaDarcysLawCacheFiller; + + void updateAdvection(const DataHandle& dataHandle) + { + handle_ = &dataHandle.advectionHandle(); + } + + const AdvectionDataHandle& dataHandle() const + { return *handle_; } + +private: + const AdvectionDataHandle* handle_; +}; + +/*! + * \ingroup CCWMpfaFlux + * \brief Specialization of the CCWMpfaDarcysLaw grids where dim=dimWorld + */ +template +class DarcysLawImplementation +{ + using Scalar = GetPropType; + using Problem = GetPropType; + using GridView = typename GetPropType::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + + using GridGeometry = GetPropType; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using ElementVolumeVariables = typename GetPropType::LocalView; + + using ElementFluxVariablesCache = typename GetPropType::LocalView; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + public: + //! state the discretization method this implementation belongs to + static const DiscretizationMethod discMethod = DiscretizationMethod::ccwmpfa; + + //! state the type for the corresponding cache + using Cache = WMpfaDarcysLawCache; + + //! Compute the advective flux + template + static Scalar flux(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + int phaseIdx, + const ElementFluxVarsCache& elemFluxVarsCache) + { + static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); + + const auto& fluxVarsCache = elemFluxVarsCache[scvf]; + const auto& dataHandle = fluxVarsCache.dataHandle(); + + // Get the inside and outside volume variables + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + const auto& insideVolVars = elemVolVars[insideScv]; + const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; + + const auto& subFluxDataIJ = dataHandle.subFluxData(); + + if (enableGravity) + { + // do averaging for the density over all neighboring elements + const auto rho = scvf.boundary() ? outsideVolVars.density(phaseIdx) + : (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5; + + const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); + + auto pseudoPotential = [&elemVolVars, phaseIdx, &g, &rho](const auto& volVars, const auto& x) + { return volVars.pressure(phaseIdx) - rho*(g*x);}; + + Scalar fluxIJ = 0.0; + Scalar wIJ = 0.5; + std::for_each(subFluxDataIJ.cbegin(), subFluxDataIJ.cend(), + [&fluxIJ, &elemVolVars, &pseudoPotential](const auto& e) + { fluxIJ += e.coefficient * pseudoPotential(elemVolVars[e.index], e.position); } ); + + Scalar fluxJI = 0.0; + Scalar wJI = 0.5; + if(!scvf.boundary()) + { + const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); + const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); + const auto& subFluxDataJI = dataHandleJ.subFluxData(); + std::for_each(subFluxDataJI.cbegin(), subFluxDataJI.cend(), + [&fluxJI, &elemVolVars, &pseudoPotential](const auto& e) + { fluxJI += e.coefficient * pseudoPotential(elemVolVars[e.index], e.position); } ); + } + else + { + wIJ = 1.0; + wJI = 0.0; + } + + return scvf.area() * (wIJ*fluxIJ - wJI*fluxJI); + } + else + { + auto pressure = [&elemVolVars, phaseIdx](const auto& volVars) + { return volVars.pressure(phaseIdx);}; + + Scalar fluxIJ = 0.0; + Scalar wIJ = 0.5; + std::for_each(subFluxDataIJ.cbegin(), subFluxDataIJ.cend(), + [&fluxIJ, &elemVolVars, &pressure](const auto& e) + { fluxIJ += e.coefficient * pressure(elemVolVars[e.index]); } ); + + Scalar fluxJI = 0.0; + Scalar wJI = 0.5; + if(!scvf.boundary()) + { + const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); + const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); + const auto& subFluxDataJI = dataHandleJ.subFluxData(); + std::for_each(subFluxDataJI.cbegin(), subFluxDataJI.cend(), + [&fluxJI, &elemVolVars, &pressure](const auto& e) + { fluxJI += e.coefficient * pressure(elemVolVars[e.index]); } ); + } + else + { + wIJ = 1.0; + wJI = 0.0; + } + + return scvf.area() * (wIJ*fluxIJ - wJI*fluxJI); + } + + return 0.0; + } + +private: + template = 0> + static decltype(auto) getPermeability_(const Problem& problem, + const VolumeVariables& volVars, + const GlobalPosition& scvfIpGlobal) + { return volVars.permeability(); } + + template = 0> + static decltype(auto) getPermeability_(const Problem& problem, + const VolumeVariables& volVars, + const GlobalPosition& scvfIpGlobal) + { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); } +}; + +} // end namespace Dumux + +#endif -- GitLab From 452452d44c43b38f7a1543235e479c5e3bf40f9e Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 15 Jul 2020 12:23:19 +0200 Subject: [PATCH 23/36] [wmpfa][neumann] Update handling of Neumann faces --- .../cellcentered/wmpfa/facedatahandle.hh | 37 +++++++++++++------ .../wmpfa/interpolationoperator.hh | 11 ++++-- .../fluxvariablescachefiller.hh | 2 +- 3 files changed, 35 insertions(+), 15 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index 62e9ff6cf2..d8698ed84e 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -90,8 +90,8 @@ public: clear(); } - template - void decompose(const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf) + template + void decompose(const Problem& problem, const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf) { boundaryFace_ = scvf.boundary(); const auto coNormal = mv(tensor(scvf.insideScvIdx()), scvf.unitOuterNormal()); @@ -102,26 +102,41 @@ public: else { WMpfaHelper::eraseZeros(coeff, indices); - updateEntries(intOp, indices, coeff, fvGeometry, scvf); + updateEntries(problem, intOp, indices, coeff, fvGeometry, scvf); } } - template - void updateEntries(const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf) + template + void updateEntries(const Problem& problem, const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf) { const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); entries_.push_back({0.0, scv.dofIndex(), scv.center()}); for(std::size_t i=0; i entries_ = {Entry(), Entry()}; }; @@ -133,11 +138,11 @@ public: + (distL*distK / (distL*tauK + distK*tauL)) * mv(tensor(insideVolVars) - tensor(outsideVolVars), scvf.unitOuterNormal()); - interpolationData_[scvf.localIndex()] = LocalInterpolationData(position, Entry(insideScvIdx, omegaK), Entry(outsideScvIdx, omegaL) ); + interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.index(), position, Entry(insideScvIdx, omegaK), Entry(outsideScvIdx, omegaL) ); } else { - interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.center(), Entry(insideScvIdx, 0.0), Entry(outsideScvIdx, 1.0) ); + interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.index(), scvf.center(), Entry(insideScvIdx, 0.0), Entry(outsideScvIdx, 1.0) ); } } isUpdated_ = true; diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh index adef8995ea..4ad7269d80 100644 --- a/dumux/porousmediumflow/fluxvariablescachefiller.hh +++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh @@ -814,7 +814,7 @@ private: void prepareHandle_(Handle& handle, const IntOP& intOp, const TensorFunction& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) { handle.prepare(); - handle.decompose(intOp, tensor, fvGeometry, scvf); + handle.decompose(problem(), intOp, tensor, fvGeometry, scvf); } //! method to fill the advective quantities -- GitLab From 58638b8d35b50d3cca40fc6dd910bf7047daedfd Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 16 Jul 2020 12:38:32 +0200 Subject: [PATCH 24/36] [flux][wmpfa] Allow for different methods in darcyslaw --- .../cellcentered/wmpfa/methods.hh | 4 +- dumux/flux/ccwmpfa/darcyslaw.hh | 113 ++++++++++-------- 2 files changed, 62 insertions(+), 55 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/methods.hh b/dumux/discretization/cellcentered/wmpfa/methods.hh index 29218bdca6..0b1c3aadf1 100644 --- a/dumux/discretization/cellcentered/wmpfa/methods.hh +++ b/dumux/discretization/cellcentered/wmpfa/methods.hh @@ -30,9 +30,9 @@ namespace Dumux { * \brief The available weighted mpfa schemes in Dumux * \ingroup CCWMpfaDiscretization */ - enum class WMpfaMethods : unsigned int + enum class WMpfaMethod { - avgMethod, nltpfa, nlmpfa + avgmpfa, nltpfa, nlmpfa }; } // end namespace Dumux diff --git a/dumux/flux/ccwmpfa/darcyslaw.hh b/dumux/flux/ccwmpfa/darcyslaw.hh index 66fec35463..977fb0cdfc 100644 --- a/dumux/flux/ccwmpfa/darcyslaw.hh +++ b/dumux/flux/ccwmpfa/darcyslaw.hh @@ -29,6 +29,7 @@ #include #include +#include #include namespace Dumux { @@ -94,6 +95,53 @@ private: const AdvectionDataHandle* handle_; }; + +template +class WMpfaDarcysLawsFluxCalculator; + +template<> +class WMpfaDarcysLawsFluxCalculator +{ +public: + + template + static void flux(Scalar& flux, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const FluxData& dJI, + const SolValues& solValues) + { + Scalar fluxIJ = 0.0; + Scalar fluxJI = 0.0; + Scalar wIJ = 0.5; + Scalar wJI = 0.5; + std::for_each(dIJ.cbegin(), dIJ.cend(), + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + + std::for_each(dJI.cbegin(), dJI.cend(), + [&fluxJI, &elemVolVars, &solValues](const auto& e) + { fluxJI += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + + flux = (wIJ*fluxIJ - wJI*fluxJI); + } + + template + static void boundaryFlux(Scalar& flux, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const SolValues& solValues) + { + Scalar fluxIJ = 0.0; + Scalar wIJ = 1.0; + std::for_each(dIJ.cbegin(), dIJ.cend(), + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + + flux = wIJ*fluxIJ; + } +}; + /*! * \ingroup CCWMpfaFlux * \brief Specialization of the CCWMpfaDarcysLaw grids where dim=dimWorld @@ -117,6 +165,8 @@ class DarcysLawImplementation using ElementFluxVariablesCache = typename GetPropType::LocalView; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using FluxCalculator = WMpfaDarcysLawsFluxCalculator; + public: //! state the discretization method this implementation belongs to static const DiscretizationMethod discMethod = DiscretizationMethod::ccwmpfa; @@ -146,6 +196,8 @@ class DarcysLawImplementation const auto& subFluxDataIJ = dataHandle.subFluxData(); + Scalar flux = 0.0; + if (enableGravity) { // do averaging for the density over all neighboring elements @@ -155,81 +207,36 @@ class DarcysLawImplementation const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); auto pseudoPotential = [&elemVolVars, phaseIdx, &g, &rho](const auto& volVars, const auto& x) - { return volVars.pressure(phaseIdx) - rho*(g*x);}; - - Scalar fluxIJ = 0.0; - Scalar wIJ = 0.5; - std::for_each(subFluxDataIJ.cbegin(), subFluxDataIJ.cend(), - [&fluxIJ, &elemVolVars, &pseudoPotential](const auto& e) - { fluxIJ += e.coefficient * pseudoPotential(elemVolVars[e.index], e.position); } ); + { return volVars.pressure(phaseIdx) - rho*(g*x);}; - Scalar fluxJI = 0.0; - Scalar wJI = 0.5; if(!scvf.boundary()) { const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); const auto& subFluxDataJI = dataHandleJ.subFluxData(); - std::for_each(subFluxDataJI.cbegin(), subFluxDataJI.cend(), - [&fluxJI, &elemVolVars, &pseudoPotential](const auto& e) - { fluxJI += e.coefficient * pseudoPotential(elemVolVars[e.index], e.position); } ); + FluxCalculator::flux(flux, elemVolVars, subFluxDataIJ, subFluxDataJI, pseudoPotential); } else - { - wIJ = 1.0; - wJI = 0.0; - } - - return scvf.area() * (wIJ*fluxIJ - wJI*fluxJI); + FluxCalculator::boundaryFlux(flux, elemVolVars, subFluxDataIJ, pseudoPotential); } else { - auto pressure = [&elemVolVars, phaseIdx](const auto& volVars) - { return volVars.pressure(phaseIdx);}; + auto pressure = [&elemVolVars, phaseIdx](const auto& volVars, const auto& x) + { return volVars.pressure(phaseIdx);}; - Scalar fluxIJ = 0.0; - Scalar wIJ = 0.5; - std::for_each(subFluxDataIJ.cbegin(), subFluxDataIJ.cend(), - [&fluxIJ, &elemVolVars, &pressure](const auto& e) - { fluxIJ += e.coefficient * pressure(elemVolVars[e.index]); } ); - - Scalar fluxJI = 0.0; - Scalar wJI = 0.5; if(!scvf.boundary()) { const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); const auto& subFluxDataJI = dataHandleJ.subFluxData(); - std::for_each(subFluxDataJI.cbegin(), subFluxDataJI.cend(), - [&fluxJI, &elemVolVars, &pressure](const auto& e) - { fluxJI += e.coefficient * pressure(elemVolVars[e.index]); } ); + FluxCalculator::flux(flux, elemVolVars, subFluxDataIJ, subFluxDataJI, pressure); } else - { - wIJ = 1.0; - wJI = 0.0; - } - - return scvf.area() * (wIJ*fluxIJ - wJI*fluxJI); + FluxCalculator::boundaryFlux(flux, elemVolVars, subFluxDataIJ, pressure); } - return 0.0; + return scvf.area() * flux; } - -private: - template = 0> - static decltype(auto) getPermeability_(const Problem& problem, - const VolumeVariables& volVars, - const GlobalPosition& scvfIpGlobal) - { return volVars.permeability(); } - - template = 0> - static decltype(auto) getPermeability_(const Problem& problem, - const VolumeVariables& volVars, - const GlobalPosition& scvfIpGlobal) - { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); } }; } // end namespace Dumux -- GitLab From aa5fb786a6ac57af10ff11155279f6a01639368f Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Tue, 21 Jul 2020 18:23:41 +0200 Subject: [PATCH 25/36] [flux][wmpfa] Add nltpfa for darcy flux discretization --- dumux/discretization/ccwmpfa.hh | 7 ++ .../cellcentered/wmpfa/facedatahandle.hh | 2 +- dumux/flux/ccwmpfa/darcyslaw.hh | 109 ++++++++++++++---- 3 files changed, 96 insertions(+), 22 deletions(-) diff --git a/dumux/discretization/ccwmpfa.hh b/dumux/discretization/ccwmpfa.hh index 86630c8e97..a3d2c3191a 100644 --- a/dumux/discretization/ccwmpfa.hh +++ b/dumux/discretization/ccwmpfa.hh @@ -41,6 +41,7 @@ #include #include #include +#include #include #include @@ -115,6 +116,12 @@ struct ElementBoundaryTypes { using type = CCElemen //! Set the BaseLocalResidual to CCLocalResidual template struct BaseLocalResidual { using type = CCLocalResidual; }; + +template +struct DiscretizationSubmethod { using type = UndefinedProperty; }; + +template +struct DiscretizationSubmethod { static constexpr WMpfaMethod value = WMpfaMethod::avgmpfa; }; } // namespace Properties namespace Detail { diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index d8698ed84e..4e8fbe77ea 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -135,7 +135,7 @@ public: auto c = coeff[i]*e.weight(); entries_[0].coefficient += c; //ToDo pos is not needed for each entry!!! - entries_.push_back({-1.0*c, e.dofIndex(), intData.position()}); + entries_.push_back({c, e.dofIndex(), intData.position()}); } } } diff --git a/dumux/flux/ccwmpfa/darcyslaw.hh b/dumux/flux/ccwmpfa/darcyslaw.hh index 977fb0cdfc..ac991f57c0 100644 --- a/dumux/flux/ccwmpfa/darcyslaw.hh +++ b/dumux/flux/ccwmpfa/darcyslaw.hh @@ -95,7 +95,6 @@ private: const AdvectionDataHandle* handle_; }; - template class WMpfaDarcysLawsFluxCalculator; @@ -104,8 +103,9 @@ class WMpfaDarcysLawsFluxCalculator { public: - template + template static void flux(Scalar& flux, + const SubControlVolumeFace& scvf, const ElementVolumeVariables& elemVolVars, const FluxData& dIJ, const FluxData& dJI, @@ -113,32 +113,97 @@ public: { Scalar fluxIJ = 0.0; Scalar fluxJI = 0.0; - Scalar wIJ = 0.5; - Scalar wJI = 0.5; - std::for_each(dIJ.cbegin(), dIJ.cend(), - [&fluxIJ, &elemVolVars, &solValues](const auto& e) - { fluxIJ += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); - std::for_each(dJI.cbegin(), dJI.cend(), + fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); + std::for_each(dIJ.cbegin()+1, dIJ.cend(), + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + fluxJI += dJI[0].coefficient * solValues(elemVolVars[dJI[0].index], dJI[0].position); + std::for_each(dJI.cbegin()+1, dJI.cend(), [&fluxJI, &elemVolVars, &solValues](const auto& e) - { fluxJI += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + { fluxJI -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); - flux = (wIJ*fluxIJ - wJI*fluxJI); + flux = scvf.area()*(0.5*fluxIJ - 0.5*fluxJI); } - template + template static void boundaryFlux(Scalar& flux, + const SubControlVolumeFace& scvf, const ElementVolumeVariables& elemVolVars, const FluxData& dIJ, const SolValues& solValues) { Scalar fluxIJ = 0.0; - Scalar wIJ = 1.0; - std::for_each(dIJ.cbegin(), dIJ.cend(), + fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); + std::for_each(dIJ.cbegin()+1, dIJ.cend(), [&fluxIJ, &elemVolVars, &solValues](const auto& e) - { fluxIJ += e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); - flux = wIJ*fluxIJ; + flux = scvf.area()*fluxIJ; + } +}; + +template<> +class WMpfaDarcysLawsFluxCalculator +{ +public: + + template + static void flux(Scalar& flux, + const SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const FluxData& dJI, + const SolValues& solValues) + { + Scalar lambdaIJ = 0.0; + Scalar lambdaJI = 0.0; + Scalar tIJ = 0.0; + Scalar tJI = 0.0; + std::for_each(dIJ.cbegin()+1, dIJ.cend(), + [&lambdaIJ, &tJI, &scvf, &elemVolVars, &solValues](const auto& e) + { + (e.index == scvf.outsideScvIdx()) + ? tJI += e.coefficient + : lambdaIJ += e.coefficient * solValues(elemVolVars[e.index], e.position); + } ); + + std::for_each(dJI.cbegin()+1, dJI.cend(), + [&lambdaJI, &tIJ, &scvf, &elemVolVars, &solValues](const auto& e) + { + (e.index == scvf.insideScvIdx()) + ? tIJ += e.coefficient + : lambdaJI += e.coefficient * solValues(elemVolVars[e.index], e.position); + } ); + + auto calcAbs = [](Scalar a){ return std::abs(a) + 1e-8; }; + + Scalar wIJ = (calcAbs(lambdaJI))/(calcAbs(lambdaIJ) + calcAbs(lambdaJI)); + Scalar wJI = (calcAbs(lambdaIJ))/(calcAbs(lambdaIJ) + calcAbs(lambdaJI)); + + tIJ *= wJI; + tIJ += wIJ*dIJ[0].coefficient; + tJI *= wIJ; + tJI += wJI*dJI[0].coefficient; + + flux = scvf.area()*( tIJ * solValues(elemVolVars[dIJ[0].index], dIJ[0].position) + -tJI * solValues(elemVolVars[dJI[0].index], dJI[0].position)); + } + + template + static void boundaryFlux(Scalar& flux, + const SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const SolValues& solValues) + { + Scalar fluxIJ = 0.0; + fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); + std::for_each(dIJ.cbegin()+1, dIJ.cend(), + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + + flux = scvf.area()*fluxIJ; } }; @@ -165,7 +230,7 @@ class DarcysLawImplementation using ElementFluxVariablesCache = typename GetPropType::LocalView; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - using FluxCalculator = WMpfaDarcysLawsFluxCalculator; + static constexpr WMpfaMethod method = getPropValue(); public: //! state the discretization method this implementation belongs to @@ -184,6 +249,8 @@ class DarcysLawImplementation int phaseIdx, const ElementFluxVarsCache& elemFluxVarsCache) { + using FluxCalculator = WMpfaDarcysLawsFluxCalculator; + static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); const auto& fluxVarsCache = elemFluxVarsCache[scvf]; @@ -214,10 +281,10 @@ class DarcysLawImplementation const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); const auto& subFluxDataJI = dataHandleJ.subFluxData(); - FluxCalculator::flux(flux, elemVolVars, subFluxDataIJ, subFluxDataJI, pseudoPotential); + FluxCalculator::flux(flux, scvf, elemVolVars, subFluxDataIJ, subFluxDataJI, pseudoPotential); } else - FluxCalculator::boundaryFlux(flux, elemVolVars, subFluxDataIJ, pseudoPotential); + FluxCalculator::boundaryFlux(flux, scvf, elemVolVars, subFluxDataIJ, pseudoPotential); } else { @@ -229,13 +296,13 @@ class DarcysLawImplementation const auto& flippedScvf = fvGeometry.flipScvf(scvf.index()); const auto& dataHandleJ = elemFluxVarsCache[flippedScvf].dataHandle(); const auto& subFluxDataJI = dataHandleJ.subFluxData(); - FluxCalculator::flux(flux, elemVolVars, subFluxDataIJ, subFluxDataJI, pressure); + FluxCalculator::flux(flux, scvf, elemVolVars, subFluxDataIJ, subFluxDataJI, pressure); } else - FluxCalculator::boundaryFlux(flux, elemVolVars, subFluxDataIJ, pressure); + FluxCalculator::boundaryFlux(flux, scvf, elemVolVars, subFluxDataIJ, pressure); } - return scvf.area() * flux; + return flux; } }; -- GitLab From 542c244bc3e07e535cc6ec5baa59ddd7314abb2e Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 22 Jul 2020 15:28:25 +0200 Subject: [PATCH 26/36] [flux][mpfa] Add possibility to add submethod for some type of mpfa schemes --- dumux/common/properties.hh | 3 ++- dumux/discretization/ccwmpfa.hh | 3 --- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh index 90f7b1ff16..e030850152 100644 --- a/dumux/common/properties.hh +++ b/dumux/common/properties.hh @@ -130,7 +130,8 @@ template struct SecondaryInteractionVolume { using type = UndefinedProperty; }; //!< The secondary interaction volume type used e.g. on the boundaries template struct DualGridNodalIndexSet { using type = UndefinedProperty; }; //!< The type used for the nodal index sets of the dual grid - +template +struct DiscretizationSubmethod { using type = UndefinedProperty; }; //!< The type used for an mpfa submethod belonging to a specific family of schemes ///////////////////////////////////////////////////////////// // Properties used by models involving flow in porous media: ///////////////////////////////////////////////////////////// diff --git a/dumux/discretization/ccwmpfa.hh b/dumux/discretization/ccwmpfa.hh index a3d2c3191a..3b606ceef5 100644 --- a/dumux/discretization/ccwmpfa.hh +++ b/dumux/discretization/ccwmpfa.hh @@ -117,9 +117,6 @@ struct ElementBoundaryTypes { using type = CCElemen template struct BaseLocalResidual { using type = CCLocalResidual; }; -template -struct DiscretizationSubmethod { using type = UndefinedProperty; }; - template struct DiscretizationSubmethod { static constexpr WMpfaMethod value = WMpfaMethod::avgmpfa; }; } // namespace Properties -- GitLab From 981f530ceecb1e58271ea5fc88a89903eeccfeef Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 22 Jul 2020 16:32:39 +0200 Subject: [PATCH 27/36] [tests][1p] Add draft for anisotropic convergence test case --- .../1p/implicit/convergence/CMakeLists.txt | 1 + .../convergence/unstructured/CMakeLists.txt | 25 ++ .../unstructured/convergencetest.py | 71 +++++ .../implicit/convergence/unstructured/main.cc | 199 ++++++++++++++ .../convergence/unstructured/params.input | 17 ++ .../convergence/unstructured/problem.hh | 251 ++++++++++++++++++ .../convergence/unstructured/spatialparams.hh | 109 ++++++++ 7 files changed, 673 insertions(+) create mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt create mode 100755 test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py create mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc create mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/params.input create mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh create mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh diff --git a/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt index fbc093277d..874247583b 100644 --- a/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt +++ b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(analyticsolution) add_subdirectory(discretesolution) +add_subdirectory(unstructured) diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt new file mode 100644 index 0000000000..5ab20d53ff --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt @@ -0,0 +1,25 @@ +add_input_file_links() +dune_symlink_to_source_files(FILES "convergencetest.py") + +dumux_add_test(NAME test_1p_convergence_structured + LABELS porousmediumflow 1p + SOURCES main.cc + COMPILE_DEFINITIONS GRIDTYPE=Dune::YaspGrid<2> + COMPILE_DEFINITIONS TYPETAG=OnePConvergence + TIMEOUT 1000 + CMAKE_GUARD HAVE_UMFPACK + COMMAND ./convergencetest.py + CMD_ARGS test_1p_convergence_structured params.input + -Problem.Name test_1p_convergence_structured) + +dumux_add_test(NAME test_1p_convergence_unstructured + LABELS porousmediumflow 1p + SOURCES main.cc + COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2> + COMPILE_DEFINITIONS TYPETAG=OnePConvergence + TIMEOUT 1000 + CMAKE_GUARD HAVE_UMFPACK + COMMAND ./convergencetest.py + CMD_ARGS test_1p_convergence_unstructured params.input + -Problem.Name test_1p_convergence_unstructured + -Grid.File ../../incompressible/grids/randomlydistorted.dgf) diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py b/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py new file mode 100755 index 0000000000..39a8d5a7c3 --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 + +from math import * +import subprocess +import sys + +if len(sys.argv) < 2: + sys.stderr.write('Please provide a single argument to the script\n') + sys.exit(1) + +executableName = str(sys.argv[1]) +testargs = [str(i) for i in sys.argv][2:] +testname = testargs[2] + +# remove the old log files +subprocess.call(['rm', testname + '.log']) +print("Removed old log file ({})!".format(testname + '.log')) + +# do the runs with different refinement +for i in [0, 1, 2, 3]: + subprocess.call(['./' + executableName] + testargs + ['-Grid.Refinement', str(i)]) + +def checkRates(): + # check the rates and append them to the log file + logfile = open(testname + '.log', "r+") + + errorP = [] + for line in logfile: + line = line.strip("\n") + line = line.strip("\[ConvergenceTest\]") + line = line.split() + errorP.append(float(line[2])) + + resultsP = [] + logfile.truncate(0) + logfile.write("n\terrorP\t\trateP\n") + logfile.write("-"*50 + "\n") + for i in range(len(errorP)-1): + if isnan(errorP[i]) or isinf(errorP[i]): + continue + if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)): + rateP = (log(errorP[i])-log(errorP[i+1]))/log(2) + message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP) + logfile.write(message) + resultsP.append(rateP) + else: + logfile.write("error: exact solution!?") + i = len(errorP)-1 + message = "{}\t{:0.4e}\n".format(i, errorP[i], "") + logfile.write(message) + + logfile.close() + print("\nComputed the following convergence rates for {}:\n".format(testname)) + + subprocess.call(['cat', testname + '_darcy.log']) + + return {"p" : resultsP} + +def checkConvRates(): + rates = checkRates() + + def mean(numbers): + return float(sum(numbers)) / len(numbers) + + # check the rates, we expect rates around 2 + if mean(rates["p"]) < 2.05 and mean(rates["p"]) < 1.95: + sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n") + sys.exit(1) + + +checkConvRates() diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc b/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc new file mode 100644 index 0000000000..4cf24a5363 --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc @@ -0,0 +1,199 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup OnePTests + * \brief Test for the one-phase CC model + */ + +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include "problem.hh" + +/*! +* \brief Creates analytical solution. +* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces +* \param problem the problem for which to evaluate the analytical solution +*/ +template +auto createAnalyticalSolution(const Problem& problem) +{ + const auto& gridGeometry = problem.gridGeometry(); + using GridView = typename std::decay_t::GridView; + + static constexpr auto dim = GridView::dimension; + static constexpr auto dimWorld = GridView::dimensionworld; + + using VelocityVector = Dune::FieldVector; + + std::vector analyticalPressure; + std::vector analyticalVelocity; + + analyticalPressure.resize(gridGeometry.numDofs()); + analyticalVelocity.resize(gridGeometry.numDofs()); + + for (const auto& element : elements(gridGeometry.gridView())) + { + auto fvGeometry = localView(gridGeometry); + fvGeometry.bindElement(element); + for (auto&& scv : scvs(fvGeometry)) + { + const auto ccDofIdx = scv.dofIndex(); + const auto ccDofPosition = scv.dofPosition(); + const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition); + analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim]; + + for (int dirIdx = 0; dirIdx < dim; ++dirIdx) + analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx]; + } + } + + return std::make_tuple(analyticalPressure, analyticalVelocity); +} + +template +void printL2Error(const Problem& problem, const SolutionVector& x) +{ + using namespace Dumux; + using Scalar = double; + + Scalar l2error = 0.0; + + for (const auto& element : elements(problem.gridGeometry().gridView())) + { + auto fvGeometry = localView(problem.gridGeometry()); + fvGeometry.bindElement(element); + + for (auto&& scv : scvs(fvGeometry)) + { + const auto dofIdx = scv.dofIndex(); + const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/]; + l2error += scv.volume()*(delta*delta); + } + } + using std::sqrt; + l2error = sqrt(l2error); + + const auto numDofs = problem.gridGeometry().numDofs(); + std::ostream tmp(std::cout.rdbuf()); + tmp << std::setprecision(8) << "** L2 error (abs) for " + << std::setw(6) << numDofs << " cc dofs " + << std::scientific + << "L2 error = " << l2error + << std::endl; + + // write the norm into a log file + std::ofstream logFile; + logFile.open(problem.name() + ".log", std::ios::app); + logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl; + logFile.close(); +} + +int main(int argc, char** argv) +{ + using namespace Dumux; + + using TypeTag = Properties::TTag::TYPETAG; + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + //////////////////////////////////////////////////////////// + // parse the command line arguments and input file + //////////////////////////////////////////////////////////// + + // parse command line arguments + Parameters::init(argc, argv); + + ////////////////////////////////////////////////////////////////////// + // try to create a grid (from the given grid file or the input file) + ///////////////////////////////////////////////////////////////////// + using Grid = GetPropType; + GridManager gridManager; + gridManager.init(); + + // we compute on the leaf grid view + const auto& leafGridView = gridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using GridGeometry = GetPropType; + auto gridGeometry = std::make_shared(leafGridView); + gridGeometry->update(); + + // the problem (boundary conditions) + using Problem = GetPropType; + auto problem = std::make_shared(gridGeometry); + + // the solution vector + using SolutionVector = GetPropType; + SolutionVector x(gridGeometry->numDofs()); + + // the grid variables + using GridVariables = GetPropType; + auto gridVariables = std::make_shared(problem, gridGeometry); + gridVariables->init(x); + + // intialize the vtk output module + VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); + using VelocityOutput = GetPropType; + vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); + using IOFields = GetPropType; + IOFields::initOutputModule(vtkWriter); // Add model specific output fields + vtkWriter.write(0.0); + + // create assembler & linear solver + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, gridGeometry, gridVariables); + + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = Dumux::NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + // linearize & solve + nonLinearSolver.solve(x); + + // output result to vtk + vtkWriter.write(1.0); + + printL2Error(*problem, x); + + if (mpiHelper.rank() == 0) + Parameters::print(); + + return 0; + +} diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input b/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input new file mode 100644 index 0000000000..6194a9cd67 --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input @@ -0,0 +1,17 @@ +[Grid] +UpperRight = 1 1 +Cells = 40 40 +Refinement = 0 + +[Problem] +Name = convergence + +[Problem] +EnableGravity = false + +[Component] +LiquidDensity = 1.0 +LiquidKinematicViscosity = 1.0 + +[Assembly] +NumericDifference.BaseEpsilon = 1e-4 diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh b/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh new file mode 100644 index 0000000000..6b4558ea3f --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh @@ -0,0 +1,251 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief The convergence test + */ + +#ifndef DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH +#define DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include "spatialparams.hh" + +namespace Dumux { +template +class ConvergenceProblem; + +namespace Properties { +// Create new type tags +namespace TTag { +struct OnePConvergence { using InheritsFrom = std::tuple; }; +} // end namespace TTag + +// Set the problem property +template +struct Problem { using type = Dumux::ConvergenceProblem; }; + +// the fluid system +template +struct FluidSystem +{ + using Scalar = GetPropType; + using type = FluidSystems::OnePLiquid > ; +}; + +// Set the grid type +template +struct Grid { using type = GRIDTYPE; }; + +template +struct SpatialParams +{ + using GridGeometry = GetPropType; + using Scalar = GetPropType; + using type = ConvergenceTestSpatialParams; +}; + +template +struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; +template +struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; +template +struct EnableGridGeometryCache { static constexpr bool value = true; }; + +template +struct DiscretizationSubmethod { static constexpr WMpfaMethod value = WMpfaMethod::nltpfa; }; + +} // end namespace Properties + +/*! + * \ingroup BoundaryTests + * \brief The convergence test + */ +template +class ConvergenceProblem : public PorousMediumFlowProblem +{ + using ParentType = PorousMediumFlowProblem; + using GridView = typename GetPropType::GridView; + using Scalar = GetPropType; + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using BoundaryTypes = Dumux::BoundaryTypes::numEq()>; + using VolumeVariables = GetPropType; + using FVElementGeometry = typename GetPropType::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using GridGeometry = GetPropType; + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + static constexpr auto velocityXIdx = 0; + static constexpr auto velocityYIdx = 1; + static constexpr auto pressureIdx = 2; + +public: + //! export the Indices + using Indices = typename GetPropType::Indices; + + ConvergenceProblem(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + {} + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the temperature within the domain in [K]. + * + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary control volume. + * + * \param element The element + * \param scvf The boundary sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + + values.setAllDirichlet(); + + return values; + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet control volume. + * + * \param element The element for which the Dirichlet boundary condition is set + * \param scvf The boundary subcontrolvolumeface + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + const auto p = analyticalSolution(scvf.center())[pressureIdx]; + return PrimaryVariables(p); + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + /*! + * \brief Evaluates the source term for all phases within a given + * sub control volume. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition& globalPos) const + { + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + using std::exp; using std::sin; using std::cos; + static constexpr Scalar omega = M_PI; + static constexpr Scalar c = 0.9; + const Scalar cosOmegaX = cos(omega*x); + static const Scalar expTwo = exp(2); + const Scalar expYPlusOne = exp(y+1); + + const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1) + + 1.5*c*expYPlusOne*cosOmegaX + + omega*omega*(expYPlusOne - expTwo + 2)) + *sin(omega*x); + return NumEqVector(result); + } + + // \} + + /*! + * \brief Evaluates the initial value for a control volume. + * + * \param element The element + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initial(const Element &element) const + { + return PrimaryVariables(0.0); + } + + /*! + * \brief Returns the analytical solution of the problem at a given position. + * + * \param globalPos The global position + */ + auto analyticalSolution(const GlobalPosition& globalPos) const + { + Dune::FieldVector sol(0.0); + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + static constexpr Scalar omega = M_PI; + static constexpr Scalar c = 0.9; + using std::exp; using std::sin; using std::cos; + const Scalar sinOmegaX = sin(omega*x); + const Scalar cosOmegaX = cos(omega*x); + static const Scalar expTwo = exp(2); + const Scalar expYPlusOne = exp(y+1); + + sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX + 10.0; + sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX + -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX; + sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX + -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX; + + return sol; + } + + // \} + +private: + static constexpr Scalar eps_ = 1e-7; +}; +} // end namespace Dumux + +#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh b/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh new file mode 100644 index 0000000000..92c6464af9 --- /dev/null +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh @@ -0,0 +1,109 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief The spatial parameters class for the test problem using the 1p cc model. + */ + +#ifndef DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH +#define DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH + +#include +#include +#include +#include +#include + +namespace Dumux { + +/*! + * \ingroup BoundaryTests + * \brief The spatial parameters class for the test problem using the + * 1p cc model. + */ +template +class ConvergenceTestSpatialParams +: public FVSpatialParamsOneP> +{ + using GridView = typename GridGeometry::GridView; + using ParentType = FVSpatialParamsOneP>; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + static constexpr int dimWorld = GridView::dimensionworld; + using DimWorldMatrix = Dune::FieldMatrix; + + +public: + // export permeability type + using PermeabilityType = DimWorldMatrix; + + ConvergenceTestSpatialParams(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + {} + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. + * + * \param element The element + * \param scv The sub control volume + * \param elemSol The element solution vector + * \return the intrinsic permeability + */ + template + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + PermeabilityType K(0.0); + + using std::cos; + using std::sin; + using std::exp; + + static constexpr Scalar c = 0.9; + static constexpr Scalar omega = M_PI; + + const Scalar x = scv.center()[0]; + K[0][0] = 1.0; + K[0][1] = -c/(2*omega) * sin(omega*x); + K[1][0] = K[0][1]; + K[1][1] = exp(-2)*(1 + c*cos(omega*x)); + + return K; + } + + /*! \brief Defines the porosity in [-]. + * + * \param globalPos The global position + */ + Scalar porosityAtPos(const GlobalPosition& globalPos) const + { return 0.4; } + +private: + Scalar permeability_; +}; + +} // end namespace Dumux + +#endif // DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH -- GitLab From ecb2026a5291f0289f7fe02a667b8af38bb2cfa4 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 30 Sep 2020 15:56:06 +0200 Subject: [PATCH 28/36] [wmpfa] Add specialization for linear solver traits --- dumux/linear/linearsolvertraits.hh | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/dumux/linear/linearsolvertraits.hh b/dumux/linear/linearsolvertraits.hh index b61ba21e14..303c508e4a 100644 --- a/dumux/linear/linearsolvertraits.hh +++ b/dumux/linear/linearsolvertraits.hh @@ -172,11 +172,16 @@ struct LinearSolverTraitsImpl { return false; } }; -//! Cell-centered mpfa: use overlapping model +//! Cell-centered mpfa: use overlapping model (like tpfa) template struct LinearSolverTraitsImpl : public LinearSolverTraitsImpl {}; +//! Weighted mpfa: use overlapping model (like tpfa) +template +struct LinearSolverTraitsImpl +: public LinearSolverTraitsImpl {}; + //! staggered: use overlapping model TODO provide staggered-specific traits, combining overlapping/non-overlapping template struct LinearSolverTraitsImpl -- GitLab From 6209146aea7ef6bc0ad637c6feaa0396dfa4affd Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 30 Sep 2020 16:04:28 +0200 Subject: [PATCH 29/36] [test][1p][wmpfa] Remove darcy suffix for log file --- .../1p/implicit/convergence/unstructured/convergencetest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py b/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py index 39a8d5a7c3..8d7c864364 100755 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py +++ b/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py @@ -52,7 +52,7 @@ def checkRates(): logfile.close() print("\nComputed the following convergence rates for {}:\n".format(testname)) - subprocess.call(['cat', testname + '_darcy.log']) + subprocess.call(['cat', testname + '.log']) return {"p" : resultsP} -- GitLab From 7de9c413652933cad1e3f4b552cbace845adac01 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 30 Sep 2020 16:04:57 +0200 Subject: [PATCH 30/36] [wmpfa] Do not add boundary volvars on non-Dirichlet boundaries --- .../cellcentered/wmpfa/elementvolumevariables.hh | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh index 7ee0572a7d..a661b2377a 100644 --- a/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/wmpfa/elementvolumevariables.hh @@ -112,17 +112,6 @@ namespace CCWMpfa { volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(scvf.outsideScvIdx()); } - else - { - VolumeVariables dirichletVolVars; - dirichletVolVars.update(elementSolution(problem.dirichlet(element, scvf)), - problem, - element, - scvI); - - volVars.emplace_back(std::move(dirichletVolVars)); - volVarIndices.push_back(scvf.outsideScvIdx()); - } } } return std::make_pair(volVars, volVarIndices); -- GitLab From c2c3f95889501bcbf26bb85d1a14e0b8464dbcb5 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 30 Sep 2020 20:23:00 +0200 Subject: [PATCH 31/36] [tests][wmpfa] Add wmpfa to analytic convergence test and delete obsolete old test --- .../1p/implicit/convergence/CMakeLists.txt | 1 - .../analyticsolution/CMakeLists.txt | 23 ++ .../analyticsolution/convergencetest.py | 2 +- .../analyticsolution/properties.hh | 6 + .../convergence/unstructured/CMakeLists.txt | 25 -- .../unstructured/convergencetest.py | 71 ----- .../implicit/convergence/unstructured/main.cc | 199 -------------- .../convergence/unstructured/params.input | 17 -- .../convergence/unstructured/problem.hh | 251 ------------------ .../convergence/unstructured/spatialparams.hh | 109 -------- 10 files changed, 30 insertions(+), 674 deletions(-) delete mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt delete mode 100755 test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py delete mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc delete mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/params.input delete mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh delete mode 100644 test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh diff --git a/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt index 874247583b..fbc093277d 100644 --- a/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt +++ b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt @@ -1,3 +1,2 @@ add_subdirectory(analyticsolution) add_subdirectory(discretesolution) -add_subdirectory(unstructured) diff --git a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/CMakeLists.txt index a50498793e..069868761b 100644 --- a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/CMakeLists.txt +++ b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/CMakeLists.txt @@ -6,15 +6,18 @@ dune_symlink_to_source_files(FILES "convergencetest.py") # the tests on unstructured grids to be skipped add_executable(test_1p_convergence_analytic_tpfa EXCLUDE_FROM_ALL main.cc) add_executable(test_1p_convergence_analytic_mpfa EXCLUDE_FROM_ALL main.cc) +add_executable(test_1p_convergence_analytic_wmpfa EXCLUDE_FROM_ALL main.cc) add_executable(test_1p_convergence_analytic_box EXCLUDE_FROM_ALL main.cc) if (dune-uggrid_FOUND) target_compile_definitions(test_1p_convergence_analytic_tpfa PUBLIC "TYPETAG=OnePConvergenceTpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>") target_compile_definitions(test_1p_convergence_analytic_mpfa PUBLIC "TYPETAG=OnePConvergenceMpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>") + target_compile_definitions(test_1p_convergence_analytic_wmpfa PUBLIC "TYPETAG=OnePConvergenceWMpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>") target_compile_definitions(test_1p_convergence_analytic_box PUBLIC "TYPETAG=OnePConvergenceBox" PUBLIC "GRIDTYPE=Dune::UGGrid<2>") else() target_compile_definitions(test_1p_convergence_analytic_tpfa PUBLIC "TYPETAG=OnePConvergenceTpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>") target_compile_definitions(test_1p_convergence_analytic_mpfa PUBLIC "TYPETAG=OnePConvergenceMpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>") + target_compile_definitions(test_1p_convergence_analytic_wmpfa PUBLIC "TYPETAG=OnePConvergenceWMpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>") target_compile_definitions(test_1p_convergence_analytic_box PUBLIC "TYPETAG=OnePConvergenceBox" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>") endif() @@ -33,6 +36,13 @@ dumux_add_test(NAME test_1p_convergence_analytic_mpfa_structured COMMAND ./convergencetest.py CMD_ARGS test_1p_convergence_analytic_mpfa mpfa_structured) +# using wmpfa and structured grid +dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_structured + TARGET test_1p_convergence_analytic_wmpfa + LABELS porousmediumflow 1p + COMMAND ./convergencetest.py + CMD_ARGS test_1p_convergence_analytic_wmpfa wmpfa_structured) + # using box and structured grid dumux_add_test(NAME test_1p_convergence_analytic_box_structured TARGET test_1p_convergence_analytic_box @@ -50,6 +60,14 @@ if (dune-uggrid_FOUND) CMD_ARGS test_1p_convergence_analytic_mpfa mpfa_unstructured params.input -Grid.File ../../incompressible/grids/randomlydistorted.dgf) + # using wmpfa and unstructured grid + dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_unstructured + TARGET test_1p_convergence_analytic_wmpfa + LABELS porousmediumflow 1p + COMMAND ./convergencetest.py + CMD_ARGS test_1p_convergence_analytic_wmpfa wmpfa_unstructured + params.input -Grid.File ../../incompressible/grids/randomlydistorted.dgf) + # using box and unstructured grid dumux_add_test(NAME test_1p_convergence_analytic_box_unstructured TARGET test_1p_convergence_analytic_box @@ -65,6 +83,11 @@ else() CMAKE_GUARD dune-uggrid_FOUND LABELS porousmediumflow 1p) + dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_unstructured + SOURCES main.cc + CMAKE_GUARD dune-uggrid_FOUND + LABELS porousmediumflow 1p) + dumux_add_test(NAME test_1p_convergence_analytic_box_unstructured SOURCES main.cc CMAKE_GUARD dune-uggrid_FOUND diff --git a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/convergencetest.py b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/convergencetest.py index b2991031c2..6f528b9a0d 100755 --- a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/convergencetest.py +++ b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/convergencetest.py @@ -69,6 +69,6 @@ def mean(numbers): return float(sum(numbers)) / len(numbers) # check the rates, we expect rates around 2 -if mean(rates["p"]) < 1.8: +if mean(rates["p"]) < 1.75: sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n") sys.exit(1) diff --git a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/properties.hh b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/properties.hh index c72ed92748..5fdaa4feba 100644 --- a/test/porousmediumflow/1p/implicit/convergence/analyticsolution/properties.hh +++ b/test/porousmediumflow/1p/implicit/convergence/analyticsolution/properties.hh @@ -31,6 +31,8 @@ #include #include +#include +#include #include #include @@ -47,6 +49,7 @@ namespace TTag { struct OnePConvergence { using InheritsFrom = std::tuple; }; struct OnePConvergenceTpfa { using InheritsFrom = std::tuple; }; struct OnePConvergenceMpfa { using InheritsFrom = std::tuple; }; +struct OnePConvergenceWMpfa { using InheritsFrom = std::tuple; }; struct OnePConvergenceBox { using InheritsFrom = std::tuple; }; } // end namespace TTag @@ -82,6 +85,9 @@ struct EnableGridFluxVariablesCache { static con template struct EnableGridGeometryCache { static constexpr bool value = true; }; +template +struct DiscretizationSubmethod { static constexpr WMpfaMethod value = WMpfaMethod::nltpfa; }; + } // end namespace Dumux::Properties #endif diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt deleted file mode 100644 index 5ab20d53ff..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt +++ /dev/null @@ -1,25 +0,0 @@ -add_input_file_links() -dune_symlink_to_source_files(FILES "convergencetest.py") - -dumux_add_test(NAME test_1p_convergence_structured - LABELS porousmediumflow 1p - SOURCES main.cc - COMPILE_DEFINITIONS GRIDTYPE=Dune::YaspGrid<2> - COMPILE_DEFINITIONS TYPETAG=OnePConvergence - TIMEOUT 1000 - CMAKE_GUARD HAVE_UMFPACK - COMMAND ./convergencetest.py - CMD_ARGS test_1p_convergence_structured params.input - -Problem.Name test_1p_convergence_structured) - -dumux_add_test(NAME test_1p_convergence_unstructured - LABELS porousmediumflow 1p - SOURCES main.cc - COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2> - COMPILE_DEFINITIONS TYPETAG=OnePConvergence - TIMEOUT 1000 - CMAKE_GUARD HAVE_UMFPACK - COMMAND ./convergencetest.py - CMD_ARGS test_1p_convergence_unstructured params.input - -Problem.Name test_1p_convergence_unstructured - -Grid.File ../../incompressible/grids/randomlydistorted.dgf) diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py b/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py deleted file mode 100755 index 8d7c864364..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py +++ /dev/null @@ -1,71 +0,0 @@ -#!/usr/bin/env python3 - -from math import * -import subprocess -import sys - -if len(sys.argv) < 2: - sys.stderr.write('Please provide a single argument to the script\n') - sys.exit(1) - -executableName = str(sys.argv[1]) -testargs = [str(i) for i in sys.argv][2:] -testname = testargs[2] - -# remove the old log files -subprocess.call(['rm', testname + '.log']) -print("Removed old log file ({})!".format(testname + '.log')) - -# do the runs with different refinement -for i in [0, 1, 2, 3]: - subprocess.call(['./' + executableName] + testargs + ['-Grid.Refinement', str(i)]) - -def checkRates(): - # check the rates and append them to the log file - logfile = open(testname + '.log', "r+") - - errorP = [] - for line in logfile: - line = line.strip("\n") - line = line.strip("\[ConvergenceTest\]") - line = line.split() - errorP.append(float(line[2])) - - resultsP = [] - logfile.truncate(0) - logfile.write("n\terrorP\t\trateP\n") - logfile.write("-"*50 + "\n") - for i in range(len(errorP)-1): - if isnan(errorP[i]) or isinf(errorP[i]): - continue - if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)): - rateP = (log(errorP[i])-log(errorP[i+1]))/log(2) - message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP) - logfile.write(message) - resultsP.append(rateP) - else: - logfile.write("error: exact solution!?") - i = len(errorP)-1 - message = "{}\t{:0.4e}\n".format(i, errorP[i], "") - logfile.write(message) - - logfile.close() - print("\nComputed the following convergence rates for {}:\n".format(testname)) - - subprocess.call(['cat', testname + '.log']) - - return {"p" : resultsP} - -def checkConvRates(): - rates = checkRates() - - def mean(numbers): - return float(sum(numbers)) / len(numbers) - - # check the rates, we expect rates around 2 - if mean(rates["p"]) < 2.05 and mean(rates["p"]) < 1.95: - sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n") - sys.exit(1) - - -checkConvRates() diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc b/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc deleted file mode 100644 index 4cf24a5363..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc +++ /dev/null @@ -1,199 +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 3 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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup OnePTests - * \brief Test for the one-phase CC model - */ - -#include - -#include - -#include -#include - -#include -#include - -#include -#include - -#include -#include -#include - -#include "problem.hh" - -/*! -* \brief Creates analytical solution. -* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces -* \param problem the problem for which to evaluate the analytical solution -*/ -template -auto createAnalyticalSolution(const Problem& problem) -{ - const auto& gridGeometry = problem.gridGeometry(); - using GridView = typename std::decay_t::GridView; - - static constexpr auto dim = GridView::dimension; - static constexpr auto dimWorld = GridView::dimensionworld; - - using VelocityVector = Dune::FieldVector; - - std::vector analyticalPressure; - std::vector analyticalVelocity; - - analyticalPressure.resize(gridGeometry.numDofs()); - analyticalVelocity.resize(gridGeometry.numDofs()); - - for (const auto& element : elements(gridGeometry.gridView())) - { - auto fvGeometry = localView(gridGeometry); - fvGeometry.bindElement(element); - for (auto&& scv : scvs(fvGeometry)) - { - const auto ccDofIdx = scv.dofIndex(); - const auto ccDofPosition = scv.dofPosition(); - const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition); - analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim]; - - for (int dirIdx = 0; dirIdx < dim; ++dirIdx) - analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx]; - } - } - - return std::make_tuple(analyticalPressure, analyticalVelocity); -} - -template -void printL2Error(const Problem& problem, const SolutionVector& x) -{ - using namespace Dumux; - using Scalar = double; - - Scalar l2error = 0.0; - - for (const auto& element : elements(problem.gridGeometry().gridView())) - { - auto fvGeometry = localView(problem.gridGeometry()); - fvGeometry.bindElement(element); - - for (auto&& scv : scvs(fvGeometry)) - { - const auto dofIdx = scv.dofIndex(); - const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/]; - l2error += scv.volume()*(delta*delta); - } - } - using std::sqrt; - l2error = sqrt(l2error); - - const auto numDofs = problem.gridGeometry().numDofs(); - std::ostream tmp(std::cout.rdbuf()); - tmp << std::setprecision(8) << "** L2 error (abs) for " - << std::setw(6) << numDofs << " cc dofs " - << std::scientific - << "L2 error = " << l2error - << std::endl; - - // write the norm into a log file - std::ofstream logFile; - logFile.open(problem.name() + ".log", std::ios::app); - logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl; - logFile.close(); -} - -int main(int argc, char** argv) -{ - using namespace Dumux; - - using TypeTag = Properties::TTag::TYPETAG; - - // initialize MPI, finalize is done automatically on exit - const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); - - //////////////////////////////////////////////////////////// - // parse the command line arguments and input file - //////////////////////////////////////////////////////////// - - // parse command line arguments - Parameters::init(argc, argv); - - ////////////////////////////////////////////////////////////////////// - // try to create a grid (from the given grid file or the input file) - ///////////////////////////////////////////////////////////////////// - using Grid = GetPropType; - GridManager gridManager; - gridManager.init(); - - // we compute on the leaf grid view - const auto& leafGridView = gridManager.grid().leafGridView(); - - // create the finite volume grid geometry - using GridGeometry = GetPropType; - auto gridGeometry = std::make_shared(leafGridView); - gridGeometry->update(); - - // the problem (boundary conditions) - using Problem = GetPropType; - auto problem = std::make_shared(gridGeometry); - - // the solution vector - using SolutionVector = GetPropType; - SolutionVector x(gridGeometry->numDofs()); - - // the grid variables - using GridVariables = GetPropType; - auto gridVariables = std::make_shared(problem, gridGeometry); - gridVariables->init(x); - - // intialize the vtk output module - VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); - using VelocityOutput = GetPropType; - vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); - using IOFields = GetPropType; - IOFields::initOutputModule(vtkWriter); // Add model specific output fields - vtkWriter.write(0.0); - - // create assembler & linear solver - using Assembler = FVAssembler; - auto assembler = std::make_shared(problem, gridGeometry, gridVariables); - - using LinearSolver = UMFPackBackend; - auto linearSolver = std::make_shared(); - - // the non-linear solver - using NewtonSolver = Dumux::NewtonSolver; - NewtonSolver nonLinearSolver(assembler, linearSolver); - - // linearize & solve - nonLinearSolver.solve(x); - - // output result to vtk - vtkWriter.write(1.0); - - printL2Error(*problem, x); - - if (mpiHelper.rank() == 0) - Parameters::print(); - - return 0; - -} diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input b/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input deleted file mode 100644 index 6194a9cd67..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/params.input +++ /dev/null @@ -1,17 +0,0 @@ -[Grid] -UpperRight = 1 1 -Cells = 40 40 -Refinement = 0 - -[Problem] -Name = convergence - -[Problem] -EnableGravity = false - -[Component] -LiquidDensity = 1.0 -LiquidKinematicViscosity = 1.0 - -[Assembly] -NumericDifference.BaseEpsilon = 1e-4 diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh b/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh deleted file mode 100644 index 6b4558ea3f..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh +++ /dev/null @@ -1,251 +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 3 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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup BoundaryTests - * \brief The convergence test - */ - -#ifndef DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH -#define DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH - -#include -#include - -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include "spatialparams.hh" - -namespace Dumux { -template -class ConvergenceProblem; - -namespace Properties { -// Create new type tags -namespace TTag { -struct OnePConvergence { using InheritsFrom = std::tuple; }; -} // end namespace TTag - -// Set the problem property -template -struct Problem { using type = Dumux::ConvergenceProblem; }; - -// the fluid system -template -struct FluidSystem -{ - using Scalar = GetPropType; - using type = FluidSystems::OnePLiquid > ; -}; - -// Set the grid type -template -struct Grid { using type = GRIDTYPE; }; - -template -struct SpatialParams -{ - using GridGeometry = GetPropType; - using Scalar = GetPropType; - using type = ConvergenceTestSpatialParams; -}; - -template -struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; -template -struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; -template -struct EnableGridGeometryCache { static constexpr bool value = true; }; - -template -struct DiscretizationSubmethod { static constexpr WMpfaMethod value = WMpfaMethod::nltpfa; }; - -} // end namespace Properties - -/*! - * \ingroup BoundaryTests - * \brief The convergence test - */ -template -class ConvergenceProblem : public PorousMediumFlowProblem -{ - using ParentType = PorousMediumFlowProblem; - using GridView = typename GetPropType::GridView; - using Scalar = GetPropType; - using PrimaryVariables = GetPropType; - using NumEqVector = GetPropType; - using BoundaryTypes = Dumux::BoundaryTypes::numEq()>; - using VolumeVariables = GetPropType; - using FVElementGeometry = typename GetPropType::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; - using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using GridGeometry = GetPropType; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - static constexpr auto velocityXIdx = 0; - static constexpr auto velocityYIdx = 1; - static constexpr auto pressureIdx = 2; - -public: - //! export the Indices - using Indices = typename GetPropType::Indices; - - ConvergenceProblem(std::shared_ptr gridGeometry) - : ParentType(gridGeometry) - {} - - /*! - * \name Problem parameters - */ - // \{ - - /*! - * \brief Returns the temperature within the domain in [K]. - * - */ - Scalar temperature() const - { return 273.15 + 10; } // 10°C - // \} - - /*! - * \name Boundary conditions - */ - // \{ - - /*! - * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary control volume. - * - * \param element The element - * \param scvf The boundary sub control volume face - */ - BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const - { - BoundaryTypes values; - - values.setAllDirichlet(); - - return values; - } - - /*! - * \brief Evaluates the boundary conditions for a Dirichlet control volume. - * - * \param element The element for which the Dirichlet boundary condition is set - * \param scvf The boundary subcontrolvolumeface - * - * For this method, the \a values parameter stores primary variables. - */ - PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const - { - const auto p = analyticalSolution(scvf.center())[pressureIdx]; - return PrimaryVariables(p); - } - - // \} - - /*! - * \name Volume terms - */ - // \{ - /*! - * \brief Evaluates the source term for all phases within a given - * sub control volume. - * - * \param globalPos The global position - */ - NumEqVector sourceAtPos(const GlobalPosition& globalPos) const - { - const Scalar x = globalPos[0]; - const Scalar y = globalPos[1]; - using std::exp; using std::sin; using std::cos; - static constexpr Scalar omega = M_PI; - static constexpr Scalar c = 0.9; - const Scalar cosOmegaX = cos(omega*x); - static const Scalar expTwo = exp(2); - const Scalar expYPlusOne = exp(y+1); - - const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1) - + 1.5*c*expYPlusOne*cosOmegaX - + omega*omega*(expYPlusOne - expTwo + 2)) - *sin(omega*x); - return NumEqVector(result); - } - - // \} - - /*! - * \brief Evaluates the initial value for a control volume. - * - * \param element The element - * - * For this method, the \a priVars parameter stores primary - * variables. - */ - PrimaryVariables initial(const Element &element) const - { - return PrimaryVariables(0.0); - } - - /*! - * \brief Returns the analytical solution of the problem at a given position. - * - * \param globalPos The global position - */ - auto analyticalSolution(const GlobalPosition& globalPos) const - { - Dune::FieldVector sol(0.0); - const Scalar x = globalPos[0]; - const Scalar y = globalPos[1]; - static constexpr Scalar omega = M_PI; - static constexpr Scalar c = 0.9; - using std::exp; using std::sin; using std::cos; - const Scalar sinOmegaX = sin(omega*x); - const Scalar cosOmegaX = cos(omega*x); - static const Scalar expTwo = exp(2); - const Scalar expYPlusOne = exp(y+1); - - sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX + 10.0; - sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX - -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX; - sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX - -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX; - - return sol; - } - - // \} - -private: - static constexpr Scalar eps_ = 1e-7; -}; -} // end namespace Dumux - -#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh b/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh deleted file mode 100644 index 92c6464af9..0000000000 --- a/test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh +++ /dev/null @@ -1,109 +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 3 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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup BoundaryTests - * \brief The spatial parameters class for the test problem using the 1p cc model. - */ - -#ifndef DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH -#define DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH - -#include -#include -#include -#include -#include - -namespace Dumux { - -/*! - * \ingroup BoundaryTests - * \brief The spatial parameters class for the test problem using the - * 1p cc model. - */ -template -class ConvergenceTestSpatialParams -: public FVSpatialParamsOneP> -{ - using GridView = typename GridGeometry::GridView; - using ParentType = FVSpatialParamsOneP>; - - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - static constexpr int dimWorld = GridView::dimensionworld; - using DimWorldMatrix = Dune::FieldMatrix; - - -public: - // export permeability type - using PermeabilityType = DimWorldMatrix; - - ConvergenceTestSpatialParams(std::shared_ptr gridGeometry) - : ParentType(gridGeometry) - {} - - /*! - * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. - * - * \param element The element - * \param scv The sub control volume - * \param elemSol The element solution vector - * \return the intrinsic permeability - */ - template - PermeabilityType permeability(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - PermeabilityType K(0.0); - - using std::cos; - using std::sin; - using std::exp; - - static constexpr Scalar c = 0.9; - static constexpr Scalar omega = M_PI; - - const Scalar x = scv.center()[0]; - K[0][0] = 1.0; - K[0][1] = -c/(2*omega) * sin(omega*x); - K[1][0] = K[0][1]; - K[1][1] = exp(-2)*(1 + c*cos(omega*x)); - - return K; - } - - /*! \brief Defines the porosity in [-]. - * - * \param globalPos The global position - */ - Scalar porosityAtPos(const GlobalPosition& globalPos) const - { return 0.4; } - -private: - Scalar permeability_; -}; - -} // end namespace Dumux - -#endif // DUMUX_CONVERGENCE_TEST_ONEP_SPATIALPARAMS_HH -- GitLab From c498aa47303e683c8c1bc0deb58cbb688a617153 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Fri, 9 Oct 2020 09:47:53 +0200 Subject: [PATCH 32/36] [fix][wmpfa][stencil] Make stencils unique --- dumux/discretization/fluxstencil.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh index 80fde0ccc5..454149cb93 100644 --- a/dumux/discretization/fluxstencil.hh +++ b/dumux/discretization/fluxstencil.hh @@ -177,7 +177,8 @@ public: stencil.push_back(scvfJ.outsideScvIdx()); } } - + std::sort(stencil.begin(), stencil.end()); + stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); return stencil; } }; -- GitLab From 060632c8a25b0cab6c5ced0af5b3a42d0f7a8cee Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Fri, 9 Oct 2020 11:53:11 +0200 Subject: [PATCH 33/36] [fix][wmpfa][indices] Safely erase entries in decomposition --- .../cellcentered/wmpfa/facedatahandle.hh | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index 4e8fbe77ea..38e83c4539 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -38,14 +38,22 @@ template void eraseZeros(Coeff& c, Stencil& s) { assert(c.size() == s.size()); + std::vector del(c.size(), false); for(std::size_t i=0; i Date: Mon, 12 Oct 2020 16:20:39 +0200 Subject: [PATCH 34/36] [wmpfa][decomposition] Allow the inclusion of the direct face neighbor --- dumux/common/vectordecomposition.hh | 62 ++++++++++++++++++- .../cellcentered/wmpfa/facedatahandle.hh | 5 +- 2 files changed, 64 insertions(+), 3 deletions(-) diff --git a/dumux/common/vectordecomposition.hh b/dumux/common/vectordecomposition.hh index 0254c55210..c6b2eb5a5d 100644 --- a/dumux/common/vectordecomposition.hh +++ b/dumux/common/vectordecomposition.hh @@ -32,6 +32,63 @@ class VectorDecomposition { public: + template + static std::tuple, std::vector, bool> calculateVectorDecomposition(const DV x, const std::vector& v, std::size_t localIdxInclusion) + { + std::size_t numVectors = v.size(); + static constexpr auto dim = DV::dimension; + if(numVectors == 0) + DUNE_THROW(Dune::InvalidStateException, "Can't perform decomposition without any given vectors!"); + + if(isAligned(x, v[localIdxInclusion])) + return {std::vector({localIdxInclusion}), calculateCoefficients(x, v[localIdxInclusion]), true}; + + using Scalar = typename DV::value_type; + std::vector coefficients(dim, std::numeric_limits::max()); + std::vector indices; + bool foundValidSolution = false; + + if constexpr (dim == 3) + { + for (std::size_t j = 0; j < numVectors - 1; ++j) + { + if(j == localIdxInclusion) + continue; + + for (std::size_t k = j+1; k < numVectors; ++k) + { + auto [coeff, valid] = calculateCoefficients(x, v[localIdxInclusion], v[j], v[k]); + using std::max_element; using std::move; + if(valid && (*max_element(coeff.begin(), coeff.end()) < *max_element(coefficients.begin(), coefficients.end()))) + { + coefficients = move(coeff); + indices = {localIdxInclusion, j, k}; + foundValidSolution = true; + } + } + } + } + else if constexpr (dim == 2) + { + for (std::size_t j = 0; j < numVectors; ++j) + { + if(j == localIdxInclusion) + continue; + + auto [coeff, valid] = calculateCoefficients(x, v[localIdxInclusion], v[j]); + using std::max_element; using std::move; + if(valid && (*max_element(coeff.begin(), coeff.end()) < *max_element(coefficients.begin(), coefficients.end()))) + { + coefficients = move(coeff); + indices = {localIdxInclusion, j}; + foundValidSolution = true; + } + } + } + + return {indices, coefficients, foundValidSolution}; + } + template static std::tuple, std::vector, bool> calculateVectorDecomposition(const DV x, const std::vector& v) { @@ -40,8 +97,9 @@ public: if(numVectors == 0) DUNE_THROW(Dune::InvalidStateException, "Can't perform decomposition without any given vectors!"); - if(isAligned(x, v[0])) - return {std::vector({0}), calculateCoefficients(x, v[0]), true}; + for (std::size_t i = 0; i < numVectors; ++i) + if(isAligned(x, v[i])) + return {std::vector({i}), calculateCoefficients(x, v[i]), true}; using Scalar = typename DV::value_type; std::vector coefficients(dim, std::numeric_limits::max()); diff --git a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh index 38e83c4539..8323170f2d 100644 --- a/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh +++ b/dumux/discretization/cellcentered/wmpfa/facedatahandle.hh @@ -101,9 +101,12 @@ public: template void decompose(const Problem& problem, const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf) { + static const bool enforceNeighborInclusion = getParamFromGroup(problem.paramGroup(), "WMPFA.EnforceNeighborInclusion", false); + boundaryFace_ = scvf.boundary(); const auto coNormal = mv(tensor(scvf.insideScvIdx()), scvf.unitOuterNormal()); - auto&& [indices, coeff, found] = VectorDecomposition::calculateVectorDecomposition(coNormal, intOp.getDistanceVectors(fvGeometry)); + auto&& [indices, coeff, found] = enforceNeighborInclusion ? VectorDecomposition::calculateVectorDecomposition(coNormal, intOp.getDistanceVectors(fvGeometry), scvf.localIndex()) + : VectorDecomposition::calculateVectorDecomposition(coNormal, intOp.getDistanceVectors(fvGeometry)); if(!found) DUNE_THROW(Dune::InvalidStateException, "CoNormal decomposition not found"); -- GitLab From 73d68e77ca768a640c2068ad93abac736a35b8b0 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Mon, 12 Oct 2020 18:00:31 +0200 Subject: [PATCH 35/36] [wmpfa][nlmpfa] First working implementation of nlmpfa --- dumux/flux/ccwmpfa/darcyslaw.hh | 76 +++++++++++++++++++++++++++++++-- 1 file changed, 72 insertions(+), 4 deletions(-) diff --git a/dumux/flux/ccwmpfa/darcyslaw.hh b/dumux/flux/ccwmpfa/darcyslaw.hh index ac991f57c0..bd680b1fa5 100644 --- a/dumux/flux/ccwmpfa/darcyslaw.hh +++ b/dumux/flux/ccwmpfa/darcyslaw.hh @@ -136,8 +136,8 @@ public: Scalar fluxIJ = 0.0; fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); std::for_each(dIJ.cbegin()+1, dIJ.cend(), - [&fluxIJ, &elemVolVars, &solValues](const auto& e) - { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); flux = scvf.area()*fluxIJ; } @@ -200,8 +200,76 @@ public: Scalar fluxIJ = 0.0; fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); std::for_each(dIJ.cbegin()+1, dIJ.cend(), - [&fluxIJ, &elemVolVars, &solValues](const auto& e) - { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); + + flux = scvf.area()*fluxIJ; + } +}; + +template<> +class WMpfaDarcysLawsFluxCalculator +{ +public: + + template + static void flux(Scalar& flux, + const SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const FluxData& dJI, + const SolValues& solValues) + { + const auto sI = solValues(elemVolVars[dIJ[0].index], dIJ[0].position); + const auto sJ = solValues(elemVolVars[dJI[0].index], dJI[0].position); + + Scalar betaIJ = 0.0; + Scalar betaJI = 0.0; + Scalar fluxIJ = 0.0; + Scalar fluxJI = 0.0; + std::for_each(dIJ.cbegin()+1, dIJ.cend(), + [&fluxIJ, &betaIJ, &sI, &scvf, &elemVolVars, &solValues](const auto& e) + { + (e.index == scvf.outsideScvIdx()) + ? betaIJ += e.coefficient + : fluxIJ += e.coefficient * (sI - solValues(elemVolVars[e.index], e.position)); + } ); + + std::for_each(dJI.cbegin()+1, dJI.cend(), + [&fluxJI, &betaJI, &sJ, &scvf, &elemVolVars, &solValues](const auto& e) + { + (e.index == scvf.insideScvIdx()) + ? betaJI += e.coefficient + : fluxJI += e.coefficient * (sJ - solValues(elemVolVars[e.index], e.position)); + } ); + + auto beta = std::max(0.0, std::min(betaIJ,betaJI)); + if(beta < 1e-30) + DUNE_THROW(Dune::InvalidStateException, "Beta coefficient must be greater than zero for NLMPFA"); + + fluxIJ += (betaIJ - beta)*(sI - sJ); + fluxJI += (betaJI - beta)*(sJ - sI); + + auto calcAbs = [](Scalar a){ return std::abs(a) + 1e-8; }; + + Scalar wIJ = (calcAbs(fluxJI))/(calcAbs(fluxIJ) + calcAbs(fluxJI)); + Scalar wJI = (calcAbs(fluxIJ))/(calcAbs(fluxIJ) + calcAbs(fluxJI)); + + flux = scvf.area()*( beta*(sI - sJ) + wIJ * fluxIJ - wJI * fluxJI); + } + + template + static void boundaryFlux(Scalar& flux, + const SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxData& dIJ, + const SolValues& solValues) + { + Scalar fluxIJ = 0.0; + fluxIJ += dIJ[0].coefficient * solValues(elemVolVars[dIJ[0].index], dIJ[0].position); + std::for_each(dIJ.cbegin()+1, dIJ.cend(), + [&fluxIJ, &elemVolVars, &solValues](const auto& e) + { fluxIJ -= e.coefficient * solValues(elemVolVars[e.index], e.position); } ); flux = scvf.area()*fluxIJ; } -- GitLab From e5f110c8d8387e3f0fd305b62cbeda84d654e839 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Mon, 12 Oct 2020 18:21:52 +0200 Subject: [PATCH 36/36] [wmpfa][nlmpfa] Update eps in flux calculation --- dumux/flux/ccwmpfa/darcyslaw.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/flux/ccwmpfa/darcyslaw.hh b/dumux/flux/ccwmpfa/darcyslaw.hh index bd680b1fa5..657ef02b1e 100644 --- a/dumux/flux/ccwmpfa/darcyslaw.hh +++ b/dumux/flux/ccwmpfa/darcyslaw.hh @@ -250,7 +250,7 @@ public: fluxIJ += (betaIJ - beta)*(sI - sJ); fluxJI += (betaJI - beta)*(sJ - sI); - auto calcAbs = [](Scalar a){ return std::abs(a) + 1e-8; }; + auto calcAbs = [](Scalar a){ return std::abs(a) + 1e-30; }; Scalar wIJ = (calcAbs(fluxJI))/(calcAbs(fluxIJ) + calcAbs(fluxJI)); Scalar wJI = (calcAbs(fluxIJ))/(calcAbs(fluxIJ) + calcAbs(fluxJI)); -- GitLab