Commit bd209117 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

box models: Allow to chose forward or backward differences instead of only...

box models: Allow to chose forward or backward differences instead of only central ones in the local jacobian

also change the 2p test problem to use forward differences

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4651 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 547ac0ee
......@@ -56,9 +56,13 @@ private:
Red = JacobianAssembler::Red,
Yellow = JacobianAssembler::Yellow,
Green = JacobianAssembler::Green
Green = JacobianAssembler::Green,
numDiffMethod = GET_PROP_VALUE(TypeTag,
PTAG(NumericDifferenceMethod))
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GridView::Grid::ctype CoordScalar;
......@@ -152,6 +156,17 @@ public:
elem_(),
fvElemGeom_,
false /* isOldSol? */);
if (numDiffMethod != 0) {
// for forward and backward differences we need the local
// residual at the evaluation point. For efficiency we
// only calculate this once.
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
evalPointResid_ = localResidual().residual();
};
// calculate the local jacobian matrix
ElementSolutionVector partialDeriv(numVertices);
......@@ -272,45 +287,60 @@ protected:
curVolVars_[scvIdx].setEvalPoint(&origVolVars);
Scalar eps = asImp_().numericEpsilon_(scvIdx, pvIdx);
Scalar delta = 0;
if (numDiffMethod >= 0) { // not backward differences
// deflect primary variables
priVars[pvIdx] += eps;
delta += eps;
// calculate the residual
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
// store the residual
dest = localResidual().residual();
}
else {
// backward differences
dest = evalPointResid_;
}
if (numDiffMethod <= 0) { // not forward differences
// deflect the primary variables
priVars[pvIdx] -= delta + eps;
delta += eps;
// calculate residual again
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
dest -= localResidual().residual();
}
else {
// forward differences
dest -= evalPointResid_;
}
// deflect primary variables
priVars[pvIdx] += eps;
// calculate the residual
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
// store the residual
dest = localResidual().residual();
// deflect the primary variables
priVars[pvIdx] -= 2*eps;
// calculate residual again
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
// central differences
dest -= localResidual().residual();
dest /= 2*eps;
dest /= delta;
// restore the orignal state of the element's secondary
// variables
......@@ -376,6 +406,7 @@ protected:
// levels
ElementVolumeVariables prevVolVars_;
ElementVolumeVariables curVolVars_;
ElementSolutionVector evalPointResid_; // residual at evaluation point
LocalResidual localResidual_;
LocalBlockMatrix A_;
......
......@@ -99,6 +99,15 @@ NEW_PROP_TAG(EnablePartialReassemble);
//! newton iterations to achive larger time step sizes
NEW_PROP_TAG(EnableTimeStepRampUp);
/*!
* \brief Specify which kind of method should be used to numerically
* calculate the partial derivatives of the residual.
*
* -1 means backward differences, 0 means central differences, 1 means
* forward differences. By default we use central differences.
*/
NEW_PROP_TAG(NumericDifferenceMethod);
// mappers from local to global indices
//! maper for vertices
......
......@@ -248,6 +248,9 @@ public:
SET_PROP(BoxModel, LocalOperator)
{ typedef typename GET_PROP(TypeTag, PTAG(GridOperatorSpace))::LocalOperator type; };
//! use central differences to calculate the jacobian by default
SET_INT_PROP(BoxModel, NumericDifferenceMethod, 0);
// disable jacobian matrix recycling by default
SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
// disable partial reassembling by default
......
......@@ -120,6 +120,9 @@ SET_BOOL_PROP(LensProblem, EnableJacobianRecycling, false);
// Write the solutions of individual newton iterations?
SET_BOOL_PROP(LensProblem, NewtonWriteConvergence, false);
// Use forward differences instead of central differences
SET_INT_PROP(LensProblem, NumericDifferenceMethod, +1);
// Enable gravity
SET_BOOL_PROP(LensProblem, EnableGravity, true);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment