From 93bf19c1e740928b0a3df7a6b5b5851d4d560193 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 13 Feb 2018 17:17:41 +0100
Subject: [PATCH] [assembly][newton] reimplement partial reassembly
Transfer the partial reassembly functionality to the current
structure.
The `NewtonSolver` is responsible for the partial reassembly.
It owns a `unique_ptr` to an object of a new class
`PartialReassembler` which is invoked in the constructor as well as
in the methods `solve`, `newtonUpdate` and `assembleLinearSystem`.
The `PartialReassembler` object is created only if the parameter
`Newton.EnablePartialReassembly` is set to true.
After each Newton iteration, the `PartialReassembler` colors the
geometrical degrees of freedom depending on the shift of the
primary variables from the last time that the derivatives have been
calculated. This shift is updated by the function
`updateDistanceFromLastLinearization_` of the `NewtonSolver` and
passed to the function `computeColors` of the `PartialReassembler`
for the actual coloring. For green-colored entities, the entries in
the Jacobian are kept, while for the other entities, they are
recalculated.
Use an enum class `EntityColor` for the coloring, enable read
access to the color of an entity via the `PartialReassembler`.
Pass a pointer to the `PartialReassembler` when calling the routine
`assembleJacobianAndResidual` of the `Assembler`. Use SFINAE and
`enable_if` so that only assemblers accepting the additional
argument are invoked. The correspnding local assemblers have to be
adapted such that no calculations are performed for green elements.
Implement the functionality for Box, TPFA and MPFA; adapt the
`FVAssembler` as well as the `Box`/`CCLocalAssembler`
correspondingly.
The `PartialReassembler` throws if constructed with a discretization
method equal to None or Staggered. All `BoxLocalAssembler`s apart
from the implicit numeric one throw if called with a non-nullptr
to a `PartialReassembler`.
Allow for a default `nullptr` parameter to the assembly functions
that take a `PartialReassembler`. This requires a new class
`DefaultPartialReassembler` to be passed as a default template
parameter for these functions. This makes it possible to use the
functions `assembleJacobianAndResidual` of the global and local
assemblers as before without the additional argument.
Allow to set the thresholds and the weight for the calculation of
the effective reassembly threshold as parameters. Try to explain
the corresponding formula in a comment.
Use the functionality in the two-phase incompressible tests.
---
dumux/assembly/CMakeLists.txt | 2 +
dumux/assembly/boxlocalassembler.hh | 73 ++-
dumux/assembly/cclocalassembler.hh | 18 +-
dumux/assembly/entitycolor.hh | 57 +++
dumux/assembly/fvassembler.hh | 21 +-
dumux/assembly/partialreassembler.hh | 481 ++++++++++++++++++
dumux/common/parameters.hh | 4 +-
dumux/nonlinear/newtonsolver.hh | 128 ++++-
.../2p/implicit/incompressible/test_2p.input | 3 +
9 files changed, 751 insertions(+), 36 deletions(-)
create mode 100644 dumux/assembly/entitycolor.hh
create mode 100644 dumux/assembly/partialreassembler.hh
diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index c8210e26b7..bd2ebb628d 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 104e218929..2601c5d90a 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 b046e732d0..13dcf1acf7 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 0000000000..4988ec347c
--- /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 032bbe8b14..efb670db12 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 0000000000..92eeee578a
--- /dev/null
+++ b/dumux/assembly/partialreassembler.hh
@@ -0,0 +1,481 @@
+// -*- 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
+ */
+ void finalizeAssembly(const FVGridGeometry& fvGridGeometry, std::ostream& outStream)
+ {
+ const auto& gridView = fvGridGeometry.gridView();
+
+ if (gridView.comm().size() > 1)
+ greenElems_ = gridView.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 7655211a20..e4a75101e8 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 0a5e1a3a02..e1c6901b72 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,11 @@ public:
*/
virtual void assembleLinearSystem(const SolutionVector& uCurrentIter)
{
- assembler_->assembleJacobianAndResidual(uCurrentIter);
+ assembleLinearSystem_(uCurrentIter);
+
+ if (enablePartialReassembly_)
+ partialReassembler_->finalizeAssembly(assembler_->fvGridGeometry(),
+ endIterMsgStream_);
}
/*!
@@ -445,9 +479,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 +749,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 +1013,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 +1053,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 +1085,6 @@ private:
Scalar residualTolerance_;
// further parameters
- bool enablePartialReassemble_;
bool useLineSearch_;
bool useChop_;
bool enableAbsoluteResidualCriterion_;
@@ -984,6 +1095,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 93f7f1f1ee..336df6cf83 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
--
GitLab