diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt index 819d019ad5f37d819bd128854c2892abba2684c9..1aae30ec2a946ece00c478cd893d618c9d53c002 100644 --- a/dumux/CMakeLists.txt +++ b/dumux/CMakeLists.txt @@ -12,4 +12,5 @@ add_subdirectory(material) add_subdirectory(multidomain) add_subdirectory(nonlinear) add_subdirectory(parallel) +add_subdirectory(porenetworkflow) add_subdirectory(porousmediumflow) diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index a43e90fb378a589d1ad6eb942091d1764ad5af9a..6fbdf151abcf8e514b6221fe8a4701f389595d22 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -1,6 +1,7 @@ add_subdirectory(box) add_subdirectory(cellcentered) add_subdirectory(fem) +add_subdirectory(porenetwork) add_subdirectory(projection) add_subdirectory(staggered) diff --git a/dumux/discretization/porenetwork/CMakeLists.txt b/dumux/discretization/porenetwork/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..27a356b4c2a9c6a29d6c5b3eb4e32e65f1390d2f --- /dev/null +++ b/dumux/discretization/porenetwork/CMakeLists.txt @@ -0,0 +1,6 @@ +install(FILES +fvelementgeometry.hh +gridgeometry.hh +subcontrolvolume.hh +subcontrolvolumeface.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/porenetwork) diff --git a/dumux/discretization/porenetwork/fvelementgeometry.hh b/dumux/discretization/porenetwork/fvelementgeometry.hh new file mode 100644 index 0000000000000000000000000000000000000000..735fe27e5c5cd7ee951f2ae41ea60b25a65e6968 --- /dev/null +++ b/dumux/discretization/porenetwork/fvelementgeometry.hh @@ -0,0 +1,330 @@ +// -*- 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 PoreNetworkDiscretization + * \brief Base class for the local geometry for porenetworks + */ +#ifndef DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH +#define DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH + +#include + +#include +#include + +namespace Dumux { + +/*! + * \ingroup PoreNetworkDiscretization + * \brief Base class for the local geometry for porenetworks + * \tparam GG the finite volume grid geometry type + * \tparam enableFVGridGeometryCache if the grid geometry is cached or not + */ +template +class PNMFVElementGeometry; + +//! specialization in case the FVElementGeometries are stored +template +class PNMFVElementGeometry +{ + using GridView = typename GG::GridView; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + using GridIndexType = typename IndexTraits::GridIndex; + using LocalIndexType = typename IndexTraits::LocalIndex; + using Element = typename GridView::template Codim<0>::Entity; + using CoordScalar = typename GridView::ctype; + using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType; +public: + //! 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; + using FVGridGeometry [[deprecated("Use GridGeometry")]] = GG; + //! the maximum number of scvs per element + static constexpr std::size_t maxNumElementScvs = 2; + + //! Constructor + PNMFVElementGeometry(const GridGeometry& gridGeometry) + : gridGeometryPtr_(&gridGeometry) {} + + //! Get a sub control volume with a local scv index + const SubControlVolume& scv(LocalIndexType scvIdx) const + { + return gridGeometry().scvs(eIdx_)[scvIdx]; + } + + //! Get a sub control volume face with a local scvf index + const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const + { + return gridGeometry().scvfs(eIdx_)[scvfIdx]; + } + + //! iterator range for sub control volumes. Iterates over + //! all scvs of the bound element. + //! 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 auto scvs(const PNMFVElementGeometry& fvGeometry) + { + const auto& g = fvGeometry.gridGeometry(); + using Iter = typename std::decay_t::const_iterator; + return Dune::IteratorRange(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end()); + } + + //! iterator range for sub control volumes faces. Iterates over + //! all scvfs of the bound element. + //! 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 auto scvfs(const PNMFVElementGeometry& fvGeometry) + { + const auto& g = fvGeometry.gridGeometry(); + using Iter = typename std::decay_t::const_iterator; + return Dune::IteratorRange(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end()); + } + + //! Get a local finite element basis + const FeLocalBasis& feLocalBasis() const + { + return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis(); + } + + //! The total number of sub control volumes + std::size_t numScv() const + { + return gridGeometry().scvs(eIdx_).size(); + } + + //! The total number of sub control volume faces + std::size_t numScvf() const + { + return gridGeometry().scvfs(eIdx_).size(); + } + + //! this function is for compatibility reasons with cc methods + //! The box stencil is always element-local so bind and bindElement + //! are identical. + void bind(const Element& element) + { + this->bindElement(element); + } + + //! Binding of an element, has to be called before using the fvgeometries + //! Prepares all the volume variables within the element + //! For compatibility reasons with the FVGeometry cache being disabled + void bindElement(const Element& element) + { + elementPtr_ = &element; + eIdx_ = gridGeometry().elementMapper().index(element); + } + + //! The global finite volume geometry we are a restriction of + [[deprecated("Use gridGeometry")]] + const GridGeometry& fvGridGeometry() const + { return *gridGeometryPtr_; } + + //! 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(eIdx_); } + +private: + const Element* elementPtr_; + const GridGeometry* gridGeometryPtr_; + + GridIndexType eIdx_; +}; + +//! specialization in case the FVElementGeometries are not stored +template +class PNMFVElementGeometry +{ + using GridView = typename GG::GridView; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + using GridIndexType = typename IndexTraits::GridIndex; + using LocalIndexType = typename IndexTraits::LocalIndex; + using Element = typename GridView::template Codim<0>::Entity; + using CoordScalar = typename GridView::ctype; + using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType; + +public: + //! 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; + using FVGridGeometry [[deprecated("Use GridGeometry")]] = GG; + //! the maximum number of scvs per element + static constexpr std::size_t maxNumElementScvs = 2; + + //! Constructor + PNMFVElementGeometry(const GridGeometry& gridGeometry) + : gridGeometryPtr_(&gridGeometry) {} + + //! Get a sub control volume with a local scv index + const SubControlVolume& scv(LocalIndexType scvIdx) const + { + return scvs_[scvIdx]; + } + + //! Get a sub control volume face with a local scvf index + const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const + { + return scvfs_[0]; + } + + //! iterator range for sub control volumes. Iterates over + //! all scvs of the bound element. + //! 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 auto scvs(const PNMFVElementGeometry& fvGeometry) + { + using Iter = typename std::decay_t::const_iterator; + return Dune::IteratorRange(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end()); + } + + //! iterator range for sub control volumes faces. Iterates over + //! all scvfs of the bound element. + //! 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 auto scvfs(const PNMFVElementGeometry& fvGeometry) + { + using Iter = typename std::decay_t::const_iterator; + return Dune::IteratorRange(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end()); + } + + //! Get a local finite element basis + const FeLocalBasis& feLocalBasis() const + { + return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis(); + } + + //! 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(); + } + + //! this function is for compatibility reasons with cc methods + //! The box stencil is always element-local so bind and bindElement + //! are identical. + void bind(const Element& element) + { + this->bindElement(element); + } + + //! Binding of an element, has to be called before using the fvgeometries + //! Prepares all the volume variables within the element + //! For compatibility reasons with the FVGeometry cache being disabled + void bindElement(const Element& element) + { + elementPtr_ = &element; + eIdx_ = gridGeometry().elementMapper().index(element); + makeElementGeometries(element); + } + + //! The global finite volume geometry we are a restriction of + [[deprecated("Use gridGeometry")]] + const GridGeometry& fvGridGeometry() const + { return *gridGeometryPtr_; } + + //! 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 hasBoundaryScvf_; } + +private: + + void makeElementGeometries(const Element& element) + { + hasBoundaryScvf_ = false; + + // get the element geometry + auto elementGeometry = element.geometry(); + + // construct the sub control volumes + for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx) + { + // get asssociated dof index + const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim); + + // get the corners + auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()}; + + // get the fractional volume asssociated with the scv + const auto volume = gridGeometry().poreVolume(dofIdxGlobal) / gridGeometry().coordinationNumber(dofIdxGlobal); + + // add scv to the local container + scvs_[scvLocalIdx] = SubControlVolume(dofIdxGlobal, + scvLocalIdx, + eIdx_, + std::move(corners), + volume); + + if (gridGeometry().poreLabel(dofIdxGlobal) > 0) + hasBoundaryScvf_ = true; + } + + // construct the inner sub control volume face + auto unitOuterNormal = elementGeometry.corner(1)-elementGeometry.corner(0); + unitOuterNormal /= unitOuterNormal.two_norm(); + LocalIndexType scvfLocalIdx = 0; + scvfs_[0] = SubControlVolumeFace(elementGeometry.center(), + std::move(unitOuterNormal), + scvfLocalIdx++, + std::vector({0, 1})); + } + + //! The bound element + const Element* elementPtr_; + GridIndexType eIdx_; + + //! The global geometry this is a restriction of + const GridGeometry* gridGeometryPtr_; + + //! vectors to store the geometries locally after binding an element + std::array scvs_; + std::array scvfs_; + + bool hasBoundaryScvf_ = false; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/porenetwork/gridgeometry.hh b/dumux/discretization/porenetwork/gridgeometry.hh new file mode 100644 index 0000000000000000000000000000000000000000..ff17acacc6ad2f9f1d74a96a0ee3ae9228969f63 --- /dev/null +++ b/dumux/discretization/porenetwork/gridgeometry.hh @@ -0,0 +1,802 @@ +// -*- 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 PoreNetworkDiscretization + * \brief Base class for the finite volume geometry for porenetwork models + */ +#ifndef DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH +#define DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace Dumux { + +/*! + * \ingroup PoreNetworkDiscretization + * \brief Base class for geometry data extraction from the grid data format + */ +template +class DefaultPNMData +{ + using GridIndex = typename IndexTraits::GridIndex; + using SmallLocalIndex = typename IndexTraits::SmallLocalIndex; + using Label = std::int_least8_t; + using Vertex = typename GridView::template Codim::Entity; + using Element = typename GridView::template Codim<0>::Entity; + + static const int dim = GridView::dimension; + +public: + + template + void update(const GridView& gridView, const GridData& gridData) + { + coordinationNumber_ = gridData.getCoordinationNumbers(); + + const auto numThroats = gridView.size(0); + throatRadius_.resize(numThroats); + throatLength_.resize(numThroats); + throatLabel_.resize(numThroats); + throatCrossSectionalArea_.resize(numThroats); + throatShapeFactor_.resize(numThroats); + + useSameGeometryForAllPores_ = true; + useSameShapeForAllThroats_ = true; + + // first check if the same geometry shall be used for all entities ... + if (hasParamInGroup(gridData.paramGroup(), "Grid.ThroatCrossSectionShape")) + { + const auto throatGeometryInput = getParamFromGroup(gridData.paramGroup(), "Grid.ThroatCrossSectionShape"); + const auto throatGeometry = Throat::shapeFromString(throatGeometryInput); + throatGeometry_.resize(1); + throatGeometry_[0] = throatGeometry; + + std::cout << "Using '" << throatGeometryInput << "' as cross-sectional shape for all throats." << std::endl; + } + else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector + { + std::cout << "Reading shape factors for throats from grid data." << std::endl; + useSameShapeForAllThroats_ = false; + throatGeometry_.resize(numThroats); + } + + // get the vertex parameters + const auto numPores = gridView.size(dim); + poreRadius_.resize(numPores); + poreLabel_.resize(numPores); + poreVolume_.resize(numPores); + + // first check if the same geometry shall be used for all entities ... + if (hasParamInGroup(gridData.paramGroup(), "Grid.PoreGeometry")) + { + const auto poreGeometryInput = getParamFromGroup(gridData.paramGroup(), "Grid.PoreGeometry"); + poreGeometry_.resize(1); + poreGeometry_[0] = Pore::shapeFromString(poreGeometryInput);; + + std::cout << "Using '" << poreGeometryInput << "' as geometry for all pores." << std::endl; + } + else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector + { + std::cout << "Reading pore shapes from grid data." << std::endl; + useSameGeometryForAllPores_ = false; + poreGeometry_.resize(numPores); + } + + + for (const auto& vertex : vertices(gridView)) + { + static const auto poreRadiusIdx = gridData.parameterIndex("PoreRadius"); + static const auto poreLabelIdx = gridData.parameterIndex("PoreLabel"); + const auto vIdx = gridView.indexSet().index(vertex); + const auto& params = gridData.parameters(vertex); + poreRadius_[vIdx] = params[poreRadiusIdx]; + assert(poreRadius_[vIdx] > 0.0); + poreLabel_[vIdx] = params[poreLabelIdx]; + + if (!useSameGeometryForAllPores()) + poreGeometry_[vIdx] = getPoreGeometry_(gridData, vertex); + + poreVolume_[vIdx] = getPoreVolume_(gridData, vertex, vIdx); + } + + for (const auto& element : elements(gridView)) + { + const int eIdx = gridView.indexSet().index(element); + const auto& params = gridData.parameters(element); + static const auto throatRadiusIdx = gridData.parameterIndex("ThroatRadius"); + static const auto throatLengthIdx = gridData.parameterIndex("ThroatLength"); + throatRadius_[eIdx] = params[throatRadiusIdx]; + throatLength_[eIdx] = params[throatLengthIdx]; + + // use a default value if no throat label is given by the grid + static const bool gridHasThroatLabel = gridData.gridHasElementParameter("ThroatLabel"); + if (gridHasThroatLabel) + { + static const auto throatLabelIdx = gridData.parameterIndex("ThroatLabel"); + throatLabel_[eIdx] = params[throatLabelIdx]; + } + else + { + const auto vIdx0 = gridView.indexSet().subIndex(element, 0, dim); + const auto vIdx1 = gridView.indexSet().subIndex(element, 1, dim); + + const auto poreLabel0 = poreLabel(vIdx0); + const auto poreLabel1 = poreLabel(vIdx1); + + if (poreLabel0 >= 0 && poreLabel1 >= 0) + { + std::cout << "\n Warning: Throat " + << eIdx << " connects two boundary pores with different pore labels. Using the greater pore label as throat label.\n" + << "Set the throat labels explicitly in your grid file, if needed." << std::endl; + } + + using std::max; + throatLabel_[eIdx] = max(poreLabel0, poreLabel1); + } + + if (!useSameShapeForAllThroats()) + { + static const auto throatShapeFactorIdx = gridData.parameterIndex("ThroatShapeFactor"); + static const auto throatAreaIdx = gridData.parameterIndex("ThroatCrossSectionalArea"); + throatShapeFactor_[eIdx] = params[throatShapeFactorIdx]; + throatGeometry_[eIdx] = Throat::shape(throatShapeFactor_[eIdx]); + throatCrossSectionalArea_[eIdx] = params[throatAreaIdx]; + } + else + { + throatCrossSectionalArea_[eIdx] = getThroatCrossSectionalArea_(gridData, element, eIdx); + throatShapeFactor_[eIdx] = getThroatShapeFactor_(gridData, element, eIdx); + } + + assert(throatRadius_[eIdx] > 0.0); + assert(throatLength_[eIdx] > 0.0); + assert(throatCrossSectionalArea_[eIdx] > 0.0); + + static const bool addThroatVolumeToPoreVolume = getParamFromGroup(gridData.paramGroup(), "Grid.AddThroatVolumeToPoreVolume", false); + if (addThroatVolumeToPoreVolume) + { + for (int vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal) + { + const auto vIdx = gridView.indexSet().subIndex(element, vIdxLocal, dim); + poreVolume_[vIdx] += 0.5 * throatCrossSectionalArea_[eIdx] * throatLength_[eIdx]; + } + } + } + + maybeResizeContainers_(); + } + + //! Returns the pore label (e.g. used for setting BCs) + Label poreLabel(const GridIndex dofIdxGlobal) const + { return poreLabel_[dofIdxGlobal]; } + + //! Returns the vector of pore labels + const std::vector