From 7e084f4a2b2d166cd9d9b9c8907e03771f331eab Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 Jan 2017 19:37:15 +0100 Subject: [PATCH] remove stencil concept The stencils are no longer used. Instead, we make use of the assembly map in cc models, where needed. For the box method the stencils concept was not necessary anyway. --- .../mpfa/elementfluxvariablescache.hh | 27 +- .../mpfa/elementvolumevariables.hh | 11 +- .../cellcentered/mpfa/fvelementgeometry.hh | 21 +- dumux/discretization/cellcentered/stencils.hh | 259 ------------------ .../tpfa/elementfluxvariablescache.hh | 30 +- .../tpfa/elementvolumevariables.hh | 31 +-- dumux/implicit/assembler.hh | 22 +- dumux/implicit/box/assembler.hh | 51 ++-- dumux/implicit/box/propertydefaults.hh | 5 - dumux/implicit/cellcentered/assembler.hh | 46 ++-- dumux/implicit/cellcentered/assemblymap.hh | 88 +++--- dumux/implicit/cellcentered/localjacobian.hh | 23 +- .../cellcentered/mpfa/propertydefaults.hh | 4 - .../implicit/cellcentered/propertydefaults.hh | 4 +- .../cellcentered/tpfa/propertydefaults.hh | 4 - dumux/implicit/model.hh | 15 - dumux/implicit/properties.hh | 3 - .../implicit/velocityoutput.hh | 21 +- 18 files changed, 160 insertions(+), 505 deletions(-) delete mode 100644 dumux/discretization/cellcentered/stencils.hh diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 9069ebeff3..01d8be6c03 100644 --- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh @@ -141,13 +141,14 @@ public: const auto& problem = globalFluxVarsCache().problem_(); const auto& globalFvGeometry = problem.model().globalFvGeometry(); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); - const auto& assemblyMap = problem.model().localJacobian().assemblyMap(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; // reserve memory - auto numNeighborScvfs = 0; - for (auto&& facesInNeighbor : assemblyMap[globalI]) numNeighborScvfs += facesInNeighbor.size(); + unsigned int numNeighborScvfs = 0; + for (auto&& dataJ : assemblyMapI) + numNeighborScvfs += dataJ.scvfsJ.size(); globalScvfIndices_.reserve(fvGeometry.numScvf() + numNeighborScvfs); // first add all the indices inside the element @@ -155,9 +156,9 @@ public: globalScvfIndices_.push_back(scvf.index()); // for the indices in the neighbors, use assembly map of the local jacobian - for (unsigned int j = 0; j < neighborStencil.size(); ++j) - for (auto fluxVarIdx : assemblyMap[globalI][j]) - globalScvfIndices_.push_back(fluxVarIdx); + for (auto&& dataJ : assemblyMapI) + for (auto scvfIdx : dataJ.scvfsJ) + globalScvfIndices_.push_back(scvfIdx); // prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class fluxVarsCache_.resize(globalScvfIndices_.size()); @@ -166,21 +167,17 @@ public: FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this); // prepare the caches in the remaining neighbors - unsigned int j = 0; - for (auto globalJ : neighborStencil) + for (auto&& dataJ : assemblyMapI) { - for (auto fluxVarIdx : assemblyMap[globalI][j]) + for (auto scvfIdx : dataJ.scvfsJ) { - const auto& scvf = fvGeometry.scvf(fluxVarIdx); + const auto& scvf = fvGeometry.scvf(scvfIdx); if (!(*this)[scvf].isUpdated()) { - auto elementJ = globalFvGeometry.element(globalJ); + auto elementJ = globalFvGeometry.element(dataJ.globalJ); FluxVariablesCacheFiller::fillFluxVarCache(problem, elementJ, fvGeometry, elemVolVars, scvf, *this); } } - - // increment counter - j++; } } diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 3d20dfa2e6..1046a21ce7 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -121,8 +121,9 @@ public: const auto& globalFvGeometry = problem.model().globalFvGeometry(); // stencil information - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); - const auto numDofs = neighborStencil.size() + 1; + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numDofs = assemblyMapI.size() + 1; // resize local containers to the required size (for internal elements) volumeVariables_.resize(numDofs); @@ -140,10 +141,10 @@ public: ++localIdx; // Update the volume variables of the neighboring elements - for (auto globalJ : neighborStencil) + for (auto&& dataJ : assemblyMapI) { - const auto& elementJ = globalFvGeometry.element(globalJ); - auto&& scvJ = fvGeometry.scv(globalJ); + const auto& elementJ = globalFvGeometry.element(dataJ.globalJ); + auto&& scvJ = fvGeometry.scv(dataJ.globalJ); volumeVariables_[localIdx].update(problem.model().elementSolution(elementJ, sol), problem, elementJ, diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index b7e74e6b68..c5bf3269a4 100644 --- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh @@ -256,30 +256,23 @@ public: // get some references for convenience const auto& problem = globalFvGeometry().problem_(); - const auto& assemblyMap = problem.model().localJacobian().assemblyMap(); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; // reserve memory - const auto numNeighbors = neighborStencil.size(); + const auto numNeighbors = assemblyMapI.size(); const auto estimate = neighborScvfEstimate_(element, numNeighbors); neighborScvs_.reserve(numNeighbors); neighborScvIndices_.reserve(numNeighbors); neighborScvfIndices_.reserve(estimate); neighborScvfs_.reserve(estimate); - // build scvfs in neighbors + // make neighbor geometries // use the assembly map to determine which faces are necessary - unsigned int j = 0; - for (auto globalJ : neighborStencil) - { - // make neighbor geometries - auto elementJ = globalFvGeometry().element(globalJ); - makeNeighborGeometries(elementJ, globalJ, assemblyMap[globalI][j]); - - // increment counter - j++; - } + for (auto&& dataJ : assemblyMapI) + makeNeighborGeometries(globalFvGeometry().element(dataJ.globalJ), + dataJ.globalJ, + dataJ.scvfsJ); // set up the scvf flip indices for network grids if (dim < dimWorld) diff --git a/dumux/discretization/cellcentered/stencils.hh b/dumux/discretization/cellcentered/stencils.hh deleted file mode 100644 index 49d70b767a..0000000000 --- a/dumux/discretization/cellcentered/stencils.hh +++ /dev/null @@ -1,259 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/***************************************************************************** - * See the file COPYING for full copying permissions. * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see . * - *****************************************************************************/ -/*! - * \file - * \brief Implements the notion of stencils for cell-centered models - */ -#ifndef DUMUX_DISCRETIZATION_CC_STENCILS_HH -#define DUMUX_DISCRETIZATION_CC_STENCILS_HH - -#include - -namespace Dumux -{ - -/*! - * \brief Element-related stencils for symmetric cc methods. - */ -template -class CCSymmetricElementStencils -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using IndexType = typename GridView::IndexSet::IndexType; - using Element = typename GridView::template Codim<0>::Entity; - // TODO a separate stencil class all stencils can derive from? - using Stencil = std::vector; -public: - //! Update the stencil. - void update(const Problem& problem, - const Element& element) - { - // restrict the FvGeometry locally and bind to the element - auto fvGeometry = localView(problem.model().globalFvGeometry()); - fvGeometry.bindElement(element); // for TPFA bind element is enough - - elementStencil_.clear(); - - // loop over sub control faces - for (auto&& scvf : scvfs(fvGeometry)) - { - FluxVariables fluxVars; - const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); - } - // make values in elementstencil unique - std::sort(elementStencil_.begin(), elementStencil_.end()); - elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); - - auto globalI = problem.elementMapper().index(element); - neighborStencil_ = elementStencil_; - - // remove the element itself and possible ghost neighbors from the neighbor stencil - neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), - [globalI](auto i){ return (i == globalI); }), - neighborStencil_.end()); - } - - //! The full element stencil (all element this element is interacting with) - const Stencil& elementStencil() const - { - return elementStencil_; - } - - //! The full element stencil without this element - const Stencil& neighborStencil() const - { - return neighborStencil_; - } - -private: - Stencil elementStencil_; - Stencil neighborStencil_; -}; - -/*! - * \ingroup CCModel - * \brief The global stencil container class for symmetric cc methods - */ -template -class CCSymmetricStencilsVector -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using StencilType = CCSymmetricElementStencils; - -public: - void update(Problem& problem) - { - problemPtr_ = &problem; - elementStencils_.resize(problem.gridView().size(0)); - for (const auto& element : elements(problem.gridView())) - { - auto eIdx = problem.elementMapper().index(element); - elementStencils_[eIdx].update(problem, element); - } - } - - //! overload for elements - template - typename std::enable_if::type - get(const Entity& entity) const - { - return elementStencils_[problemPtr_->elementMapper().index(entity)]; - } - -private: - std::vector elementStencils_; - const Problem* problemPtr_; -}; - -//! Forward declaration of the global stencil class for nonsymmetric cc methods -template class CCNonSymmetricStencilsVector; - -/*! - * \brief Element-related stencils for non-symmetric cc methods - */ -template -class CCNonSymmetricElementStencils -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using IndexType = typename GridView::IndexSet::IndexType; - using Element = typename GridView::template Codim<0>::Entity; - // TODO a separate stencil class all stencils can derive from? - using Stencil = std::vector; -public: - //! Update the stencil. We expect a bound fvGeometry - void update(const Problem& problem, - const Element& element, - CCNonSymmetricStencilsVector& stencilsVector) - { - // the element index of our own here - auto eIdx = problem.elementMapper().index(element); - elementStencil_.clear(); - - // restrict the FvGeometry locally and bind to the element - auto fvGeometry = localView(problem.model().globalFvGeometry()); - fvGeometry.bindElement(element); - - // loop over sub control faces - for (auto&& scvf : scvfs(fvGeometry)) - { - FluxVariables fluxVars; - const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - - // insert stencil into the element stencil - elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); - - // insert our index in the neighbor stencils of the elements in the flux stencil - for (auto globalI : stencil) - { - if (globalI != eIdx) - stencilsVector[globalI].addNeighbor(eIdx); - } - } - - // make values in element stencil unique - std::sort(elementStencil_.begin(), elementStencil_.end()); - elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); - } - - //! The full element stencil (all element this element is interacting with) - const Stencil& elementStencil() const - { - return elementStencil_; - } - - //! The full element stencil without this element - const Stencil& neighborStencil() const - { - return neighborStencil_; - } - - void addNeighbor(const IndexType nIdx) - { - neighborStencil_.push_back(nIdx); - } - - void makeNeighborStencilUnique() - { - // make values in neighbor stencil unique - std::sort(neighborStencil_.begin(), neighborStencil_.end()); - neighborStencil_.erase(std::unique(neighborStencil_.begin(), neighborStencil_.end()), neighborStencil_.end()); - } - -private: - Stencil elementStencil_; - Stencil neighborStencil_; -}; - -/*! - * \ingroup CCModel - * \brief The global stencil container class for non-symmetric cc methods - */ -template -class CCNonSymmetricStencilsVector -{ - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; - using StencilType = CCNonSymmetricElementStencils; - - // The stencil type of an element has to access bracket operator - friend StencilType; - -public: - void update(Problem& problem) - { - problemPtr_ = &problem; - elementStencils_.resize(problem.gridView().size(0)); - for (const auto& element : elements(problem.gridView())) - { - auto eIdx = problem.elementMapper().index(element); - elementStencils_[eIdx].update(problem, element, *this); - } - - for (auto&& stencil : elementStencils_) - stencil.makeNeighborStencilUnique(); - } - - //! overload for elements - template - typename std::enable_if::type - get(const Entity& entity) const - { - return elementStencils_[problemPtr_->elementMapper().index(entity)]; - } - -private: - StencilType& operator[] (const IndexType globalIdx) - { return elementStencils_[globalIdx]; } - - std::vector elementStencils_; - const Problem* problemPtr_; -}; - -} // end namespace - -#endif diff --git a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh index 0636ba0478..343e7d13e7 100644 --- a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh @@ -145,39 +145,35 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - const auto globalI = globalFluxVarsCache().problem_().elementMapper().index(element); - const auto& neighborStencil = globalFluxVarsCache().problem_().model().stencils(element).neighborStencil(); - const auto numNeighbors = neighborStencil.size(); + const auto& problem = globalFluxVarsCache().problem_(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numNeighbors = assemblyMapI.size(); // find the number of scv faces that need to be prepared auto numScvf = fvGeometry.numScvf(); - for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) - { - const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - numScvf += fluxVarIndicesJ.size(); - } + for (unsigned int localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) + numScvf += assemblyMapI[localIdxJ].scvfsJ.size(); // fill the containers with the data on the scv faces inside the actual element fluxVarsCache_.resize(numScvf); globalScvfIndices_.resize(numScvf); - IndexType localScvfIdx = 0; + unsigned int localScvfIdx = 0; for (auto&& scvf : scvfs(fvGeometry)) { - fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); + fluxVarsCache_[localScvfIdx].update(problem, element, fvGeometry, elemVolVars, scvf); globalScvfIndices_[localScvfIdx] = scvf.index(); localScvfIdx++; } // add required data on the scv faces in the neighboring elements - for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) + for (unsigned int localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) { - const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - const auto elementJ = fvGeometry.globalFvGeometry().element(neighborStencil[localIdxJ]); - - for (auto fluxVarIdx : fluxVarIndicesJ) + const auto elementJ = fvGeometry.globalFvGeometry().element(assemblyMapI[localIdxJ].globalJ); + for (auto scvfIdx : assemblyMapI[localIdxJ].scvfsJ) { - auto&& scvfJ = fvGeometry.scvf(fluxVarIdx); - fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), elementJ, fvGeometry, elemVolVars, scvfJ); + auto&& scvfJ = fvGeometry.scvf(scvfIdx); + fluxVarsCache_[localScvfIdx].update(problem, elementJ, fvGeometry, elemVolVars, scvfJ); globalScvfIndices_[localScvfIdx] = scvfJ.index(); localScvfIdx++; } diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index 9080810de0..6cab6287d8 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -117,11 +117,10 @@ public: const FVElementGeometry& fvGeometry, const SolutionVector& sol) { - auto eIdx = globalVolVars().problem_().elementMapper().index(element); - - // stencil information - const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil(); - const auto numDofs = neighborStencil.size() + 1; + const auto& problem = globalVolVars().problem_(); + const auto globalI = problem.elementMapper().index(element); + const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI]; + const auto numDofs = assemblyMapI.size() + 1; // resize local containers to the required size (for internal elements) volumeVariables_.resize(numDofs); @@ -129,21 +128,21 @@ public: int localIdx = 0; // update the volume variables of the element at hand - auto&& scvI = fvGeometry.scv(eIdx); - volumeVariables_[localIdx].update(globalVolVars().problem_().model().elementSolution(element, sol), - globalVolVars().problem_(), + auto&& scvI = fvGeometry.scv(globalI); + volumeVariables_[localIdx].update(problem.model().elementSolution(element, sol), + problem, element, scvI); volVarIndices_[localIdx] = scvI.index(); ++localIdx; // Update the volume variables of the neighboring elements - for (auto globalJ : neighborStencil) + for (const auto& dataJ : assemblyMapI) { - const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ); - auto&& scvJ = fvGeometry.scv(globalJ); - volumeVariables_[localIdx].update(globalVolVars().problem_().model().elementSolution(elementJ, sol), - globalVolVars().problem_(), + const auto& elementJ = fvGeometry.globalFvGeometry().element(dataJ.globalJ); + auto&& scvJ = fvGeometry.scv(dataJ.globalJ); + volumeVariables_[localIdx].update(problem.model().elementSolution(elementJ, sol), + problem, elementJ, scvJ); volVarIndices_[localIdx] = scvJ.index(); @@ -158,15 +157,15 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf); + const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - const ElementSolution dirichletPriVars({globalVolVars().problem_().dirichlet(element, scvf)}); + const ElementSolution dirichletPriVars({problem.dirichlet(element, scvf)}); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); volumeVariables_[localIdx].update(dirichletPriVars, - globalVolVars().problem_(), + problem, element, scvI); volVarIndices_[localIdx] = scvf.outsideScvIdx(); diff --git a/dumux/implicit/assembler.hh b/dumux/implicit/assembler.hh index f54a127bbe..4fc310e2bc 100644 --- a/dumux/implicit/assembler.hh +++ b/dumux/implicit/assembler.hh @@ -176,27 +176,7 @@ private: // Construct the BCRS matrix for the global jacobian void createMatrix_() { - auto numDofs = problem_().model().numDofs(); - - // allocate raw matrix - matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - - // set the row sizes - asImp_().setRowSizes_(); - - // set the indices - asImp_().addIndices_(); - } - - //! Set the row sizes - void setRowSizes_() - { - DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a setRowSizes_() method!"); - } - - void addIndices_() - { - DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a addIndices_() method!"); + DUNE_THROW(Dune::NotImplemented, "Actual implementation does not provide a createMatrix_() method!"); } Implementation &asImp_() diff --git a/dumux/implicit/box/assembler.hh b/dumux/implicit/box/assembler.hh index 33a18bc256..fa34084202 100644 --- a/dumux/implicit/box/assembler.hh +++ b/dumux/implicit/box/assembler.hh @@ -23,6 +23,8 @@ #ifndef DUMUX_BOX_ASSEMBLER_HH #define DUMUX_BOX_ASSEMBLER_HH +#include + #include #include @@ -35,40 +37,37 @@ namespace Dumux { template class BoxAssembler : public ImplicitAssembler { - typedef ImplicitAssembler ParentType; friend class ImplicitAssembler; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - - typedef typename GridView::IndexSet::IndexType IndexType; + using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + static const int dim = GridView::dimension; - void setRowSizes_() + void createMatrix_() { - for (const auto& vertex : vertices(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->vertexMapper_().index(vertex); - const auto& stencil = this->model_().stencils(vertex).vertexStencil(); + const auto numDofs = this->problem_().model().numDofs(); - this->matrix().setrowsize(globalI, stencil.size()); - } - this->matrix().endrowsizes(); + // allocate raw matrix + this->matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - } + // get occupation pattern of the matrix + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numDofs, numDofs); - void addIndices_() - { - for (const auto& vertex : vertices(this->gridView_())) + for (const auto& element : elements(this->problem_().gridView())) { - // the global index of the element at hand - const auto globalI = this->vertexMapper_().index(vertex); - const auto& stencil = this->model_().stencils(vertex).vertexStencil(); - - for (auto&& globalJ : stencil) - this->matrix().addindex(globalI, globalJ); + for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) + { + const auto globalI = this->problem_().vertexMapper().subIndex(element, vIdx, dim); + for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2) + { + const auto globalJ = this->problem_().vertexMapper().subIndex(element, vIdx2, dim); + occupationPattern.add(globalI, globalJ); + occupationPattern.add(globalJ, globalI); + } + } } - this->matrix().endindices(); + // export pattern to matrix + occupationPattern.exportIdx(*this->matrix_); } }; diff --git a/dumux/implicit/box/propertydefaults.hh b/dumux/implicit/box/propertydefaults.hh index bc26b2565a..c95a65d75b 100644 --- a/dumux/implicit/box/propertydefaults.hh +++ b/dumux/implicit/box/propertydefaults.hh @@ -35,7 +35,6 @@ #include #include #include -#include #include #include @@ -53,7 +52,6 @@ namespace Dumux { // forward declarations template class BoxLocalResidual; template class BoxElementBoundaryTypes; -template class BoxStencilsVector; namespace Properties { //! Set the corresponding discretization method property @@ -101,9 +99,6 @@ SET_TYPE_PROP(BoxModel, ElementBoundaryTypes, BoxElementBoundaryTypes); //! Mapper for the degrees of freedoms. SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper)); -//! The stencil container -SET_TYPE_PROP(BoxModel, StencilsVector, BoxStencilsVector); - //! The global volume variables vector class SET_TYPE_PROP(BoxModel, GlobalVolumeVariables, BoxGlobalVolumeVariables); diff --git a/dumux/implicit/cellcentered/assembler.hh b/dumux/implicit/cellcentered/assembler.hh index a2430d7cbd..d376a026a5 100644 --- a/dumux/implicit/cellcentered/assembler.hh +++ b/dumux/implicit/cellcentered/assembler.hh @@ -23,6 +23,8 @@ #ifndef DUMUX_CC_ASSEMBLER_HH #define DUMUX_CC_ASSEMBLER_HH +#include + #include #include @@ -35,40 +37,30 @@ namespace Dumux { template class CCAssembler : public ImplicitAssembler { - typedef ImplicitAssembler ParentType; friend class ImplicitAssembler; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::IndexSet::IndexType IndexType; + using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix); - void setRowSizes_() + void createMatrix_() { - for (const auto& element : elements(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->elementMapper_().index(element); - const auto& stencil = this->model_().stencils(element).elementStencil(); + const auto numDofs = this->problem_().model().numDofs(); - this->matrix().setrowsize(globalI, stencil.size()); - } - this->matrix().endrowsizes(); - } + // allocate raw matrix + this->matrix_ = std::make_shared(numDofs, numDofs, JacobianMatrix::random); - void addIndices_() - { - for (const auto& element : elements(this->gridView_())) - { - // the global index of the element at hand - const auto globalI = this->elementMapper_().index(element); - const auto& stencil = this->model_().stencils(element).elementStencil(); + // get occupation pattern of the matrix + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numDofs, numDofs); - for (auto&& globalJ : stencil) - this->matrix().addindex(globalI, globalJ); + const auto& assemblyMap = this->problem_().model().localJacobian().assemblyMap(); + for (unsigned int globalI = 0; globalI < numDofs; ++globalI) + { + occupationPattern.add(globalI, globalI); + for (const auto& dataJ : assemblyMap[globalI]) + occupationPattern.add(dataJ.globalJ, globalI); } - this->matrix().endindices(); + + // export pattern to matrix + occupationPattern.exportIdx(*this->matrix_); } }; diff --git a/dumux/implicit/cellcentered/assemblymap.hh b/dumux/implicit/cellcentered/assemblymap.hh index 4bce9c22cb..6e8e1d78ca 100644 --- a/dumux/implicit/cellcentered/assemblymap.hh +++ b/dumux/implicit/cellcentered/assemblymap.hh @@ -19,26 +19,35 @@ /*! * \file * \brief Stores the face indices corresponding to the neighbors of an element - * that contribute to the derivative calculation + * that contribute to the derivative calculation. This is used for + * finite-volume schemes with symmetric sparsity pattern in the global matrix. */ -#ifndef DUMUX_CC_ASSEMBLY_MAP_HH -#define DUMUX_CC_ASSEMBLY_MAP_HH +#ifndef DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH +#define DUMUX_CC_SYMMETRIC_ASSEMBLY_MAP_HH -#include +#include +#include namespace Dumux { template -class CCAssemblyMap +class CCSymmetricAssemblyMap { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); using IndexType = typename GridView::IndexSet::IndexType; - using Map = std::vector< std::vector< std::vector > >; + + struct DataJ + { + IndexType globalJ; + std::vector scvfsJ; + }; + + using Map = std::vector>; public: @@ -52,58 +61,51 @@ public: map_.resize(problem.gridView().size(0)); for (const auto& element : elements(problem.gridView())) { - // get a local finite volume geometry object that is bindable - auto fvGeometryJ = localView(problem.model().globalFvGeometry()); + // We are looking for the elements I, for which this element J is in the flux stencil + auto globalJ = problem.elementMapper().index(element); - auto globalI = problem.elementMapper().index(element); - const auto& neighborStencil = problem.model().stencils(element).neighborStencil(); + auto fvGeometry = localView(problem.model().globalFvGeometry()); + fvGeometry.bindElement(element); - map_[globalI].reserve(neighborStencil.size()); - for (auto globalJ : neighborStencil) - { - const auto& elementJ = fvGeometryJ.globalFvGeometry().element(globalJ); + // obtain the data of J in elements I + std::vector> dataJForI; - // find the flux vars needed for the calculation of the flux into element - std::vector fluxVarIndices; + // loop over sub control faces + for (auto&& scvf : scvfs(fvGeometry)) + { + FluxVariables fluxVars; + const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf); - // only non-ghost neighbors (J) have to be considered, derivatives from non-ghost to ghost dofs - // are assembled when assembling the ghost element (I) - if (elementJ.partitionType() != Dune::GhostEntity) + // insert our index in the neighbor stencils of the elements in the flux stencil + for (auto globalI : stencil) { - fvGeometryJ.bindElement(elementJ); - for (auto&& scvFaceJ : scvfs(fvGeometryJ)) - { - auto fluxVarsIdx = scvFaceJ.index(); - - // if globalI is in flux var stencil, add to list - FluxVariables fluxVars; - const auto fluxStencil = fluxVars.computeStencil(problem, elementJ, fvGeometryJ, scvFaceJ); - - for (auto globalIdx : fluxStencil) - { - if (globalIdx == globalI) - { - fluxVarIndices.push_back(fluxVarsIdx); - break; - } - } - } + 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()) + it->second.scvfsJ.push_back(scvf.index()); + else + dataJForI.emplace_back(std::make_pair(globalI, + DataJ({globalJ, std::vector({scvf.index()})}))); } - map_[globalI].emplace_back(std::move(fluxVarIndices)); } + + for (auto&& pair : dataJForI) + map_[pair.first].emplace_back(std::move(pair.second)); } } - //! Some implementations of the assembly map can be solution-dependent and need to be updated - void update() {} - - const std::vector>& operator [] (const IndexType globalI) const + const std::vector& operator[] (const IndexType globalI) const { return map_[globalI]; } private: Map map_; }; -} +} // end namespace Dumux #endif diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 6ac50bf3f9..7a5fafd63b 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -211,8 +211,7 @@ private: const bool isGhost) { // get stencil informations - const auto& neighborStencil = this->model_().stencils(element).neighborStencil(); - const auto numNeighbors = neighborStencil.size(); + const auto numNeighbors = assemblyMap_[globalI_].size(); // the localresidual class used for the flux calculations LocalResidual localRes; @@ -227,14 +226,14 @@ private: origFlux = 0.0; unsigned int j = 0; - for (auto globalJ : neighborStencil) + for (const auto& dataJ : assemblyMap_[globalI_]) { - auto elementJ = fvGeometry.globalFvGeometry().element(globalJ); + auto elementJ = fvGeometry.globalFvGeometry().element(dataJ.globalJ); neighborElements.push_back(elementJ); - for (auto fluxVarIdx : assemblyMap_[globalI_][j]) + for (auto scvfIdx : dataJ.scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, curElemVolVars, @@ -295,9 +294,9 @@ private: // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) { - for (auto fluxVarIdx : assemblyMap_[globalI_][k]) + for (auto scvfIdx : assemblyMap_[globalI_][k].scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, @@ -346,9 +345,9 @@ private: // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) { - for (auto fluxVarIdx : assemblyMap_[globalI_][k]) + for (auto scvfIdx : assemblyMap_[globalI_][k].scvfsJ) { - auto&& scvf = fvGeometry.scvf(fluxVarIdx); + auto&& scvf = fvGeometry.scvf(scvfIdx); neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, @@ -383,8 +382,8 @@ private: this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv); j = 0; - for (auto globalJ : neighborStencil) - this->updateGlobalJacobian_(matrix, globalJ, globalI_, pvIdx, neighborDeriv[j++]); + for (const auto& dataJ : assemblyMap_[globalI_]) + this->updateGlobalJacobian_(matrix, dataJ.globalJ, globalI_, pvIdx, neighborDeriv[j++]); } } diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh index 9ff21ee167..3278dafaef 100644 --- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh @@ -39,7 +39,6 @@ #include #include #include -#include #include #include #include @@ -93,9 +92,6 @@ SET_TYPE_PROP(CCMpfaModel, ElementVolumeVariables, Dumux::CCMpfaElementVolumeVar //! The local flux variables cache vector class SET_TYPE_PROP(CCMpfaModel, ElementFluxVariablesCache, Dumux::CCMpfaElementFluxVariablesCache); -//! The stencil container, some mpfa methods might lead to non-symmetric global systems -SET_TYPE_PROP(CCMpfaModel, StencilsVector, Dumux::CCMpfaStencilsVector); - SET_PROP(CCMpfaModel, SubControlVolumeFace) { private: diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh index 03c46d81af..20cadb848c 100644 --- a/dumux/implicit/cellcentered/propertydefaults.hh +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -76,8 +76,8 @@ SET_TYPE_PROP(CCModel, GlobalVolumeVariables, CCGlobalVolumeVariables); -//! Set the AssemblyMap to the CCAssemblyMap -SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCAssemblyMap); +//! Set the AssemblyMap to the CCAssemblyMap (default: symmetric occupation pattern) +SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSymmetricAssemblyMap); //! indicate that this is no box discretization SET_BOOL_PROP(CCModel, ImplicitIsBox, false); diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index f029f862d6..1422045377 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -35,7 +35,6 @@ #include #include #include -#include #include #include @@ -66,9 +65,6 @@ SET_TYPE_PROP(CCTpfaModel, ElementVolumeVariables, Dumux::CCTpfaElementVolumeVar //! The local flux variables cache vector class SET_TYPE_PROP(CCTpfaModel, ElementFluxVariablesCache, Dumux::CCTpfaElementFluxVariablesCache); -//! The stencil container, tpfa leads to a symmetric global matrix -SET_TYPE_PROP(CCTpfaModel, StencilsVector, Dumux::CCSymmetricStencilsVector); - SET_PROP(CCTpfaModel, SubControlVolumeFace) { private: diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 5d846855c3..a067d55a69 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -71,7 +71,6 @@ class ImplicitModel using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using GlobalFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GlobalFluxVariablesCache); - using StencilsVector = typename GET_PROP_TYPE(TypeTag, StencilsVector); enum { numEq = GET_PROP_VALUE(TypeTag, NumEq), @@ -129,9 +128,6 @@ public: // resize and update the volVars with the initial solution curGlobalVolVars_.update(problem, curSol()); - // update stencils - stencilsVector_.update(problem); - // initialize assembler and create matrix localJacobian_.init(problem); jacAsm_ = std::make_shared(); @@ -845,23 +841,12 @@ protected: // the flux variables cache vector vector GlobalFluxVariablesCache globalfluxVarsCache_; - // the stencils vector - StencilsVector stencilsVector_; - // the finite volume element geometries std::shared_ptr globalFvGeometryPtr_; // container to store the box volumes Dune::BlockVector > boxVolume_; -public: - //! Get stencils related to an entity of specified codimension - // Cell-centered discretizations typically only implement element stencils - // Vertex-centered discretizations both vertex and element stencils - template - auto stencils(const Entity& entity) const -> decltype(stencilsVector_.get(entity)) - { return stencilsVector_.get(entity); } - private: /*! * \brief Returns whether messages should be printed diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index d081d15030..44b5c023c8 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -103,9 +103,6 @@ NEW_PROP_TAG(SolutionDependentHeatConduction); //!< specifies if the parameters NEW_PROP_TAG(VtkAddVelocity); //!< specifies if an element velocity it reconstructed for the output NEW_PROP_TAG(VtkAddProcessRank); //!< specifies if the process rank should be added the output -// stencils -NEW_PROP_TAG(StencilsVector); //! The type of the global vector of stencils per element - // high level simulation control NEW_PROP_TAG(TimeManager); //!< Manages the simulation time NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh index b9539a510a..3e75f61866 100644 --- a/dumux/porousmediumflow/implicit/velocityoutput.hh +++ b/dumux/porousmediumflow/implicit/velocityoutput.hh @@ -65,7 +65,6 @@ class ImplicitVelocityOutput using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; using CoordScalar = typename GridView::ctype; - using Stencil = std::vector; using GlobalPosition = Dune::FieldVector; using ReferenceElements = Dune::ReferenceElements; @@ -89,8 +88,9 @@ public: // resize to the number of vertices of the grid cellNum_.assign(problem.gridView().size(dim), 0); - for (const auto& vertex : vertices(problem.gridView())) - cellNum_[problem.vertexMapper().index(vertex)] = getStencil(vertex).size(); + for (const auto& element : elements(problem.gridView())) + for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) + ++cellNum_[problem.vertexMapper().subIndex(element, vIdx, dim)]; } } } @@ -111,8 +111,7 @@ public: // compiler error. However, that call is needed for calculating velocities // if the cell-centered discretization is used. By proceeding as in the // following lines, that call will only be compiled if cell-centered - // actually is used. For the same reason we also provide a isBox-specific - // implementation of the getStencil method below. + // actually is used. template typename std::enable_if::type problemBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const @@ -124,18 +123,6 @@ public: problemBoundaryTypes(const Element& element, const SubControlVolume& scv) const { return BoundaryTypes(); } - //! returns the elements connected to a vertex - template - const typename std::enable_if::type& - getStencil(const Vertex& vertex) const - { return problem_.model().stencils(vertex).elementIndices(); } - - //! we should never call this method for cc models - template - const typename std::enable_if::type - getStencil(const Vertex& vertex) const - { return Stencil(0); } - //! Calculate the velocities for the scvs in the element //! We assume the local containers to be bound to the complete stencil template -- GitLab