diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index f7ea8bb601bd42ca26054537cdfdf28f24bcaf70..c8210e26b71d2cb8abd48c67ee76b5e36395c934 100644
--- a/dumux/assembly/CMakeLists.txt
+++ b/dumux/assembly/CMakeLists.txt
@@ -6,7 +6,9 @@ cclocalassembler.hh
cclocalresidual.hh
diffmethod.hh
fvassembler.hh
+fvlocalassemblerbase.hh
fvlocalresidual.hh
+jacobianpattern.hh
staggeredfvassembler.hh
staggeredlocalassembler.hh
staggeredlocalresidual.hh
diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index 9487428f95ef5381fd0ae9118cd6f782b80f256e..9c9a3df0b383a9af57d7d4ad96b8513d05f70b7a 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -30,6 +30,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/common/numericdifferentiation.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/assembly/fvlocalassemblerbase.hh>
@@ -205,6 +206,8 @@ class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/tr
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { dim = GET_PROP_TYPE(TypeTag, GridView)::dimension };
@@ -221,7 +224,7 @@ public:
*
* \return The element residual at the current solution.
*/
- auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
// get some aliases for convenience
const auto& element = this->element();
@@ -243,7 +246,7 @@ public:
// create the element solution
ElementSolutionVector elemSol(element, curSol, fvGeometry);
- using ElementResidualVector = std::decay_t<decltype(origResiduals)>;
+ // create the vector storing the partial derivatives
ElementResidualVector partialDerivs(element.subEntities(dim));
// calculation of the derivatives
@@ -268,7 +271,8 @@ public:
};
// derive the residuals numerically
- NumericDifferentiation::partialDerivative(evalResiduals, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals);
+ static const int numDiffMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals, numDiffMethod);
// update the global stiffness matrix with the current partial derivatives
for (auto&& scvJ : scvs(fvGeometry))
@@ -313,6 +317,8 @@ class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/fa
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { dim = GET_PROP_TYPE(TypeTag, GridView)::dimension };
@@ -327,7 +333,7 @@ public:
*
* \return The element residual at the current solution.
*/
- auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
// get some aliases for convenience
const auto& element = this->element();
@@ -349,7 +355,7 @@ public:
// create the element solution
ElementSolutionVector elemSol(element, curSol, fvGeometry);
- using ElementResidualVector = std::decay_t<decltype(origResiduals)>;
+ // create the vector storing the partial derivatives
ElementResidualVector partialDerivs(element.subEntities(dim));
// calculation of the derivatives
@@ -374,7 +380,8 @@ public:
};
// derive the residuals numerically
- NumericDifferentiation::partialDerivative(evalStorage, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals);
+ static const int numDiffMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalStorage, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals, numDiffMethod);
// update the global stiffness matrix with the current partial derivatives
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
@@ -412,6 +419,8 @@ class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/t
using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
public:
@@ -423,7 +432,7 @@ public:
*
* \return The element residual at the current solution.
*/
- auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
// get some aliases for convenience
const auto& element = this->element();
@@ -524,6 +533,8 @@ class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/f
using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
public:
@@ -535,7 +546,7 @@ public:
*
* \return The element residual at the current solution.
*/
- auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
// get some aliases for convenience
const auto& element = this->element();
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 00d9c51e2d82bf7b0d11df075241cce916f99710..49774a6c97302b5086ad947d1576bc762bc3294e 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -32,8 +32,8 @@
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/common/numericdifferentiation.hh>
#include <dumux/assembly/diffmethod.hh>
-#include <dumux/assembly/numericdifferentiation.hh>
#include <dumux/assembly/fvlocalassemblerbase.hh>
#include <dumux/discretization/fluxstencil.hh>
@@ -245,7 +245,8 @@ public:
};
// derive the residuals numerically
- NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals);
+ static const int numDiffMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals, numDiffMethod);
// add the current partial derivatives to the global jacobian matrix
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
@@ -363,7 +364,10 @@ public:
// for non-ghosts compute the derivative numerically
if (!this->elementIsGhost())
- NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, residual);
+ {
+ static const int numDiffMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, residual, numDiffMethod);
+ }
// for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
// as we always solve for a delta of the solution with repect to the initial solution this
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index 4310071b2a7d97005dc21e045ac3cd38a330c8a6..abeeb48e597ae5335e52e7f881cee5f92451dfa3 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -33,6 +33,7 @@
#include <dumux/discretization/methods.hh>
#include <dumux/parallel/vertexhandles.hh>
+#include "jacobianpattern.hh"
#include "diffmethod.hh"
#include "boxlocalassembler.hh"
#include "cclocalassembler.hh"
@@ -227,11 +228,7 @@ public:
jacobian_->setSize(numDofs, numDofs);
// create occupation pattern of the jacobian
- Dune::MatrixIndexSet occupationPattern;
- occupationPattern.resize(numDofs, numDofs);
-
- // set the jacobian pattern depending on space and time discretization
- setJacobianPattern_(occupationPattern, numDofs);
+ const auto occupationPattern = getJacobianPattern<isImplicit>(fvGridGeometry());
// export pattern to jacobian
occupationPattern.exportIdx(*jacobian_);
@@ -336,65 +333,6 @@ private:
DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
}
- //! Implicit jacobian pattern for cell-centered fv schemes
- template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) != DiscretizationMethods::Box, int> = 0>
- void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
- {
- // matrix pattern for implicit Jacobians
- if (isImplicit)
- {
- for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
- {
- pattern.add(globalI, globalI);
- for (const auto& dataJ : fvGridGeometry().connectivityMap()[globalI])
- pattern.add(dataJ.globalJ, globalI);
-
- // reserve index for additional user defined DOF dependencies
- // const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
- // for (auto globalJ : additionalDofDependencies)
- // pattern.add(globalI, globalJ);
- }
- }
-
- // matrix pattern for explicit Jacobians -> diagonal matrix
- else
- {
- for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
- pattern.add(globalI, globalI);
- }
- }
-
- //! Implicit jacobian pattern for vertex-centered fv schemes
- template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Box, int> = 0>
- void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
- {
- // matrix pattern for implicit Jacobians
- if (isImplicit)
- {
- static constexpr int dim = GridView::dimension;
- for (const auto& element : elements(fvGridGeometry().gridView()))
- {
- for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
- {
- const auto globalI = fvGridGeometry().vertexMapper().subIndex(element, vIdx, dim);
- for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
- {
- const auto globalJ = fvGridGeometry().vertexMapper().subIndex(element, vIdx2, dim);
- pattern.add(globalI, globalJ);
- pattern.add(globalJ, globalI);
- }
- }
- }
- }
-
- // matrix pattern for explicit Jacobians -> diagonal matrix
- else
- {
- for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
- pattern.add(globalI, globalI);
- }
- }
-
/*!
* \brief A method assembling something per element
* \note Handles exceptions for parallel runs
diff --git a/dumux/assembly/fvlocalassemblerbase.hh b/dumux/assembly/fvlocalassemblerbase.hh
index dbae670f0b5ca590c0bb6867d820e56afccde073..7758e8af9b06a12c8e160bb2f2422d474d6ea149 100644
--- a/dumux/assembly/fvlocalassemblerbase.hh
+++ b/dumux/assembly/fvlocalassemblerbase.hh
@@ -32,7 +32,6 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/assembly/diffmethod.hh>
-#include <dumux/assembly/numericdifferentiation.hh>
namespace Dumux {
@@ -41,13 +40,14 @@ namespace Dumux {
* \brief A base class for all local assemblers
* \tparam TypeTag The TypeTag
* \tparam Assembler The assembler type
- * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
+ * \tparam isImplicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
*/
-template<class TypeTag, class Assembler, class Implementation, bool implicit>
+template<class TypeTag, class Assembler, class Implementation, bool isImplicit>
class FVLocalAssemblerBase
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
@@ -86,30 +86,28 @@ public:
* \brief Convenience function to evaluate the complete local residual for the current element. Automatically chooses the the appropriate
* element volume variables.
*/
- auto evalLocalResidual() const
+ ElementResidualVector evalLocalResidual() const
{
- if(this->assembler().isStationaryProblem() && !isImplicit())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+ if (!isImplicit)
+ if (this->assembler().isStationaryProblem())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
- if(elementIsGhost())
- {
- using ResdiualType = decltype(evalLocalResidual(std::declval<ElementVolumeVariables>()));
- return ResdiualType(0.0);
- }
+ if (elementIsGhost())
+ return ElementResidualVector(0.0);
- return isImplicit() ? evalLocalResidual(curElemVolVars())
- : evalLocalResidual(prevElemVolVars());
+ return isImplicit ? evalLocalResidual(curElemVolVars())
+ : evalLocalResidual(prevElemVolVars());
}
/*!
* \brief Evaluates the complete local residual for the current element.
* \param elemVolVars The element volume variables
*/
- auto evalLocalResidual(const ElementVolumeVariables& elemVolVars) const
+ ElementResidualVector evalLocalResidual(const ElementVolumeVariables& elemVolVars) const
{
if (!assembler().isStationaryProblem())
{
- auto residual = evalLocalFluxAndSourceResidual(elemVolVars);
+ ElementResidualVector residual = evalLocalFluxAndSourceResidual(elemVolVars);
residual += evalLocalStorageResidual();
return residual;
}
@@ -122,10 +120,10 @@ public:
* of the local residual for the current element. Automatically chooses the the appropriate
* element volume variables.
*/
- auto evalLocalFluxAndSourceResidual() const
+ ElementResidualVector evalLocalFluxAndSourceResidual() const
{
- return isImplicit() ? evalLocalFluxAndSourceResidual(curElemVolVars())
- : evalLocalFluxAndSourceResidual(prevElemVolVars());
+ return isImplicit ? evalLocalFluxAndSourceResidual(curElemVolVars())
+ : evalLocalFluxAndSourceResidual(prevElemVolVars());
}
/*!
@@ -134,9 +132,9 @@ public:
*
* \param elemVolVars The element volume variables
*/
- auto evalLocalFluxAndSourceResidual(const ElementVolumeVariables& elemVolVars) const
+ ElementResidualVector evalLocalFluxAndSourceResidual(const ElementVolumeVariables& elemVolVars) const
{
- return localResidual_.evalFluxSource(element_, fvGeometry_, elemVolVars, elemFluxVarsCache_, elemBcTypes_);
+ return localResidual_.evalFluxAndSource(element_, fvGeometry_, elemVolVars, elemFluxVarsCache_, elemBcTypes_);
}
/*!
@@ -144,7 +142,7 @@ public:
* of the local residual for the current element. Automatically chooses the the appropriate
* element volume variables.
*/
- auto evalLocalStorageResidual() const
+ ElementResidualVector evalLocalStorageResidual() const
{
return localResidual_.evalStorage(element_, fvGeometry_, prevElemVolVars_, curElemVolVars_);
}
@@ -167,7 +165,7 @@ public:
// bind the caches
fvGeometry.bind(element);
- if(isImplicit())
+ if (isImplicit)
{
curElemVolVars.bind(element, fvGeometry, curSol);
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
@@ -250,9 +248,6 @@ public:
const LocalResidual& localResidual() const
{ return localResidual_; }
- constexpr bool isImplicit() const
- { return implicit; }
-
protected:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
diff --git a/dumux/assembly/fvlocalresidual.hh b/dumux/assembly/fvlocalresidual.hh
index 3e06c94d388f9e506e94027724775678c5d5f0d1..a26e90e227ffbccf9e57a3481c359aa802629e12 100644
--- a/dumux/assembly/fvlocalresidual.hh
+++ b/dumux/assembly/fvlocalresidual.hh
@@ -144,7 +144,7 @@ public:
* \param bcTypes The types of the boundary conditions for all boundary entities of an element
* \param elemFluxVarsCache The flux variable caches for the element stencil
*/
- DUNE_DEPRECATED_MSG("eval is deprecated because it doesn't allow to specify on which time level to evaluate. Use evalFluxSource, and evalStorage instead!")
+ DUNE_DEPRECATED_MSG("eval is deprecated because it doesn't allow to specify on which time level to evaluate. Use evalFluxAndSource, and evalStorage instead!")
ElementResidualVector eval(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -219,7 +219,7 @@ public:
* \param bcTypes The types of the boundary conditions for all boundary entities of an element
* \param elemFluxVarsCache The flux variable caches for the element stencil
*/
- DUNE_DEPRECATED_MSG("Use evalFluxSource instead!")
+ DUNE_DEPRECATED_MSG("Use evalFluxAndSource instead!")
ElementResidualVector eval(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -227,14 +227,14 @@ public:
const ElementBoundaryTypes &bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- return evalFluxSource(element, fvGeometry, curElemVolVars, elemFluxVarsCache, bcTypes);
+ return evalFluxAndSource(element, fvGeometry, curElemVolVars, elemFluxVarsCache, bcTypes);
}
- ElementResidualVector evalFluxSource(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes &bcTypes) const
+ ElementResidualVector evalFluxAndSource(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes &bcTypes) const
{
// initialize the residual vector for all scvs in this element
ElementResidualVector residual(fvGeometry.numScv());
diff --git a/dumux/assembly/jacobianpattern.hh b/dumux/assembly/jacobianpattern.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1847c116213aafef95713ded2028032301323ec7
--- /dev/null
+++ b/dumux/assembly/jacobianpattern.hh
@@ -0,0 +1,110 @@
+// -*- 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 Helper function to generate Jacobian pattern for different discretization methods
+ */
+#ifndef DUMUX_JACOBIAN_PATTERN_HH
+#define DUMUX_JACOBIAN_PATTERN_HH
+
+#include <type_traits>
+#include <dune/istl/matrixindexset.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for the box method
+ */
+template<bool isImplicit, class GridGeometry,
+ typename std::enable_if_t<(GridGeometry::discretizationMethod == DiscretizationMethods::Box), int> = 0>
+Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
+{
+ const auto numDofs = gridGeometry.numDofs();
+ Dune::MatrixIndexSet pattern;
+ pattern.resize(numDofs, numDofs);
+
+ // matrix pattern for implicit Jacobians
+ if (isImplicit)
+ {
+ static constexpr int dim = std::decay_t<decltype(gridGeometry.gridView())>::dimension;
+ for (const auto& element : elements(gridGeometry.gridView()))
+ {
+ for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
+ {
+ const auto globalI = gridGeometry.vertexMapper().subIndex(element, vIdx, dim);
+ for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
+ {
+ const auto globalJ = gridGeometry.vertexMapper().subIndex(element, vIdx2, dim);
+ pattern.add(globalI, globalJ);
+ pattern.add(globalJ, globalI);
+ }
+ }
+ }
+ }
+
+ // matrix pattern for explicit Jacobians -> diagonal matrix
+ else
+ {
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ pattern.add(globalI, globalI);
+ }
+
+ return pattern;
+}
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for cell-centered methods
+ */
+template<bool isImplicit, class GridGeometry,
+ typename std::enable_if_t<( (GridGeometry::discretizationMethod == DiscretizationMethods::CCTpfa)
+ || (GridGeometry::discretizationMethod == DiscretizationMethods::CCMpfa) ), int> = 0>
+Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
+{
+ const auto numDofs = gridGeometry.numDofs();
+ Dune::MatrixIndexSet pattern;
+ pattern.resize(numDofs, numDofs);
+
+ // matrix pattern for implicit Jacobians
+ if (isImplicit)
+ {
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ {
+ pattern.add(globalI, globalI);
+ for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
+ pattern.add(dataJ.globalJ, globalI);
+ }
+ }
+
+ // matrix pattern for explicit Jacobians -> diagonal matrix
+ else
+ {
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ pattern.add(globalI, globalI);
+ }
+
+ return pattern;
+}
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index 445c4932df251ead8bd512527abc8b4a9c52093c..8addc554443d474c48b85f7377b9512f369fb8fd 100644
--- a/dumux/common/CMakeLists.txt
+++ b/dumux/common/CMakeLists.txt
@@ -22,6 +22,7 @@ loggingparametertree.hh
math.hh
matrixvectorhelper.hh
optional.hh
+numericdifferentiation.hh
parameters.hh
pointsource.hh
properties.hh
diff --git a/dumux/assembly/numericdifferentiation.hh b/dumux/common/numericdifferentiation.hh
similarity index 91%
rename from dumux/assembly/numericdifferentiation.hh
rename to dumux/common/numericdifferentiation.hh
index cdf0ee3f1bb810e5a3d1d4f0196e2baad6c87332..fdc34ffa6b3c81e5f1ab4fd1e5e74018adc66531 100644
--- a/dumux/assembly/numericdifferentiation.hh
+++ b/dumux/common/numericdifferentiation.hh
@@ -27,7 +27,6 @@
#include <cmath>
#include <limits>
-#include <dumux/common/parameters.hh>
namespace Dumux {
@@ -61,8 +60,9 @@ public:
template<class Function, class Scalar, class FunctionEvalType>
static void partialDerivative(const Function& function, Scalar x0,
FunctionEvalType& derivative,
- const FunctionEvalType& fx0)
- { partialDerivative(function, x0, derivative, fx0, epsilon(x0)); }
+ const FunctionEvalType& fx0,
+ const int numericDifferenceMethod = 1)
+ { partialDerivative(function, x0, derivative, fx0, epsilon(x0), numericDifferenceMethod); }
/*!
* \brief Computes the derivative of a function with repect to a function parameter
@@ -71,15 +71,16 @@ public:
* \param derivative The partial derivative (output)
* \param fx0 The result of the function evaluated at x0
* \param eps The numeric epsilon used in the differentiation
+ * \param numericDifferenceMethod The numeric difference method
+ * (1: foward differences (default), 0: central differences, -1: backward differences)
*/
template<class Function, class Scalar, class FunctionEvalType>
static void partialDerivative(const Function& function, Scalar x0,
FunctionEvalType& derivative,
const FunctionEvalType& fx0,
- const Scalar eps)
+ const Scalar eps,
+ const int numericDifferenceMethod = 1)
{
- static const int numericDifferenceMethod = numDiffMethod();
-
Scalar delta = 0.0;
// we are using forward or central differences, i.e. we need to calculate f(x + \epsilon)
if (numericDifferenceMethod >= 0)
@@ -112,13 +113,6 @@ public:
// deflections between the two function evaluation
derivative /= delta;
}
-
- static int numDiffMethod()
- {
- static const int numericDifferenceMethod
- = getParam<int>("Assembly.NumericDifferenceMethod", 1);
- return numericDifferenceMethod;
- }
};
} // end namespace Dumux
diff --git a/dumux/common/parameters.hh b/dumux/common/parameters.hh
index d81b5fd1b209fe4528a3ddbc8183933c8561e70e..7ffe4836518eb806fc3818320609556b316a0f73 100644
--- a/dumux/common/parameters.hh
+++ b/dumux/common/parameters.hh
@@ -354,7 +354,9 @@ private:
params["Implicit.UpwindWeight"] = "1.0";
params["Implicit.EnablePartialReassemble"] = "false";
params["Implicit.EnableJacobianRecycling"] = "false";
- params["Implicit.NumericDifferenceMethod"] = "1";
+
+ // parameters in the assembly group
+ params["Assembly.NumericDifferenceMethod"] = "1";
// parameters in the linear solver group
params["LinearSolver.GMResRestart"] = "10";
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index e3b457bebea5fa0230ab8c757d40cc9884fdf58f..71ea05f596b0285cb446132659b122a2beb95e32 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -28,6 +28,7 @@
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dumux/discretization/methods.hh>
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
#include <dumux/discretization/box/fvelementgeometry.hh>
@@ -71,6 +72,9 @@ class BoxFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using GeometryHelper = BoxGeometryHelper<GridView, dim, SubControlVolume, SubControlVolumeFace>;
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::Box;
+
//! Constructor
BoxFVGridGeometry(const GridView gridView)
: ParentType(gridView) {}
@@ -263,6 +267,9 @@ class BoxFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::Box;
+
//! Constructor
BoxFVGridGeometry(const GridView gridView)
: ParentType(gridView) {}
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index dfba3592e209700e45d4642e397f517b5433f180..248c0129d9739c625ab49c728ab6f79c0414d3da 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -32,6 +32,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/discretization/methods.hh>
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/connectivitymap.hh>
@@ -87,6 +88,9 @@ class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCMpfa;
+
using SecondaryIvIndicatorType = std::function<bool(const Element&, const Intersection&, bool)>;
//! Constructor without indicator function for secondary interaction volumes
@@ -429,6 +433,9 @@ class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCMpfa;
+
using SecondaryIvIndicator = std::function<bool(const Element&, const Intersection&, bool)>;
//! Constructor without indicator function for secondary interaction volumes
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index 3a20ad88449e6cfb1d6809240a930c123f73a2a1..ce805e818a37a5fa7f8c5395f284b6fd4d15a26f 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -81,7 +81,7 @@ private:
public:
// Per default, we allow for 8 neighbors on network/surface grids
- static const std::size_t value = dim < dimWorld ? 9 : 2;
+ static constexpr std::size_t value = dim < dimWorld ? 9 : 2;
};
//! Set the index set type used on the dual grid nodes
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index f71a30ea60989698801ff931a22dad27b616e5a4..9d38e45f97d540192eec3a8dcf194fdddbce4ba2 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
@@ -28,6 +28,7 @@
#include <dune/common/version.hh>
+#include <dumux/discretization/methods.hh>
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/connectivitymap.hh>
@@ -74,6 +75,9 @@ class CCTpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCTpfa;
+
//! Constructor
CCTpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
@@ -331,6 +335,9 @@ class CCTpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
public:
+ //! export discretization method
+ static constexpr DiscretizationMethods discretizationMethod = DiscretizationMethods::CCTpfa;
+
//! Constructor
CCTpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
diff --git a/dumux/discretization/cellcentered/tpfa/properties.hh b/dumux/discretization/cellcentered/tpfa/properties.hh
index 3710cf0fa7f339df01da81fc79ae8b570dac3f81..c9262c4b141561ec9b65bdec6de69c0019461f51 100644
--- a/dumux/discretization/cellcentered/tpfa/properties.hh
+++ b/dumux/discretization/cellcentered/tpfa/properties.hh
@@ -92,7 +92,7 @@ private:
public:
// Per default, we allow for 8 neighbors on network/surface grids
- static const std::size_t value = dim < dimWorld ? 9 : 2;
+ static constexpr std::size_t value = dim < dimWorld ? 9 : 2;
};
//! The sub control volume
diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
index dfdc173f082d9f9aced2fd6b416775cca63f27c4..1a0f1177ca8b88f7634b6c7458b54117124d8759 100644
--- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
@@ -44,7 +44,6 @@ class OnePIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
index 60ee8808597d7ca8c116c615a4035266cf236e04..73affdffac48d9ffdedaa9c69635e0dffd5f8760 100644
--- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
@@ -47,7 +47,7 @@ class TwoPIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -188,10 +188,10 @@ public:
const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
insideScv,
- ElementResidualVector({insideVolVars.priVars()}));
+ ElementSolutionVector({insideVolVars.priVars()}));
const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement,
outsideScv,
- ElementResidualVector({outsideVolVars.priVars()}));
+ ElementSolutionVector({outsideVolVars.priVars()}));
// get references to the two participating derivative matrices
auto& dI_dI = derivativeMatrices[insideScvIdx];
@@ -303,7 +303,7 @@ public:
const auto& outsideVolVars = curElemVolVars[outsideScv];
// we need the element solution for the material parameters
- ElementResidualVector elemSol(fvGeometry.numScv());
+ ElementSolutionVector elemSol(fvGeometry.numScv());
for (const auto& scv : scvs(fvGeometry)) elemSol[scv.indexInElement()] = curElemVolVars[scv].priVars();
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, elemSol);
const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(element, outsideScv, elemSol);
@@ -449,7 +449,7 @@ public:
const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
insideScv,
- ElementResidualVector({insideVolVars.priVars()}));
+ ElementSolutionVector({insideVolVars.priVars()}));
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
diff --git a/dumux/porousmediumflow/mpnc/localresidual.hh b/dumux/porousmediumflow/mpnc/localresidual.hh
index d982de82f85d55f6c8ae25fa820ff196f0967036..737b507b71e7c23231995d98e480a2b8192ab4a7 100644
--- a/dumux/porousmediumflow/mpnc/localresidual.hh
+++ b/dumux/porousmediumflow/mpnc/localresidual.hh
@@ -60,13 +60,13 @@ public:
using typename ParentType::ElementResidualVector;
- ElementResidualVector evalFluxSource(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementFluxVariablesCache& elemFluxVarsCache,
- const ElementBoundaryTypes &bcTypes) const
+ ElementResidualVector evalFluxAndSource(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes &bcTypes) const
{
- ElementResidualVector residual = ParentType::evalFluxSource(element, fvGeometry, elemVolVars, elemFluxVarsCache, bcTypes);
+ ElementResidualVector residual = ParentType::evalFluxAndSource(element, fvGeometry, elemVolVars, elemFluxVarsCache, bcTypes);
for (auto&& scv : scvs(fvGeometry))
{
diff --git a/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input
index ef445b1a852b728898f0ae79c34506e0d9d8c446..d910fd45112ccb392383ce3737ace10045a15e58 100644
--- a/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input
+++ b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input
@@ -9,7 +9,7 @@ Cells = 40 4
[Problem]
Name = infiltration
-[Implicit]
+[Assembly]
NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences
[Newton]
diff --git a/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input b/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input
index c1848acdb19e54f7bc697e71ff5daf79569d8ea3..80efdf4037028cba0518be7849c1d6a4d7d0d54f 100644
--- a/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input
+++ b/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input
@@ -10,5 +10,5 @@ Cells = 1 120
[Problem]
Name = columnxylol # name passed to the output routines
-[Implicit]
+[Assembly]
NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences
diff --git a/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input b/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input
index 1738b0ae022d312f8d7a7c9cc4c9e6077951bd03..e48af1cd0e96d072b820eb0b3925fb5cf1ee53dc 100644
--- a/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input
+++ b/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input
@@ -11,7 +11,7 @@ Cells = 15 8
[Problem]
Name = kuevette # name passed to the output routines
-[Implicit]
+[Assembly]
NumericDifferenceMethod = 0 # use central differences (backward -1, forward +1)
[Newton]