From 6547ce08dc9f629636e31f869476c8544b1303ee Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 16 Feb 2017 19:43:09 +0100
Subject: [PATCH] [staggeredGrid][localJacobian] Use function to determine
numerical eps
* Do not use fixed values anymore
---
dumux/implicit/staggered/localjacobian.hh | 30 ++++++++++++++++-------
1 file changed, 21 insertions(+), 9 deletions(-)
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 0d6dca42ae..027747031f 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -266,11 +266,10 @@ private:
for(auto pvIdx : PriVarIndices(cellCenterIdx))
{
- const Scalar eps = 1e-4; // TODO: do properly
-
PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().curSol()[cellCenterIdx][globalJ]),
FacePrimaryVariables(0.0));
+ const Scalar eps = numericEpsilon(priVars[pvIdx]);
priVars[pvIdx] += eps;
ElementSolutionVector elemSol{std::move(priVars)};
curVolVars.update(elemSol, this->problem_(), elementJ, scvJ);
@@ -318,9 +317,8 @@ private:
for(auto pvIdx : PriVarIndices(faceIdx))
{
PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
- const Scalar eps = 1e-4; // TODO: do properly
+ const Scalar eps = numericEpsilon(priVars[pvIdx]);
priVars[pvIdx] += eps;
-
curFaceVars.update(priVars[faceIdx]);
this->localResidual().evalCellCenter(element, fvGeometry, scvI,
@@ -370,11 +368,10 @@ private:
for(auto pvIdx : PriVarIndices(cellCenterIdx))
{
- const Scalar eps = 1e-4; // TODO: do properly
-
PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().curSol()[cellCenterIdx][globalJ]),
FacePrimaryVariables(0.0));
+ const Scalar eps = numericEpsilon(priVars[pvIdx]);
priVars[pvIdx] += eps;
ElementSolutionVector elemSol{std::move(priVars)};
curVolVars.update(elemSol, this->problem_(), elementJ, scvJ);
@@ -424,11 +421,10 @@ private:
for(auto pvIdx : PriVarIndices(faceIdx))
{
- const Scalar eps = 1e-10; // TODO: do properly
-
PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
- priVars[pvIdx] += eps;
+ const Scalar eps = numericEpsilon(priVars[pvIdx]);
+ priVars[pvIdx] += eps;
curFaceVars.update(priVars[faceIdx]);
this->localResidual().evalFace(element, fvGeometry, scvf,
@@ -491,6 +487,22 @@ private:
const JacobianAssembler &jacAsm_() const
{ return model_().jacobianAssembler(); }
+ Scalar numericEpsilon(const Scalar priVar) const
+ {
+ // define the base epsilon as the geometric mean of 1 and the
+ // resolution of the scalar type. E.g. for standard 64 bit
+ // floating point values, the resolution is about 10^-16 and
+ // the base epsilon is thus approximately 10^-8.
+ /*
+ static const Scalar baseEps
+ = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
+ */
+ static const Scalar baseEps = 1e-10;
+ assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
+ // the epsilon value used for the numeric differentiation is
+ // now scaled by the absolute value of the primary variable...
+ return baseEps*(std::abs(priVar) + 1.0);
+ }
/*!
* \brief Updates the current global Jacobian matrix with the
* partial derivatives of all equations in regard to the
--
GitLab