diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index c8210e26b71d2cb8abd48c67ee76b5e36395c934..bd2ebb628d681a0c1f904bbd7bd1e3b195791f25 100644
--- a/dumux/assembly/CMakeLists.txt
+++ b/dumux/assembly/CMakeLists.txt
@@ -5,10 +5,12 @@ boxlocalresidual.hh
cclocalassembler.hh
cclocalresidual.hh
diffmethod.hh
+entitycolor.hh
fvassembler.hh
fvlocalassemblerbase.hh
fvlocalresidual.hh
jacobianpattern.hh
+partialreassembler.hh
staggeredfvassembler.hh
staggeredlocalassembler.hh
staggeredlocalresidual.hh
diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index 104e2189294622b223061866396c55485bc4f534..2601c5d90a88216d347584828a9e162612ec562e 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -34,6 +34,8 @@
#include <dumux/assembly/numericepsilon.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/assembly/fvlocalassemblerbase.hh>
+#include <dumux/assembly/partialreassembler.hh>
+#include <dumux/assembly/entitycolor.hh>
namespace Dumux {
@@ -72,14 +74,25 @@ public:
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix. The element residual is written into the right hand side.
*/
- void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
{
this->asImp_().bindLocalViews();
-
- const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables);
-
- for (const auto& scv : scvs(this->fvGeometry()))
- res[scv.dofIndex()] += residual[scv.indexInElement()];
+ const auto eIdxGlobal = this->assembler().fvGridGeometry().elementMapper().index(this->element());
+ if (partialReassembler
+ && partialReassembler->elementColor(eIdxGlobal) == EntityColor::green)
+ {
+ const auto residual = this->asImp_().evalLocalResidual(); // forward to the internal implementation
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.indexInElement()];
+ }
+ else
+ {
+ const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables, partialReassembler); // forward to the internal implementation
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.indexInElement()];
+ }
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
@@ -225,7 +238,9 @@ public:
*
* \return The element residual at the current solution.
*/
- ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
{
// get some aliases for convenience
const auto& element = this->element();
@@ -233,7 +248,7 @@ public:
const auto& curSol = this->curSol();
auto&& curElemVolVars = this->curElemVolVars();
- // get the vecor of the acutal element residuals
+ // get the vector of the actual element residuals
const auto origResiduals = this->evalLocalResidual();
//////////////////////////////////////////////////////////////////////////////////////////////////
@@ -280,14 +295,19 @@ public:
// update the global stiffness matrix with the current partial derivatives
for (auto&& scvJ : scvs(fvGeometry))
{
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // A[i][col][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at dof 'i'
- // depending on the primary variable 'pvIdx' at dof
- // 'col'.
- A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
- }
+ // don't add derivatives for green vertices
+ if (!partialReassembler
+ || partialReassembler->vertexColor(scvJ.dofIndex()) != EntityColor::green)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+ A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
+ }
+ }
}
// restore the original state of the scv's volume variables
@@ -336,8 +356,13 @@ public:
*
* \return The element residual at the current solution.
*/
- ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
{
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for explicit time discretization");
+
// get some aliases for convenience
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
@@ -437,8 +462,13 @@ public:
*
* \return The element residual at the current solution.
*/
- ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
{
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for analytic differentiation");
+
// get some aliases for convenience
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
@@ -551,8 +581,13 @@ public:
*
* \return The element residual at the current solution.
*/
- ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
{
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for explicit time discretization");
+
// get some aliases for convenience
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index b046e732d03db6554f5673fa99c0971d98cb4c7f..13dcf1acf746f4b0dd11529521a0c3fbaea28800 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -36,6 +36,8 @@
#include <dumux/assembly/numericepsilon.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/assembly/fvlocalassemblerbase.hh>
+#include <dumux/assembly/entitycolor.hh>
+#include <dumux/assembly/partialreassembler.hh>
#include <dumux/discretization/fluxstencil.hh>
namespace Dumux {
@@ -47,7 +49,7 @@ namespace Dumux {
* \tparam TypeTag The TypeTag
* \tparam Assembler The assembler type
* \tparam Implementation The actual implementation
- * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
+ * \tparam implicit Specifies whether the time discretization is implicit or not (i.e. explicit)
*/
template<class TypeTag, class Assembler, class Implementation, bool implicit>
class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
@@ -68,11 +70,21 @@ public:
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix. The element residual is written into the right hand side.
*/
- void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables)
+ template <class PartialReassembler = DefaultPartialReassembler>
+ void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler)
{
this->asImp_().bindLocalViews();
const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
- res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
+ if (partialReassembler
+ && partialReassembler->elementColor(globalI) == EntityColor::green)
+ {
+ res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
+ }
+ else
+ {
+ res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
+ }
}
/*!
diff --git a/dumux/assembly/entitycolor.hh b/dumux/assembly/entitycolor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4988ec347c0d17e25b97cace2e1e6922f31c8f6c
--- /dev/null
+++ b/dumux/assembly/entitycolor.hh
@@ -0,0 +1,57 @@
+// -*- 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
+ * \ingroup Assembly
+ * \brief An enum class to define the colors of elements and vertices required
+ * for partial Jacobian reassembly.
+ */
+#ifndef DUMUX_ENTITY_COLOR_HH
+#define DUMUX_ENTITY_COLOR_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief The colors of elements and vertices required for partial
+ * Jacobian reassembly.
+ */
+enum class EntityColor {
+ //! distance from last linearization is above the tolerance
+ red,
+
+ //! neighboring entity is red
+ yellow,
+
+ /*!
+ * A yellow entity that has only non-green neighbor elements.
+ *
+ * This means that its relative error is below the tolerance,
+ * but its defect can be linearized without any additional
+ * cost.
+ */
+ orange,
+
+ //! does not need to be reassembled
+ green
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index 032bbe8b14bcfed75fedcd9a2e1e7b1023937f87..efb670db12bdeabfbb04e087651527ed34495757 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -58,7 +58,8 @@ class FVAssembler
using TimeLoop = TimeLoopBase<typename GET_PROP_TYPE(TypeTag, Scalar)>;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+ static constexpr DiscretizationMethods discMethod = GET_PROP_VALUE(TypeTag, DiscretizationMethod);
+ static constexpr bool isBox = discMethod == DiscretizationMethods::Box;
using ThisType = FVAssembler<TypeTag, diffMethod, isImplicit>;
using LocalAssembler = std::conditional_t<isBox, BoxLocalAssembler<TypeTag, ThisType, diffMethod, isImplicit>,
@@ -68,6 +69,7 @@ public:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
using ResidualType = SolutionVector;
/*!
@@ -107,16 +109,17 @@ public:
* \brief Assembles the global Jacobian of the residual
* and the residual for the current solution.
*/
- void assembleJacobianAndResidual(const SolutionVector& curSol)
+ template<class PartialReassembler = DefaultPartialReassembler>
+ void assembleJacobianAndResidual(const SolutionVector& curSol, const PartialReassembler* partialReassembler = nullptr)
{
checkAssemblerState_();
- resetJacobian_();
+ resetJacobian_(partialReassembler);
resetResidual_();
assemble_([&](const Element& element)
{
LocalAssembler localAssembler(*this, element, curSol);
- localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_);
+ localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
});
}
@@ -329,8 +332,9 @@ private:
(*residual_) = 0.0;
}
- // reset the jacobian vector to 0.0
- void resetJacobian_()
+ // reset the Jacobian matrix to 0.0
+ template <class PartialReassembler = DefaultPartialReassembler>
+ void resetJacobian_(const PartialReassembler *partialReassembler = nullptr)
{
if(!jacobian_)
{
@@ -339,7 +343,10 @@ private:
setJacobianPattern();
}
- (*jacobian_) = 0.0;
+ if (partialReassembler)
+ partialReassembler->resetJacobian(*jacobian_, fvGridGeometry());
+ else
+ *jacobian_ = 0.0;
}
// check if the assembler is in a correct state for assembly
diff --git a/dumux/assembly/partialreassembler.hh b/dumux/assembly/partialreassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0daad5cb67e5b75bbcbf14f5fd0e4233e3dccbec
--- /dev/null
+++ b/dumux/assembly/partialreassembler.hh
@@ -0,0 +1,480 @@
+// -*- 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
+ * \ingroup Assembly
+ * \brief detects which entries in the Jacobian have to be recomputed
+ */
+#ifndef DUMUX_PARTIAL_REASSEMBLER_HH
+#define DUMUX_PARTIAL_REASSEMBLER_HH
+
+#include <algorithm>
+#include <vector>
+
+#include <dune/grid/common/gridenums.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+#include <dumux/common/typetraits/isvalid.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/parallel/vertexhandles.hh>
+
+#include "entitycolor.hh"
+
+namespace Dumux {
+
+class DefaultPartialReassembler
+{
+public:
+ DefaultPartialReassembler() = delete;
+
+ template<typename... Args>
+ void resetJacobian(Args&&... args) const {}
+
+ EntityColor elementColor(size_t idx) const
+ { return EntityColor::red; }
+
+ EntityColor vertexColor(size_t idx) const
+ { return EntityColor::red; }
+};
+
+//! the partial reassembler engine specialized for discretization methods
+template<class Assembler, DiscretizationMethods discMethod>
+class PartialReassemblerEngine
+{
+public:
+ PartialReassemblerEngine(const typename Assembler::FVGridGeometry&)
+ { DUNE_THROW(Dune::NotImplemented, "PartialReassembler for this discretization method!"); }
+
+ EntityColor elementColor(size_t idx) const
+ { return EntityColor::red; }
+
+ EntityColor dofColor(size_t idx) const
+ { return EntityColor::red; }
+
+ template<typename... Args>
+ std::size_t computeColors(Args&&... args) { return 0; }
+
+ template<typename... Args>
+ void resetJacobian(Args&&... args) const {}
+
+ template<typename... Args>
+ void resetColors(Args&&... args) {}
+};
+
+//! the partial reassembler engine specialized for the box method
+template<class Assembler>
+class PartialReassemblerEngine<Assembler, DiscretizationMethods::Box>
+{
+ using Scalar = typename Assembler::Scalar;
+ using FVGridGeometry = typename Assembler::FVGridGeometry;
+ using JacobianMatrix = typename Assembler::JacobianMatrix;
+ using VertexMapper = typename FVGridGeometry::VertexMapper;
+ static constexpr int dim = FVGridGeometry::GridView::dimension;
+
+public:
+ PartialReassemblerEngine(const FVGridGeometry& fvGridGeometry)
+ : elementColor_(fvGridGeometry.elementMapper().size(), EntityColor::red)
+ , vertexColor_(fvGridGeometry.vertexMapper().size(), EntityColor::red)
+ {}
+
+ // returns number of green elements
+ std::size_t computeColors(Scalar threshold,
+ const FVGridGeometry& fvGridGeometry,
+ const std::vector<Scalar>& distanceFromLastLinearization)
+ {
+ const auto& gridView = fvGridGeometry.gridView();
+ const auto& elementMapper = fvGridGeometry.elementMapper();
+ const auto& vertexMapper = fvGridGeometry.vertexMapper();
+
+ // set all vertices to green
+ vertexColor_.assign(vertexColor_.size(), EntityColor::green);
+
+ // mark the red vertices
+ for (unsigned int i = 0; i < vertexColor_.size(); ++i)
+ {
+ using std::max;
+ if (distanceFromLastLinearization[i] > threshold)
+ // mark vertex as red if discrepancy is larger than
+ // the relative tolerance
+ vertexColor_[i] = EntityColor::red;
+ }
+
+ // Mark all red elements
+ for (const auto& element : elements(gridView))
+ {
+ // find out whether the current element features a red vertex
+ bool isRed = false;
+
+ int numVertices = element.subEntities(dim);
+
+ for (int i = 0; i < numVertices; ++i) {
+ int globalI = vertexMapper.subIndex(element, i, dim);
+
+ if (vertexColor_[globalI] == EntityColor::red) {
+ isRed = true;
+ break;
+ }
+ }
+
+ int eIdx = elementMapper.index(element);
+ // if a vertex is red, the element color is also red, otherwise green
+ if (isRed)
+ elementColor_[eIdx] = EntityColor::red;
+ else
+ elementColor_[eIdx] = EntityColor::green;
+ }
+
+ // mark orange vertices
+ for (const auto& element : elements(gridView))
+ {
+ int eIdx = elementMapper.index(element);
+
+ // only red elements tint vertices yellow
+ if (elementColor_[eIdx] == EntityColor::red)
+ {
+ int numVertices = element.subEntities(dim);
+
+ for (int i = 0; i < numVertices; ++i) {
+ int globalI = vertexMapper.subIndex(element, i, dim);
+
+ // red vertices don't become orange
+ if (vertexColor_[globalI] != EntityColor::red)
+ vertexColor_[globalI] = EntityColor::orange;
+ }
+ }
+ }
+
+ // at this point we communicate the yellow vertices to the
+ // neighboring processes because a neigbor process may not see
+ // the red vertex for yellow border vertices
+ VertexHandleMin<EntityColor, std::vector<EntityColor>, VertexMapper>
+ minHandle(vertexColor_, vertexMapper);
+ gridView.communicate(minHandle,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
+
+ // mark yellow elements
+ for (const auto& element : elements(gridView))
+ {
+ int eIdx = elementMapper.index(element);
+
+ // only treat non-red elements
+ if (elementColor_[eIdx] != EntityColor::red)
+ {
+ // check whether the element features a orange vertex
+ bool isOrange = false;
+ int numVertices = element.subEntities(dim);
+
+ for (int i = 0; i < numVertices; ++i) {
+ int globalI = vertexMapper.subIndex(element, i, dim);
+
+ if (vertexColor_[globalI] == EntityColor::orange) {
+ isOrange = true;
+ break;
+ }
+ }
+
+ if (isOrange)
+ elementColor_[eIdx] = EntityColor::yellow;
+ }
+ }
+
+ // change orange vertices to yellow ones if it has at least
+ // one green element as a neighbor
+ for (const auto& element : elements(gridView))
+ {
+ int eIdx = elementMapper.index(element);
+
+ // only green elements are considered
+ if (elementColor_[eIdx] == EntityColor::green)
+ {
+ int numVertices = element.subEntities(dim);
+
+ for (int i = 0; i < numVertices; ++i) {
+ int globalI = vertexMapper.subIndex(element, i, dim);
+
+ // if a vertex is orange, recolor it to yellow
+ if (vertexColor_[globalI] == EntityColor::orange)
+ vertexColor_[globalI] = EntityColor::yellow;
+ }
+ }
+ }
+
+ // demote the border orange vertices
+ VertexHandleMax<EntityColor, std::vector<EntityColor>, VertexMapper>
+ maxHandle(vertexColor_, vertexMapper);
+ gridView.communicate(maxHandle,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
+
+ // promote the remaining orange vertices to red
+ for (unsigned int i=0; i < vertexColor_.size(); ++i) {
+ // if a vertex is green or yellow don't do anything
+ if (vertexColor_[i] == EntityColor::green || vertexColor_[i] == EntityColor::yellow)
+ continue;
+
+ // set the vertex to red
+ vertexColor_[i] = EntityColor::red;
+ }
+
+ // count green elements
+ return std::count_if(elementColor_.begin(), elementColor_.end(),
+ [](EntityColor c){ return c == EntityColor::green; });
+ }
+
+ void resetJacobian(JacobianMatrix& jacobian,
+ const FVGridGeometry& fvGridGeometry) const
+ {
+ // loop over all dofs
+ for (unsigned int rowIdx = 0; rowIdx < jacobian.N(); ++rowIdx)
+ {
+ // reset all entries corrosponding to a non-green vertex
+ if (vertexColor_[rowIdx] != EntityColor::green)
+ {
+ // set all matrix entries in the row to 0
+ auto colIt = jacobian[rowIdx].begin();
+ const auto& colEndIt = jacobian[rowIdx].end();
+ for (; colIt != colEndIt; ++colIt) {
+ *colIt = 0.0;
+ }
+ }
+ }
+ }
+
+ void resetColors()
+ {
+ elementColor_.assign(elementColor_.size(), EntityColor::red);
+ vertexColor_.assign(vertexColor_.size(), EntityColor::red);
+ }
+
+ EntityColor elementColor(size_t idx) const
+ { return elementColor_[idx]; }
+
+ EntityColor vertexColor(size_t idx) const
+ { return vertexColor_[idx]; }
+
+ EntityColor dofColor(size_t idx) const
+ { return vertexColor_[idx]; }
+
+private:
+ //! entity colors for partial reassembly
+ std::vector<EntityColor> elementColor_;
+ std::vector<EntityColor> vertexColor_;
+};
+
+//! the partial reassembler engine specialized for the box method
+template<class Assembler>
+class PartialReassemblerEngine<Assembler, DiscretizationMethods::CCTpfa>
+{
+ using Scalar = typename Assembler::Scalar;
+ using FVGridGeometry = typename Assembler::FVGridGeometry;
+ using JacobianMatrix = typename Assembler::JacobianMatrix;
+
+public:
+ PartialReassemblerEngine(const FVGridGeometry& fvGridGeometry)
+ : elementColor_(fvGridGeometry.elementMapper().size(), EntityColor::red)
+ {}
+
+ // returns number of green elements
+ std::size_t computeColors(Scalar threshold,
+ const FVGridGeometry& fvGridGeometry,
+ const std::vector<Scalar>& distanceFromLastLinearization)
+ {
+ const auto& gridView = fvGridGeometry.gridView();
+ const auto& elementMapper = fvGridGeometry.elementMapper();
+
+ // mark the red elements
+ for (const auto& element : elements(gridView))
+ {
+ int eIdx = elementMapper.index(element);
+
+ if (distanceFromLastLinearization[eIdx] > threshold)
+ {
+ // mark element as red if discrepancy is larger than
+ // the relative tolerance
+ elementColor_[eIdx] = EntityColor::red;
+ }
+ else
+ {
+ elementColor_[eIdx] = EntityColor::green;
+ }
+ }
+
+ // mark the neighbors also red
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ for (unsigned eIdx = 0; eIdx < elementColor_.size(); ++eIdx)
+ {
+ if (elementColor_[eIdx] == EntityColor::red)
+ continue; // element is red already!
+
+ if (distanceFromLastLinearization[eIdx] > threshold)
+ {
+ for (const auto& connectedDof : connectivityMap[eIdx])
+ elementColor_[connectedDof.globalJ] = EntityColor::red;
+ }
+ }
+
+ // count green elements
+ return std::count_if(elementColor_.begin(), elementColor_.end(),
+ [](EntityColor c){return c == EntityColor::green;});
+
+ }
+
+ void resetJacobian(JacobianMatrix& jacobian,
+ const FVGridGeometry& fvGridGeometry) const
+ {
+ // loop over all dofs
+ for (unsigned int colIdx = 0; colIdx < jacobian.M(); ++colIdx)
+ {
+ // reset all entries corresponding to a non-green element
+ if (elementColor_[colIdx] != EntityColor::green)
+ {
+ // set all matrix entries in the column to 0
+ jacobian[colIdx][colIdx] = 0;
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ for (const auto& dataJ : connectivityMap[colIdx])
+ jacobian[dataJ.globalJ][colIdx] = 0;
+ }
+ }
+ }
+
+ void resetColors()
+ {
+ elementColor_.assign(elementColor_.size(), EntityColor::red);
+ }
+
+ EntityColor elementColor(size_t idx) const
+ { return elementColor_[idx]; }
+
+ EntityColor dofColor(size_t idx) const
+ { return elementColor_[idx]; }
+
+private:
+ //! entity colors for partial reassembly
+ std::vector<EntityColor> elementColor_;
+};
+
+//! the partial reassembler engine specialized for the mpfa method
+template<class Assembler>
+class PartialReassemblerEngine<Assembler, DiscretizationMethods::CCMpfa>
+: public PartialReassemblerEngine<Assembler, DiscretizationMethods::CCTpfa>
+{
+ using ParentType = PartialReassemblerEngine<Assembler, DiscretizationMethods::CCTpfa>;
+public:
+ using ParentType::ParentType;
+};
+
+/*!
+ * \ingroup Assembly
+ * \brief detects which entries in the Jacobian have to be recomputed
+ * \tparam TypeTag The TypeTag
+ */
+template<class Assembler>
+class PartialReassembler
+{
+ using Scalar = typename Assembler::Scalar;
+ using FVGridGeometry = typename Assembler::FVGridGeometry;
+ using JacobianMatrix = typename Assembler::JacobianMatrix;
+ using VertexMapper = typename FVGridGeometry::VertexMapper;
+
+ static constexpr DiscretizationMethods discMethod = FVGridGeometry::discretizationMethod;
+ using Engine = PartialReassemblerEngine<Assembler, discMethod>;
+
+ static constexpr auto hasVertexColor = Dumux::isValid([](auto&& a) -> decltype(a.vertexColor(0)) {});
+ static constexpr bool engineHasVertexColor = decltype(hasVertexColor(std::declval<Engine>()))::value;
+
+public:
+
+ /*!
+ * \brief constructor
+ * \param fvGridGeometry the employed FVGridGeometry
+ */
+ PartialReassembler(const FVGridGeometry& fvGridGeometry)
+ : engine_(fvGridGeometry)
+ , greenElems_(0)
+ {
+ totalElems_ = fvGridGeometry.elementMapper().size();
+ totalElems_ = fvGridGeometry.gridView().comm().sum(totalElems_);
+ }
+
+ /*!
+ * \brief Determine the colors of entities for partial reassembly.
+ *
+ * The following approach is used:
+ *
+ * - Set all elements to 'green'
+ * - Mark all elements as 'red' which exhibit an relative error above
+ * the tolerance
+ * - Mark all neighbors of 'red' elements also 'red'
+ *
+ * \param threshold Reassemble only if the distance from the last
+ * linearization is above this value.
+ */
+ void computeColors(Scalar threshold,
+ const FVGridGeometry& fvGridGeometry,
+ const std::vector<Scalar>& distanceFromLastLinearization)
+ {
+ greenElems_ = engine_.computeColors(threshold, fvGridGeometry, distanceFromLastLinearization);
+ }
+
+ void resetColors()
+ {
+ engine_.resetColors();
+ }
+
+ void resetJacobian(JacobianMatrix& jacobian, const FVGridGeometry& fvGridGeometry) const
+ {
+ engine_.resetJacobian(jacobian, fvGridGeometry);
+ }
+
+ /*!
+ * \brief called by the assembler after successful assembly
+ */
+ template <class Communication>
+ void report(const Communication& comm, std::ostream& outStream)
+ {
+ if (comm.size() > 1)
+ greenElems_ = comm.sum(greenElems_);
+
+ auto reassembledElems = totalElems_ - greenElems_;
+ auto width = std::to_string(totalElems_).size();
+ outStream << ", reassembled " << std::setw(width)
+ << reassembledElems << " (" << std::setw(3)
+ << 100*reassembledElems/totalElems_ << "%) elements";
+ }
+
+ EntityColor elementColor(size_t idx) const
+ { return engine_.elementColor(idx); }
+
+ EntityColor dofColor(size_t idx) const
+ { return engine_.dofColor(idx); }
+
+ template<bool enable = engineHasVertexColor, typename std::enable_if_t<enable, int> = 0>
+ EntityColor vertexColor(size_t idx) const
+ { return engine_.vertexColor(idx); }
+
+private:
+ Engine engine_;
+ size_t totalElems_;
+ size_t greenElems_;
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/common/parameters.hh b/dumux/common/parameters.hh
index 7655211a20516297f6d507b40e4229ea473d9bef..e4a75101e831fcdb2524ce11b892f91a6ae694f3 100644
--- a/dumux/common/parameters.hh
+++ b/dumux/common/parameters.hh
@@ -352,7 +352,6 @@ private:
{
// parameters in the implicit group
params["Implicit.UpwindWeight"] = "1.0";
- params["Implicit.EnablePartialReassemble"] = "false";
params["Implicit.EnableJacobianRecycling"] = "false";
// parameters in the assembly group
@@ -369,7 +368,7 @@ private:
// parameters in the problem group
params["Problem.EnableGravity"] = "true";
- // parameters in the newton group
+ // parameters in the Newton group
params["Newton.MaxSteps"] = "18";
params["Newton.TargetSteps"] = "10";
params["Newton.UseLineSearch"] = "false";
@@ -381,6 +380,7 @@ private:
params["Newton.EnableAbsoluteResidualCriterion"] = "false";
params["Newton.MaxAbsoluteResidual"] = "1e-5";
params["Newton.SatisfyResidualAndShiftCriterion"] = "false";
+ params["Newton.EnablePartialReassembly"] = "false";
// parameters in the time loop group
params["TimeLoop.MaxTimeStepSize"] = "1e300";
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 0a5e1a3a02caf5690df1c2f6c84afce3fc3327ce..aa7dfb5eb60fcd2142f6f2717565cc4fc2017bc4 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -41,8 +41,10 @@
#include <dumux/common/exceptions.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/typetraits/vector.hh>
+#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
+#include <dumux/assembly/partialreassembler.hh>
#include "newtonconvergencewriter.hh"
@@ -65,7 +67,15 @@ class NewtonSolver
using JacobianMatrix = typename Assembler::JacobianMatrix;
using SolutionVector = typename Assembler::ResidualType;
using ConvergenceWriter = ConvergenceWriterInterface<SolutionVector>;
-
+ using Reassembler = PartialReassembler<Assembler>;
+
+ // TODO: Require Reassembler argument from the standard interface
+ static constexpr auto supportsPartialReassembly =
+ Dumux::isValid([](auto&& a) -> decltype(
+ a.assembleJacobianAndResidual(std::declval<Assembler>(),
+ std::declval<const SolutionVector&>(),
+ std::declval<const Reassembler*>())
+ ){});
public:
using Communication = Comm;
@@ -91,6 +101,13 @@ public:
// set a different default for the linear solver residual reduction
// within the Newton the linear solver doesn't need to solve too exact
linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(paramGroup, "LinearSolver.ResidualReduction", 1e-6));
+
+ // initialize the partial reassembler
+ if (enablePartialReassembly_)
+ {
+ partialReassembler_ = std::make_unique<Reassembler>(assembler_->fvGridGeometry());
+ distanceFromLastLinearization_.resize(assembler_->fvGridGeometry().numDofs(), 0);
+ }
}
/*!
@@ -116,6 +133,13 @@ public:
// set a different default for the linear solver residual reduction
// within the Newton the linear solver doesn't need to solve too exact
linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(paramGroup, "LinearSolver.ResidualReduction", 1e-6));
+
+ // initialize the partial reassembler
+ if (enablePartialReassembly_)
+ {
+ partialReassembler_ = std::make_unique<Reassembler>(assembler_->fvGridGeometry());
+ distanceFromLastLinearization_.resize(assembler_->fvGridGeometry().numDofs(), 0);
+ }
}
//! the communicator for parallel runs
@@ -188,6 +212,12 @@ public:
Dune::Timer solveTimer(false);
Dune::Timer updateTimer(false);
+ if (enablePartialReassembly_)
+ {
+ partialReassembler_->resetColors();
+ std::fill(distanceFromLastLinearization_.begin(),
+ distanceFromLastLinearization_.end(), 0.0);
+ }
newtonBegin(uCurrentIter);
// execute the method as long as the controller thinks
@@ -360,7 +390,10 @@ public:
*/
virtual void assembleLinearSystem(const SolutionVector& uCurrentIter)
{
- assembler_->assembleJacobianAndResidual(uCurrentIter);
+ assembleLinearSystem_(uCurrentIter);
+
+ if (enablePartialReassembly_)
+ partialReassembler_->report(comm_, endIterMsgStream_);
}
/*!
@@ -445,9 +478,44 @@ public:
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
- if (enableShiftCriterion_)
+ if (enableShiftCriterion_ || enablePartialReassembly_)
newtonUpdateShift_(uLastIter, deltaU);
+ if (enablePartialReassembly_) {
+ // Determine the threshold 'eps' that is used for the partial reassembly.
+ // Every entity where the primary variables exhibit a relative shift
+ // summed up since the last linearization above 'eps' will be colored
+ // red yielding a reassembly.
+ // The user can provide three parameters to influence the threshold:
+ // 'minEps' by 'Newton.ReassemblyMinThreshold' (1e-1*shiftTolerance_ default)
+ // 'maxEps' by 'Newton.ReassemblyMaxThreshold' (1e2*shiftTolerance_ default)
+ // 'omega' by 'Newton.ReassemblyShiftWeight' (1e-3 default)
+ // The threshold is calculated from the currently achieved maximum
+ // relative shift according to the formula
+ // eps = max( minEps, min(maxEps, omega*shift) ).
+ // Increasing/decreasing 'minEps' leads to less/more reassembly if
+ // 'omega*shift' is small, i.e., for the last Newton iterations.
+ // Increasing/decreasing 'maxEps' leads to less/more reassembly if
+ // 'omega*shift' is large, i.e., for the first Newton iterations.
+ // Increasing/decreasing 'omega' leads to more/less first and last
+ // iterations in this sense.
+ using std::max;
+ using std::min;
+ auto reassemblyThreshold = max(reassemblyMinThreshold_,
+ min(reassemblyMaxThreshold_,
+ shift_*reassemblyShiftWeight_));
+
+ updateDistanceFromLastLinearization_(uLastIter, deltaU);
+ partialReassembler_->computeColors(reassemblyThreshold,
+ assembler_->fvGridGeometry(),
+ distanceFromLastLinearization_);
+
+ // set the discrepancy of the red entities to zero
+ for (unsigned int i = 0; i < distanceFromLastLinearization_.size(); i++)
+ if (partialReassembler_->dofColor(i) == EntityColor::red)
+ distanceFromLastLinearization_[i] = 0;
+ }
+
if (useLineSearch_)
lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
@@ -680,6 +748,22 @@ protected:
private:
+ //! assembleLinearSystem_ for assemblers that support partial reassembly
+ template<class A = Assembler>
+ auto assembleLinearSystem_(const SolutionVector& uCurrentIter)
+ -> typename std::enable_if_t<decltype(supportsPartialReassembly(std::declval<A>()))::value, void>
+ {
+ assembler_->assembleJacobianAndResidual(uCurrentIter, partialReassembler_.get());
+ }
+
+ //! assembleLinearSystem_ for assemblers that don't support partial reassembly
+ template<class A = Assembler>
+ auto assembleLinearSystem_(const SolutionVector& uCurrentIter)
+ -> typename std::enable_if_t<!decltype(supportsPartialReassembly(std::declval<A>()))::value, void>
+ {
+ assembler_->assembleJacobianAndResidual(uCurrentIter);
+ }
+
/*!
* \brief Update the maximum relative shift of the solution compared to
* the previous iteration. Overload for "normal" solution vectors.
@@ -928,10 +1012,37 @@ private:
setTargetSteps(getParamFromGroup<int>(group, "Newton.TargetSteps"));
setMaxSteps(getParamFromGroup<int>(group, "Newton.MaxSteps"));
+ enablePartialReassembly_ = getParamFromGroup<bool>(group, "Newton.EnablePartialReassembly");
+ reassemblyMinThreshold_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyMinThreshold", 1e-1*shiftTolerance_);
+ reassemblyMaxThreshold_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyMaxThreshold", 1e2*shiftTolerance_);
+ reassemblyShiftWeight_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyShiftWeight", 1e-3);
+
verbose_ = comm_.rank() == 0;
numSteps_ = 0;
}
+ template<class SolVec>
+ void updateDistanceFromLastLinearization_(const SolVec &u,
+ const SolVec &uDelta)
+ {
+ for (size_t i = 0; i < u.size(); ++i) {
+ const auto& currentPriVars(u[i]);
+ auto nextPriVars(currentPriVars);
+ nextPriVars -= uDelta[i];
+
+ // add the current relative shift for this degree of freedom
+ auto shift = relativeShiftAtDof_(currentPriVars, nextPriVars);
+ distanceFromLastLinearization_[i] += shift;
+ }
+ }
+
+ template<class ...Args>
+ void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<Args...> &uLastIter,
+ const Dune::MultiTypeBlockVector<Args...> &deltaU)
+ {
+ //! \todo implement the function for MultiTypeBlockVectors
+ }
+
/*!
* \brief Returns the maximum relative shift between two vectors of
* primary variables.
@@ -941,7 +1052,7 @@ private:
*/
template<class PrimaryVariables>
Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
- const PrimaryVariables &priVars2)
+ const PrimaryVariables &priVars2) const
{
Scalar result = 0.0;
using std::abs;
@@ -973,7 +1084,6 @@ private:
Scalar residualTolerance_;
// further parameters
- bool enablePartialReassemble_;
bool useLineSearch_;
bool useChop_;
bool enableAbsoluteResidualCriterion_;
@@ -984,6 +1094,13 @@ private:
//! the parameter group for getting parameters from the parameter tree
std::string paramGroup_;
+ // infrastructure for partial reassembly
+ bool enablePartialReassembly_;
+ std::unique_ptr<Reassembler> partialReassembler_;
+ std::vector<Scalar> distanceFromLastLinearization_;
+ Scalar reassemblyMinThreshold_;
+ Scalar reassemblyMaxThreshold_;
+ Scalar reassemblyShiftWeight_;
};
} // end namespace Dumux
diff --git a/test/porousmediumflow/2p/implicit/incompressible/test_2p.input b/test/porousmediumflow/2p/implicit/incompressible/test_2p.input
index 93f7f1f1ee922e699bd27f853eda26be4521f705..336df6cf83b0f396fee77f0b57e0fabb6bcb1eaf 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/test_2p.input
+++ b/test/porousmediumflow/2p/implicit/incompressible/test_2p.input
@@ -14,3 +14,6 @@ LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
[Problem]
Name = 2p
EnableGravity = true
+
+[Newton]
+EnablePartialReassembly = true