diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index f9326416f21caa5c19f970f35aedbed3d920b3a7..ae1da06f025141e543e0d76e4e6d90f3fc5c9bc7 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -1,3 +1,77 @@
-#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
diff --git a/dumux/boxmodels/common/porousmediaboxproblem.hh b/dumux/boxmodels/common/porousmediaboxproblem.hh
index 63cd95cb08977333ad9fdab9814b3b4b3f695f80..c86c92e6d4a1fcd3c0f5c00a362823a30dbabe59 100644
--- a/dumux/boxmodels/common/porousmediaboxproblem.hh
+++ b/dumux/boxmodels/common/porousmediaboxproblem.hh
@@ -1,7 +1,6 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * 2012 by Bernd Flemisch *
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
diff --git a/dumux/freeflow/stokes/stokeslocaljacobian.hh b/dumux/freeflow/stokes/stokeslocaljacobian.hh
index 9dee8d50b57c54f969b6c4709d4280468033e448..074443ef7c815683daec28b1ed10c7daf6237b22 100644
--- a/dumux/freeflow/stokes/stokeslocaljacobian.hh
+++ b/dumux/freeflow/stokes/stokeslocaljacobian.hh
@@ -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
{
diff --git a/dumux/implicit/box/boxfvelementgeometry.hh b/dumux/implicit/box/boxfvelementgeometry.hh
index 2a8dcf4c97767c2b14ca59d88a5527e7f3f180ca..8f0811b4401dede22e891806ae014892943397ae 100644
--- a/dumux/implicit/box/boxfvelementgeometry.hh
+++ b/dumux/implicit/box/boxfvelementgeometry.hh
@@ -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);
diff --git a/dumux/implicit/box/boxlocaljacobian.hh b/dumux/implicit/box/boxlocaljacobian.hh
deleted file mode 100644
index 7fc5140d2545e0ba467331e73bb39df09999ead6..0000000000000000000000000000000000000000
--- a/dumux/implicit/box/boxlocaljacobian.hh
+++ /dev/null
@@ -1,534 +0,0 @@
-// -*- 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];
- }
-
- for (int i = 0; i < fvElemGeom_.numVertices; i++)
- {
- // Green vertices are not to be changed!
- if (jacAsm_().vertexColor(element_(), i) != Green) {
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // A[i][scvIdx][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at vertex 'i'
- // depending on the primary variable 'pvIdx' at vertex
- // 'scvIdx'.
- this->A_[i][scvIdx][eqIdx][pvIdx] = partialDeriv[i][eqIdx];
- Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
- }
- }
- }
- }
-
- const Element *elemPtr_;
- FVElementGeometry fvElemGeom_;
-
- ElementBoundaryTypes bcTypes_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- // secondary variables at the previous and at the current time
- // levels
- ElementVolumeVariables prevVolVars_;
- ElementVolumeVariables curVolVars_;
-
- LocalResidual localResidual_;
-
- LocalBlockMatrix A_;
- std::vector<MatrixBlock> storageJacobian_;
-
- ElementSolutionVector residual_;
- ElementSolutionVector storageTerm_;
-
- int numericDifferenceMethod_;
-};
-}
-
-#endif
diff --git a/dumux/implicit/box/boxlocalresidual.hh b/dumux/implicit/box/boxlocalresidual.hh
index 8c02754133eec4a5cde1c8afdbff9bbc3d93ece1..3b8fb1a63a8578a11a886719f474c48ca09f53e2 100644
--- a/dumux/implicit/box/boxlocalresidual.hh
+++ b/dumux/implicit/box/boxlocalresidual.hh
@@ -27,6 +27,7 @@
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/common/implicitlocalresidual.hh>
#include "boxproperties.hh"
@@ -35,18 +36,16 @@ namespace Dumux
/*!
* \ingroup BoxModel
* \ingroup BoxLocalResidual
- * \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * \brief Element-wise calculation of the residual for models
+ * based on the fully implicit box scheme.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class BoxLocalResidual
+class BoxLocalResidual : public ImplicitLocalResidual<TypeTag>
{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef ImplicitLocalResidual<TypeTag> ParentType;
+ friend class ImplicitLocalResidual<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -63,273 +62,18 @@ private:
typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
- 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, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
BoxLocalResidual(const BoxLocalResidual &);
public:
- BoxLocalResidual()
+ BoxLocalResidual() : ParentType()
{ }
- ~BoxLocalResidual()
- { }
-
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarily initialized completely.
- *
- * \param problem The representation of the physical problem to be
- * solved.
- */
- void init(Problem &problem)
- { problemPtr_ = &problem; }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- */
- void eval(const Element &element)
- {
- FVElementGeometry fvGeometry;
-
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementVolumeVariables volVarsPrev, volVarsCur;
- // update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
-
- volVarsPrev.update(problem_(),
- element,
- fvGeometry_(),
- true /* oldSol? */);
- volVarsCur.update(problem_(),
- element,
- fvGeometry_(),
- false /* oldSol? */);
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- // 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);
-
- asImp_().eval(element, fvGeometry_(), volVarsPrev, volVarsCur, bcTypes);
- }
-
- /*!
- * \brief Compute the storage term for the current solution.
- *
- * This can be used to figure out how much of each conservation
- * quantity is inside the element.
- *
- * \param element The DUNE Codim<0> entity for which the storage
- * term ought to be calculated
- */
- void evalStorage(const Element &element)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
-
- // no previous volume variables!
- prevVolVarsPtr_ = 0;
-
- ElementVolumeVariables volVars;
-
- // update the hints
- model_().setHints(element, volVars);
-
- // calculate volume current variables
- volVars.update(problem_(), element, fvGeometry_(), false);
- curVolVarsPtr_ = &volVars;
-
- asImp_().evalStorage_();
- }
-
- /*!
- * \brief Compute the flux term for the current solution.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element
- */
- void evalFluxes(const Element &element,
- const ElementVolumeVariables &curVolVars)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- residual_.resize(fvGeometry_().numVertices);
- residual_ = 0;
-
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = 0;
- curVolVarsPtr_ = &curVolVars;
- asImp_().evalFluxes_();
- }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvGeometry The finite-volume geometry of the element
- * \param prevVolVars The volume averaged variables for all
- * sub-control volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- void eval(const Element &element,
- const FVElementGeometry &fvGeometry,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
- const ElementBoundaryTypes &bcTypes)
- {
- Valgrind::CheckDefined(prevVolVars);
- Valgrind::CheckDefined(curVolVars);
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry.numVertices; i++) {
- prevVolVars[i].checkDefined();
- curVolVars[i].checkDefined();
- }
-#endif // HAVE_VALGRIND
-
- elemPtr_ = &element;
- fvElemGeomPtr_ = &fvGeometry;
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = &prevVolVars;
- curVolVarsPtr_ = &curVolVars;
-
- // resize the vectors for all terms
- int numVerts = fvGeometry_().numVertices;
- residual_.resize(numVerts);
- storageTerm_.resize(numVerts);
-
- residual_ = 0.0;
- storageTerm_ = 0.0;
-
- asImp_().evalFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().evalVolumeTerms_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++) {
- Valgrind::CheckDefined(residual_[i]);
- }
-#endif // HAVE_VALGRIND
-
- // evaluate the boundary conditions
- asImp_().evalBoundary_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
- }
-
- /*!
- * \brief Returns the local residual for all sub-control
- * volumes of the element.
- */
- const ElementSolutionVector &residual() const
- { return residual_; }
-
- /*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
- *
- * \param scvIdx The local index of the sub-control volume
- * (i.e. the element's local vertex index)
- */
- const PrimaryVariables &residual(const int scvIdx) const
- { return residual_[scvIdx]; }
-
- /*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
- */
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
-
- /*!
- * \brief Returns the storage term for a given sub-control volumes
- * of the element.
- */
- const PrimaryVariables &storageTerm(const int scvIdx) const
- { return storageTerm_[scvIdx]; }
-
protected:
- Implementation &asImp_()
- {
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
-
- const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
- }
-
- /*!
- * \brief Evaluate the boundary conditions
- * of the current element.
- */
- void evalBoundary_()
- {
- if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
- asImp_().evalBoundaryFluxes_();
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
- if (bcTypes_().hasDirichlet())
- asImp_().evalDirichlet_();
- }
-
/*!
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
@@ -337,14 +81,14 @@ protected:
void evalDirichlet_()
{
PrimaryVariables dirichletValues(0);
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) {
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+ for (int scvIdx = 0; scvIdx < this->fvGeometry_().numVertices; ++scvIdx) {
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
if (bcTypes.hasDirichlet()) {
// ask the problem for the dirichlet values
- const VertexPointer vPtr = element_().template subEntity<dim>(scvIdx);
+ const VertexPointer vPtr = this->element_().template subEntity<dim>(scvIdx);
Valgrind::SetUndefined(dirichletValues);
- asImp_().problem_().dirichlet(dirichletValues, *vPtr);
+ this->asImp_().problem_().dirichlet(dirichletValues, *vPtr);
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
@@ -353,10 +97,10 @@ protected:
assert(0 <= pvIdx && pvIdx < numEq);
Valgrind::CheckDefined(dirichletValues[pvIdx]);
- residual_[scvIdx][eqIdx] =
- curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
+ this->residual_[scvIdx][eqIdx] =
+ this->curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
- storageTerm_[scvIdx][eqIdx] = 0.0;
+ this->storageTerm_[scvIdx][eqIdx] = 0.0;
}
}
}
@@ -369,11 +113,11 @@ protected:
*/
void evalBoundaryFluxes_()
{
- Dune::GeometryType geoType = element_().geometry().type();
+ Dune::GeometryType geoType = this->element_().geometry().type();
const ReferenceElement &refElement = ReferenceElements::general(geoType);
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -392,19 +136,17 @@ protected:
dim);
int boundaryFaceIdx =
- fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
- asImp_().evalNeumannSegment_(isIt,
- scvIdx,
- boundaryFaceIdx);
+ this->asImp_().evalNeumannSegment_(isIt,
+ scvIdx,
+ boundaryFaceIdx);
// evaluate the outflow conditions at the boundary face
- // ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model
- // available and not thoroughly tested.
- asImp_().evalOutflowSegment_(isIt,
- scvIdx,
- boundaryFaceIdx);
+ this->asImp_().evalOutflowSegment_(isIt,
+ scvIdx,
+ boundaryFaceIdx);
}
}
}
@@ -420,28 +162,28 @@ protected:
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables neumannFlux(0.0);
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
Valgrind::SetUndefined(neumannFlux);
- problem_().boxSDNeumann(neumannFlux,
- element_(),
- fvGeometry_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- curVolVars_());
+ this->problem_().boxSDNeumann(neumannFlux,
+ this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx,
+ this->curVolVars_());
neumannFlux *=
- fvGeometry_().boundaryFace[boundaryFaceIdx].area
- * curVolVars_(scvIdx).extrusionFactor();
+ this->fvGeometry_().boundaryFace[boundaryFaceIdx].area
+ * this->curVolVars_(scvIdx).extrusionFactor();
Valgrind::CheckDefined(neumannFlux);
// set the neumann conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isNeumann(eqIdx))
continue;
- residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
+ this->residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
}
}
}
@@ -464,7 +206,7 @@ protected:
{
//calculate outflow fluxes
PrimaryVariables values(0.0);
- asImp_().computeFlux(values, boundaryFaceIdx, true);
+ this->asImp_().computeFlux(values, boundaryFaceIdx, true);
Valgrind::CheckDefined(values);
for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
@@ -485,20 +227,20 @@ protected:
{
// calculate the mass flux over the faces and subtract
// it from the local rates
- for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++)
+ for (int scvfIdx = 0; scvfIdx < this->fvGeometry_().numEdges; scvfIdx++)
{
- int i = fvGeometry_().subContVolFace[scvfIdx].i;
- int j = fvGeometry_().subContVolFace[scvfIdx].j;
+ int i = this->fvGeometry_().subContVolFace[scvfIdx].i;
+ int j = this->fvGeometry_().subContVolFace[scvfIdx].j;
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, scvfIdx);
+ this->asImp_().computeFlux(flux, scvfIdx);
Valgrind::CheckDefined(flux);
Scalar extrusionFactor =
- (curVolVars_(i).extrusionFactor()
- + curVolVars_(j).extrusionFactor())
+ (this->curVolVars_(i).extrusionFactor()
+ + this->curVolVars_(j).extrusionFactor())
/ 2;
flux *= extrusionFactor;
@@ -516,225 +258,13 @@ protected:
// of sub-control volume i and _INTO_ sub-control volume
// j, we need to add the flux to finite volume i and
// subtract it from finite volume j
- residual_[i] += flux;
- residual_[j] -= flux;
+ this->residual_[i] += flux;
+ this->residual_[j] -= flux;
}
}
-
- /*!
- * \brief Set the local residual to the storage terms of all
- * sub-control volumes of the current element.
- */
- void evalStorage_()
- {
- storageTerm_.resize(fvGeometry_().numVertices);
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) {
- Valgrind::SetUndefined(storageTerm_[scvIdx]);
- asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, /*isOldSol=*/false);
- storageTerm_[scvIdx] *=
- fvGeometry_().subContVol[scvIdx].volume
- * curVolVars_(scvIdx).extrusionFactor();
- Valgrind::CheckDefined(storageTerm_[scvIdx]);
- }
- }
-
- /*!
- * \brief Add the change in the storage terms and the source term
- * to the local residual of all sub-control volumes of the
- * current element.
- */
- void evalVolumeTerms_()
- {
- // evaluate the volume terms (storage + source terms)
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++)
- {
- Scalar extrusionFactor =
- curVolVars_(scvIdx).extrusionFactor();
-
- PrimaryVariables values(0.0);
-
- // mass balance within the element. this is the
- // \f$\frac{m}{\partial t}\f$ term if using implicit
- // euler as time discretization.
- //
- // TODO (?): we might need a more explicit way for
- // doing the time discretization...
- Valgrind::SetUndefined(storageTerm_[scvIdx]);
- Valgrind::SetUndefined(values);
- asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, false);
- asImp_().computeStorage(values, scvIdx, true);
- Valgrind::CheckDefined(storageTerm_[scvIdx]);
- Valgrind::CheckDefined(values);
-
- storageTerm_[scvIdx] -= values;
- storageTerm_[scvIdx] *=
- fvGeometry_().subContVol[scvIdx].volume
- / problem_().timeManager().timeStepSize()
- * extrusionFactor;
- residual_[scvIdx] += storageTerm_[scvIdx];
-
- // subtract the source term from the local rate
- Valgrind::SetUndefined(values);
- asImp_().computeSource(values, scvIdx);
- Valgrind::CheckDefined(values);
- values *= fvGeometry_().subContVol[scvIdx].volume * extrusionFactor;
- residual_[scvIdx] -= values;
-
- // make sure that only defined quantities were used
- // to calculate the residual.
- Valgrind::CheckDefined(residual_[scvIdx]);
- }
- }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return *problemPtr_; };
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); };
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the current element.
- */
- const Element &element_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry &fvGeometry_() const
- {
- Valgrind::CheckDefined(fvElemGeomPtr_);
- return *fvElemGeomPtr_;
- }
-
- /*!
- * \brief Returns a reference to the primary variables of
- * the last time step of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &prevPriVars_(const int i) const
- {
- return prevVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns a reference to the primary variables of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &curPriVars_(const int i) const
- {
- return curVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns the j'th primary of the i'th sub-control volume
- * of the current element.
- */
- Scalar curPriVar_(const int i, const int j) const
- {
- return curVolVars_(i).priVar(j);
- }
-
- /*!
- * \brief Returns a reference to the current volume variables of
- * all sub-control volumes of the current element.
- */
- const ElementVolumeVariables &curVolVars_() const
- {
- Valgrind::CheckDefined(curVolVarsPtr_);
- return *curVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the volume variables of the i-th
- * sub-control volume of the current element.
- */
- const VolumeVariables &curVolVars_(const int i) const
- {
- return curVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of all sub-control volumes of the current
- * element.
- */
- const ElementVolumeVariables &prevVolVars_() const
- {
- Valgrind::CheckDefined(prevVolVarsPtr_);
- return *prevVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of the i-th sub-control volume of the current
- * element.
- */
- const VolumeVariables &prevVolVars_(const int i) const
- {
- return prevVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the boundary types of all
- * sub-control volumes of the current element.
- */
- const ElementBoundaryTypes &bcTypes_() const
- {
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
- }
-
- /*!
- * \brief Returns a reference to the boundary types of the i-th
- * sub-control volume of the current element.
- */
- const BoundaryTypes &bcTypes_(const int i) const
- {
- return bcTypes_()[i];
- }
-
-protected:
- ElementSolutionVector storageTerm_;
- ElementSolutionVector residual_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- const Element *elemPtr_;
- const FVElementGeometry *fvElemGeomPtr_;
-
- // current and previous secondary variables for the element
- const ElementVolumeVariables *prevVolVarsPtr_;
- const ElementVolumeVariables *curVolVarsPtr_;
-
- const ElementBoundaryTypes *bcTypesPtr_;
+ const VertexMapper &vertexMapper_() const
+ { return this->problem_().vertexMapper(); };
};
}
diff --git a/dumux/implicit/box/boxmodel.hh b/dumux/implicit/box/boxmodel.hh
index 567d706dd7fafb52f8c4a9115904b634eefc243a..a73f4741e5d390771c31a6b9fb1733d031957448 100644
--- a/dumux/implicit/box/boxmodel.hh
+++ b/dumux/implicit/box/boxmodel.hh
@@ -28,7 +28,6 @@
#include "boxpropertydefaults.hh"
#include "boxelementvolumevariables.hh"
-#include "boxlocaljacobian.hh"
#include "boxlocalresidual.hh"
#include <dumux/parallel/vertexhandles.hh>
diff --git a/dumux/implicit/box/boxpropertydefaults.hh b/dumux/implicit/box/boxpropertydefaults.hh
index 49b96ba12db6d10b6c358677b627c877a04ede05..6e592e0b4519823b1b4a73f420ee25e681b6cb08 100644
--- a/dumux/implicit/box/boxpropertydefaults.hh
+++ b/dumux/implicit/box/boxpropertydefaults.hh
@@ -32,7 +32,6 @@
#include "boxmodel.hh"
#include "boxfvelementgeometry.hh"
#include "boxelementboundarytypes.hh"
-#include "boxlocaljacobian.hh"
#include "boxlocalresidual.hh"
#include "boxelementvolumevariables.hh"
#include "boxproperties.hh"
@@ -64,9 +63,6 @@ SET_TYPE_PROP(BoxModel, BaseLocalResidual, Dumux::BoxLocalResidual<TypeTag>);
//! Set the BaseModel to BoxModel
SET_TYPE_PROP(BoxModel, BaseModel, Dumux::BoxModel<TypeTag>);
-//! The local jacobian operator for the box scheme
-SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);
-
//! An array of secondary variable containers
SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
diff --git a/dumux/implicit/cellcentered/cclocaljacobian.hh b/dumux/implicit/cellcentered/cclocaljacobian.hh
deleted file mode 100644
index 1aa4d7ff50bf071ac15e734397c70fbcb1b64403..0000000000000000000000000000000000000000
--- a/dumux/implicit/cellcentered/cclocaljacobian.hh
+++ /dev/null
@@ -1,541 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Copyright (C) 2007-2009 by Bernd Flemisch *
- * Institute for Modelling Hydraulic and Environmental Systems *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * 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_CC_LOCAL_JACOBIAN_HH
-#define DUMUX_CC_LOCAL_JACOBIAN_HH
-
-#include <dune/istl/matrix.hh>
-
-#include <dumux/common/math.hh>
-#include "ccelementboundarytypes.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup CCModel
- * \ingroup CCLocalJacobian
- * \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 CCLocalJacobian
-{
-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
- };
-
- 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
- CCLocalJacobian(const CCLocalJacobian &);
-
-public:
- CCLocalJacobian()
- {
- numericDifferenceMethod_ = GET_PARAM(TypeTag, int, ImplicitNumericDifferenceMethod);
- 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(1, 2<<dim);
- storageJacobian_.resize(1);
- }
-
- /*!
- * \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_(), elem_(), 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(elem_());
-
- // 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_(),
- elem_(),
- fvElemGeom_,
- true /* isOldSol? */);
-
- curVolVars_.update(problem_(),
- elem_(),
- fvElemGeom_,
- false /* isOldSol? */);
-
- // update the hints of the model
- model_().updateCurHints(element, curVolVars_);
-
- // calculate the local residual
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
- residual_ = localResidual().residual();
- storageTerm_ = localResidual().storageTerm();
-
- model_().updatePVWeights(elem_(), curVolVars_);
-
- // calculate the local jacobian matrix
- int numNeighbors = fvElemGeom_.numNeighbors;
- ElementSolutionVector partialDeriv(1);
- PrimaryVariables storageDeriv(0.0);
- for (int j = 0; j < numNeighbors; 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(int i, 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(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(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(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(int scvIdx,
- 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 pv = this->curVolVars_[scvIdx].priVar(pvIdx);
- return baseEps*(std::abs(pv) + 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 &elem_() 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_()
- {
- for (int i = 0; i < 1; ++ i) {
- storageJacobian_[i] = 0.0;
- for (int j = 0; j < A_.M(); ++ 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 dest The vector storing the partial derivatives of all
- * equations
- * \param destStorage 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 &dest,
- PrimaryVariables &destStorage,
- int neighborIdx,
- int pvIdx)
- {
- const Element& neighbor = *(fvElemGeom_.neighbors[neighborIdx]);
- FVElementGeometry neighborFVGeom;
- neighborFVGeom.updateInner(neighbor);
-
- int globalIdx = problemPtr_->elementMapper().map(neighbor);
-
- PrimaryVariables priVars(model_().curSol()[globalIdx]);
- VolumeVariables origVolVars(curVolVars_[neighborIdx]);
-
- curVolVars_[neighborIdx].setEvalPoint(&origVolVars);
- Scalar eps = asImp_().numericEpsilon(neighborIdx, 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_[neighborIdx].update(priVars,
- problem_(),
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- false);
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
-
- // store the residual and the storage term
- dest = localResidual().residual();
- if (neighborIdx == 0)
- destStorage = localResidual().storageTerm()[neighborIdx];
- }
- 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)
- dest = residual_;
- if (neighborIdx == 0)
- destStorage = storageTerm_[neighborIdx];
- }
-
-
- 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_[neighborIdx].update(priVars,
- problem_(),
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- false);
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
- dest -= localResidual().residual();
- if (neighborIdx == 0)
- destStorage -= localResidual().storageTerm()[neighborIdx];
- }
- 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)
- dest -= residual_;
- if (neighborIdx == 0)
- destStorage -= storageTerm_[neighborIdx];
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- dest /= delta;
- destStorage /= delta;
-
- // restore the original state of the element's volume variables
- curVolVars_[neighborIdx] = origVolVars;
-
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < dest.size(); ++i)
- Valgrind::CheckDefined(dest[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_(int neighborIdx,
- int pvIdx,
- const ElementSolutionVector &deriv,
- const PrimaryVariables &storageDeriv)
- {
- // store the derivative of the storage term
- if (neighborIdx == 0)
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- storageJacobian_[neighborIdx][eqIdx][pvIdx] = storageDeriv[eqIdx];
- }
-
- for (int i = 0; i < 1; i++)
- {
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // A[i][scvIdx][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at vertex 'i'
- // depending on the primary variable 'pvIdx' at vertex
- // 'scvIdx'.
- this->A_[i][neighborIdx][eqIdx][pvIdx] = deriv[i][eqIdx];
- Valgrind::CheckDefined(this->A_[i][neighborIdx][eqIdx][pvIdx]);
- }
- }
- }
-
- const Element *elemPtr_;
- FVElementGeometry fvElemGeom_;
-
- ElementBoundaryTypes bcTypes_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- // secondary variables at the previous and at the current time
- // levels
- ElementVolumeVariables prevVolVars_;
- ElementVolumeVariables curVolVars_;
-
- LocalResidual localResidual_;
-
- LocalBlockMatrix A_;
- std::vector<MatrixBlock> storageJacobian_;
-
- ElementSolutionVector residual_;
- ElementSolutionVector storageTerm_;
-
- int numericDifferenceMethod_;
-};
-}
-
-#endif
diff --git a/dumux/implicit/cellcentered/cclocalresidual.hh b/dumux/implicit/cellcentered/cclocalresidual.hh
index cdc914e8b2d0b9c446f6486eb2e8a92dcfca1881..abe2f9f0b53319fc00e1f2c0c5806eb0672132cc 100644
--- a/dumux/implicit/cellcentered/cclocalresidual.hh
+++ b/dumux/implicit/cellcentered/cclocalresidual.hh
@@ -31,6 +31,7 @@
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/common/implicitlocalresidual.hh>
#include "ccproperties.hh"
@@ -39,18 +40,16 @@ namespace Dumux
/*!
* \ingroup CCModel
* \ingroup CCLocalResidual
- * \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * \brief Element-wise calculation of the residual for models
+ * based on the fully implicit cell-centered scheme.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class CCLocalResidual
+class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef ImplicitLocalResidual<TypeTag> ParentType;
+ friend class ImplicitLocalResidual<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -62,302 +61,17 @@ private:
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- 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, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
CCLocalResidual(const CCLocalResidual &);
public:
- CCLocalResidual()
- { }
-
- ~CCLocalResidual()
+ CCLocalResidual() : ParentType()
{ }
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarily initialized completely.
- *
- * \param prob The representation of the physical problem to be
- * solved.
- */
- void init(Problem &prob)
- { problemPtr_ = &prob; }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- */
- void eval(const Element &element)
- {
- FVElementGeometry fvElemGeom;
-
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementVolumeVariables volVarsPrev, volVarsCur;
- // update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
-
- volVarsPrev.update(problem_(),
- element,
- fvGeometry_(),
- true /* oldSol? */);
- volVarsCur.update(problem_(),
- element,
- fvGeometry_(),
- false /* oldSol? */);
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- // 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);
-
- asImp_().eval(element, fvGeometry_(), volVarsPrev, volVarsCur, bcTypes);
- }
-
- /*!
- * \brief Compute the storage term for the current solution.
- *
- * This can be used to figure out how much of each conservation
- * quantity is inside the element.
- *
- * \param element The DUNE Codim<0> entity for which the storage
- * term ought to be calculated
- */
- void evalStorage(const Element &element)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvElemGeom;
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
-
- // no previous volume variables!
- prevVolVarsPtr_ = 0;
-
- ElementVolumeVariables volVars;
-
- // update the hints
- model_().setHints(element, volVars);
-
- // calculate volume current variables
- volVars.update(problem_(), element, fvGeometry_(), false);
- curVolVarsPtr_ = &volVars;
-
- asImp_().evalStorage_();
- }
-
- /*!
- * \brief Compute the flux term for the current solution.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element
- */
- void evalFluxes(const Element &element,
- const ElementVolumeVariables &curVolVars)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvElemGeom;
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- residual_.resize(1);
- residual_ = 0;
-
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = 0;
- curVolVarsPtr_ = &curVolVars;
- asImp_().evalFluxes_();
- }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvElemGeom The finite-volume geometry of the element
- * \param prevVolVars The volume averaged variables for all
- * sub-control volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- void eval(const Element &element,
- const FVElementGeometry &fvElemGeom,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
- const ElementBoundaryTypes &bcTypes)
- {
- Valgrind::CheckDefined(prevVolVars);
- Valgrind::CheckDefined(curVolVars);
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom.numNeighbors; i++) {
- prevVolVars[i].checkDefined();
- curVolVars[i].checkDefined();
- }
-#endif // HAVE_VALGRIND
-
- elemPtr_ = &element;
- fvElemGeomPtr_ = &fvElemGeom;
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = &prevVolVars;
- curVolVarsPtr_ = &curVolVars;
-
- // resize the vectors for all terms
- residual_.resize(1);
- storageTerm_.resize(1);
-
- residual_ = 0.0;
- storageTerm_ = 0.0;
-
- asImp_().evalFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- asImp_().evalVolumeTerms_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- // evaluate the boundary conditions
- asImp_().evalBoundary_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
- }
-
- /*!
- * \brief Add Outflow boundary conditions for a single sub-control
- * volume face to the local residual.
- *
- * \note This is so far an empty method doing not more than
- * nothing. A beta version is available for the 2p2c and
- * 2p2cni model. There you can find a sample implementation.
- */
- void evalOutflowSegment(const IntersectionIterator &isIt,
- int scvIdx,
- int boundaryFaceIdx)
- {}
-
- /*!
- * \brief Returns the local residual for all sub-control
- * volumes of the element.
- */
- const ElementSolutionVector &residual() const
- { return residual_; }
-
- /*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
- *
- * \param scvIdx The local index of the sub-control volume
- * (i.e. the element's local vertex index)
- */
- const PrimaryVariables &residual(int scvIdx) const
- { return residual_[scvIdx]; }
-
- /*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
- */
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
-
- /*!
- * \brief Returns the storage term for a given sub-control volumes
- * of the element.
- */
- const PrimaryVariables &storageTerm(int scvIdx) const
- { return storageTerm_[scvIdx]; }
-
- void evalOutflowSegment(const IntersectionIterator &isIt)
- {
- const BoundaryTypes &bcTypes = this->bcTypes_(0);
-
- // deal with outflow boundaries
- if (bcTypes.hasOutflow())
- {
- PrimaryVariables values(0.0);
- asImp_().computeFlux(values,
- /*boundaryFaceIdx=*/isIt->indexInInside(),
- true);
- Valgrind::CheckDefined(values);
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (!bcTypes.isOutflow(eqIdx) )
- continue;
- this->residual_[0][eqIdx] += values[eqIdx];
- }
- }
- }
protected:
- Implementation &asImp_()
- {
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
-
- const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
- }
-
- /*!
- * \brief Evaluate the boundary conditions
- * of the current element.
- */
- void evalBoundary_()
- {
- if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
- asImp_().evalBoundaryFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- if (bcTypes_().hasDirichlet())
- asImp_().evalDirichlet_();
- }
/*!
* \brief Set the values of the Dirichlet boundary control volumes
@@ -365,15 +79,15 @@ protected:
*/
void evalDirichlet_()
{
- PrimaryVariables tmp(0);
+ PrimaryVariables dirichletValues(0);
- const BoundaryTypes &bcTypes = bcTypes_(0);
+ const BoundaryTypes &bcTypes = this->bcTypes_(0);
if (! bcTypes.hasDirichlet())
return;
- Valgrind::SetUndefined(tmp);
- // HACK: ask for Dirichlet value at element center
- asImp_().problem_().dirichletAtPos(tmp, element_().geometry().center());
+ Valgrind::SetUndefined(dirichletValues);
+ // HACK: ask for Dirichlet value at element center
+ this->asImp_().problem_().dirichletAtPos(dirichletValues, this->element_().geometry().center());
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
@@ -381,12 +95,12 @@ protected:
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
assert(0 <= pvIdx && pvIdx < numEq);
- Valgrind::CheckDefined(tmp[pvIdx]);
+ Valgrind::CheckDefined(dirichletValues[pvIdx]);
- residual_[0][eqIdx] =
- (curPriVar_(0, pvIdx) - tmp[pvIdx]);
+ this->residual_[0][eqIdx] =
+ (this->curPriVar_(0, pvIdx) - dirichletValues[pvIdx]);
- storageTerm_[0][eqIdx] = 0.0;
+ this->storageTerm_[0][eqIdx] = 0.0;
}
}
@@ -396,8 +110,8 @@ protected:
*/
void evalBoundaryFluxes_()
{
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -406,18 +120,15 @@ protected:
// add the residual of all vertices of the boundary
// segment
- asImp_().evalNeumannSegment_(isIt);
-
+ this->asImp_().evalNeumannSegment_(isIt);
+
// evaluate the outflow conditions at the boundary face
- // ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model
- // available and not thoroughly tested.
- asImp_().evalOutflowSegment(isIt);
+ this->asImp_().evalOutflowSegment_(isIt);
}
}
/*!
- * \brief Add Neumann boundary conditions for a single sub-control
- * volume face to the local residual.
+ * \brief Add Neumann boundary conditions for a single intersection
*/
void evalNeumannSegment_(const IntersectionIterator &isIt)
{
@@ -425,42 +136,66 @@ protected:
PrimaryVariables values(0.0);
BoundaryTypes bcTypes;
- problem_().boundaryTypes(bcTypes, *isIt);
+ this->problem_().boundaryTypes(bcTypes, *isIt);
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
Valgrind::SetUndefined(values);
- problem_().boxSDNeumann(values,
- element_(),
- fvGeometry_(),
- *isIt,
- /*scvIdx=*/0,
- /*boundaryFaceIdx=*/isIt->indexInInside(),
- curVolVars_());
+ this->problem_().boxSDNeumann(values,
+ this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ /*scvIdx=*/0,
+ /*boundaryFaceIdx=*/isIt->indexInInside(),
+ this->curVolVars_());
values *= isIt->geometry().volume()
- * curVolVars_(0).extrusionFactor();
+ * this->curVolVars_(0).extrusionFactor();
Valgrind::CheckDefined(values);
// set the neumann conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isNeumann(eqIdx))
continue;
- residual_[0][eqIdx] += values[eqIdx];
+ this->residual_[0][eqIdx] += values[eqIdx];
+ }
+ }
+ }
+
+ /*!
+ * \brief Add outflow boundary conditions for a single intersection
+ */
+ void evalOutflowSegment_(const IntersectionIterator &isIt)
+ {
+ const BoundaryTypes &bcTypes = this->bcTypes_(0);
+
+ // deal with outflow boundaries
+ if (bcTypes.hasOutflow())
+ {
+ PrimaryVariables values(0.0);
+ this->asImp_().computeFlux(values,
+ /*boundaryFaceIdx=*/isIt->indexInInside(),
+ true);
+ Valgrind::CheckDefined(values);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (!bcTypes.isOutflow(eqIdx) )
+ continue;
+ this->residual_[0][eqIdx] += values[eqIdx];
}
}
}
/*!
- * \brief Add the flux terms to the local residual of all
- * sub-control volumes of the current element.
+ * \brief Add the flux terms to the local residual of the current element
*/
void evalFluxes_()
{
// calculate the mass flux over the faces and subtract
// it from the local rates
int faceIdx = -1;
- IntersectionIterator isIt = gridView_().ibegin(element_());
- IntersectionIterator isEndIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ IntersectionIterator isEndIt = this->gridView_().iend(this->element_());
for (; isIt != isEndIt; ++isIt) {
if (!isIt->neighbor())
continue;
@@ -469,240 +204,14 @@ protected:
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, faceIdx);
+ this->asImp_().computeFlux(flux, faceIdx);
Valgrind::CheckDefined(flux);
- flux *= curVolVars_(0).extrusionFactor();
-
- // The balance equation for a finite volume is:
- //
- // dStorage/dt = Flux + Source
- //
- // where the 'Flux' and the 'Source' terms represent the
- // mass per second which _ENTER_ the finite
- // volume. Re-arranging this, we get
- //
- // dStorage/dt - Source - Flux = 0
- //
- // Since the flux calculated by computeFlux() goes _OUT_
- // of sub-control volume i and _INTO_ sub-control volume
- // j, we need to add the flux to finite volume i and
- // subtract it from finite volume j
- residual_[0] += flux;
- }
- }
-
- /*!
- * \brief Set the local residual to the storage terms of all
- * sub-control volumes of the current element.
- */
- void evalStorage_()
- {
- storageTerm_.resize(1);
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- Valgrind::SetUndefined(storageTerm_[0]);
- asImp_().computeStorage(storageTerm_[0], 0, /*isOldSol=*/false);
- storageTerm_[0] *= fvGeometry_().subContVol[0].volume
- * curVolVars_(0).extrusionFactor();
- Valgrind::CheckDefined(storageTerm_[0]);
- }
-
- /*!
- * \brief Add the change in the storage terms and the source term
- * to the local residual of all sub-control volumes of the
- * current element.
- */
- void evalVolumeTerms_()
- {
- // evaluate the volume terms (storage + source terms)
- for (int i=0; i < 1; i++)
- {
- Scalar extrusionFactor =
- curVolVars_(i).extrusionFactor();
-
- PrimaryVariables tmp(0.);
-
- // mass balance within the element. this is the
- // \f$\frac{m}{\partial t}\f$ term if using implicit
- // euler as time discretization.
- //
- // TODO (?): we might need a more explicit way for
- // doing the time discretization...
- Valgrind::SetUndefined(storageTerm_[i]);
- Valgrind::SetUndefined(tmp);
- asImp_().computeStorage(storageTerm_[i], i, false);
- asImp_().computeStorage(tmp, i, true);
- Valgrind::CheckDefined(storageTerm_[i]);
- Valgrind::CheckDefined(tmp);
+ flux *= this->curVolVars_(0).extrusionFactor();
- storageTerm_[i] -= tmp;
- storageTerm_[i] *=
- fvGeometry_().subContVol[i].volume
- / problem_().timeManager().timeStepSize()
- * extrusionFactor;
- residual_[i] += storageTerm_[i];
-
- // subtract the source term from the local rate
- Valgrind::SetUndefined(tmp);
- asImp_().computeSource(tmp, i);
- Valgrind::CheckDefined(tmp);
- tmp *= fvGeometry_().subContVol[i].volume * extrusionFactor;
- residual_[i] -= tmp;
-
- // make sure that only defined quantities were used
- // to calculate the residual.
- Valgrind::CheckDefined(residual_[i]);
+ this->residual_[0] += flux;
}
}
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return *problemPtr_; };
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); };
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the current element.
- */
- const Element &element_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry &fvGeometry_() const
- {
- Valgrind::CheckDefined(fvElemGeomPtr_);
- return *fvElemGeomPtr_;
- }
-
- /*!
- * \brief Returns a reference to the primary variables of
- * the last time step of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &prevPriVars_(int i) const
- {
- return prevVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns a reference to the primary variables of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &curPriVars_(int i) const
- {
- return curVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns the j'th primary of the i'th sub-control volume
- * of the current element.
- */
- Scalar curPriVar_(int i, int j) const
- {
- return curVolVars_(i).priVar(j);
- }
-
- /*!
- * \brief Returns a reference to the current volume variables of
- * all sub-control volumes of the current element.
- */
- const ElementVolumeVariables &curVolVars_() const
- {
- Valgrind::CheckDefined(curVolVarsPtr_);
- return *curVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the volume variables of the i-th
- * sub-control volume of the current element.
- */
- const VolumeVariables &curVolVars_(int i) const
- {
- return curVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of all sub-control volumes of the current
- * element.
- */
- const ElementVolumeVariables &prevVolVars_() const
- {
- Valgrind::CheckDefined(prevVolVarsPtr_);
- return *prevVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of the i-th sub-control volume of the current
- * element.
- */
- const VolumeVariables &prevVolVars_(int i) const
- {
- return prevVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the boundary types of all
- * sub-control volumes of the current element.
- */
- const ElementBoundaryTypes &bcTypes_() const
- {
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
- }
-
- /*!
- * \brief Returns a reference to the boundary types of the i-th
- * sub-control volume of the current element.
- */
- const BoundaryTypes &bcTypes_(int i) const
- {
- return bcTypes_()[i];
- }
-
-protected:
- ElementSolutionVector storageTerm_;
- ElementSolutionVector residual_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- const Element *elemPtr_;
- const FVElementGeometry *fvElemGeomPtr_;
-
- // current and previous secondary variables for the element
- const ElementVolumeVariables *prevVolVarsPtr_;
- const ElementVolumeVariables *curVolVarsPtr_;
-
- const ElementBoundaryTypes *bcTypesPtr_;
};
}
diff --git a/dumux/implicit/cellcentered/ccpropertydefaults.hh b/dumux/implicit/cellcentered/ccpropertydefaults.hh
index 9e930ec6d0185f47b8dd2c84d74a97ff3073ae4f..57904da7ed7594bc6fa170f5647ba99aa81e61ff 100644
--- a/dumux/implicit/cellcentered/ccpropertydefaults.hh
+++ b/dumux/implicit/cellcentered/ccpropertydefaults.hh
@@ -36,7 +36,6 @@
#include "ccmodel.hh"
#include "ccfvelementgeometry.hh"
#include "ccelementboundarytypes.hh"
-#include "cclocaljacobian.hh"
#include "cclocalresidual.hh"
#include "ccelementvolumevariables.hh"
#include "ccproperties.hh"
@@ -63,9 +62,6 @@ SET_TYPE_PROP(CCModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>);
//! Set the BaseModel to CCModel
SET_TYPE_PROP(CCModel, BaseModel, Dumux::CCModel<TypeTag>);
-//! The local jacobian operator for the box scheme
-SET_TYPE_PROP(CCModel, LocalJacobian, Dumux::CCLocalJacobian<TypeTag>);
-
//! An array of secondary variable containers
SET_TYPE_PROP(CCModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag>);
diff --git a/dumux/implicit/common/implicitlocaljacobian.hh b/dumux/implicit/common/implicitlocaljacobian.hh
index 7449afbcaf2a2f9854378c76b98db595379b0574..a285c9f2844e9a2eeee49537a1c28a614ffcd4a2 100644
--- a/dumux/implicit/common/implicitlocaljacobian.hh
+++ b/dumux/implicit/common/implicitlocaljacobian.hh
@@ -20,19 +20,17 @@
* \file
* \brief Caculates the Jacobian of the local residual for box models
*/
-#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
-#define DUMUX_BOX_LOCAL_JACOBIAN_HH
+#ifndef DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
+#define DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
#include <dune/istl/matrix.hh>
-
#include <dumux/common/math.hh>
-#include "implicitelementboundarytypes.hh"
namespace Dumux
{
/*!
- * \ingroup BoxModel
- * \ingroup BoxLocalJacobian
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitLocalJacobian
* \brief Calculates the Jacobian of the local residual for box models
*
* The default behavior is to use numeric differentiation, i.e.
@@ -63,7 +61,7 @@ namespace Dumux
* evaluation point and \f$\epsilon\f$ is a small value larger than 0.
*/
template<class TypeTag>
-class BoxLocalJacobian
+class ImplicitLocalJacobian
{
private:
typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) Implementation;
@@ -72,6 +70,7 @@ private:
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 GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
enum {
@@ -94,11 +93,13 @@ private:
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
// copying a local jacobian is not a good idea
- BoxLocalJacobian(const BoxLocalJacobian &);
+ ImplicitLocalJacobian(const ImplicitLocalJacobian &);
public:
- BoxLocalJacobian()
+ ImplicitLocalJacobian()
{
numericDifferenceMethod_ = GET_PARAM_FROM_GROUP(TypeTag, int, Implicit, NumericDifferenceMethod);
Valgrind::SetUndefined(problemPtr_);
@@ -119,8 +120,16 @@ public:
localResidual_.init(prob);
// assume quadrilinears as elements with most vertices
- A_.setSize(2<<dim, 2<<dim);
- storageJacobian_.resize(2<<dim);
+ if (isBox)
+ {
+ A_.setSize(2<<dim, 2<<dim);
+ storageJacobian_.resize(2<<dim);
+ }
+ else
+ {
+ A_.setSize(1, 2<<dim);
+ storageJacobian_.resize(1);
+ }
}
/*!
@@ -176,19 +185,28 @@ public:
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++) {
+ int numRows, numCols;
+ if (isBox)
+ {
+ numRows = numCols = fvElemGeom_.numVertices;
+ }
+ else
+ {
+ numRows = 1;
+ numCols = fvElemGeom_.numNeighbors;
+ }
+ ElementSolutionVector partialDeriv(numRows);
+ PrimaryVariables storageDeriv(0.0);
+ for (int col = 0; col < numCols; col++) {
for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
asImp_().evalPartialDerivative_(partialDeriv,
storageDeriv,
- j,
+ col,
pvIdx);
// update the local stiffness matrix with the current partial
// derivatives
- updateLocalJacobian_(j,
+ updateLocalJacobian_(col,
pvIdx,
partialDeriv,
storageDeriv);
@@ -328,10 +346,9 @@ protected:
*/
void reset_()
{
- int n = element_().template count<dim>();
- for (int i = 0; i < n; ++ i) {
+ for (int i = 0; i < A_.N(); ++ i) {
storageJacobian_[i] = 0.0;
- for (int j = 0; j < n; ++ j) {
+ for (int j = 0; j < A_.M(); ++ j) {
A_[i][j] = 0.0;
}
}
@@ -374,26 +391,37 @@ protected:
* \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
+ * \param col The block column index of the degree of freedom
+ * for which the partial derivative is calculated.
+ * Box: a sub-control volume index.
+ * Cell centered: a neighbor index.
+ * \param pvIdx The index of the primary variable
+ * for which the partial derivative is calculated
*/
void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
PrimaryVariables &storageDeriv,
- const int scvIdx,
+ const int col,
const int pvIdx)
{
- int globalIdx = vertexMapper_().map(element_(), scvIdx, dim);
+ int globalIdx;
+ FVElementGeometry neighborFVGeom;
+ ElementPointer neighbor(element_());
+ if (isBox)
+ {
+ globalIdx = vertexMapper_().map(element_(), col, dim);
+ }
+ else
+ {
+ neighbor = fvElemGeom_.neighbors[col];
+ neighborFVGeom.updateInner(*neighbor);
+ globalIdx = problemPtr_->elementMapper().map(*neighbor);
+ }
PrimaryVariables priVars(model_().curSol()[globalIdx]);
- VolumeVariables origVolVars(curVolVars_[scvIdx]);
+ VolumeVariables origVolVars(curVolVars_[col]);
- curVolVars_[scvIdx].setEvalPoint(&origVolVars);
- Scalar eps = asImp_().numericEpsilon(scvIdx, pvIdx);
+ curVolVars_[col].setEvalPoint(&origVolVars);
+ Scalar eps = asImp_().numericEpsilon(col, pvIdx);
Scalar delta = 0;
if (numericDifferenceMethod_ >= 0) {
@@ -405,12 +433,21 @@ protected:
delta += eps;
// calculate the residual
- curVolVars_[scvIdx].update(priVars,
- problem_(),
- element_(),
- fvElemGeom_,
- scvIdx,
- false);
+ if (isBox)
+ curVolVars_[col].update(priVars,
+ problem_(),
+ element_(),
+ fvElemGeom_,
+ col,
+ false);
+ else
+ curVolVars_[col].update(priVars,
+ problem_(),
+ *neighbor,
+ neighborFVGeom,
+ /*scvIdx=*/0,
+ false);
+
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
@@ -419,14 +456,15 @@ protected:
// store the residual and the storage term
partialDeriv = localResidual().residual();
- storageDeriv = localResidual().storageTerm()[scvIdx];
+ if (isBox || col == 0)
+ storageDeriv = localResidual().storageTerm()[col];
}
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];
+ storageDeriv = storageTerm_[col];
}
@@ -439,26 +477,37 @@ protected:
delta += eps;
// calculate residual again
- curVolVars_[scvIdx].update(priVars,
- problem_(),
- element_(),
- fvElemGeom_,
- scvIdx,
- false);
+ if (isBox)
+ curVolVars_[col].update(priVars,
+ problem_(),
+ element_(),
+ fvElemGeom_,
+ col,
+ false);
+ else
+ curVolVars_[col].update(priVars,
+ problem_(),
+ *neighbor,
+ neighborFVGeom,
+ /*scvIdx=*/0,
+ false);
+
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
partialDeriv -= localResidual().residual();
- storageDeriv -= localResidual().storageTerm()[scvIdx];
+ if (isBox || col == 0)
+ storageDeriv -= localResidual().storageTerm()[col];
}
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];
+ if (isBox || col == 0)
+ storageDeriv -= storageTerm_[col];
}
// divide difference in residuals by the magnitude of the
@@ -467,7 +516,7 @@ protected:
storageDeriv /= delta;
// restore the original state of the element's volume variables
- curVolVars_[scvIdx] = origVolVars;
+ curVolVars_[col] = origVolVars;
#if HAVE_VALGRIND
for (unsigned i = 0; i < partialDeriv.size(); ++i)
@@ -478,29 +527,32 @@ protected:
/*!
* \brief Updates the current local Jacobian matrix with the
* partial derivatives of all equations in regard to the
- * primary variable 'pvIdx' at vertex 'scvIdx' .
+ * primary variable 'pvIdx' at dof 'col' .
*/
- void updateLocalJacobian_(const int scvIdx,
+ void updateLocalJacobian_(const int col,
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];
+ if (isBox || col == 0)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
+ storageJacobian_[col][eqIdx][pvIdx] = storageDeriv[eqIdx];
+ }
}
- for (int i = 0; i < fvElemGeom_.numVertices; i++)
+ for (int i = 0; i < fvElemGeom_.numSCV; i++)
{
// Green vertices are not to be changed!
- if (jacAsm_().vertexColor(element_(), i) != Green) {
+ if (!isBox || jacAsm_().vertexColor(element_(), i) != Green) {
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // A[i][scvIdx][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at vertex 'i'
- // depending on the primary variable 'pvIdx' at vertex
- // 'scvIdx'.
- this->A_[i][scvIdx][eqIdx][pvIdx] = partialDeriv[i][eqIdx];
- Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+ this->A_[i][col][eqIdx][pvIdx] = partialDeriv[i][eqIdx];
+ Valgrind::CheckDefined(this->A_[i][col][eqIdx][pvIdx]);
}
}
}
diff --git a/dumux/implicit/common/implicitlocalresidual.hh b/dumux/implicit/common/implicitlocalresidual.hh
index 53eaa0bf554990078effe06e67ccab58387bf4e9..f020e8bbd1a32fd1a387bc4cfda69a1a7a56a0fd 100644
--- a/dumux/implicit/common/implicitlocalresidual.hh
+++ b/dumux/implicit/common/implicitlocalresidual.hh
@@ -20,8 +20,8 @@
* \file
* \brief Calculates the residual of models based on the box scheme element-wise.
*/
-#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
-#define DUMUX_BOX_LOCAL_RESIDUAL_HH
+#ifndef DUMUX_IMPLICIT_LOCAL_RESIDUAL_HH
+#define DUMUX_IMPLICIT_LOCAL_RESIDUAL_HH
#include <dune/istl/matrix.hh>
#include <dune/grid/common/geometry.hh>
@@ -33,15 +33,15 @@
namespace Dumux
{
/*!
- * \ingroup BoxModel
- * \ingroup BoxLocalResidual
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * using a fully implicit discretization.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class BoxLocalResidual
+class ImplicitLocalResidual
{
private:
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
@@ -49,22 +49,10 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
- dim = GridView::dimension
- };
-
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename GridView::Grid::ctype CoordScalar;
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-
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, BoundaryTypes) BoundaryTypes;
@@ -73,13 +61,13 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
- BoxLocalResidual(const BoxLocalResidual &);
+ ImplicitLocalResidual(const ImplicitLocalResidual &);
public:
- BoxLocalResidual()
+ ImplicitLocalResidual()
{ }
- ~BoxLocalResidual()
+ ~ImplicitLocalResidual()
{ }
/*!
@@ -190,7 +178,7 @@ public:
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
- residual_.resize(fvGeometry_().numVertices);
+ residual_.resize(fvGeometry_().numSCV);
residual_ = 0;
bcTypesPtr_ = &bcTypes;
@@ -225,7 +213,7 @@ public:
Valgrind::CheckDefined(curVolVars);
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry.numVertices; i++) {
+ for (int i = 0; i < prevVolVars.size(); i++) {
prevVolVars[i].checkDefined();
curVolVars[i].checkDefined();
}
@@ -238,9 +226,9 @@ public:
curVolVarsPtr_ = &curVolVars;
// resize the vectors for all terms
- int numVerts = fvGeometry_().numVertices;
- residual_.resize(numVerts);
- storageTerm_.resize(numVerts);
+ int numSCV = fvGeometry_().numSCV;
+ residual_.resize(numSCV);
+ storageTerm_.resize(numSCV);
residual_ = 0.0;
storageTerm_ = 0.0;
@@ -248,14 +236,14 @@ public:
asImp_().evalFluxes_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
asImp_().evalVolumeTerms_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++) {
+ for (int i=0; i < fvGeometry_().numSCV; i++) {
Valgrind::CheckDefined(residual_[i]);
}
#endif // HAVE_VALGRIND
@@ -264,7 +252,7 @@ public:
asImp_().evalBoundary_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
}
@@ -322,7 +310,7 @@ protected:
if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
asImp_().evalBoundaryFluxes_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
@@ -330,209 +318,18 @@ protected:
asImp_().evalDirichlet_();
}
- /*!
- * \brief Set the values of the Dirichlet boundary control volumes
- * of the current element.
- */
- void evalDirichlet_()
- {
- PrimaryVariables dirichletValues(0);
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) {
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
-
- if (bcTypes.hasDirichlet()) {
- // ask the problem for the dirichlet values
- const VertexPointer vPtr = element_().template subEntity<dim>(scvIdx);
- Valgrind::SetUndefined(dirichletValues);
- asImp_().problem_().dirichlet(dirichletValues, *vPtr);
-
- // set the dirichlet conditions
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (bcTypes.isDirichlet(eqIdx)) {
- int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- assert(0 <= pvIdx && pvIdx < numEq);
- Valgrind::CheckDefined(dirichletValues[pvIdx]);
-
- residual_[scvIdx][eqIdx] =
- curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
-
- storageTerm_[scvIdx][eqIdx] = 0.0;
- }
- }
- }
- }
- }
-
- /*!
- * \brief Add all Neumann and outflow boundary conditions to the local
- * residual.
- */
- void evalBoundaryFluxes_()
- {
- Dune::GeometryType geoType = element_().geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geoType);
-
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
- for (; isIt != endIt; ++isIt)
- {
- // handle only faces on the boundary
- if (isIt->boundary()) {
- // Assemble the boundary for all vertices of the current
- // face
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElement.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- ++faceVertIdx)
- {
- int scvIdx = refElement.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
-
- int boundaryFaceIdx =
- fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
-
- // add the residual of all vertices of the boundary
- // segment
- asImp_().evalNeumannSegment_(isIt,
- scvIdx,
- boundaryFaceIdx);
- // evaluate the outflow conditions at the boundary face
- // ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model
- // available and not thoroughly tested.
- asImp_().evalOutflowSegment_(isIt,
- scvIdx,
- boundaryFaceIdx);
- }
- }
- }
- }
-
- /*!
- * \brief Add Neumann boundary conditions for a single sub-control
- * volume face to the local residual.
- */
- void evalNeumannSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
- {
- // temporary vector to store the neumann boundary fluxes
- PrimaryVariables neumannFlux(0.0);
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
-
- // deal with neumann boundaries
- if (bcTypes.hasNeumann()) {
- Valgrind::SetUndefined(neumannFlux);
- problem_().boxSDNeumann(neumannFlux,
- element_(),
- fvGeometry_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- curVolVars_());
- neumannFlux *=
- fvGeometry_().boundaryFace[boundaryFaceIdx].area
- * curVolVars_(scvIdx).extrusionFactor();
- Valgrind::CheckDefined(neumannFlux);
-
- // set the neumann conditions
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!bcTypes.isNeumann(eqIdx))
- continue;
- residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
- }
- }
- }
-
- /*!
- * \brief Add outflow boundary conditions for a single sub-control
- * volume face to the local residual.
- *
- * \param isIt The intersection iterator of current element
- * \param scvIdx The index of the considered face of the sub-control volume
- * \param boundaryFaceIdx The index of the considered boundary face of the sub control volume
- */
- void evalOutflowSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
- {
- const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
- // deal with outflow boundaries
- if (bcTypes.hasOutflow())
- {
- //calculate outflow fluxes
- PrimaryVariables values(0.0);
- asImp_().computeFlux(values, boundaryFaceIdx, true);
- Valgrind::CheckDefined(values);
-
- for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
- {
- if (!bcTypes.isOutflow(equationIdx) )
- continue;
- // deduce outflow
- this->residual_[scvIdx][equationIdx] += values[equationIdx];
- }
- }
- }
-
- /*!
- * \brief Add the flux terms to the local residual of all
- * sub-control volumes of the current element.
- */
- void evalFluxes_()
- {
- // calculate the mass flux over the faces and subtract
- // it from the local rates
- for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++)
- {
- int i = fvGeometry_().subContVolFace[scvfIdx].i;
- int j = fvGeometry_().subContVolFace[scvfIdx].j;
-
- PrimaryVariables flux;
-
- Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, scvfIdx);
- Valgrind::CheckDefined(flux);
-
- Scalar extrusionFactor =
- (curVolVars_(i).extrusionFactor()
- + curVolVars_(j).extrusionFactor())
- / 2;
- flux *= extrusionFactor;
-
- // The balance equation for a finite volume is:
- //
- // dStorage/dt = Flux + Source
- //
- // where the 'Flux' and the 'Source' terms represent the
- // mass per second which _ENTER_ the finite
- // volume. Re-arranging this, we get
- //
- // dStorage/dt - Source - Flux = 0
- //
- // Since the flux calculated by computeFlux() goes _OUT_
- // of sub-control volume i and _INTO_ sub-control volume
- // j, we need to add the flux to finite volume i and
- // subtract it from finite volume j
- residual_[i] += flux;
- residual_[j] -= flux;
- }
- }
-
/*!
* \brief Set the local residual to the storage terms of all
* sub-control volumes of the current element.
*/
void evalStorage_()
{
- storageTerm_.resize(fvGeometry_().numVertices);
+ storageTerm_.resize(fvGeometry_().numSCV);
storageTerm_ = 0;
// calculate the amount of conservation each quantity inside
// all sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) {
+ for (int scvIdx = 0; scvIdx < fvGeometry_().numSCV; scvIdx++) {
Valgrind::SetUndefined(storageTerm_[scvIdx]);
asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, /*isOldSol=*/false);
storageTerm_[scvIdx] *=
@@ -550,7 +347,7 @@ protected:
void evalVolumeTerms_()
{
// evaluate the volume terms (storage + source terms)
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++)
+ for (int scvIdx = 0; scvIdx < fvGeometry_().numSCV; scvIdx++)
{
Scalar extrusionFactor =
curVolVars_(scvIdx).extrusionFactor();
@@ -602,12 +399,6 @@ protected:
const Model &model_() const
{ return problem_().model(); };
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
/*!
* \brief Returns a reference to the grid view.
*/
diff --git a/dumux/implicit/common/implicitporousmediaproblem.hh b/dumux/implicit/common/implicitporousmediaproblem.hh
index 0fbde7b6edbbef7299712e892a35a8ce6c6128b3..fab231a7b6486b762d62ef2c0ca6cfd9410c0872 100644
--- a/dumux/implicit/common/implicitporousmediaproblem.hh
+++ b/dumux/implicit/common/implicitporousmediaproblem.hh
@@ -1,7 +1,6 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * 2012 by Bernd Flemisch *
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
diff --git a/dumux/implicit/common/implicitpropertydefaults.hh b/dumux/implicit/common/implicitpropertydefaults.hh
index 69bdeb04a1475ade88b8acd1ee26cb210073dc7e..032ab4faa5b98cfe0ae9cad92fc650a2338cf58d 100644
--- a/dumux/implicit/common/implicitpropertydefaults.hh
+++ b/dumux/implicit/common/implicitpropertydefaults.hh
@@ -33,6 +33,7 @@
#include <dumux/common/timemanager.hh>
#include "implicitproperties.hh"
+#include "implicitlocaljacobian.hh"
#include "implicitvolumevariables.hh"
namespace Dumux {
@@ -73,7 +74,10 @@ SET_TYPE_PROP(ImplicitBase,
Dune::MCMGElementLayout>);
//! The volume variable class, to be overloaded by the model
-SET_TYPE_PROP(BoxModel, VolumeVariables, Dumux::ImplicitVolumeVariables<TypeTag>);
+SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables<TypeTag>);
+
+//! The local jacobian operator for the box scheme
+SET_TYPE_PROP(ImplicitBase, LocalJacobian, Dumux::ImplicitLocalJacobian<TypeTag>);
//! The type of a solution for the whole grid at a fixed time
SET_TYPE_PROP(ImplicitBase,