Commit d6c8b766 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

implicit branch: unify Box- and CC-LocalJacobian to ImplicitLocalJacobian....

implicit branch: unify Box- and CC-LocalJacobian to ImplicitLocalJacobian. Delete cclocaljacobian.hh. Deprecate BoxLocalJacobian in favor of ImplicitLocalJacobian. Adapt includes and names.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9749 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 99fa9709
#warning This file is deprecated. Include dumux/implicit/box/boxlocaljacobian.hh instead.
// -*- 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
* \brief Caculates the Jacobian of the local residual for box models
*/
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH
#include <dumux/implicit/box/boxlocaljacobian.hh>
#include <dumux/implicit/common/implicitlocaljacobian.hh>
namespace Dumux
{
/*!
* \ingroup BoxModel
* \ingroup BoxLocalJacobian
* \brief Calculates the Jacobian of the local residual for box models
*
* The default behavior is to use numeric differentiation, i.e.
* forward or backward differences (2nd order), or central
* differences (3rd order). The method used is determined by the
* "NumericDifferenceMethod" property:
*
* - if the value of this property is smaller than 0, backward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
* \f]
*
* - if the value of this property is 0, central
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
* \f]
*
* - if the value of this property is larger than 0, forward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
* \f]
*
* Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
* is the value of a sub-control volume's primary variable at the
* evaluation point and \f$\epsilon\f$ is a small value larger than 0.
*/
template<class TypeTag>
class BoxLocalJacobian : public ImplicitLocalJacobian<TypeTag>
{
typedef ImplicitLocalJacobian<TypeTag> ParentType;
// copying a local jacobian is not a good idea
BoxLocalJacobian(const BoxLocalJacobian &);
public:
DUNE_DEPRECATED_MSG("Use ImplicitLocalJacobian from dumux/implicit/common/implicitlocaljacobian.hh.")
BoxLocalJacobian() : ParentType()
{ }
};
}
#endif
......@@ -25,7 +25,7 @@
#include <dune/istl/matrix.hh>
#include "stokesproperties.hh"
#include <dumux/implicit/box/boxlocaljacobian.hh>
#include <dumux/implicit/common/implicitlocaljacobian.hh>
namespace Dumux
{
......@@ -38,13 +38,13 @@ namespace Dumux
* The momentum balance equation uses larger epsilons than the rest.
*/
template<class TypeTag>
class StokesLocalJacobian : public BoxLocalJacobian<TypeTag>
class StokesLocalJacobian : public ImplicitLocalJacobian<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
public:
//! \copydoc BoxLocalJacobian::numericEpsilon()
//! \copydoc ImplicitLocalJacobian::numericEpsilon()
Scalar numericEpsilon(const int scvIdx,
const int pvIdx) const
{
......
......@@ -587,6 +587,7 @@ public:
int numEdges; //!< number of edges
int numFaces; //!< number of faces (0 in < 3D)
int numSCV; //!< number of subcontrol volumes
int numNeighbors; //!< needed for compatibility with cc models
int numFAP; //!< number of flux approximation points
std::vector<ElementPointer> neighbors; //!< needed for compatibility with cc models
......@@ -614,6 +615,7 @@ public:
numEdges = referenceElement.size(dim-1);
numFaces = (dim<3)?0:referenceElement.size(1);
numSCV = numVertices;
numNeighbors = 0;
bool useTwoPointFlux
= GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, UseTwoPointFlux);
......
// -*- 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
* \brief Caculates the Jacobian of the local residual for box models
*/
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH
#include <dune/istl/matrix.hh>
#include <dumux/common/math.hh>
#include "boxelementboundarytypes.hh"
namespace Dumux
{
/*!
* \ingroup BoxModel
* \ingroup BoxLocalJacobian
* \brief Calculates the Jacobian of the local residual for box models
*
* The default behavior is to use numeric differentiation, i.e.
* forward or backward differences (2nd order), or central
* differences (3rd order). The method used is determined by the
* "NumericDifferenceMethod" property:
*
* - if the value of this property is smaller than 0, backward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
* \f]
*
* - if the value of this property is 0, central
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
* \f]
*
* - if the value of this property is larger than 0, forward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
* \f]
*
* Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
* is the value of a sub-control volume's primary variable at the
* evaluation point and \f$\epsilon\f$ is a small value larger than 0.
*/
template<class TypeTag>
class BoxLocalJacobian
{
private:
typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
enum {
numEq = GET_PROP_VALUE(TypeTag, NumEq),
dim = GridView::dimension,
Green = JacobianAssembler::Green
};
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
// copying a local jacobian is not a good idea
BoxLocalJacobian(const BoxLocalJacobian &);
public:
BoxLocalJacobian()
{
numericDifferenceMethod_ = GET_PARAM_FROM_GROUP(TypeTag, int, Implicit, NumericDifferenceMethod);
Valgrind::SetUndefined(problemPtr_);
}
/*!
* \brief Initialize the local Jacobian object.
*
* At this point we can assume that everything has been allocated,
* although some objects may not yet be completely initialized.
*
* \param prob The problem which we want to simulate.
*/
void init(Problem &prob)
{
problemPtr_ = &prob;
localResidual_.init(prob);
// assume quadrilinears as elements with most vertices
A_.setSize(2<<dim, 2<<dim);
storageJacobian_.resize(2<<dim);
}
/*!
* \brief Assemble an element's local Jacobian matrix of the
* defect.
*
* \param element The DUNE Codim<0> entity which we look at.
*/
void assemble(const Element &element)
{
// set the current grid element and update the element's
// finite volume geometry
elemPtr_ = &element;
fvElemGeom_.update(gridView_(), element);
reset_();
bcTypes_.update(problem_(), element_(), fvElemGeom_);
// this is pretty much a HACK because the internal state of
// the problem is not supposed to be changed during the
// evaluation of the residual. (Reasons: It is a violation of
// abstraction, makes everything more prone to errors and is
// not thread save.) The real solution are context objects!
problem_().updateCouplingParams(element_());
// set the hints for the volume variables
model_().setHints(element, prevVolVars_, curVolVars_);
// update the secondary variables for the element at the last
// and the current time levels
prevVolVars_.update(problem_(),
element_(),
fvElemGeom_,
true /* isOldSol? */);
curVolVars_.update(problem_(),
element_(),
fvElemGeom_,
false /* isOldSol? */);
// update the hints of the model
model_().updateCurHints(element, curVolVars_);
// calculate the local residual
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
residual_ = localResidual().residual();
storageTerm_ = localResidual().storageTerm();
model_().updatePVWeights(element_(), curVolVars_);
// calculate the local jacobian matrix
int numVertices = fvElemGeom_.numVertices;
ElementSolutionVector partialDeriv(numVertices);
PrimaryVariables storageDeriv;
for (int j = 0; j < numVertices; j++) {
for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
asImp_().evalPartialDerivative_(partialDeriv,
storageDeriv,
j,
pvIdx);
// update the local stiffness matrix with the current partial
// derivatives
updateLocalJacobian_(j,
pvIdx,
partialDeriv,
storageDeriv);
}
}
}
/*!
* \brief Returns a reference to the object which calculates the
* local residual.
*/
const LocalResidual &localResidual() const
{ return localResidual_; }
/*!
* \brief Returns a reference to the object which calculates the
* local residual.
*/
LocalResidual &localResidual()
{ return localResidual_; }
/*!
* \brief Returns the Jacobian of the equations at vertex i to the
* primary variables at vertex j.
*
* \param i The local vertex (or sub-control volume) index on which
* the equations are defined
* \param j The local vertex (or sub-control volume) index which holds
* primary variables
*/
const MatrixBlock &mat(const int i, const int j) const
{ return A_[i][j]; }
/*!
* \brief Returns the Jacobian of the storage term at vertex i.
*
* \param i The local vertex (or sub-control volume) index
*/
const MatrixBlock &storageJacobian(const int i) const
{ return storageJacobian_[i]; }
/*!
* \brief Returns the residual of the equations at vertex i.
*
* \param i The local vertex (or sub-control volume) index on which
* the equations are defined
*/
const PrimaryVariables &residual(const int i) const
{ return residual_[i]; }
/*!
* \brief Returns the storage term for vertex i.
*
* \param i The local vertex (or sub-control volume) index on which
* the equations are defined
*/
const PrimaryVariables &storageTerm(const int i) const
{ return storageTerm_[i]; }
/*!
* \brief Returns the epsilon value which is added and removed
* from the current solution.
*
* \param scvIdx The local index of the element's vertex for
* which the local derivative ought to be calculated.
* \param pvIdx The index of the primary variable which gets varied
*/
Scalar numericEpsilon(const int scvIdx,
const int pvIdx) 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...
Scalar priVar = this->curVolVars_[scvIdx].priVar(pvIdx);
return baseEps*(std::abs(priVar) + 1.0);
}
protected:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
/*!
* \brief Returns a reference to the problem.
*/
const Problem &problem_() const
{
Valgrind::CheckDefined(problemPtr_);
return *problemPtr_;
};
/*!
* \brief Returns a reference to the grid view.
*/
const GridView &gridView_() const
{ return problem_().gridView(); }
/*!
* \brief Returns a reference to the element.
*/
const Element &element_() const
{
Valgrind::CheckDefined(elemPtr_);
return *elemPtr_;
};
/*!
* \brief Returns a reference to the model.
*/
const Model &model_() const
{ return problem_().model(); };
/*!
* \brief Returns a reference to the jacobian assembler.
*/
const JacobianAssembler &jacAsm_() const
{ return model_().jacobianAssembler(); }
/*!
* \brief Returns a reference to the vertex mapper.
*/
const VertexMapper &vertexMapper_() const
{ return problem_().vertexMapper(); };
/*!
* \brief Reset the local jacobian matrix to 0
*/
void reset_()
{
int n = element_().template count<dim>();
for (int i = 0; i < n; ++ i) {
storageJacobian_[i] = 0.0;
for (int j = 0; j < n; ++ j) {
A_[i][j] = 0.0;
}
}
}
/*!
* \brief Compute the partial derivatives to a primary variable at
* an degree of freedom.
*
* This method can be overwritten by the implementation if a
* better scheme than numerical differentiation is available.
*
* The default implementation of this method uses numeric
* differentiation, i.e. forward or backward differences (2nd
* order), or central differences (3rd order). The method used is
* determined by the "NumericDifferenceMethod" property:
*
* - if the value of this property is smaller than 0, backward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
* \f]
*
* - if the value of this property is 0, central
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
* \f]
*
* - if the value of this property is larger than 0, forward
* differences are used, i.e.:
* \f[
\frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
* \f]
*
* Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
* is the value of a sub-control volume's primary variable at the
* evaluation point and \f$\epsilon\f$ is a small value larger than 0.
*
* \param partialDeriv The vector storing the partial derivatives of all
* equations
* \param storageDeriv the mass matrix contributions
* \param scvIdx The sub-control volume index of the current
* finite element for which the partial derivative
* ought to be calculated
* \param pvIdx The index of the primary variable at the scvIdx'
* sub-control volume of the current finite element
* for which the partial derivative ought to be
* calculated
*/
void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
PrimaryVariables &storageDeriv,
const int scvIdx,
const int pvIdx)
{
int globalIdx = vertexMapper_().map(element_(), scvIdx, dim);
PrimaryVariables priVars(model_().curSol()[globalIdx]);
VolumeVariables origVolVars(curVolVars_[scvIdx]);
curVolVars_[scvIdx].setEvalPoint(&origVolVars);
Scalar eps = asImp_().numericEpsilon(scvIdx, pvIdx);
Scalar delta = 0;
if (numericDifferenceMethod_ >= 0) {
// we are not using backward differences, i.e. we need to
// calculate f(x + \epsilon)
// deflect primary variables
priVars[pvIdx] += eps;
delta += eps;
// calculate the residual
curVolVars_[scvIdx].update(priVars,
problem_(),
element_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
// store the residual and the storage term
partialDeriv = localResidual().residual();
storageDeriv = localResidual().storageTerm()[scvIdx];
}
else {
// we are using backward differences, i.e. we don't need
// to calculate f(x + \epsilon) and we can recycle the
// (already calculated) residual f(x)
partialDeriv = residual_;
storageDeriv = storageTerm_[scvIdx];
}
if (numericDifferenceMethod_ <= 0) {
// we are not using forward differences, i.e. we don't
// need to calculate f(x - \epsilon)
// deflect the primary variables
priVars[pvIdx] -= delta + eps;
delta += eps;
// calculate residual again
curVolVars_[scvIdx].update(priVars,
problem_(),
element_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
partialDeriv -= localResidual().residual();
storageDeriv -= localResidual().storageTerm()[scvIdx];
}
else {
// we are using forward differences, i.e. we don't need to
// calculate f(x - \epsilon) and we can recycle the
// (already calculated) residual f(x)
partialDeriv -= residual_;
storageDeriv -= storageTerm_[scvIdx];
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
partialDeriv /= delta;
storageDeriv /= delta;
// restore the original state of the element's volume variables
curVolVars_[scvIdx] = origVolVars;
#if HAVE_VALGRIND
for (unsigned i = 0; i < partialDeriv.size(); ++i)
Valgrind::CheckDefined(partialDeriv[i]);
#endif
}
/*!
* \brief Updates the current local Jacobian matrix with the
* partial derivatives of all equations in regard to the
* primary variable 'pvIdx' at vertex 'scvIdx' .
*/
void updateLocalJacobian_(const int scvIdx,
const int pvIdx,
const ElementSolutionVector &partialDeriv,
const PrimaryVariables &storageDeriv)
{
// store the derivative of the storage term
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
storageJacobian_[scvIdx][eqIdx][pvIdx] = storageDeriv[eqIdx];
}