From 61479b08be2d64bd19e388f01459554979bbdf95 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Tue, 5 Jul 2016 11:32:16 +0200
Subject: [PATCH] [implicit][fluxVarCacheVec] specialize flux var cache
The flux var cache is now specialized for the box and cc methods
and furthermore permits the binding of an element prior to flux calculations
in order to recycle data later in the flux calculations when the global
caching is disabled.
---
.../cellcentered/fluxvariablescachevector.hh | 194 ++++++++++++++++++
dumux/implicit/cellcentered/localjacobian.hh | 7 +-
dumux/implicit/fluxvariablescachevector.hh | 81 --------
.../constitutivelaws/darcyslaw.hh | 14 +-
4 files changed, 200 insertions(+), 96 deletions(-)
create mode 100644 dumux/implicit/cellcentered/fluxvariablescachevector.hh
delete mode 100644 dumux/implicit/fluxvariablescachevector.hh
diff --git a/dumux/implicit/cellcentered/fluxvariablescachevector.hh b/dumux/implicit/cellcentered/fluxvariablescachevector.hh
new file mode 100644
index 0000000000..ffae3dbf96
--- /dev/null
+++ b/dumux/implicit/cellcentered/fluxvariablescachevector.hh
@@ -0,0 +1,194 @@
+// -*- 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 Base class for the volume variables vector
+ */
+#ifndef DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH
+#define DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH
+
+#include <dumux/implicit/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the flux variables cache vector, we store one cache per face
+ */
+template<class TypeTag>
+class CCFluxVariablesCacheVector
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+public:
+ // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
+ template <typename T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ {
+ problemPtr_ = &problem;
+
+ fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
+ for (const auto& element : elements(problem.gridView()))
+ {
+ // Prepare the geometries within the elements of the stencil
+ problem.model().fvGeometries_().bind(element);
+
+ const auto& fvGeometry = problem.model().fvGeometries(element);
+ for (const auto& scvf : fvGeometry.scvfs())
+ {
+ (*this)[scvf].update(problem, element, scvf);
+ }
+ }
+ }
+
+ // When global flux variables caching is disabled, we don't need to update the cache
+ template <typename T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ {
+ problemPtr_ = &problem;
+ }
+
+ // 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.
+ template <typename T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bindElement(const Element& element)
+ {
+ const auto& fvGeometry = problem_().model().fvGeometries(element);
+
+ // resizing of the cache
+ const auto numScvf = fvGeometry.numScvf();
+ fluxVarsCache_.resize(numScvf);
+ globalScvfIndices_.resize(numScvf);
+ IndexType localScvfIdx = 0;
+
+ // fill the containers
+ for (const auto& scvf : fvGeometry.scvfs())
+ {
+ fluxVarsCache_[localScvfIdx].update(problem_(), element, scvf);
+ globalScvfIndices_[localScvfIdx] = scvf.index();
+ localScvfIdx++;
+ }
+ }
+
+ // Specialization for the global caching being enabled - do nothing here
+ template <typename T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bindElement(const Element& element) {}
+
+ // This function is called by the CCLocalResidual before flux calculations during assembly.
+ // Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound.
+ template <typename T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bind(const Element& element)
+ {
+ const auto globalI = problem_().elementMapper().index(element);
+ const auto& fvGeometry = problem_().model().fvGeometries(element);
+ const auto& neighborStencil = problem_().model().stencils(element).neighborStencil();
+ const auto numNeighbors = neighborStencil.size();
+
+ // find the number of scv faces that need to be prepared
+ IndexType numScvf = fvGeometry.numScvf();
+ for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
+ {
+ const auto& fluxVarIndicesJ = problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
+ numScvf += fluxVarIndicesJ.size();
+ }
+
+ // fill the containers with the data on the scv faces inside the actual element
+ fluxVarsCache_.resize(numScvf);
+ globalScvfIndices_.resize(numScvf);
+ IndexType localScvfIdx = 0;
+ for (const auto& scvf : fvGeometry.scvfs())
+ {
+ fluxVarsCache_[localScvfIdx].update(problem_(), element, scvf);
+ globalScvfIndices_[localScvfIdx] = scvf.index();
+ localScvfIdx++;
+ }
+
+ // add required data on the scv faces in the neighboring elements
+ for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
+ {
+ const auto& fluxVarIndicesJ = problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
+ const auto elementJ = problem_().model().fvGeometries().element(neighborStencil[localIdxJ]);
+ for (auto fluxVarIdx : fluxVarIndicesJ)
+ {
+ const auto& scvfJ = problem_().model().fvGeometries().subControlVolumeFace(fluxVarIdx);
+ fluxVarsCache_[localScvfIdx].update(problem_(), elementJ, scvfJ);
+ globalScvfIndices_[localScvfIdx] = scvfJ.index();
+ localScvfIdx++;
+ }
+ }
+ }
+
+ // Specialization for the global caching being enabled - do nothing here
+ template <typename T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bind(const Element& element) {}
+
+ // access operators in the case of caching
+ template <typename T = TypeTag>
+ const typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
+ operator [](const SubControlVolumeFace& scvf) const
+ { return fluxVarsCache_[scvf.index()]; }
+
+ template <typename T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
+ operator [](const SubControlVolumeFace& scvf)
+ { return fluxVarsCache_[scvf.index()]; }
+
+ // access operators in the case of no caching
+ template <typename T = TypeTag>
+ const typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
+ operator [](const SubControlVolumeFace& scvf) const
+ { return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; }
+
+ template <typename T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
+ operator [](const SubControlVolumeFace& scvf)
+ { return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; }
+
+private:
+ // get index of scvf in the local container
+ template <typename T = TypeTag>
+ const typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), int>::type
+ getLocalScvfIdx_(const int scvfIdx) const
+ {
+ auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
+
+ if (it != globalScvfIndices_.end())
+ return std::distance(globalScvfIndices_.begin(), it);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Could not find the flux vars cache for scvfIdx = " << scvfIdx);
+ }
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+
+ std::vector<FluxVariablesCache> fluxVarsCache_;
+ std::vector<IndexType> globalScvfIndices_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh
index 9767424060..961f828d10 100644
--- a/dumux/implicit/cellcentered/localjacobian.hh
+++ b/dumux/implicit/cellcentered/localjacobian.hh
@@ -159,6 +159,7 @@ public:
this->model_().fvGeometries_().bind(element);
this->model_().curVolVars_().bind(element);
this->model_().prevVolVars_().bindElement(element);
+ this->model_().fluxVariablesCache_().bind(element);
// set the current grid element and update the element's
// finite volume geometry
@@ -259,8 +260,9 @@ private:
priVars[pvIdx] += eps;
delta += eps;
- // update the volume variables
+ // update the volume variables and bind the flux var cache again
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
+ this->model_().fluxVariablesCache_().bind(element);
if (!isGhost)
{
@@ -297,8 +299,9 @@ private:
priVars[pvIdx] -= delta + eps;
delta += eps;
- // update the volume variables
+ // update the volume variables and bind the flux var cache again
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
+ this->model_().fluxVariablesCache_().bind(element);
if (!isGhost)
{
diff --git a/dumux/implicit/fluxvariablescachevector.hh b/dumux/implicit/fluxvariablescachevector.hh
deleted file mode 100644
index e426ac6dd2..0000000000
--- a/dumux/implicit/fluxvariablescachevector.hh
+++ /dev/null
@@ -1,81 +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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Base class for the volume variables vector
- */
-#ifndef DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH
-#define DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH
-
-#include <dumux/implicit/properties.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ImplicitModel
- * \brief Base class for the flux variables cache vector, we store one cache per face
- */
-template<class TypeTag>
-class FluxVariablesCacheVector
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-
-public:
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
- {
- fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
- for (const auto& element : elements(problem.gridView()))
- {
- // bind the geometries and volume variables to the element (all the elements in stencil)
- problem.model().fvGeometries_().bind(element);
- problem.model().curVolVars_().bind(element);
-
- const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scvf : fvGeometry.scvfs())
- {
- (*this)[scvf.index()].update(problem, element, scvf);
- }
- }
- }
-
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
- {}
-
- template <typename T = TypeTag>
- const typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), FluxVariablesCache>::type& operator [](IndexType scvfIdx) const
- { return fluxVarsCache_[scvfIdx]; }
-
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), FluxVariablesCache>::type& operator [](IndexType scvfIdx)
- { return fluxVarsCache_[scvfIdx]; }
-
-private:
- std::vector<FluxVariablesCache> fluxVarsCache_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh b/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh
index bdd907a423..91cb2e6869 100644
--- a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh
+++ b/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh
@@ -82,7 +82,7 @@ public:
const SubControlVolumeFace& scvFace,
const IndexType phaseIdx)
{
- const auto& tij = getTransmissibilities(problem, scvFace);
+ const auto& tij = problem.model().fluxVarsCache(scvFace).tij();
// Get the inside volume variables
const auto insideScvIdx = scvFace.insideScvIdx();
@@ -143,18 +143,6 @@ public:
return stencil;
}
- template <typename T = TypeTag>
- static const typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), Scalar>::type& getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- return problem.model().fluxVarsCache(scvFace).tij();
- }
-
- template <typename T = TypeTag>
- static const typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache), Scalar>::type getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- return calculateTransmissibilities(problem, scvFace);
- }
-
static Scalar calculateTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
{
Scalar tij;
--
GitLab