From 2a6840cc4726417f5e7dcd6abbf19e2e01ac53a1 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de> Date: Fri, 27 Jul 2018 11:54:24 +0200 Subject: [PATCH] [staggered] Fix caching of fluxVariablesCache * empty cache filler used * contains several todos * fluxVarsCache not used anywehere, but might be required in future --- .../staggered/elementfluxvariablescache.hh | 112 ++++++++++++------ .../staggered/gridfluxvariablescache.hh | 84 ++++++++----- dumux/discretization/staggered/properties.hh | 4 +- 3 files changed, 136 insertions(+), 64 deletions(-) diff --git a/dumux/discretization/staggered/elementfluxvariablescache.hh b/dumux/discretization/staggered/elementfluxvariablescache.hh index 7be5655c75..18e85b92b8 100644 --- a/dumux/discretization/staggered/elementfluxvariablescache.hh +++ b/dumux/discretization/staggered/elementfluxvariablescache.hh @@ -53,6 +53,9 @@ public: //! export the type of the flux variables cache using FluxVariablesCache = typename GFVC::FluxVariablesCache; + //! export the type of the flux variables cache filler + using FluxVariablesCacheFiller = typename GFVC::FluxVariablesCacheFiller; + StaggeredElementFluxVariablesCache(const GridFluxVariablesCache& global) : gridFluxVarsCachePtr_(&global) {} @@ -75,6 +78,15 @@ public: const ElementVolumeVariables& elemVolVars, const typename FVElementGeometry::SubControlVolumeFace& scvf) {} + //! Specialization for the global caching being enabled - do nothing here + template<class FVElementGeometry, class ElementVolumeVariables> + void update(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) + { + DUNE_THROW(Dune::InvalidStateException, "In case of enabled caching, the grid flux variables cache has to be updated"); + } + //! operators in the case of caching template<class SubControlVolumeFace> const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const @@ -103,11 +115,17 @@ public: //! export the type of the flux variables cache using FluxVariablesCache = typename GFVC::FluxVariablesCache; + //! export the type of the flux variables cache filler + using FluxVariablesCacheFiller = typename GFVC::FluxVariablesCacheFiller; + StaggeredElementFluxVariablesCache(const GridFluxVariablesCache& global) : gridFluxVarsCachePtr_(&global) {} - //! This function has to be called prior to flux calculations on the element. - //! Prepares the transmissibilities of the scv faces in an element. The FvGeometry is assumed to be bound. + /*! + * \brief Prepares the transmissibilities of the scv faces in an element + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on the element. + */ template<class FVElementGeometry, class ElementVolumeVariables> void bindElement(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, @@ -117,63 +135,55 @@ public: const auto numScvf = fvGeometry.numScvf(); fluxVarsCache_.resize(numScvf); globalScvfIndices_.resize(numScvf); - std::size_t localScvfIdx = 0; + // instantiate helper class to fill the caches + // FluxVariablesCacheFiller filler(gridFluxVarsCache().problem()); TODO: use proper ctor + FluxVariablesCacheFiller filler; + + std::size_t localScvfIdx = 0; // fill the containers for (auto&& scvf : scvfs(fvGeometry)) { - fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); + filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf, true); globalScvfIndices_[localScvfIdx] = scvf.index(); localScvfIdx++; } } - //! This function is called by the StaggeredLocalResidual before flux calculations during assembly. - //! Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound. + /*! + * \brief Prepares the transmissibilities of the scv faces in the stencil of an element + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on the element. + */ template<class FVElementGeometry, class ElementVolumeVariables> void bind(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - const auto globalI = gridFluxVarsCache().problem_().elementMapper().index(element); - const auto& neighborStencil = gridFluxVarsCache().problem_().model().stencils(element).neighborStencil(); - const auto numNeighbors = neighborStencil.size(); + // instantiate helper class to fill the caches + // FluxVariablesCacheFiller filler(problem); TODO: use proper ctor + FluxVariablesCacheFiller filler; // find the number of scv faces that need to be prepared - auto numScvf = fvGeometry.numScvf(); - for (std::size_t localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) - { - const auto& fluxVarIndicesJ = gridFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - numScvf += fluxVarIndicesJ.size(); - } + const auto numScvf = fvGeometry.numScvf(); // fill the containers with the data on the scv faces inside the actual element fluxVarsCache_.resize(numScvf); globalScvfIndices_.resize(numScvf); - std::size_t localScvfIdx = 0; + unsigned int localScvfIdx = 0; for (auto&& scvf : scvfs(fvGeometry)) { - fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); + filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf, true); globalScvfIndices_[localScvfIdx] = scvf.index(); localScvfIdx++; } - - // add required data on the scv faces in the neighboring elements - for (std::size_t localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ) - { - const auto& fluxVarIndicesJ = gridFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ]; - const auto elementJ = fvGeometry.fvGridGeometry().element(neighborStencil[localIdxJ]); - - for (auto fluxVarIdx : fluxVarIndicesJ) - { - auto&& scvfJ = fvGeometry.scvf(fluxVarIdx); - fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), elementJ, fvGeometry, elemVolVars, scvfJ); - globalScvfIndices_[localScvfIdx] = scvfJ.index(); - localScvfIdx++; - } - } } + /*! + * \brief Prepares the transmissibilities of a single scv face + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on the element. + */ template<class FVElementGeometry, class ElementVolumeVariables> void bindScvf(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, @@ -183,10 +193,46 @@ public: fluxVarsCache_.resize(1); globalScvfIndices_.resize(1); - fluxVarsCache_[0].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); + // instantiate helper class to fill the caches + // FluxVariablesCacheFiller filler(gridFluxVarsCache().problem()); + FluxVariablesCacheFiller filler; // TODO: use proper ctor + + filler.fill(*this, fluxVarsCache_[0], element, fvGeometry, elemVolVars, scvf, true); globalScvfIndices_[0] = scvf.index(); } + /*! + * \brief Update the transmissibilities if the volume variables have changed + * \note Results in undefined behaviour if called before bind() or with a different element + */ + template<class FVElementGeometry, class ElementVolumeVariables> + void update(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) + { + // if (FluxVariablesCacheFiller::isSolDependent) TODO + // { + // const auto& problem = gridFluxVarsCache().problem(); + // const auto globalI = fvGeometry.fvGridGeometry().elementMapper().index(element); + // + // // instantiate filler class + // FluxVariablesCacheFiller filler(problem); + // + // // let the filler class update the caches + // for (unsigned int localScvfIdx = 0; localScvfIdx < fluxVarsCache_.size(); ++localScvfIdx) + // { + // const auto& scvf = fvGeometry.scvf(globalScvfIndices_[localScvfIdx]); + // + // const auto scvfInsideScvIdx = scvf.insideScvIdx(); + // const auto& insideElement = scvfInsideScvIdx == globalI ? + // element : + // fvGeometry.fvGridGeometry().element(scvfInsideScvIdx); + // + // filler.fill(*this, fluxVarsCache_[localScvfIdx], insideElement, fvGeometry, elemVolVars, scvf); + // } + // } + } + //! access operators in the case of no caching template<class SubControlVolumeFace> const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const diff --git a/dumux/discretization/staggered/gridfluxvariablescache.hh b/dumux/discretization/staggered/gridfluxvariablescache.hh index 4d02b339f1..17af99801a 100644 --- a/dumux/discretization/staggered/gridfluxvariablescache.hh +++ b/dumux/discretization/staggered/gridfluxvariablescache.hh @@ -36,11 +36,12 @@ namespace Dumux { * \tparam P The problem type * \tparam FVC The flux variables cache type */ -template<class P, class FVC> +template<class P, class FVC, class FVCF> struct StaggeredDefaultGridFluxVariablesCacheTraits { using Problem = P; using FluxVariablesCache = FVC; + using FluxVariablesCacheFiller = FVCF; template<class GridFluxVariablesCache, bool cachingEnabled> using LocalView = StaggeredElementFluxVariablesCache<GridFluxVariablesCache, cachingEnabled>; @@ -52,8 +53,9 @@ struct StaggeredDefaultGridFluxVariablesCacheTraits */ template<class Problem, class FluxVariablesCache, - bool cachingEnabled = false, - class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache> > + class FluxVariablesCacheFiller, + bool EnableGridFluxVariablesCache = false, + class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache, FluxVariablesCacheFiller> > class StaggeredGridFluxVariablesCache; /*! @@ -61,16 +63,19 @@ class StaggeredGridFluxVariablesCache; * \brief Flux variables cache class for staggered models. Specialization in case of storing the flux cache. */ -template<class P, class FVC, class Traits> -class StaggeredGridFluxVariablesCache<P, FVC, true, Traits> +template<class P, class FVC, class FVCF, class Traits> +class StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits> { - using ThisType = StaggeredGridFluxVariablesCache<P, FVC, true, Traits>; using Problem = typename Traits::Problem; + using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>; public: //! export the flux variable cache type using FluxVariablesCache = typename Traits::FluxVariablesCache; + //! export the flux variable cache filler type + using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller; + //! make it possible to query if caching is enabled static constexpr bool cachingEnabled = true; @@ -86,21 +91,31 @@ public: const SolutionVector& sol, bool forceUpdate = false) { - // fluxVarsCache_.resize(fvGridGeometry.numScvf()); - // for (const auto& element : elements(fvGridGeometry.gridView())) - // { - // // Prepare the geometries within the elements of the stencil - // auto fvGeometry = localView(fvGridGeometry); - // fvGeometry.bind(element); - // - // auto elemVolVars = localView(gridVolVars); - // elemVolVars.bind(element, fvGeometry, sol); - // - // for (auto&& scvf : scvfs(fvGeometry)) - // { - // fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf); - // } - // } + // only do the update if fluxes are solution dependent or if update is forced + // TODO: so far, the staggered models do not use any fluxVar caches, therefore an empty cache filler + // is used which does not implement isSolDependent + if (/*FluxVariablesCacheFiller::isSolDependent ||*/ forceUpdate) + { + // instantiate helper class to fill the caches + // FluxVariablesCacheFiller filler(problem()); TODO: use proper ctor + FluxVariablesCacheFiller filler; + + fluxVarsCache_.resize(fvGridGeometry.numScvf()); + for (const auto& element : elements(fvGridGeometry.gridView())) + { + // Prepare the geometries within the elements of the stencil + auto fvGeometry = localView(fvGridGeometry); + fvGeometry.bind(element); + + auto elemVolVars = localView(gridVolVars); + elemVolVars.bind(element, fvGeometry, sol); + + for (auto&& scvf : scvfs(fvGeometry)) + { + filler.fill(*this, fluxVarsCache_[scvf.index()], element, fvGeometry, elemVolVars, scvf, forceUpdate); + } + } + } } const Problem& problem() const @@ -125,31 +140,40 @@ private: * \brief Flux variables cache class for staggered models. Specialization in case of not storing the flux cache. */ -template<class P, class FVC, class Traits> -class StaggeredGridFluxVariablesCache<P, FVC, false, Traits> +template<class P, class FVC, class FVCF, class Traits> +class StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits> { - using ThisType = StaggeredGridFluxVariablesCache<P, FVC, false, Traits>; using Problem = typename Traits::Problem; + using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>; public: //! export the flux variable cache type using FluxVariablesCache = typename Traits::FluxVariablesCache; + //! export the flux variable cache filler type + using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller; + //! make it possible to query if caching is enabled - static constexpr bool cachingEnabled = true; + static constexpr bool cachingEnabled = false; //! export the type of the local view using LocalView = typename Traits::template LocalView<ThisType, cachingEnabled>; - // When global flux variables caching is disabled, we don't need to update the cache - void update(Problem& problem) - { problemPtr_ = &problem; } + StaggeredGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {} -private: + // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces + template<class FVGridGeometry, class GridVolumeVariables, class SolutionVector> + void update(const FVGridGeometry& fvGridGeometry, + const GridVolumeVariables& gridVolVars, + const SolutionVector& sol, + bool forceUpdate = false) {} - const Problem& problem_() const + const Problem& problem() const { return *problemPtr_; } +private: + + const Problem* problemPtr_; }; diff --git a/dumux/discretization/staggered/properties.hh b/dumux/discretization/staggered/properties.hh index 51f0661255..1671641215 100644 --- a/dumux/discretization/staggered/properties.hh +++ b/dumux/discretization/staggered/properties.hh @@ -32,6 +32,7 @@ #include <dumux/discretization/methods.hh> #include <dumux/discretization/fvproperties.hh> +#include <dumux/discretization/fluxvariablescaching.hh> #include <dumux/discretization/cellcentered/elementboundarytypes.hh> #include <dumux/assembly/staggeredlocalresidual.hh> @@ -78,8 +79,9 @@ private: static constexpr auto enableCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); + using FluxVariablesCacheFiller = FluxVariablesCaching::EmptyCacheFiller; public: - using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, enableCache>; + using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, FluxVariablesCacheFiller, enableCache>; }; //! Set the face solution type -- GitLab