diff --git a/dumux/implicit/box/localjacobian.hh b/dumux/implicit/box/localjacobian.hh
index 7ce07c5560c2bfed6ee187336c2dd7f7c8cefbb3..14041d5f88cf53376bed8a39596b63e0b0a6ee0f 100644
--- a/dumux/implicit/box/localjacobian.hh
+++ b/dumux/implicit/box/localjacobian.hh
@@ -18,10 +18,10 @@
*****************************************************************************/
/*!
* \file
- * \brief Caculates the Jacobian of the local residual for fully-implicit models
+ * \brief Caculates the Jacobian of the local residual for fully-implicit box models
*/
-#ifndef DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
-#define DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
+#ifndef DUMUX_IMPLICIT_BOX_LOCAL_JACOBIAN_HH
+#define DUMUX_IMPLICIT_BOX_LOCAL_JACOBIAN_HH
#include <dune/istl/io.hh>
#include <dune/istl/matrix.hh>
@@ -82,18 +82,16 @@ class BoxLocalJacobian : public ImplicitLocalJacobian<TypeTag>
dim = GridView::dimension,
};
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
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, ElementBoundaryTypes) ElementBoundaryTypes;
-
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
- typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
public:
BoxLocalJacobian()
@@ -101,62 +99,56 @@ public:
numericDifferenceMethod_ = GET_PARAM_FROM_GROUP(TypeTag, int, Implicit, NumericDifferenceMethod);
}
-
/*!
* \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, JacobianMatrix& matrix)
+ void assemble(const Element& element, JacobianMatrix& matrix, SolutionVector& residual)
{
- // set the current grid element and update the element's
- // finite volume geometry
- elemPtr_ = &element;
- reset_();
+ // prepare the volvars/fvGeometries for the case when caching is disabled
+ this->model_().fvGeometries_().bind(element);
+ this->model_().curVolVars_().bind(element);
+ this->model_().prevVolVars_().bindElement(element);
+ this->model_().fluxVariablesCache_().bind(element);
- bcTypes_.update(problem_(), element, fvElemGeom_());
+ // the finite volume element geometry
+ const auto& fvGeometry = this->model_().fvGeometries(element);
- // calculate the local residual
+ // the element boundary types
+ bcTypes_.update(this->problem_(), element, fvGeometry);
+
+ // calculate the actual element residual
this->localResidual().eval(element, bcTypes_);
this->residual_ = this->localResidual().residual();
- this->model_().updatePVWeights(fvElemGeom_());
-
- // get stencil informations
- const auto& elementStencil = this->model_().stencils(element).elementStencil();
+ this->model_().updatePVWeights(fvGeometry);
- // calculate derivatives for the dofs inside the element
- ElementSolutionVector partialDeriv(numRows);
- for (auto&& scv : fvElemGeom_().scvs())
+ // calculation of the derivatives
+ for (const auto& scv : fvGeometry.scvs())
{
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ const auto localScvIdx = scv.indexInElement();
+ VolumeVariables origVolVars(this->model_().curVolVars(scv));
+
+ // add precalculated residual for this scv into the global container
+ residual[dofIdx] += this->residual_[localScvIdx];
+
+ // calculate derivatives w.r.t to the privars at the dof at hand
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
{
- asImp_().evalPartialDerivative_(partialDeriv, element, scv, pvIdx);
+ evalPartialDerivative_(matrix, element, scv, pvIdx);
- // update the local stiffness matrix with the current partial derivatives
- for (auto&& scvJ : fvElemGeom_().scvs())
- this->updateGlobalJacobian_(matrix, scv.dofIndex(), scvJ.dofIndex(), pvIdx, partialDeriv[scvJ.indexInElement()]);
+ // restore the original state of the scv's volume variables
+ this->model_().curVolVars_(scv) = origVolVars;
}
- }
- // TODO: calculate derivatives in the case of an extended source stencil
- // const auto& extendedSourceStencil = model_().stencils(element).extendedSourceStencil();
- // for (auto&& globalJ : extendedSourceStencil)
- // {
- // for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- // {
- // evalPartialDerivativeSource_(partialDeriv, globalJ, pvIdx, fluxVarIndices);
-
- // // update the local stiffness matrix with the partial derivatives
- // updateLocalJacobian_(j, pvIdx, partialDeriv);
- // }
- // ++j;
- // }
-
- // for cellcentered methods, calculate the derivatives w.r.t cells in stencil
+ // TODO: what if we have an extended source stencil????
+ }
}
-
+protected:
/*!
* \brief Compute the partial derivatives to a primary variable at
* an degree of freedom.
@@ -193,6 +185,7 @@ public:
*
* \param partialDeriv The vector storing the partial derivatives of all
* equations
+ * \param storageDeriv the mass matrix contributions
* \param col The block column index of the degree of freedom
* for which the partial derivative is calculated.
* Box: a sub-control volume index.
@@ -200,19 +193,19 @@ public:
* \param pvIdx The index of the primary variable
* for which the partial derivative is calculated
*/
- void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
+ void evalPartialDerivative_(JacobianMatrix& matrix,
const Element& element,
- const SubControlVolume &scv,
+ const SubControlVolume& scv,
const int pvIdx)
{
- auto dofIdxGlobal = scv.dofIndex();
+ auto dofIdx = scv.dofIndex();
- PrimaryVariables priVars(model_().curSol()[dofIdxGlobal]);
- VolumeVariables origVolVars(model_().curVolVars(scv));
-
- Scalar eps = asImp_().numericEpsilon(scv, pvIdx);
+ ElementSolutionVector partialDeriv(element.subEntities(dim));
+ PrimaryVariables priVars(this->model_().curSol()[dofIdx]);
+ Scalar eps = this->numericEpsilon(scv, pvIdx);
Scalar delta = 0;
+ // calculate the residual with the forward deflected primary variables
if (numericDifferenceMethod_ >= 0)
{
// we are not using backward differences, i.e. we need to
@@ -222,21 +215,21 @@ public:
priVars[pvIdx] += eps;
delta += eps;
- // update the volume variables
- model_().curVolVars_(scv).update(priVars, problem_(), element, scv);
+ // update the volume variables connected to the dof
+ this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
- // calculate the residual with the deflected primary variables
- localResidual().eval(element, bcTypes_);
+ // calculate the deflected residual
+ this->localResidual().eval(element, bcTypes_);
- // store the residual and the storage term
- partialDeriv = localResidual().residual();
+ // store the residual
+ partialDeriv = this->localResidual().residual();
}
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_;
+ partialDeriv = this->residual_;
}
if (numericDifferenceMethod_ <= 0)
@@ -248,45 +241,38 @@ public:
priVars[pvIdx] -= delta + eps;
delta += eps;
- // update the volume variables
- model_().curVolVars_(scv).update(priVars, problem_(), element, scv);
+ // update the volume variables connected to the dof
+ this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
- // calculate the residual with the deflected primary variables
- localResidual().eval(element, bcTypes_);
+ // calculate the deflected residual
+ this->localResidual().eval(element, bcTypes_);
// subtract the residual from the derivative storage
- partialDeriv -= localResidual().residual();
+ partialDeriv -= this->localResidual().residual();
}
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_;
+ partialDeriv -= this->residual_;
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
partialDeriv /= delta;
- // restore the original state of the scv's volume variables
- model_().curVolVars_(scv) = origVolVars;
-
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < partialDeriv.size(); ++i)
- Valgrind::CheckDefined(partialDeriv[i]);
-#endif
+ // update the global stiffness matrix with the current partial derivatives
+ const auto& fvGeometry = this->model_().fvGeometries(element);
+ for (const auto& scvJ : fvGeometry.scvs())
+ this->updateGlobalJacobian_(matrix, scvJ.dofIndex(), dofIdx, pvIdx, partialDeriv[scvJ.indexInElement()]);
}
- const FVElementGeometry& fvElemGeom_() const
- {
- return model_().fvGeometries(eIdxGlobal_);
- }
-
- IndexType eIdxGlobal_;
- ElementBoundaryTypes bcTypes_;
+private:
int numericDifferenceMethod_;
+
+ ElementBoundaryTypes bcTypes_;
};
}
diff --git a/dumux/implicit/box/localresidual.hh b/dumux/implicit/box/localresidual.hh
index 1eacceb3b1711de1024d04f5ab68b0717caa6627..a43ee2c14a495d4c7a3686f813f31d855ec440a9 100644
--- a/dumux/implicit/box/localresidual.hh
+++ b/dumux/implicit/box/localresidual.hh
@@ -55,14 +55,10 @@ class BoxLocalResidual : public ImplicitLocalResidual<TypeTag>
};
typedef typename GridView::template Codim<0>::Entity Element;
-
- typedef typename GridView::Grid::ctype CoordScalar;
- typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
-
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+ typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
// copying the local residual class is not a good idea
BoxLocalResidual(const BoxLocalResidual &);
@@ -72,116 +68,86 @@ public:
{ }
protected:
+
+ void evalFluxes_()
+ {
+ // calculate the mass flux over the scv faces and subtract
+ const auto& fvGeometry = this->fvGeometry_();
+ for (const auto& scvf : fvGeometry.scvfs())
+ {
+ if (!scvf.boundary())
+ {
+ auto flux = this->asImp_().computeFlux(scvf);
+ const auto& insideScv = this->model_().fvGeometries().subControlVolume(scvf.insideScvIdx());
+ const auto& outsideScv = this->model_().fvGeometries().subControlVolume(scvf.outsideScvIdx());
+
+ this->residual_[insideScv.indexInElement()] += flux;
+ this->residual_[outsideScv.indexInElement()] -= flux;
+ }
+ }
+ }
+
/*!
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
*/
- void evalDirichlet_()
+ void evalDirichlet_(const SubControlVolume& scv, const BoundaryTypes& bcTypes)
{
- PrimaryVariables dirichletValues(0);
- for (int scvIdx = 0; scvIdx < this->fvGeometry_().numScv; ++scvIdx) {
- const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
-
- if (bcTypes.hasDirichlet()) {
- // ask the problem for the dirichlet values
- Valgrind::SetUndefined(dirichletValues);
- this->asImp_().problem_().dirichlet(dirichletValues, this->element_().template subEntity<dim>(scvIdx));
+ PrimaryVariables dirichletValues = this->asImp_().problem_().dirichlet(this->element_(), scv);
- // 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]);
+ // 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]);
- this->residual_[scvIdx][eqIdx] =
- this->curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
+ // get the primary variables
+ const auto& priVars = this->model_().curVolVars(scv).priVars();
- this->storageTerm_[scvIdx][eqIdx] = 0.0;
- }
- }
+ this->residual_[scv.indexInElement()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
}
}
}
/*!
- * \brief Add all Neumann and outflow boundary conditions to the local
- * residual.
+ * \brief Add all fluxes resulting from Neumann and outflow boundary conditions to the local residual.
*/
- void evalBoundaryFluxes_()
+ void evalBoundaryFluxes_(const SubControlVolumeFace &scvf, const SubControlVolume& insideScv, const BoundaryTypes& bcTypes)
{
- Dune::GeometryType geoType = this->element_().geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geoType);
- for (const auto& intersection : intersections(this->gridView_(), this->element_()))
- {
- // handle only faces on the boundary
- if (intersection.boundary()) {
- // Assemble the boundary for all vertices of the current
- // face
- int fIdx = intersection.indexInInside();
- int numFaceVerts = refElement.size(fIdx, 1, dim);
- for (int faceVertexIdx = 0;
- faceVertexIdx < numFaceVerts;
- ++faceVertexIdx)
- {
- int scvIdx = refElement.subEntity(fIdx,
- 1,
- faceVertexIdx,
- dim);
-
- int boundaryFaceIdx =
- this->fvGeometry_().boundaryFaceIndex(fIdx, faceVertexIdx);
+ // evaluate the Neumann conditions at the boundary face
+ if (bcTypes.hasNeumann())
+ this->residual_[insideScv.indexInElement()] += this->asImp_().evalNeumannSegment_(scvf, insideScv, bcTypes);
- // add the residual of all vertices of the boundary
- // segment
- this->asImp_().evalNeumannSegment_(&intersection,
- scvIdx,
- boundaryFaceIdx);
- // evaluate the outflow conditions at the boundary face
- this->asImp_().evalOutflowSegment_(&intersection,
- scvIdx,
- boundaryFaceIdx);
- }
- }
- }
+ // TODO: evaluate the outflow conditions at the boundary face
+ //if (bcTypes.hasOutflow())
+ // flux += this->asImp_().evalOutflowSegment_(&intersection, bcTypes);
}
/*!
- * \brief Add Neumann boundary conditions for a single sub-control
- * volume face to the local residual.
+ * \brief Add Neumann boundary conditions for a single scv face
*/
- template <class IntersectionIterator>
- void evalNeumannSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
+ PrimaryVariables evalNeumannSegment_(const SubControlVolumeFace &scvf,
+ const SubControlVolume& insideScv,
+ const BoundaryTypes &bcTypes)
{
// temporary vector to store the neumann boundary fluxes
- PrimaryVariables neumannFlux(0.0);
- const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+ PrimaryVariables flux(0);
- // deal with neumann boundaries
- if (bcTypes.hasNeumann()) {
- Valgrind::SetUndefined(neumannFlux);
- this->problem_().solDependentNeumann(neumannFlux,
- this->element_(),
- this->fvGeometry_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- this->curVolVars_());
- neumannFlux *=
- this->fvGeometry_().boundaryFace[boundaryFaceIdx].area
- * this->curVolVars_(scvIdx).extrusionFactor();
- Valgrind::CheckDefined(neumannFlux);
+ auto neumannFluxes = this->problem_().neumann(this->element_(), scvf);
- // set the neumann conditions
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!bcTypes.isNeumann(eqIdx))
- continue;
- this->residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
- }
- }
+ // multiply neumann fluxes with the area and the extrusion factor
+ neumannFluxes *= scvf.area()*this->problem_().model().curVolVars(insideScv).extrusionFactor();
+
+ // add fluxes to the temporary vector
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isNeumann(eqIdx))
+ flux[eqIdx] += neumannFluxes[eqIdx];
+
+ return flux;
}
/*!
@@ -192,77 +158,59 @@ protected:
* \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
*/
- template <class IntersectionIterator>
- void evalOutflowSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
+ // template <class IntersectionIterator>
+ // 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);
+ // this->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];
+ // }
+ // }
+ // }
+
+ void evalBoundary_()
{
- const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
- // deal with outflow boundaries
- if (bcTypes.hasOutflow())
+ const auto& fvGeometry = this->fvGeometry_();
+ if (this->bcTypes_().hasNeumann() || this->bcTypes_().hasOutflow())
{
- //calculate outflow fluxes
- PrimaryVariables values(0.0);
- this->asImp_().computeFlux(values, boundaryFaceIdx, true);
- values *= this->curVolVars_(scvIdx).extrusionFactor();
- Valgrind::CheckDefined(values);
- for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
+ for (const auto& scvf : fvGeometry.scvfs())
{
- if (!bcTypes.isOutflow(equationIdx) )
- continue;
- // deduce outflow
- this->residual_[scvIdx][equationIdx] += values[equationIdx];
+ if (scvf.boundary())
+ {
+ const auto& scv = this->problem_().model().fvGeometries().subControlVolume(scvf.insideScvIdx());
+ this->asImp_().evalBoundaryFluxes_(scvf, scv, this->bcTypes_(scv.indexInElement()));
+ }
}
}
- }
- /*!
- * \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 < this->fvGeometry_().numScvf; scvfIdx++)
+ // additionally treat mixed D/N conditions in a strong sense
+ if (this->bcTypes_().hasDirichlet())
{
- int i = this->fvGeometry_().subContVolFace[scvfIdx].i;
- int j = this->fvGeometry_().subContVolFace[scvfIdx].j;
-
- PrimaryVariables flux;
-
- Valgrind::SetUndefined(flux);
- this->asImp_().computeFlux(flux, scvfIdx);
- Valgrind::CheckDefined(flux);
-
- Scalar extrusionFactor =
- (this->curVolVars_(i).extrusionFactor()
- + this->curVolVars_(j).extrusionFactor())
- / 2;
- flux *= extrusionFactor;
+ for (const auto& scv : fvGeometry.scvs())
+ {
+ BoundaryTypes bcTypes = this->bcTypes_(scv.indexInElement());
+ if (!bcTypes.hasDirichlet())
+ continue;
- // 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
- this->residual_[i] += flux;
- this->residual_[j] -= flux;
+ this->asImp_().evalDirichlet_(scv, bcTypes);
+ }
}
}
-
- const VertexMapper &vertexMapper_() const
- { return this->problem_().vertexMapper(); }
};
}