From 32e0397dac724974ad2a50afa00e4e33b785eff2 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sat, 13 Jan 2018 17:14:03 +0100
Subject: [PATCH] [assembly] Hand in numeric diff method to the partial
derivative function
---
dumux/assembly/boxlocalassembler.hh | 7 +++++--
dumux/assembly/cclocalassembler.hh | 8 ++++++--
dumux/assembly/fvlocalassemblerbase.hh | 1 -
dumux/assembly/numericdifferentiation.hh | 20 +++++++------------
dumux/common/parameters.hh | 4 +++-
.../3p3c/implicit/test_3p3c_fv.input | 2 +-
.../3p3c/implicit/test_columnxylol_fv.input | 2 +-
.../3p3c/implicit/test_kuvette_fv.input | 2 +-
8 files changed, 24 insertions(+), 22 deletions(-)
diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index 86644925ff..4679e60e35 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -32,6 +32,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/assembly/fvlocalassemblerbase.hh>
+#include <dumux/assembly/numericdifferentiation.hh>
namespace Dumux {
@@ -270,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))
@@ -378,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++)
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 00d9c51e2d..b7d5d3ea74 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.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/fvlocalassemblerbase.hh b/dumux/assembly/fvlocalassemblerbase.hh
index 2659946d05..7758e8af9b 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 {
diff --git a/dumux/assembly/numericdifferentiation.hh b/dumux/assembly/numericdifferentiation.hh
index cdf0ee3f1b..fdc34ffa6b 100644
--- a/dumux/assembly/numericdifferentiation.hh
+++ b/dumux/assembly/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 d81b5fd1b2..7ffe483651 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/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input
index ef445b1a85..d910fd4511 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 c1848acdb1..80efdf4037 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 1738b0ae02..e48af1cd0e 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]
--
GitLab