diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index e79a41db5eca80335665f9585cb48af9f002a811..3be824293ef4db33eb2fb33d195895e6209847e9 100644 --- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh @@ -261,18 +261,19 @@ public: // reserve memory const auto numNeighbors = assemblyMapI.size(); - const auto estimate = neighborScvfEstimate_(element, numNeighbors); + const auto numNeighborScvfs = numNeighborScvfs_(assemblyMapI); neighborScvs_.reserve(numNeighbors); neighborScvIndices_.reserve(numNeighbors); - neighborScvfIndices_.reserve(estimate); - neighborScvfs_.reserve(estimate); + neighborScvfIndices_.reserve(numNeighborScvfs); + neighborScvfs_.reserve(numNeighborScvfs); // make neighbor geometries // use the assembly map to determine which faces are necessary for (auto&& dataJ : assemblyMapI) makeNeighborGeometries(globalFvGeometry().element(dataJ.globalJ), dataJ.globalJ, - dataJ.scvfsJ); + dataJ.scvfsJ, + dataJ.additionalScvfs); // set up the scvf flip indices for network grids if (dim < dimWorld) @@ -309,31 +310,13 @@ private: //! Estimates the number of neighboring scvfs that have to be prepared //! The estimate holds for inner elements, overestimats on the boundary - int neighborScvfEstimate_(const Element& element, int numNeighbors) + template<class DataJ> + unsigned int numNeighborScvfs_(const std::vector<DataJ>& dataJVector) { - auto type = element.geometry().type(); - - if (type == Dune::GeometryType(Dune::GeometryType::simplex, dim)) - { - if (dim == 2) - return 6 + ( (numNeighbors-3 > 0) ? (numNeighbors-3)*2 : 0 ); - if (dim == 3) - return 12 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 ); - } - else if (type == Dune::GeometryType(Dune::GeometryType::cube, dim)) - { - if (dim == 2) - return 16; - if (dim == 3) - return 120; - } - else if (type == Dune::GeometryType(Dune::GeometryType::pyramid, dim)) - return 16 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 ); - else if (type == Dune::GeometryType(Dune::GeometryType::prism, dim)) - return 18 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 ); - else - DUNE_THROW(Dune::NotImplemented, "Mpfa for given geometry type."); - + unsigned int numNeighborScvfs = 0; + for (const auto& dataJ : dataJVector) + numNeighborScvfs += dataJ.scvfsJ.size() + dataJ.additionalScvfs.size(); + return numNeighborScvfs; } //! create scvs and scvfs of the bound element @@ -424,7 +407,10 @@ private: //! create the scv and necessary scvfs of a neighboring element template<typename IndexVector> - void makeNeighborGeometries(const Element& element, IndexType eIdxGlobal, const IndexVector& scvfIndices) + void makeNeighborGeometries(const Element& element, + IndexType eIdxGlobal, + const IndexVector& scvfIndices, + const IndexVector& additionalScvfs) { // create the neighbor scv if it doesn't exist yet neighborScvs_.emplace_back(element.geometry(), eIdxGlobal); @@ -483,7 +469,8 @@ private: continue; // only build the scvf if it is in the list of necessary indices - if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter])) + if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter]) + || MpfaHelper::contains(additionalScvfs, scvFaceIndices[scvfCounter])) { neighborScvfs_.emplace_back(MpfaHelper(), MpfaHelper::getScvfCorners(isCorners, c), diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh index 1e2e74c3c46f5d5ab492eb56ef3125d4eb04d9cd..d9fe0f67c7a7c67b03f904fe58100e5f5c151d66 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh @@ -57,10 +57,6 @@ public: SeedVector& seeds, BoundarySeedVector& boundarySeeds) { - // print warning - std::cout << "\nWARNING: Note that the mpfa-l method currently only works for unstructured grids" - << " when the global volume variables and global fv geometry caches are enabled!\n" << std::endl; - seeds.clear(); boundarySeeds.clear(); scvfIndexMap.clear(); diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh index e56cb825c51042d6d3666643ac1281760895774e..6ab3030efd62216e261933b8b7cb2e4d56fe7e6c 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh @@ -277,43 +277,26 @@ private: { const auto& globalFvGeometry = problem_().model().globalFvGeometry(); - if (regionUnique_) + interactionRegions_.reserve(2); + auto&& scvSeed1 = seed.scvSeed(0); + auto&& scvSeed2 = seed.scvSeed(1); + auto&& outerScvSeed1 = seed.outerScvSeed(0); + auto&& outerScvSeed2 = seed.outerScvSeed(1); + auto e1 = globalFvGeometry.element(scvSeed1.globalIndex()); + auto e2 = globalFvGeometry.element(scvSeed2.globalIndex()); + auto e3 = globalFvGeometry.element(outerScvSeed1.globalIndex()); + auto e4 = globalFvGeometry.element(outerScvSeed2.globalIndex()); + + // scvSeed1 is the one the seed construction began at + if (scvSeed1.contiFaceLocalIdx() == 0) { - interactionRegions_.reserve(1); - auto&& scvSeed1 = seed.scvSeed(0); - auto&& scvSeed2 = seed.outerScvSeed(0); - auto&& scvSeed3 = seed.outerScvSeed(1); - auto e1 = globalFvGeometry.element(scvSeed1.globalIndex()); - auto e2 = globalFvGeometry.element(scvSeed2.globalIndex()); - auto e3 = globalFvGeometry.element(scvSeed3.globalIndex()); - if (scvSeed1.contiFaceLocalIdx() == 0) - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed2, scvSeed3, e1, e2, e3); - else - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed3, scvSeed2, e1, e3, e2); + interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, OuterScvSeedType(scvSeed2), outerScvSeed1, e1, e2, e3); + interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, outerScvSeed2, OuterScvSeedType(scvSeed1), e2, e4, e1); } else { - interactionRegions_.reserve(2); - auto&& scvSeed1 = seed.scvSeed(0); - auto&& scvSeed2 = seed.scvSeed(1); - auto&& outerScvSeed1 = seed.outerScvSeed(0); - auto&& outerScvSeed2 = seed.outerScvSeed(1); - auto e1 = globalFvGeometry.element(scvSeed1.globalIndex()); - auto e2 = globalFvGeometry.element(scvSeed2.globalIndex()); - auto e3 = globalFvGeometry.element(outerScvSeed1.globalIndex()); - auto e4 = globalFvGeometry.element(outerScvSeed2.globalIndex()); - - // scvSeed1 is the one the seed construction began at - if (scvSeed1.contiFaceLocalIdx() == 0) - { - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, OuterScvSeedType(scvSeed2), outerScvSeed1, e1, e2, e3); - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, outerScvSeed2, OuterScvSeedType(scvSeed1), e2, e4, e1); - } - else - { - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, outerScvSeed1, OuterScvSeedType(scvSeed2), e1, e3, e2); - interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, OuterScvSeedType(scvSeed1), outerScvSeed2, e2, e1, e4); - } + interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, outerScvSeed1, OuterScvSeedType(scvSeed2), e1, e3, e2); + interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, OuterScvSeedType(scvSeed1), outerScvSeed2, e2, e1, e4); } // as an initial guess we set the first interaction region as the one used diff --git a/dumux/implicit/cellcentered/assemblymap.hh b/dumux/implicit/cellcentered/assemblymap.hh index 58c40166350bd8dbb09184e0f71891be73101beb..eaf1ee205d8ffba74960161ea98ad8f6cdfdc363 100644 --- a/dumux/implicit/cellcentered/assemblymap.hh +++ b/dumux/implicit/cellcentered/assemblymap.hh @@ -22,8 +22,8 @@ * that contribute to the derivative calculation. This is used for * finite-volume schemes with symmetric sparsity pattern in the global matrix. */ -#ifndef DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH -#define DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH +#ifndef DUMUX_CC_ASSEMBLY_MAP_HH +#define DUMUX_CC_ASSEMBLY_MAP_HH #include <dune/istl/bcrsmatrix.hh> @@ -32,8 +32,18 @@ namespace Dumux { +/*! + * \ingroup CellCentered + * \brief A simple version of the assembly map for cellcentered schemes. + * This implementation works for schemes in which for a given cell I only + * those cells J have to be prepared in whose stencil the cell I appears. + * This means that for the flux calculations in the cells J (in order to compute + * the derivatives with respect to cell I), we do not need data on any additional cells J + * to compute these fluxes. The same holds for scvfs in the cells J, i.e. we need only those + * scvfs in the cells J in which the cell I is in the stencil. + */ template<class TypeTag> -class CCSymmetricAssemblyMap +class CCSimpleAssemblyMap { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); @@ -45,6 +55,9 @@ class CCSymmetricAssemblyMap { IndexType globalJ; std::vector<IndexType> scvfsJ; + // A list of additional scvfs is needed for compatibility + // reasons with more complex assembly maps (see mpfa) + std::vector<IndexType> additionalScvfs; }; using Map = std::vector<std::vector<DataJ>>; diff --git a/dumux/implicit/cellcentered/mpfa/assemblymap.hh b/dumux/implicit/cellcentered/mpfa/assemblymap.hh index cf9b374fba30f5b29447b9a43a31bcd151e5011d..76b6bca57a3928ff6ec4e86903acff330dae3132 100644 --- a/dumux/implicit/cellcentered/mpfa/assemblymap.hh +++ b/dumux/implicit/cellcentered/mpfa/assemblymap.hh @@ -25,6 +25,7 @@ #define DUMUX_CC_MPFA_ASSEMBLY_MAP_HH #include <dumux/implicit/cellcentered/assemblymap.hh> +#include <dumux/implicit/cellcentered/mpfa/generalassemblymap.hh> namespace Dumux { @@ -36,9 +37,13 @@ class CCMpfaAssemblyMapImplementation; template<class TypeTag> using CCMpfaAssemblyMap = CCMpfaAssemblyMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>; -//! The default is simply the CCAssemblyMap +//! The default is the general assembly map for mpfa schemes template<class TypeTag, MpfaMethods method> -class CCMpfaAssemblyMapImplementation : public CCAssemblyMap<TypeTag> {}; +class CCMpfaAssemblyMapImplementation : public CCMpfaGeneralAssemblyMap<TypeTag> {}; + +//! The o-method can use the simple assembly map +template<class TypeTag> +class CCMpfaAssemblyMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleAssemblyMap<TypeTag> {}; } #endif diff --git a/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh b/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh new file mode 100644 index 0000000000000000000000000000000000000000..c69911e08de49ca1dc55e01f785ea6bc8f3421e7 --- /dev/null +++ b/dumux/implicit/cellcentered/mpfa/generalassemblymap.hh @@ -0,0 +1,201 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \brief Stores the face indices corresponding to the neighbors of an element + * that contribute to the derivative calculation + */ +#ifndef DUMUX_CC_MPFA_GENERAL_ASSEMBLY_MAP_HH +#define DUMUX_CC_MPFA_GENERAL_ASSEMBLY_MAP_HH + +#include <dumux/implicit/cellcentered/assemblymap.hh> + +namespace Dumux +{ + +/*! + * \ingroup CellCentered + * \brief General version of the assembly map for cellcentered schemes. To each + * cell I we store a list of cells J that are needed to compute the fluxes + * in these cells J that depend on cell I. Furthermore, we store for each cell J + * a list of scvfs in which cell I is in the stencil, as well as additional scvfs + * that are also required to set up the transmissibilities. + */ +template<class TypeTag> +class CCMpfaGeneralAssemblyMap +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + + using IndexType = typename GridView::IndexSet::IndexType; + + // To each cell "globalI" there will be a list of "globalJ", in which globalI is part + // of the stencil. We save the scvfs over which fluxes depend on globalI and a list of + // additional scvfs which are needed temporarily to set up the transmissibilities of the scvfsJ + struct DataJ + { + IndexType globalJ; + std::vector<IndexType> scvfsJ; + std::vector<IndexType> additionalScvfs; + }; + + using Map = std::vector<std::vector<DataJ>>; + +public: + + /*! + * \brief Initialize the AssemblyMap object. + * + * \param problem The problem which we want to simulate. + */ + void init(const Problem& problem) + { + map_.resize(problem.gridView().size(0)); + for (const auto& element : elements(problem.gridView())) + { + // We are looking for the elements I, for which this element J is in the flux stencil + auto globalJ = problem.elementMapper().index(element); + + auto fvGeometry = localView(problem.model().globalFvGeometry()); + fvGeometry.bindElement(element); + + // obtain the data of J in elements I + std::vector<std::pair<IndexType, std::vector<DataJ>>> dataJForI; + + // loop over sub control faces + for (auto&& scvf : scvfs(fvGeometry)) + { + FluxVariables fluxVars; + const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); + + // insert our index in the neighbor stencils of the elements in the flux stencil + for (auto globalI : stencil) + { + if (globalI == globalJ) + continue; + + auto it = std::find_if(dataJForI.begin(), + dataJForI.end(), + [globalI](const auto& pair) { return pair.first == globalI; }); + + if (it != dataJForI.end()) + { + // get the data J which corresponds to the actual global J + // This will be the first entry, as we do so in the else statement (see below) + auto& globalJDataJ = it->second[0]; + + // insert actual scvf index in the list of scvfs which couple I and J + globalJDataJ.scvfsJ.push_back(scvf.index()); + + // Also, all scvfs connected to a vertex together with the actual scvf + // land in the list of additional scvfs. Of that list we will delete those + // that are already in the list of scvfsJ later... + const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); + std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); + for (unsigned int i = 0; i < scvfVectorAtVertex.size(); ++i) + scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index(); + globalJDataJ.additionalScvfs.insert(globalJDataJ.additionalScvfs.end(), + scvfIndicesAtVertex.begin(), + scvfIndicesAtVertex.end()); + + // all the other dofs in the stencil have to appear as "globalJ" to globalI as well + for (auto globalJ2 : stencil) + { + if (globalJ2 == globalJ || globalJ2 == globalI) + continue; + + auto it2 = std::find_if(it->second.begin(), + it->second.end(), + [globalJ2](const auto& dataJ) { return dataJ.globalJ == globalJ2; }); + + // if entry for globalJ2 does not exist yet, add globalJ2 to the J-data of globalI + // with an empty set of scvfs over which I and J are coupled (i.e. they aren't coupled) + if (it2 == it->second.end()) + it->second.push_back(DataJ({globalJ2, std::vector<IndexType>()})); + } + } + else + { + // No DataJ for globalI exists yet. Make it and insert data on the actual + // global J as first entry in the vector of DataJs belonging to globalI + dataJForI.emplace_back(std::make_pair(globalI, + std::vector<DataJ>({DataJ({globalJ, std::vector<IndexType>({scvf.index()})})}))); + + // Also, all scvfs connected to a vertex together with the actual scvf + // land in the list of additional scvfs. Of that list we will delete those + // that are already in the list of scvfsJ later... + const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); + std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); + for (unsigned int i = 0; i < scvfVectorAtVertex.size(); ++i) + scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index(); + dataJForI.back().second[0].additionalScvfs.insert(dataJForI.back().second[0].additionalScvfs.end(), + scvfIndicesAtVertex.begin(), + scvfIndicesAtVertex.end()); + + // all the other dofs in the stencil will be "globalJ" to globalI as well + for (auto globalJ2 : stencil) + if (globalJ2 != globalJ && globalJ2 != globalI) + dataJForI.back().second.push_back(DataJ({globalJ2, std::vector<IndexType>()})); + } + } + } + + // Insert the data into the global map + for (auto&& pair : dataJForI) + { + // obtain the corresponding entry in the map + auto& dataJVector = map_[pair.first]; + for (auto&& dataJ : pair.second) + { + // delete those additionalScvfs indices that are already in the list of scvfs + dataJ.additionalScvfs.erase(std::remove_if(dataJ.additionalScvfs.begin(), + dataJ.additionalScvfs.end(), + [&dataJ] (const auto& idx) { return MpfaHelper::contains(dataJ.scvfsJ, idx); }), + dataJ.additionalScvfs.end()); + + // if entry for j exists in the map already add scvf and additional scvf indices, create otherwise + auto it = std::find_if(dataJVector.begin(), + dataJVector.end(), + [&dataJ](const auto& dataJofMap) { return dataJofMap.globalJ == dataJ.globalJ; }); + + if (it != dataJVector.end()) + { + it->scvfsJ.insert(it->scvfsJ.end(), dataJ.scvfsJ.begin(), dataJ.scvfsJ.end()); + it->additionalScvfs.insert(it->additionalScvfs.end(), dataJ.additionalScvfs.begin(), dataJ.additionalScvfs.end()); + } + else + dataJVector.emplace_back(std::move(dataJ)); + } + } + } + } + + const std::vector<DataJ>& operator[] (const IndexType globalI) const + { return map_[globalI]; } + +private: + Map map_; +}; +} + +#endif diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh index 483f831e8367bbe97157c1a43e74944574674f9d..0f0a0dbc26d5b7728dc8aa0605630c2f8dfca79d 100644 --- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh @@ -43,6 +43,7 @@ #include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh> #include <dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh> #include <dumux/implicit/cellcentered/mpfa/localresidual.hh> +#include <dumux/implicit/cellcentered/mpfa/assemblymap.hh> #include <dumux/implicit/cellcentered/properties.hh> namespace Dumux { @@ -66,6 +67,9 @@ SET_PROP(CCMpfaModel, MpfaMethod) //! The mpfa helper class SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>); +//! The assembly map for mpfa schemes +SET_TYPE_PROP(CCMpfaModel, AssemblyMap, CCMpfaAssemblyMap<TypeTag>); + //! The interaction volume class SET_TYPE_PROP(CCMpfaModel, InteractionVolume, CCMpfaInteractionVolume<TypeTag>); diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh index 702259c31e8b32641980917887ad55f61c822ec8..99fbb8925c9fadc81dfe84ac82ad9ae3929d5e99 100644 --- a/dumux/implicit/cellcentered/propertydefaults.hh +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -76,8 +76,8 @@ SET_TYPE_PROP(CCModel, GlobalVolumeVariables, CCGlobalVolumeVariables<TypeTag, G //! Set the BaseLocalResidual to CCLocalResidual SET_TYPE_PROP(CCModel, BaseLocalResidual, CCLocalResidual<TypeTag>); -//! Set the AssemblyMap to the CCAssemblyMap (default: symmetric occupation pattern) -SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSymmetricAssemblyMap<TypeTag>); +//! Set the AssemblyMap to the CCAssemblyMap per default +SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSimpleAssemblyMap<TypeTag>); // By default, we disable interior boundaries SET_BOOL_PROP(CCModel, EnableInteriorBoundaries, false);