diff --git a/dumux/boxmodels/2p/2pfluxvariables.hh b/dumux/boxmodels/2p/2pfluxvariables.hh
index 4ea7096c07d4c285f68fd0dc19ac8e1bcdabc179..082d62298db5175332c4ea8370ac2e9e309516cd 100644
--- a/dumux/boxmodels/2p/2pfluxvariables.hh
+++ b/dumux/boxmodels/2p/2pfluxvariables.hh
@@ -58,7 +58,7 @@ class TwoPFluxVariables
typedef typename GridView::ctype CoordScalar;
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
enum {
dim = GridView::dimension,
@@ -77,100 +77,90 @@ class TwoPFluxVariables
public:
/*
* \brief The constructor
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemGeom The finite-volume geometry in the box scheme
- * \param faceIdx The local index of the SCV (sub-control-volume) face
- * \param elemDat The volume variables of the current element
*/
- TwoPFluxVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int faceIdx,
- const ElementVolumeVariables &elemDat)
- : fvElemGeom_(elemGeom)
- {
- scvfIdx_ = faceIdx;
+ TwoPFluxVariables()
+ {}
- for (int phase = 0; phase < numPhases; ++phase) {
- potentialGrad_[phase] = Scalar(0);
- }
+#warning Docme
+ void update(const ElementVariables &elemVars, int scvfIdx)
+ {
+ insideScvIdx_ = elemVars.fvElemGeom().subContVolFace[scvfIdx].i;
+ outsideScvIdx_ = elemVars.fvElemGeom().subContVolFace[scvfIdx].j;
- calculateGradients_(problem, element, elemDat);
- calculateK_(problem, element, elemDat);
+ calculateGradients_(elemVars, scvfIdx);
+ calculateNormalFlux_(elemVars, scvfIdx);
};
-public:
- /*
- * \brief Return the intrinsic permeability.
+ /*!
+ * \brief Return the extrusion factor of the SCVF.
*/
- const Tensor &intrinsicPermeability() const
- { return K_; }
+ Scalar extrusionFactor() const
+ { return 1.0; }
/*!
- * \brief Return the pressure potential gradient.
+ * \brief Return a phase's pressure potential gradient.
*
* \param phaseIdx The index of the fluid phase
*/
const Vector &potentialGrad(int phaseIdx) const
{ return potentialGrad_[phaseIdx]; }
+ /*!
+ * \brief Return a phase's pressure potential gradient times
+ * intrinsic permeability times the normal of the sub
+ * control volume face times the area of the SCVF.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar normalFlux(int phaseIdx) const
+ { return normalFlux_[phaseIdx]; }
+
/*!
* \brief Return the local index of the downstream control volume
* for a given phase as a function of the normal flux.
*
- * \param normalFlux The normal flux i.e. the given intrinsic permeability
- * times the pressure potential gradient and SCV face normal.
+ * \param phaseIdx The index of the fluid phase for which the downstream
+ * direction is requested.
*/
- int downstreamIdx(Scalar normalFlux) const
- { return (normalFlux >= 0)?face().j:face().i; }
+ int downstreamIdx(int phaseIdx) const
+ { return (normalFlux_[phaseIdx] > 0)?outsideScvIdx_:insideScvIdx_; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase as a function of the normal flux.
*
- * \param normalFlux The normal flux i.e. the given intrinsic permeability
- * times the pressure potential gradient and SCV face normal.
+ * \param phaseIdx The index of the fluid phase for which the upstream
+ * direction is requested.
*/
- int upstreamIdx(Scalar normalFlux) const
- { return (normalFlux > 0)?face().i:face().j; }
-
- /*!
- * \brief Return the SCV (sub-control-volume) face
- */
- const SCVFace &face() const
- { return fvElemGeom_.subContVolFace[scvfIdx_]; }
+ int upstreamIdx(int phaseIdx) const
+ { return (normalFlux_[phaseIdx] > 0)?insideScvIdx_:outsideScvIdx_; }
protected:
- const FVElementGeometry &fvElemGeom_;
- int scvfIdx_;
-
- // gradients
- Vector potentialGrad_[numPhases];
-
- // intrinsic permeability
- Tensor K_;
-
-private:
- void calculateGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
+ void calculateGradients_(const ElementVariables &elemVars,
+ int scvfIdx)
{
+ // reset all gradients to 0
+ for (int phase = 0; phase < numPhases; ++phase) {
+ potentialGrad_[phase] = Scalar(0);
+ }
+
+ typedef typename FVElementGeometry::SubControlVolumeFace Scvf;
+ const Scvf &scvf = elemVars.fvElemGeom().subContVolFace[scvfIdx];
+
// calculate gradients
- for (int idx = 0;
- idx < fvElemGeom_.numVertices;
- idx++) // loop over adjacent vertices
+ for (int scvIdx = 0;
+ scvIdx < elemVars.numScv();
+ scvIdx ++) // loop over adjacent vertices
{
// FE gradient at vertex idx
- const Vector &feGrad = face().grad[idx];
+ const Vector &feGrad = scvf.grad[scvIdx];
// compute sum of pressure gradients for each phase
for (int phase = 0; phase < numPhases; phase++)
{
// the pressure gradient
Vector tmp(feGrad);
- tmp *= elemVolVars[idx].pressure(phase);
+ tmp *= elemVars.volVars(scvIdx, /*historyIdx=*/0).pressure(phase);
potentialGrad_[phase] += tmp;
}
}
@@ -182,18 +172,22 @@ private:
{
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
- Vector g(problem.boxGravity(element, fvElemGeom_, face().i));
- g += problem.boxGravity(element, fvElemGeom_, face().j);
+ Vector g(elemVars.problem().boxGravity(elemVars.element(),
+ elemVars.fvElemGeom(),
+ insideScvIdx_));
+ g += elemVars.problem().boxGravity(elemVars.element(),
+ elemVars.fvElemGeom(),
+ outsideScvIdx_);
g /= 2;
-
+
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
// calculate the phase density at the integration point. we
// only do this if the wetting phase is present in both cells
- Scalar SI = elemVolVars[face().i].saturation(phaseIdx);
- Scalar SJ = elemVolVars[face().j].saturation(phaseIdx);
- Scalar rhoI = elemVolVars[face().i].density(phaseIdx);
- Scalar rhoJ = elemVolVars[face().j].density(phaseIdx);
+ Scalar SI = elemVars.volVars(insideScvIdx_, /*historyIdx=*/0).saturation(phaseIdx);
+ Scalar SJ = elemVars.volVars(outsideScvIdx_, /*historyIdx=*/0).saturation(phaseIdx);
+ Scalar rhoI = elemVars.volVars(insideScvIdx_, /*historyIdx=*/0).density(phaseIdx);
+ Scalar rhoJ = elemVars.volVars(outsideScvIdx_, /*historyIdx=*/0).density(phaseIdx);
Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
if (fI + fJ == 0)
@@ -212,20 +206,52 @@ private:
}
}
- void calculateK_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
+ void calculateNormalFlux_(const ElementVariables &elemVars,
+ int scvfIdx)
{
- const SpatialParameters &spatialParams = problem.spatialParameters();
+ const SpatialParameters &spatialParams = elemVars.problem().spatialParameters();
+
// calculate the intrinsic permeability
- spatialParams.meanK(K_,
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().j));
+ Tensor K;
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(elemVars,
+ insideScvIdx_),
+ spatialParams.intrinsicPermeability(elemVars,
+ outsideScvIdx_));
+
+ const Vector &normal = elemVars.fvElemGeom().subContVolFace[scvfIdx].normal;
+
+ // calculate the flux in the normal direction of the
+ // current sub control volume face:
+ //
+ // v = - (K grad p) * n
+ //
+ // (the minus comes from the Darcy law which states that
+ // the flux is from high to low pressure potentials.)
+ Vector tmpVec;
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ K.mv(potentialGrad(phaseIdx), tmpVec);
+
+ // scalar product with the face normal
+ normalFlux_[phaseIdx] = 0.0;
+ for (int i = 0; i < Vector::size; ++i)
+ normalFlux_[phaseIdx] += tmpVec[i]*normal[i];
+
+ // flux is along negative pressure gradients
+ normalFlux_[phaseIdx] *= -1;
+ }
}
+
+ // local indices of the inside and the outside sub-control volumes
+ int insideScvIdx_;
+ int outsideScvIdx_;
+
+ // gradients
+ Vector potentialGrad_[numPhases];
+
+ // normal fluxes
+ Scalar normalFlux_[numPhases];
};
} // end namepace
diff --git a/dumux/boxmodels/2p/2plocalresidual.hh b/dumux/boxmodels/2p/2plocalresidual.hh
index 0b554f26fdb2beb559e25dff62c57e86428bf911..8c45c76ba8bc1a42b0893da20ed72b049a0c8ba4 100644
--- a/dumux/boxmodels/2p/2plocalresidual.hh
+++ b/dumux/boxmodels/2p/2plocalresidual.hh
@@ -82,7 +82,7 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
@@ -107,23 +107,26 @@ public:
* \param scvIdx The SCV (sub-control-volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
- void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ void computeStorage(PrimaryVariables &result,
+ const ElementVariables &elemVars,
+ int scvIdx,
+ int historyIdx) const
{
- // if flag usePrevSol is set, the solution from the previous
- // time step is used, otherwise the current solution is
- // used. The secondary variables are used accordingly. This
- // is required to compute the derivative of the storage term
- // using the implicit euler method.
- const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
- const VolumeVariables &vertDat = elemDat[scvIdx];
+ // retrieve the volume variables for the SCV at the specified
+ // point in time
+ const VolumeVariables &volVars = elemVars.volVars(scvIdx, historyIdx);
// wetting phase mass
- result[contiWEqIdx] = vertDat.density(wPhaseIdx) * vertDat.porosity()
- * vertDat.saturation(wPhaseIdx);
+ result[contiWEqIdx] =
+ volVars.density(wPhaseIdx)
+ * volVars.porosity()
+ * volVars.saturation(wPhaseIdx);
// non-wetting phase mass
- result[contiNEqIdx] = vertDat.density(nPhaseIdx) * vertDat.porosity()
- * vertDat.saturation(nPhaseIdx);
+ result[contiNEqIdx] =
+ volVars.density(nPhaseIdx)
+ * volVars.porosity()
+ * volVars.saturation(nPhaseIdx);
}
/*!
@@ -133,16 +136,13 @@ public:
* \param flux The flux over the SCV (sub-control-volume) face for each phase
* \param faceIdx The index of the SCV face
*/
- void computeFlux(PrimaryVariables &flux, int faceIdx) const
+ void computeFlux(PrimaryVariables &flux,
+ const ElementVariables &elemVars,
+ int scvfIdx) const
{
- FluxVariables vars(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
- faceIdx,
- this->curVolVars_());
flux = 0;
- asImp_()->computeAdvectiveFlux(flux, vars);
- asImp_()->computeDiffusiveFlux(flux, vars);
+ asImp_()->computeAdvectiveFlux(flux, elemVars, scvfIdx);
+ asImp_()->computeDiffusiveFlux(flux, elemVars, scvfIdx);
}
/*!
@@ -155,31 +155,40 @@ public:
* This method is called by compute flux and is mainly there for
* derived models to ease adding equations selectively.
*/
- void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const ElementVariables &elemVars,
+ int scvfIdx) const
{
+ const FluxVariables &fluxVars = elemVars.fluxVars(scvfIdx);
+ const FluxVariables &evalPointFluxVars = elemVars.evalPointFluxVars(scvfIdx);
+
////////
// advective fluxes of all components in all phases
////////
Vector tmpVec;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- // calculate the flux in the normal direction of the
+ // data attached to upstream and the downstream vertices
+ // of the current phase. The upstream decision has to be
+ // made using the evaluation point and *not* the current
+ // solution. (although the actual secondary variables must
+ // obviously come from the current solution.)
+ int upIdx = evalPointFluxVars.upstreamIdx(phaseIdx);
+ int dnIdx = evalPointFluxVars.downstreamIdx(phaseIdx);
+ const VolumeVariables &up = elemVars.volVars(upIdx, /*historyIdx=*/0);
+ const VolumeVariables &dn = elemVars.volVars(dnIdx, /*historyIdx=*/0);
+
+ // retrieve the 'velocity' in the normal direction of the
// current sub control volume face:
//
- // v = - (K grad p) * n
+ // v_\alpha := - (K grad p) * n
//
- // (the minus comes from the Darcy law which states that
- // the flux is from high to low pressure potentials.)
- fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx), tmpVec);
- Scalar normalFlux = 0;
- for (int i = 0; i < Vector::size; ++i)
- normalFlux += tmpVec[i]*fluxVars.face().normal[i];
- normalFlux *= -1;
-
- // data attached to upstream and the downstream vertices
- // of the current phase
- const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
- const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
+ // note that v_\alpha is *not* equivalent to the Darcy
+ // velocity because the relative permebility is not
+ // included. Its purpose is to make the upwind decision,
+ // i.e. to decide whether the flow goes from sub control
+ // volume i to j or vice versa.
+ Scalar normalFlux = fluxVars.normalFlux(phaseIdx);
// add advective flux of current component in current
// phase
@@ -204,10 +213,11 @@ public:
* non-isothermal two-phase models to calculate diffusive heat
* fluxes
*/
- void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxData) const
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const ElementVariables &elemVars,
+ int scvfIdx) const
{
- // diffusive fluxes
- flux += 0.0;
+ // no diffusive fluxes for immiscible models
}
/*!
@@ -217,14 +227,12 @@ public:
* \param localVertexIdx The index of the SCV
*
*/
- void computeSource(PrimaryVariables &q, int localVertexIdx) const
+ void computeSource(PrimaryVariables &values,
+ const ElementVariables &elemVars,
+ int scvIdx) const
{
// retrieve the source term intrinsic to the problem
- this->problem_().boxSDSource(q,
- this->elem_(),
- this->fvElemGeom_(),
- localVertexIdx,
- this->curVolVars_());
+ elemVars.problem().source(values, elemVars, scvIdx);
}
diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 604d538d0d87294433bc5962c4fa453089356518..8c2126b9a3f41d834be180f831b9c798f75174f9 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/dumux/boxmodels/2p/2pmodel.hh
@@ -93,7 +93,7 @@ class TwoPModel : public BoxModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
enum {
@@ -122,8 +122,11 @@ public:
*/
Scalar primaryVarWeight(int globalVertexIdx, int pvIdx) const
{
- if (Indices::pressureIdx == pvIdx)
- return std::min(10.0/this->prevSol()[globalVertexIdx][pvIdx], 1.0);
+ if (Indices::pressureIdx == pvIdx) {
+ Scalar absPv =
+ std::abs(this->solution(/*historyIdx=*/1)[globalVertexIdx][pvIdx]);
+ return std::min(10.0/absPv, 1.0);
+ }
return 1;
}
@@ -173,9 +176,7 @@ public:
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
- FVElementGeometry fvElemGeom;
- VolumeVariables volVars;
- ElementVolumeVariables elemVolVars;
+ ElementVariables elemVars(this->problem_());
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
@@ -188,19 +189,14 @@ public:
int idx = this->elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
- fvElemGeom.update(this->gridView_(), *elemIt);
+ elemVars.updateFVElemGeom(*elemIt);
+ elemVars.updateScvVars(/*historyIdx=*/0);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
- volVars.update(sol[globalIdx],
- this->problem_(),
- *elemIt,
- fvElemGeom,
- i,
- false);
+ int globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ const VolumeVariables &volVars = elemVars.volVars(scvIdx, /*historyIdx=*/0);
(*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
(*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
(*pC)[globalIdx] = volVars.capillaryPressure();
@@ -218,6 +214,7 @@ public:
}
};
+#if 0
if(velocityOutput)
{
// calculate vertex velocities
@@ -273,6 +270,7 @@ public:
// Transformation of the global normal vector to normal vector in the reference element
const Dune::FieldMatrix<CoordScalar, dim, dim> jacobianT1 = elemIt->geometry().jacobianTransposed(localPosIP);
+ elemIt->geometry().jacobianTransposed(localPosIP);
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
@@ -331,7 +329,9 @@ public:
(*velocityN)[globalIdx] += scvVelocity;
}
}
+#endif
}
+#if 0
if(velocityOutput)
{
// divide the vertex velocities by the number of adjacent scvs i.e. cells
@@ -340,7 +340,7 @@ public:
(*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
}
}
-
+#endif
writer.attachVertexData(*Sn, "Sn");
writer.attachVertexData(*Sw, "Sw");
writer.attachVertexData(*pN, "pn");
@@ -352,11 +352,13 @@ public:
writer.attachVertexData(*mobN, "mobN");
writer.attachVertexData(*poro, "porosity");
writer.attachVertexData(*Te, "temperature");
+#if 0
if(velocityOutput) // check if velocity output is demanded
{
writer.attachVertexData(*velocityW, "velocityW", dim);
writer.attachVertexData(*velocityN, "velocityN", dim);
}
+#endif
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/boxmodels/2p/2pproblem.hh b/dumux/boxmodels/2p/2pproblem.hh
index 6981be5acc2f130a392677a9fa65c720929c1003..bb30c5d275166b2197d61e3679a6a4d0ef24a764 100644
--- a/dumux/boxmodels/2p/2pproblem.hh
+++ b/dumux/boxmodels/2p/2pproblem.hh
@@ -111,6 +111,26 @@ public:
*/
// \{
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume.
+ *
+ * \param context Container for the volume variables, element,
+ * fvElementGeometry, etc
+ * \param localIdx The local index of the sub control volume inside
+ * the element
+ */
+ template <class Context>
+ // if you get an deprecated warning here, please use context
+ // objects to specify your problem!
+ DUNE_DEPRECATED
+ Scalar temperature(const Context &context,
+ int localIdx) const
+ {
+ return asImp_().boxTemperature(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+ };
+
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume.
*
@@ -125,6 +145,7 @@ public:
Scalar boxTemperature(const Element &element,
const FVElementGeometry fvGeom,
int scvIdx) const
+ DUNE_DEPRECATED // use context objects!
{ return asImp_().temperatureAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
@@ -136,6 +157,7 @@ public:
* \param pos The position in global coordinates where the temperature should be specified.
*/
Scalar temperatureAtPos(const GlobalPosition &pos) const
+ DUNE_DEPRECATED // use context objects!
{ return asImp_().temperature(); }
/*!
@@ -146,8 +168,30 @@ public:
* no energy equation is used.
*/
Scalar temperature() const
+ DUNE_DEPRECATED // use context objects!
{ DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem"); };
+
+ /*!
+ * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
+ *
+ * \param context Container for the volume variables, element,
+ * fvElementGeometry, etc
+ * \param localIdx The local index of the sub control volume inside
+ * the element
+ */
+ template <class Context>
+ // if you get an deprecated warning here, please use context
+ // objects to specify your problem!
+ DUNE_DEPRECATED
+ const Vector &gravity(const Context &context,
+ int localIdx) const
+ {
+ return asImp_().boxGravity(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+ };
+
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
@@ -155,8 +199,9 @@ public:
* it just calls gravityAtPos().
*/
const Vector &boxGravity(const Element &element,
- const FVElementGeometry &fvGeom,
- int scvIdx) const
+ const FVElementGeometry &fvGeom,
+ int scvIdx) const
+ DUNE_DEPRECATED // use context objects!
{ return asImp_().gravityAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index 8612933a7e4ad1a204ccd5deb2d18fe7f2a60262..52391ab80b82605d9c05af62eae5289ddd974dea 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -56,6 +56,8 @@ class TwoPVolumeVariables : public BoxVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
enum {
@@ -89,28 +91,26 @@ public:
* \param isOldSol Evaluate function with solution of current or previous time step
*/
void update(const PrimaryVariables &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
+ const ElementVariables &elemVars,
int scvIdx,
- bool isOldSol)
+ int historyIdx)
{
ParentType::update(priVars,
- problem,
- element,
- elemGeom,
+ elemVars,
scvIdx,
- isOldSol);
+ historyIdx);
asImp().updateTemperature_(priVars,
- element,
- elemGeom,
+ elemVars,
scvIdx,
- problem);
+ historyIdx);
- // material law parameters
+ const SpatialParameters &spatialParams =
+ elemVars.problem().spatialParameters();
+
+ // material law parameters
const MaterialLawParams &materialParams =
- problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+ spatialParams.materialLawParams(elemVars, scvIdx);
Scalar p[numPhases];
Scalar Sn;
@@ -147,9 +147,7 @@ public:
fluidState_);
// porosity
- porosity_ = problem.spatialParameters().porosity(element,
- elemGeom,
- scvIdx);
+ porosity_ = spatialParams.porosity(elemVars, scvIdx);
}
/*!
@@ -218,12 +216,12 @@ public:
protected:
void updateTemperature_(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
+ const ElementVariables &elemVars,
int scvIdx,
- const Problem &problem)
+ int historyIdx)
{
- temperature_ = problem.boxTemperature(element, elemGeom, scvIdx);
+ temperature_ =
+ elemVars.problem().temperature(elemVars, scvIdx);
}
FluidState fluidState_;
diff --git a/dumux/boxmodels/common/boxassembler.hh b/dumux/boxmodels/common/boxassembler.hh
index 045ad46f77fd5b9033efc95b492d44fb60eb403e..ef31386e032977419c8bd9438fe49eddef70a7b5 100644
--- a/dumux/boxmodels/common/boxassembler.hh
+++ b/dumux/boxmodels/common/boxassembler.hh
@@ -113,9 +113,6 @@ public:
BoxAssembler()
{
- enableJacobianRecycling_ = GET_PARAM(TypeTag, bool, EnableJacobianRecycling);
- enablePartialReassemble_ = GET_PARAM(TypeTag, bool, EnablePartialReassemble);
-
problemPtr_ = 0;
matrix_ = 0;
@@ -159,7 +156,7 @@ public:
// initialize the storage part of the Jacobian matrix. Since
// we only need this if Jacobian matrix recycling is enabled,
// we do not waste space if it is disabled
- if (enableJacobianRecycling_) {
+ if (enableJacobianRecycling_()) {
storageJacobian_.resize(numVerts);
storageTerm_.resize(numVerts);
}
@@ -167,7 +164,7 @@ public:
totalElems_ = gridView_().comm().sum(numElems);
// initialize data needed for partial reassembly
- if (enablePartialReassemble_) {
+ if (enablePartialReassemble_()) {
vertexColor_.resize(numVerts);
vertexDelta_.resize(numVerts);
elementColor_.resize(numElems);
@@ -176,14 +173,14 @@ public:
}
/*!
- * \brief Assemble the local jacobian of the problem.
+ * \brief Assemble the global Jacobian of the residual and the residual for the current solution.
*
* The current state of affairs (esp. the previous and the current
* solutions) is represented by the model object.
*/
void assemble()
{
- bool printReassembleStatistics = enablePartialReassemble_ && !reuseMatrix_;
+ bool printReassembleStatistics = enablePartialReassemble_() && !reuseMatrix_;
int succeeded;
try {
assemble_();
@@ -231,7 +228,7 @@ public:
*/
void setMatrixReuseable(bool yesno = true)
{
- if (enableJacobianRecycling_)
+ if (enableJacobianRecycling_())
reuseMatrix_ = yesno;
}
@@ -247,7 +244,7 @@ public:
// do not use partial reassembly for the next iteration
nextReassembleAccuracy_ = 0.0;
- if (enablePartialReassemble_) {
+ if (enablePartialReassemble_()) {
std::fill(vertexColor_.begin(),
vertexColor_.end(),
Red);
@@ -283,7 +280,7 @@ public:
void updateDiscrepancy(const SolutionVector &u,
const SolutionVector &uDelta)
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return;
// update the vector with the distances of the current
@@ -342,7 +339,7 @@ public:
*/
void computeColors(Scalar relTol)
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return;
ElementIterator elemIt = gridView_().template begin<0>();
@@ -486,7 +483,7 @@ public:
*/
int vertexColor(const Element &element, int vertIdx) const
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return Red; // reassemble unconditionally!
int globalIdx = vertexMapper_().map(element, vertIdx, dim);
@@ -500,7 +497,7 @@ public:
*/
int vertexColor(int globalVertIdx) const
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return Red; // reassemble unconditionally!
return vertexColor_[globalVertIdx];
}
@@ -512,7 +509,7 @@ public:
*/
int elementColor(const Element &element) const
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return Red; // reassemble unconditionally!
int globalIdx = elementMapper_().map(element);
@@ -526,7 +523,7 @@ public:
*/
int elementColor(int globalElementIdx) const
{
- if (!enablePartialReassemble_)
+ if (!enablePartialReassemble_())
return Red; // reassemble unconditionally!
return elementColor_[globalElementIdx];
}
@@ -543,8 +540,27 @@ public:
const SolutionVector& residual() const
{ return residual_; }
-
private:
+ static bool enableJacobianRecycling_()
+ { return GET_PARAM(TypeTag, bool, EnableJacobianRecycling); }
+ static bool enablePartialReassemble_()
+ { return GET_PARAM(TypeTag, bool, EnablePartialReassemble); }
+
+ Problem &problem_()
+ { return *problemPtr_; }
+ const Problem &problem_() const
+ { return *problemPtr_; }
+ const Model &model_() const
+ { return problem_().model(); }
+ Model &model_()
+ { return problem_().model(); }
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
+ const VertexMapper &vertexMapper_() const
+ { return problem_().vertexMapper(); }
+ const ElementMapper &elementMapper_() const
+ { return problem_().elementMapper(); }
+
// Construct the BCRS matrix for the global jacobian
void createMatrix_()
{
@@ -623,7 +639,7 @@ private:
// always reset the right hand side.
residual_ = 0.0;
- if (!enablePartialReassemble_) {
+ if (!enablePartialReassemble_()) {
// If partial reassembly of the jacobian is not enabled,
// we can just reset everything!
(*matrix_) = 0;
@@ -740,21 +756,20 @@ private:
residual_[globI] += model_().localJacobian().residual(i);
if (enableJacobianRecycling_) {
- // save the flux term and the jacobian of the
- // storage term in case we can reuse the current
- // linearization...
+ // save storage term and its Jacobian in case we can
+ // reuse the current linearization...
storageJacobian_[globI] +=
- model_().localJacobian().storageJacobian(i);
+ model_().localJacobian().jacobianStorage(i);
storageTerm_[globI] +=
- model_().localJacobian().storageTerm(i);
+ model_().localJacobian().residualStorage(i);
}
// update the jacobian matrix
for (int j=0; j < numVertices; ++ j) {
int globJ = vertexMapper_().map(elem, j, dim);
(*matrix_)[globI][globJ] +=
- model_().localJacobian().mat(i,j);
+ model_().localJacobian().jacobian(i, j);
}
}
}
@@ -785,22 +800,6 @@ private:
}
}
-
- Problem &problem_()
- { return *problemPtr_; }
- const Problem &problem_() const
- { return *problemPtr_; }
- const Model &model_() const
- { return problem_().model(); }
- Model &model_()
- { return problem_().model(); }
- const GridView &gridView_() const
- { return problem_().gridView(); }
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); }
- const ElementMapper &elementMapper_() const
- { return problem_().elementMapper(); }
-
Problem *problemPtr_;
// the jacobian matrix
@@ -827,9 +826,6 @@ private:
Scalar nextReassembleAccuracy_;
Scalar reassembleAccuracy_;
-
- bool enableJacobianRecycling_;
- bool enablePartialReassemble_;
};
} // namespace Dumux
diff --git a/dumux/boxmodels/common/boxelementvolumevariables.hh b/dumux/boxmodels/common/boxelementvolumevariables.hh
index 8e692ab0d072b8e8e4d6446f7f3ac915241082a6..63ac9daf88d3e2406c314304cf51db6e79da28a6 100644
--- a/dumux/boxmodels/common/boxelementvolumevariables.hh
+++ b/dumux/boxmodels/common/boxelementvolumevariables.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2010 by Andreas Lauser *
+ * Copyright (C) 2010-2011 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
@@ -22,8 +22,8 @@
*
* \brief Volume variables gathered on an element
*/
-#ifndef DUMUX_BOX_ELEMENT_VOLUME_VARIABLES_HH
-#define DUMUX_BOX_ELEMENT_VOLUME_VARIABLES_HH
+#ifndef DUMUX_BOX_ELEMENT_VARIABLES_HH
+#define DUMUX_BOX_ELEMENT_VARIABLES_HH
#include "boxproperties.hh"
@@ -38,16 +38,28 @@ namespace Dumux
* volume variables object for each of the element's vertices
*/
template<class TypeTag>
-class BoxElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) >
+class BoxElementVariables
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- typedef std::vector<VolumeVariables> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+
+ // the history size of the time discretization in number of steps
+ enum { timeDiscHistorySize = GET_PROP_VALUE(TypeTag, PTAG(TimeDiscHistorySize)) };
+
+ struct ScvStore_ {
+ VolumeVariables volVars[timeDiscHistorySize];
+ VolumeVariables bcTypes;
+ };
+ typedef std::vector<ScvStore_> ScvVarsVector;
+ typedef std::vector<FluxVariables> ScvfVarsVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
@@ -57,47 +69,107 @@ class BoxElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(Type
enum { dim = GridView::dimension };
enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+
public:
/*!
* \brief The constructor.
*/
- BoxElementVolumeVariables()
- { }
+ explicit BoxElementVariables(const Problem &problem)
+ : gridView_(problem.gridView())
+ {
+ // remember the problem object
+ problemPtr_ = &problem;
+ modelPtr_ = &problem.model();
+ }
/*!
- * \brief Construct the volume variables for all of vertices of an element.
+ * \brief Construct the volume variables of an element from scratch.
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvElemGeom The finite volume geometry of the element
- * \param oldSol Tells whether the model's previous or current solution should be used.
*/
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- bool oldSol)
+ void updateAll(const Element &elem)
+ {
+ updateFVElemGeom(elem);
+ updateBoundaryTypes();
+ updateAllScvVars();
+ updateAllScvfVars();
+ }
+
+ void updateFVElemGeom(const Element &elem)
+ {
+ // remember the current element
+ elemPtr_ = &elem;
+
+ // update the finite element geometry
+ fvElemGeom_.update(gridView_, elem);
+ }
+
+ void updateBoundaryTypes()
{
- const SolutionVector &globalSol =
- oldSol?
- problem.model().prevSol():
- problem.model().curSol();
- const VertexMapper &vertexMapper = problem.vertexMapper();
- // we assert that the i-th shape function is
- // associated to the i-th vert of the element.
- int n = element.template count<dim>();
- this->resize(n);
- for (int i = 0; i < n; i++) {
- const PrimaryVariables &solI
- = globalSol[vertexMapper.map(element, i, dim)];
- (*this)[i].update(solI,
- problem,
- element,
- fvElemGeom,
- i,
- oldSol);
+ // resize the SCV and the SCVF arrays
+ if (fvElemGeom_.numVertices > scvVars_.size()) {
+ scvVars_.resize(fvElemGeom_.numVertices);
+ boundaryTypes_.resize(fvElemGeom_.numVertices);
}
+ if (fvElemGeom_.numEdges > scvfVars_.size())
+ scvfVars_.resize(fvElemGeom_.numEdges);
+
+ // fill the boundary types stuff
+ hasNeumann_ = false;
+ hasDirichlet_ = false;
+ for (int scvIdx = 0; scvIdx < numScv(); ++scvIdx) {
+ int globalIdx = model().vertexMapper().map(element(), scvIdx, dim);
+ boundaryTypes_[scvIdx] = model().boundaryTypes(globalIdx);
+ hasDirichlet_ = hasDirichlet_ || boundaryTypes_[scvIdx]->hasDirichlet();
+ hasNeumann_ = hasNeumann_ || boundaryTypes_[scvIdx]->hasNeumann();
+ }
+
+ // set the default evaluation points
+ scvIdxSaved_ = -1;
+ scvfVarsEval_ = &scvfVars_;
+ }
+
+ void updateAllScvVars()
+ {
+ for (int historyIdx = 0; historyIdx < timeDiscHistorySize; ++ historyIdx)
+ updateScvVars(historyIdx);
};
+ void updateScvVars(int historyIdx)
+ {
+ // update the volume variables for the whole history
+ const VertexMapper &vertexMapper = problem().vertexMapper();
+ const SolutionVector &globalSol = model().solution(historyIdx);
+
+ // set the hints for the volume variables
+ model().setHints(*this, historyIdx);
+
+ int nScv = numScv();
+ for (int scvIdx = 0; scvIdx < nScv; scvIdx++) {
+ const PrimaryVariables &scvSol
+ = globalSol[vertexMapper.map(element(), scvIdx, dim)];
+ updateScvVars(scvSol, scvIdx, historyIdx);
+ }
+ };
+
+ void updateScvVars(const PrimaryVariables &priVars, int scvIdx, int historyIdx)
+ {
+ scvVars_[scvIdx].volVars[historyIdx].update(priVars,
+ /*context=*/*this,
+ scvIdx,
+ historyIdx);
+ }
+
+ void updateAllScvfVars()
+ {
+ for (int scvfIdx = 0; scvfIdx < numScvf(); scvfIdx++) {
+ scvfVars_[scvfIdx].update(/*context=*/ *this,
+ /*localIndex=*/scvfIdx);
+ }
+ }
+
/*!
* \brief Construct the volume variables for all of vertices of an
* element given a solution vector computed by PDELab.
@@ -112,27 +184,216 @@ public:
* \param elementSolVector The local solution for the element using PDELab ordering
*/
template<typename ElemSolVectorType>
- void updatePDELab(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- const ElemSolVectorType& elementSolVector)
+ void updatePDELab(const Element &element,
+ const ElemSolVectorType &elementSolVector)
{
- int n = element.template count<dim>();
- this->resize(n);
- for (int vertexIdx = 0; vertexIdx < n; vertexIdx++)
+ updateFVElemGeom(element);
+ updateBoundaryTypes(element);
+
+ // update the current time step's volume variables
+ PrimaryVariables scvSol;
+ for (int scvIdx = 0; scvIdx < numScv(); scvIdx++)
{
- PrimaryVariables solI(0);
- for (int eqnIdx=0; eqnIdx<numEq; eqnIdx++)
- solI[eqnIdx] = elementSolVector[vertexIdx + eqnIdx*n];
- (*this)[vertexIdx].update(solI,
- problem,
- element,
- fvElemGeom,
- vertexIdx,
- false);
+ // reorder the solution
+ for (int eqnIdx = 0; eqnIdx < numEq; eqnIdx++)
+ scvSol[eqnIdx] = elementSolVector[scvIdx + eqnIdx*numScv()];
+ // update the volume variables for the newest history index
+ updateScvVars_(scvSol, /*historyIdx=*/0, scvIdx);
}
};
+
+ /*!
+ * \brief Return a reference to the problem.
+ */
+ const Problem &problem() const
+ { return *problemPtr_; }
+
+ /*!
+ * \brief Return a reference to the model.
+ */
+ const Model &model() const
+ { return *modelPtr_; }
+
+ /*!
+ * \brief Return a reference to the grid view.
+ */
+ const GridView &gridView() const
+ { return gridView_; }
+
+ /*!
+ * \brief Return the current element.
+ */
+ const Element &element() const
+ { return *elemPtr_; }
+
+ /*!
+ * \brief Return the number of sub-control volumes of the current element.
+ */
+ int numScv() const
+ { return fvElemGeom_.numVertices; }
+
+ /*!
+ * \brief Return the number of sub-control volume faces of the current element.
+ */
+ int numScvf() const
+ { return fvElemGeom_.numEdges; }
+
+ /*!
+ * \brief Return the current finite element geometry.
+ */
+ const FVElementGeometry &fvElemGeom() const
+ { return fvElemGeom_; }
+
+ /*!
+ * \brief Return the position of a local entities in global coordinates
+ */
+ const GlobalPosition &pos(int scvIdx) const
+ { return fvElemGeom_.subContVol[scvIdx].global; }
+
+ /*!
+ * \brief Returns whether the current element is on the domain's
+ * boundary.
+ */
+ bool onBoundary() const
+ { return hasNeumann_ || hasDirichlet_; };
+
+ /*!
+ * \brief Returns whether the current element has a Neumann boundary segment.
+ */
+ bool hasNeumann() const
+ { return hasNeumann_; };
+
+ /*!
+ * \brief Returns whether the current element has a Dirichlet vertex
+ */
+ bool hasDirichlet() const
+ { return hasDirichlet_; };
+
+ /*!
+ * \brief Returns the boundary types for a given vertex
+ */
+ const BoundaryTypes &boundaryTypes(int scvIdx) const
+ {
+ return *boundaryTypes_[scvIdx];
+ }
+
+ /*!
+ * \brief Save the current flux variables and use them as the
+ * evaluation point.
+ */
+ void saveScvfVars()
+ {
+ scvfVarsSaved_ = scvfVars_;
+
+ // change evaluation point
+ scvfVarsEval_ = &scvfVarsSaved_;
+ }
+
+ /*!
+ * \brief Restore current flux variables from the saved ones.
+ */
+ void restoreScvfVars()
+ {
+ //scvfVarsSaved_ = scvfVars_; // not needed
+
+ // change evaluation point
+ scvfVarsEval_ = &scvfVars_;
+ }
+
+ /*!
+ * \brief Return a reference to the volume variables of a
+ * sub-control volume at a given time.
+ *
+ * If the time step index is not given, return the volume
+ * variables for the current time.
+ *
+ * \param scvIdx The local index of the sub-control volume for
+ * which the volume variables are requested
+ * \param historyIdx The index of the time step for which the
+ * volume variables are requested. 0 means
+ * current time step, 1 previous time step,
+ * 2 next-to-previous, etc.
+ */
+ const VolumeVariables &volVars(int scvIdx, int historyIdx = 0) const
+ { return scvVars_[scvIdx].volVars[historyIdx]; }
+
+ /*!
+ * \copydoc volVars()
+ */
+ VolumeVariables &volVars(int scvIdx, int historyIdx = 0)
+ { return scvVars_[scvIdx].volVars[historyIdx]; }
+
+ /*!
+ * \brief Returns the volume variables at the evaluation point.
+ */
+ void saveScvVars(int scvIdx)
+ {
+ scvIdxSaved_ = scvIdx;
+ scvVarsSaved_ = scvVars_[scvIdx].volVars[/*historyIdx=*/0];
+ }
+
+ /*!
+ * \brief Restores the volume variables at the evaluation point.
+ */
+ void restoreScvVars(int scvIdx)
+ {
+ scvIdxSaved_ = -1;
+ scvVars_[scvIdx].volVars[/*historyIdx=*/0] = scvVarsSaved_;
+ }
+
+ /*!
+ * \brief Return a reference to the flux variables of a
+ * sub-control volume face.
+ *
+ * \param scvfIdx The local index of the sub-control volume face for
+ * which the flux variables are requested
+ */
+ const FluxVariables &fluxVars(int scvIdx) const
+ { return scvfVars_[scvIdx]; }
+
+ /*!
+ * \brief Return a reference to the flux variables of a
+ * sub-control volume face for the evaluation point.
+ *
+ * \param scvfIdx The local index of the sub-control volume face for
+ * which the flux variables are requested
+ */
+ const FluxVariables &evalPointFluxVars(int scvfIdx) const
+ {
+ return (*scvfVarsEval_)[scvfIdx];
+ }
+
+ /*!
+ * \brief Returns the volume variables for history index 0 at the
+ * evaluation point.
+ */
+ const VolumeVariables &evalPointVolVars(int scvIdx) const
+ {
+ if (scvIdxSaved_ == scvIdx)
+ return scvVarsSaved_;
+ return volVars(scvIdx, /*historyIdx=*/0);
+ }
+
+protected:
+ ScvVarsVector scvVars_;
+
+ int scvIdxSaved_;
+ VolumeVariables scvVarsSaved_;
+
+ ScvfVarsVector scvfVars_;
+ ScvfVarsVector scvfVarsSaved_;
+
+ ScvfVarsVector *scvfVarsEval_;
+
+ const Problem *problemPtr_;
+ const Model *modelPtr_;
+ const Element *elemPtr_;
+ const GridView gridView_;
+ FVElementGeometry fvElemGeom_;
+ bool hasNeumann_;
+ bool hasDirichlet_;
+ std::vector<const BoundaryTypes*> boundaryTypes_;
};
} // namespace Dumux
diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index 9ccf1ffed81e285d28b74674d5e9d9fcff0eeab8..05fee5964a1c63a6dcf78e22c55fb9cbb0720c4b 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -30,6 +30,7 @@
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH
+#include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh>
#include "boxelementboundarytypes.hh"
@@ -77,17 +78,12 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
-
- Red = JacobianAssembler::Red,
- Yellow = JacobianAssembler::Yellow,
- Green = JacobianAssembler::Green,
};
@@ -113,22 +109,29 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+ typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> LocalBlockVector;
+ typedef Dune::BlockVector<MatrixBlock> LocalStorageMatrix;
+
// copying a local jacobian is not a good idea
BoxLocalJacobian(const BoxLocalJacobian &);
public:
BoxLocalJacobian()
- {
- numericDifferenceMethod_ = GET_PARAM(TypeTag, int, NumericDifferenceMethod);
- Valgrind::SetUndefined(problemPtr_);
+ {
+ internalElemVars_ = 0;
}
+ ~BoxLocalJacobian()
+ {
+ delete internalElemVars_;
+ }
/*!
* \brief Initialize the local Jacobian object.
@@ -139,139 +142,66 @@ public:
* \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);
+ modelPtr_ = &prob.model();
+ internalElemVars_ = new ElementVariables(prob);
}
/*!
* \brief Assemble an element's local Jacobian matrix of the
* defect.
*
- * \param element The DUNE Codim<0> entity which we look at.
+ * This assembles the 'grad f(x^k)' and 'f(x^k)' part of the newton update
*/
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_);
+ internalElemVars_->updateAll(element);
+ assemble(*internalElemVars_);
+ }
+
+ /*!
+ * \brief Assemble an element's local Jacobian matrix of the
+ * defect, given all secondary variables for the element.
+ *
+ * After calling this method the ElementVariables are in undefined
+ * state, so do not use it anymore!
+ */
+ void assemble(ElementVariables &elemVars)
+ {
+ // update the weights of the primary variables using the
+ // current element variables
+ model_().updatePVWeights(elemVars);
+
+ resize_(elemVars);
+ reset_(elemVars);
+
// calculate the local residual
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
- residual_ = localResidual().residual();
- storageTerm_ = localResidual().storageTerm();
+ localResidual_.eval(residual_, residualStorage_, elemVars);
- model_().updatePVWeights(elem_(), curVolVars_);
+ // save all flux variables calculated using the unmodified
+ // primary variables. This automatically makes these flux
+ // variables the evaluation point.
+ elemVars.saveScvfVars();
// calculate the local jacobian matrix
- int numVertices = fvElemGeom_.numVertices;
- ElementSolutionVector partialDeriv(numVertices);
- PrimaryVariables storageDeriv;
- for (int j = 0; j < numVertices; j++) {
+ int numScv = elemVars.numScv();
+ for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
- asImp_().evalPartialDerivative_(partialDeriv,
- storageDeriv,
- j,
+ asImp_().evalPartialDerivative_(elemVars,
+ scvIdx,
pvIdx);
- // update the local stiffness matrix with the current partial
- // derivatives
- updateLocalJacobian_(j,
- pvIdx,
- partialDeriv,
- storageDeriv);
+ // update the local stiffness matrix with the current
+ // partial derivatives
+ updateLocalJacobian_(elemVars, scvIdx, pvIdx);
}
}
- }
-
- /*!
- * \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]; }
+ // restore flux variables.
+ //elemVars.restoreScvfVars(); // not necessary
+ }
/*!
* \brief Returns the epsilon value which is added and removed
@@ -281,7 +211,8 @@ public:
* which the local derivative ought to be calculated.
* \param pvIdx The index of the primary variable which gets varied
*/
- Scalar numericEpsilon(int scvIdx,
+ Scalar numericEpsilon(const ElementVariables &elemVars,
+ int scvIdx,
int pvIdx) const
{
// define the base epsilon as the geometric mean of 1 and the
@@ -291,71 +222,102 @@ public:
static const Scalar baseEps
= Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(),
1.0);
-
+
// the epsilon value used for the numeric differentiation is
// now scaled by the absolute value of the primary variable...
- Scalar pv = this->curVolVars_[scvIdx].primaryVar(pvIdx);
+ Scalar pv = elemVars.volVars(scvIdx, /*historyIdx=*/0).primaryVar(pvIdx);
return baseEps*(std::abs(pv) + 1);
}
-protected:
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
+ /*!
+ * \brief Return reference to the local residual.
+ */
+ LocalResidual &localResidual()
+ { return localResidual_; }
+ const LocalResidual &localResidual() const
+ { return localResidual_; }
/*!
- * \brief Returns a reference to the problem.
+ * \brief Returns the local Jacobian matrix of the residual of a sub-control volume.
+ *
+ * \param domainScvIdx The local index of the sub control volume which contains the independents
+ * \param rangeScvIdx The local index of the sub control volume which contains the local residual
*/
- const Problem &problem_() const
- {
- Valgrind::CheckDefined(problemPtr_);
- return *problemPtr_;
- };
+ const MatrixBlock &jacobian(int domainScvIdx, int rangeScvIdx) const
+ { return jacobian_[domainScvIdx][rangeScvIdx]; }
/*!
- * \brief Returns a reference to the grid view.
+ * \brief Returns the local Jacobian matrix the storage term of a sub-control volume.
+ *
+ * \param scvIdx The local index of sub control volume
*/
- const GridView &gridView_() const
- { return problem_().gridView(); }
+ const MatrixBlock &jacobianStorage(int scvIdx) const
+ { return jacobianStorage_[scvIdx]; }
/*!
- * \brief Returns a reference to the element.
+ * \brief Returns the local residual of a sub-control volume.
+ *
+ * \param scvIdx The local index of the sub control volume
*/
- const Element &elem_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- };
+ const VectorBlock &residual(int scvIdx) const
+ { return residual_[scvIdx]; }
/*!
- * \brief Returns a reference to the model.
+ * \brief Returns the local storage term of a sub-control volume.
+ *
+ * \param scvIdx The local index of the sub control volume
*/
+ const VectorBlock &residualStorage(int scvIdx) const
+ { return residualStorage_[scvIdx]; }
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ const Problem &problem_() const
+ { return *problemPtr_; }
const Model &model_() const
- { return problem_().model(); };
+ { return *modelPtr_; }
/*!
- * \brief Returns a reference to the jacobian assembler.
+ * \brief Returns the numeric difference method which is applied.
*/
- const JacobianAssembler &jacAsm_() const
- { return model_().jacobianAssembler(); }
+ static int numericDifferenceMethod_()
+ { return GET_PARAM(TypeTag, int, NumericDifferenceMethod); }
/*!
- * \brief Returns a reference to the vertex mapper.
+ * \brief Resize all internal attributes to the size of the
+ * element.
*/
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
+ void resize_(const ElementVariables &elemVars)
+ {
+ int n = elemVars.numScv();
+
+ jacobian_.setSize(n, n);
+ jacobianStorage_.resize(n);
+
+ residual_.resize(n);
+ residualStorage_.resize(n);
+
+ derivResidual_.resize(n);
+ derivStorage_.resize(n);
+ };
/*!
- * \brief Reset the local jacobian matrix to 0
+ * \brief Reset the all relevant internal attributes to 0
*/
- void reset_()
+ void reset_(const ElementVariables &elemVars)
{
- int n = elem_().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;
+ int numScv = elemVars.numScv();
+ for (int i = 0; i < numScv; ++ i) {
+ residual_[i] = 0.0;
+ residualStorage_[i] = 0.0;
+
+ jacobianStorage_[i] = 0.0;
+ for (int j = 0; j < numScv; ++ j) {
+ jacobian_[i][j] = 0.0;
}
}
}
@@ -404,21 +366,19 @@ protected:
* for which the partial derivative ought to be
* calculated
*/
- void evalPartialDerivative_(ElementSolutionVector &dest,
- PrimaryVariables &destStorage,
+ void evalPartialDerivative_(ElementVariables &elemVars,
int scvIdx,
int pvIdx)
{
- int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
-
- PrimaryVariables priVars(model_().curSol()[globalIdx]);
- VolumeVariables origVolVars(curVolVars_[scvIdx]);
+ // save all quantities which depend on the specified primary
+ // variable at the given sub control volume
+ elemVars.saveScvVars(scvIdx);
- curVolVars_[scvIdx].setEvalPoint(&origVolVars);
- Scalar eps = asImp_().numericEpsilon(scvIdx, pvIdx);
+ PrimaryVariables priVars(elemVars.volVars(scvIdx, /*historyIdx=*/0).primaryVars());
+ Scalar eps = asImp_().numericEpsilon(elemVars, scvIdx, pvIdx);
Scalar delta = 0;
- if (numericDifferenceMethod_ >= 0) {
+ if (numericDifferenceMethod_() >= 0) {
// we are not using backward differences, i.e. we need to
// calculate f(x + \epsilon)
@@ -427,32 +387,19 @@ protected:
delta += eps;
// calculate the residual
- curVolVars_[scvIdx].update(priVars,
- problem_(),
- elem_(),
- fvElemGeom_,
- scvIdx,
- false);
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
-
- // store the residual and the storage term
- dest = localResidual().residual();
- destStorage = localResidual().storageTerm()[scvIdx];
+ elemVars.updateScvVars(priVars, scvIdx, /*historyIdx=*/0);
+ elemVars.updateAllScvfVars();
+ localResidual_.eval(derivResidual_, derivStorage_, elemVars);
}
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_;
- destStorage = storageTerm_[scvIdx];
+ derivResidual_ = residual_;
+ derivStorage_[scvIdx] = residualStorage_[scvIdx];
}
-
- if (numericDifferenceMethod_ <= 0) {
+ if (numericDifferenceMethod_() <= 0) {
// we are not using forward differences, i.e. we don't
// need to calculate f(x - \epsilon)
@@ -460,40 +407,35 @@ protected:
priVars[pvIdx] -= delta + eps;
delta += eps;
- // calculate residual again
- curVolVars_[scvIdx].update(priVars,
- problem_(),
- elem_(),
- fvElemGeom_,
- scvIdx,
- false);
- localResidual().eval(elem_(),
- fvElemGeom_,
- prevVolVars_,
- curVolVars_,
- bcTypes_);
- dest -= localResidual().residual();
- destStorage -= localResidual().storageTerm()[scvIdx];
+ // calculate residual again, this time we use the local
+ // residual's internal storage.
+ elemVars.updateScvVars(priVars, scvIdx, /*historyIdx=*/0);
+ elemVars.updateAllScvfVars();
+ localResidual_.eval(elemVars);
+
+ derivResidual_ -= localResidual_.residual();
+ derivStorage_[scvIdx] -= 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)
- dest -= residual_;
- destStorage -= storageTerm_[scvIdx];
+ derivResidual_ -= residual_;
+ derivStorage_[scvIdx] -= residualStorage_[scvIdx];
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
- dest /= delta;
- destStorage /= delta;
+ derivResidual_ /= delta;
+ derivStorage_[scvIdx] /= delta;
- // restore the original state of the element's volume variables
- curVolVars_[scvIdx] = origVolVars;
+ // restore the original state of the element's volume
+ // variables
+ elemVars.restoreScvVars(scvIdx);
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < dest.size(); ++i)
- Valgrind::CheckDefined(dest[i]);
+#ifndef NDEBUG
+ for (unsigned i = 0; i < derivResidual_.size(); ++i)
+ Valgrind::CheckDefined(derivResidual_[i]);
#endif
}
@@ -502,56 +444,44 @@ protected:
* partial derivatives of all equations in regard to the
* primary variable 'pvIdx' at vertex 'scvIdx' .
*/
- void updateLocalJacobian_(int scvIdx,
- int pvIdx,
- const ElementSolutionVector &deriv,
- const PrimaryVariables &storageDeriv)
+ void updateLocalJacobian_(const ElementVariables &elemVars,
+ int scvIdx,
+ int pvIdx)
{
// store the derivative of the storage term
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- storageJacobian_[scvIdx][eqIdx][pvIdx] = storageDeriv[eqIdx];
+ jacobianStorage_[scvIdx][eqIdx][pvIdx] = derivStorage_[scvIdx][eqIdx];
}
- for (int i = 0; i < fvElemGeom_.numVertices; i++)
+ int numScv = elemVars.numScv();
+ for (int eqScvIdx = 0; eqScvIdx < numScv; eqScvIdx++)
{
- if (jacAsm_().vertexColor(elem_(), i) == Green) {
- // Green vertices are not to be changed!
- continue;
- }
-
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] = deriv[i][eqIdx];
- Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
+ // A[eqScvIdx][scvIdx][eqIdx][pvIdx] is the rate of
+ // change of the residual of equation 'eqIdx' at
+ // vertex 'eqScvIdx' depending on the primary variable
+ // 'pvIdx' at vertex 'scvIdx'.
+ jacobian_[eqScvIdx][scvIdx][eqIdx][pvIdx] = derivResidual_[eqScvIdx][eqIdx];
+ Valgrind::CheckDefined(jacobian_[eqScvIdx][scvIdx][eqIdx][pvIdx]);
}
}
}
- const Element *elemPtr_;
- FVElementGeometry fvElemGeom_;
-
- ElementBoundaryTypes bcTypes_;
-
- // The problem we would like to solve
Problem *problemPtr_;
+ Model *modelPtr_;
- // secondary variables at the previous and at the current time
- // levels
- ElementVolumeVariables prevVolVars_;
- ElementVolumeVariables curVolVars_;
+ ElementVariables *internalElemVars_;
- LocalResidual localResidual_;
+ LocalBlockMatrix jacobian_;
+ LocalStorageMatrix jacobianStorage_;
- LocalBlockMatrix A_;
- std::vector<MatrixBlock> storageJacobian_;
+ LocalBlockVector residual_;
+ LocalBlockVector residualStorage_;
- ElementSolutionVector residual_;
- ElementSolutionVector storageTerm_;
+ LocalBlockVector derivResidual_;
+ LocalBlockVector derivStorage_;
- int numericDifferenceMethod_;
+ LocalResidual localResidual_;
};
}
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index d22a826c8fa3dc6143c212761b96aac1e6918a3c..b3176627521054f4dd00e5e3a041a1efc27dda12 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -25,7 +25,7 @@
#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
#define DUMUX_BOX_LOCAL_RESIDUAL_HH
-#include <dune/istl/matrix.hh>
+#include <dune/istl/bvector.hh>
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
@@ -47,48 +47,33 @@ class BoxLocalResidual
private:
typedef BoxLocalResidual<TypeTag> ThisType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
-
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GridView::Grid::ctype CoordScalar;
-
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum { dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
+ typedef typename GridView::ctype CoordScalar;
typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
+ typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
-
- typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
- typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+ enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
+ typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> LocalBlockVector;
+
// copying the local residual class is not a good idea
- BoxLocalResidual(const BoxLocalResidual &);
+ BoxLocalResidual(const BoxLocalResidual &)
+ {};
public:
BoxLocalResidual()
@@ -97,312 +82,217 @@ public:
~BoxLocalResidual()
{ }
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarrily 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.
+ * conservation equations from zero and store the results
+ * internally.
*
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
+ * The results can be requested afterwards using the residual()
+ * and storageTerm() methods.
*/
- void eval(const Element &element)
+ void eval(const ElementVariables &elemVars)
{
- FVElementGeometry fvElemGeom;
-
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementVolumeVariables volVarsPrev, volVarsCur;
- // update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
-
- volVarsPrev.update(problem_(),
- element,
- fvElemGeom_(),
- true /* oldSol? */);
- volVarsCur.update(problem_(),
- element,
- fvElemGeom_(),
- false /* oldSol? */);
-
- ElementBoundaryTypes bcTypes;
- 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);
-
- asImp_().eval(element, fvElemGeom_(), volVarsPrev, volVarsCur, bcTypes);
- }
+ int numScv = elemVars.numScv();
+ internalResidual_.resize(numScv);
+ internalStorageTerm_.resize(numScv);
+ eval(internalResidual_, internalStorageTerm_, elemVars);
+ };
+
/*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero. Gets a solution vector computed by PDELab
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param elementSolVector The local solution for the element using PDELab ordering
+ * \brief Return the result of the eval() call using internal
+ * storage.
*/
- template<typename ElemSolVectorType>
- void evalPDELab(const Element &element, const ElemSolVectorType& elementSolVector)
- {
- FVElementGeometry fvGeom;
- fvGeom.update(gridView_(), element);
- ElementVolumeVariables volVarsPrev, volVarsCur;
- volVarsPrev.update(problem_(),
- element,
- fvGeom,
- true /* oldSol? */);
- volVarsCur.updatePDELab(problem_(),
- element,
- fvGeom,
- elementSolVector);
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeom);
-
- // 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, fvGeom, volVarsPrev, volVarsCur, bcTypes);
- }
+ const LocalBlockVector &residual() const
+ { return internalResidual_; }
+ const VectorBlock &residual(int scvIdx) const
+ { return internalResidual_[scvIdx]; }
/*!
- * \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
+ * \brief Return the storage term calculated using the last call
+ * to eval() using internal storage.
*/
- void evalStorage(const Element &element)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvElemGeom;
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvElemGeom_());
- 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, fvElemGeom_(), 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, fvElemGeom_());
-
- residual_.resize(fvElemGeom_().numVertices);
- residual_ = 0;
-
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = 0;
- curVolVarsPtr_ = &curVolVars;
- asImp_().evalFluxes_();
- }
-
+ const LocalBlockVector &storageTerm() const
+ { return internalStorageTerm_; }
+ const VectorBlock &storageTerm(int scvIdx) const
+ { return internalStorageTerm_[scvIdx]; }
+
/*!
* \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
+ * conservation equations from zero.
*
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvGeom The finite-volume geometry of the element
- * \param prevVolVars The volume averaged variables for all
- * sub-contol volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
+ * \param residual Stores the residual vector
+ * \param storageTerm Stores the derivative of the storage term to time
+ * \param elemVars All the element's secondary variables required to calculate the local residual
*/
- void eval(const Element &element,
- const FVElementGeometry &fvGeom,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
- const ElementBoundaryTypes &bcTypes)
+ void eval(LocalBlockVector &residual,
+ LocalBlockVector &storageTerm,
+ const ElementVariables &elemVars) const
{
- //Valgrind::CheckDefined(fvGeom);
- Valgrind::CheckDefined(prevVolVars);
- Valgrind::CheckDefined(curVolVars);
+ residual = 0.0;
+ storageTerm = 0.0;
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeom.numVertices; i++) {
- prevVolVars[i].checkDefined();
- curVolVars[i].checkDefined();
- }
-#endif // HAVE_VALGRIND
-
- elemPtr_ = &element;
- fvElemGeomPtr_ = &fvGeom;
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = &prevVolVars;
- curVolVarsPtr_ = &curVolVars;
-
- // resize the vectors for all terms
- int numVerts = fvElemGeom_().numVertices;
- residual_.resize(numVerts);
- storageTerm_.resize(numVerts);
-
- residual_ = 0.0;
- storageTerm_ = 0.0;
-
- asImp_().evalFluxes_();
+ // evaluate the flux terms
+ asImp_().evalFluxes(residual, elemVars);
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
+ for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
+ Valgrind::CheckDefined(residual[i]);
#endif // HAVE_VALGRIND
- asImp_().evalVolumeTerms_();
+ // evaluate the storage and the source terms
+ asImp_().evalVolumeTerms_(residual, storageTerm, elemVars);
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom_().numVertices; i++) {
- Valgrind::CheckDefined(residual_[i]);
- }
-#endif // HAVE_VALGRIND
+ for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
+ Valgrind::CheckDefined(residual[i]);
+#endif // !defined NDEBUG && HAVE_VALGRIND
// evaluate the boundary conditions
- asImp_().evalBoundary_();
+ asImp_().evalBoundary_(residual, storageTerm, elemVars);
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
+ for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
+ Valgrind::CheckDefined(residual[i]);
#endif // HAVE_VALGRIND
}
/*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
+ * \brief Calculate the amount of all conservation quantities
+ * stored in all element's sub-control volumes for a given
+ * history index.
+ *
+ * This is used to figure out how much of each conservation
+ * quantity is inside the element.
*/
- const ElementSolutionVector &residual() const
- { return residual_; }
+ void evalStorage(const ElementVariables &elemVars,
+ int historyIdx)
+ {
+ int numScv = elemVars.numScv();
+ internalStorageTerm_.resize(numScv);
+ evalStorage(internalStorageTerm_,
+ elemVars,
+ historyIdx);
+ }
/*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
+ * \brief Calculate the amount of all conservation quantities
+ * stored in all element's sub-control volumes for a given
+ * history index.
*
- * \param scvIdx The local index of the sub-control volume
- * (i.e. the element's local vertex index)
+ * This is used to figure out how much of each conservation
+ * quantity is inside the element.
*/
- const PrimaryVariables &residual(int scvIdx) const
- { return residual_[scvIdx]; }
+ void evalStorage(LocalBlockVector &storage,
+ const ElementVariables &elemVars,
+ int historyIdx) const
+ {
+ // calculate the amount of conservation each quantity inside
+ // all sub control volumes
+ for (int scvIdx=0; scvIdx < elemVars.numScv(); scvIdx++)
+ {
+ storage[scvIdx] = 0.0;
+ asImp_().computeStorage(storage[scvIdx],
+ elemVars,
+ scvIdx,
+ historyIdx);
+ storage[scvIdx] *=
+ elemVars.fvElemGeom().subContVol[scvIdx].volume
+ * elemVars.volVars(scvIdx).extrusionFactor();
+ }
+ }
/*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
+ * \brief Add the flux term to a local residual.
*/
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
-
-protected:
- Implementation &asImp_()
+ void evalFluxes(LocalBlockVector &residual,
+ const ElementVariables &elemVars) const
{
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
+ PrimaryVariables flux;
+
+ // calculate the mass flux over the sub-control volume faces
+ for (int scvfIdx = 0;
+ scvfIdx < elemVars.fvElemGeom().numEdges;
+ scvfIdx++)
+ {
+ int i = elemVars.fvElemGeom().subContVolFace[scvfIdx].i;
+ int j = elemVars.fvElemGeom().subContVolFace[scvfIdx].j;
- const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
+ Valgrind::SetUndefined(flux);
+ asImp_().computeFlux(flux, /*context=*/elemVars, scvfIdx);
+ flux *= elemVars.fluxVars(scvfIdx).extrusionFactor();
+ Valgrind::CheckDefined(flux);
+
+ // The balance equation for a finite volume is given by
+ //
+ // 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 mass flux as 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;
+ }
}
+protected:
/*!
- * \brief Evaluate the boundary conditions
- * of the current element.
+ * \brief Evaluate the boundary conditions of an element.
*/
- void evalBoundary_()
+ void evalBoundary_(LocalBlockVector &residual,
+ LocalBlockVector &storageTerm,
+ const ElementVariables &elemVars) const
{
- if (bcTypes_().hasNeumann())
- asImp_().evalNeumann_();
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
- if (bcTypes_().hasDirichlet())
- asImp_().evalDirichlet_();
+ if (!elemVars.onBoundary()) {
+ return;
+ }
+
+ if (elemVars.hasNeumann())
+ asImp_().evalNeumann_(residual, elemVars);
+
+ if (elemVars.hasDirichlet())
+ asImp_().evalDirichlet_(residual, storageTerm, elemVars);
}
/*!
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
*/
- void evalDirichlet_()
+ void evalDirichlet_(LocalBlockVector &residual,
+ LocalBlockVector &storageTerm,
+ const ElementVariables &elemVars) const
{
PrimaryVariables tmp(0);
- for (int i = 0; i < fvElemGeom_().numVertices; ++i) {
- const BoundaryTypes &bcTypes = bcTypes_(i);
- if (! bcTypes.hasDirichlet())
+ for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
+ const BoundaryTypes &bcTypes = elemVars.boundaryTypes(scvIdx);
+ if (!bcTypes.hasDirichlet())
continue;
-
+
// ask the problem for the dirichlet values
- const VertexPointer vPtr = elem_().template subEntity<dim>(i);
- asImp_().problem_().dirichlet(tmp, *vPtr);
+ asImp_().computeDirichlet_(tmp, elemVars, scvIdx);
+
+ const PrimaryVariables &priVars =
+ elemVars.volVars(scvIdx).primaryVars();
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isDirichlet(eqIdx))
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+
assert(0 <= pvIdx && pvIdx < numEq);
Valgrind::CheckDefined(pvIdx);
- Valgrind::CheckDefined(curPrimaryVar_(i, pvIdx));
+ Valgrind::CheckDefined(priVars);
Valgrind::CheckDefined(tmp[pvIdx]);
- this->residual_[i][eqIdx] =
- curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
-
- this->storageTerm_[i][eqIdx] = 0.0;
+ residual[scvIdx][eqIdx] = priVars[pvIdx] - tmp[pvIdx];
+ storageTerm[scvIdx][eqIdx] = 0.0;
};
};
}
@@ -411,13 +301,16 @@ protected:
* \brief Add all Neumann boundary conditions to the local
* residual.
*/
- void evalNeumann_()
+ void evalNeumann_(LocalBlockVector &residual,
+ const ElementVariables &elemVars) const
{
- Dune::GeometryType geoType = elem_().geometry().type();
+ const Element &elem = elemVars.element();
+ Dune::GeometryType geoType = elem.geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
- IntersectionIterator isIt = gridView_().ibegin(elem_());
- const IntersectionIterator &endIt = gridView_().iend(elem_());
+ const GridView &gridView = elemVars.gridView();
+ IntersectionIterator isIt = gridView.ibegin(elem);
+ const IntersectionIterator &endIt = gridView.iend(elem);
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -432,112 +325,59 @@ protected:
faceVertIdx < numFaceVerts;
++faceVertIdx)
{
- int elemVertIdx = refElem.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
-
+ int scvIdx = refElem.subEntity(/*entityIdx=*/faceIdx,
+ /*entityCodim=*/1,
+ /*subEntityIdx=*/faceVertIdx,
+ /*subEntityCodim=*/dim);
+
int boundaryFaceIdx =
- fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ elemVars.fvElemGeom().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
- evalNeumannSegment_(isIt,
- elemVertIdx,
+ evalNeumannSegment_(residual,
+ elemVars,
+ *isIt,
+ scvIdx,
boundaryFaceIdx);
}
}
}
+ void computeDirichlet_(PrimaryVariables &values,
+ const ElementVariables &elemVars,
+ int scvIdx) const
+ { elemVars.problem().dirichlet(values, elemVars, scvIdx); }
+
+
/*!
* \brief Add Neumann boundary conditions for a single sub-control
* volume face to the local residual.
*/
- void evalNeumannSegment_(const IntersectionIterator &isIt,
+ void evalNeumannSegment_(LocalBlockVector &residual,
+ const ElementVariables &elemVars,
+ const Intersection &is,
int scvIdx,
- int boundaryFaceIdx)
+ int boundaryFaceIdx) const
{
// temporary vector to store the neumann boundary fluxes
- PrimaryVariables values(0.0);
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
-
+ const BoundaryTypes &bcTypes = elemVars.boundaryTypes(scvIdx);
+ PrimaryVariables values;
+
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
- problem_().boxSDNeumann(values,
- elem_(),
- fvElemGeom_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- curVolVars_());
- if (!Valgrind::CheckDefined(values))
- std::cerr << "undefined: " << values << "\n";
- values *=
- fvElemGeom_().boundaryFace[boundaryFaceIdx].area
- * curVolVars_(scvIdx).extrusionFactor();
+ elemVars.problem().neumann(values,
+ elemVars,
+ is,
+ scvIdx,
+ boundaryFaceIdx);
Valgrind::CheckDefined(values);
- residual_[scvIdx] += values;
- }
- }
-
- /*!
- * \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 k = 0; k < fvElemGeom_().numEdges; k++)
- {
- int i = fvElemGeom_().subContVolFace[k].i;
- int j = fvElemGeom_().subContVolFace[k].j;
- PrimaryVariables flux;
- Valgrind::SetUndefined(flux);
- this->asImp_().computeFlux(flux, k);
- Scalar extrusionFactor =
- (curVolVars_(i).extrusionFactor()
- + curVolVars_(j).extrusionFactor())
- / 2;
- flux *= extrusionFactor;
- Valgrind::CheckDefined(flux);
-
- // 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(fvElemGeom_().numVertices);
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- for (int i=0; i < fvElemGeom_().numVertices; i++) {
- this->asImp_().computeStorage(storageTerm_[i], i, /*isOldSol=*/false);
- storageTerm_[i] *=
- fvElemGeom_().subContVol[i].volume
- * curVolVars_(i).extrusionFactor();
+ values *=
+ elemVars.fvElemGeom().boundaryFace[boundaryFaceIdx].area
+ * elemVars.volVars(scvIdx, /*historyIdx*/0).extrusionFactor();
+ Valgrind::CheckDefined(values);
+ residual[scvIdx] += values;
}
}
@@ -546,15 +386,17 @@ protected:
* to the local residual of all sub-control volumes of the
* current element.
*/
- void evalVolumeTerms_()
+ void evalVolumeTerms_(LocalBlockVector &residual,
+ LocalBlockVector &storageTerm,
+ const ElementVariables &elemVars) const
{
+ PrimaryVariables tmp, tmp2;
+
// evaluate the volume terms (storage + source terms)
- for (int i=0; i < fvElemGeom_().numVertices; i++)
+ for (int scvIdx=0; scvIdx < elemVars.numScv(); scvIdx++)
{
Scalar extrusionFactor =
- curVolVars_(i).extrusionFactor();
-
- PrimaryVariables tmp(0);
+ elemVars.volVars(scvIdx, /*historyIdx=*/0).extrusionFactor();
// mass balance within the element. this is the
// $\frac{m}{\partial t}$ term if using implicit
@@ -562,162 +404,51 @@ protected:
//
// TODO (?): we might need a more explicit way for
// doing the time discretization...
- this->asImp_().computeStorage(storageTerm_[i], i, false);
- this->asImp_().computeStorage(tmp, i, true);
+ asImp_().computeStorage(tmp2,
+ elemVars,
+ scvIdx,
+ /*historyIdx=*/1);
+ asImp_().computeStorage(tmp,
+ elemVars,
+ scvIdx,
+ /*historyIdx=*/0);
- storageTerm_[i] -= tmp;
- storageTerm_[i] *=
- fvElemGeom_().subContVol[i].volume
- / problem_().timeManager().timeStepSize()
+ tmp -= tmp2;
+ tmp *=
+ elemVars.fvElemGeom().subContVol[scvIdx].volume
+ / elemVars.problem().timeManager().timeStepSize()
* extrusionFactor;
- residual_[i] += storageTerm_[i];
- // subtract the source term from the local rate
- this->asImp_().computeSource(tmp, i);
- tmp *= fvElemGeom_().subContVol[i].volume * extrusionFactor;
- residual_[i] -= tmp;
+ storageTerm[scvIdx] += tmp;
+ residual[scvIdx] += tmp;
+
+ // subtract the source term from the residual
+ asImp_().computeSource(tmp2, elemVars, scvIdx);
+ tmp2 *= elemVars.fvElemGeom().subContVol[scvIdx].volume * extrusionFactor;
+ residual[scvIdx] -= tmp2;
// make sure that only defined quantities were used
// to calculate the residual.
- Valgrind::CheckDefined(residual_[i]);
+ Valgrind::CheckDefined(storageTerm[scvIdx]);
+ 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 &elem_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry &fvElemGeom_() const
- {
- Valgrind::CheckDefined(fvElemGeomPtr_);
- return *fvElemGeomPtr_;
- }
-
- /*!
- * \brief Returns a reference to the primary variables of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &curPrimaryVars_(int i) const
- {
- return curVolVars_(i).primaryVars();
- }
-
- /*!
- * \brief Returns the j'th primary of the i'th sub-control volume
- * of the current element.
- */
- Scalar curPrimaryVar_(int i, int j) const
- {
- return curVolVars_(i).primaryVar(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
+private:
+ Implementation &asImp_()
{
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
+ assert(static_cast<Implementation*>(this) != 0);
+ return *static_cast<Implementation*>(this);
}
- /*!
- * \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
+ const Implementation &asImp_() const
{
- return bcTypes_()[i];
+ assert(static_cast<const Implementation*>(this) != 0);
+ return *static_cast<const Implementation*>(this);
}
-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_;
+ LocalBlockVector internalResidual_;
+ LocalBlockVector internalStorageTerm_;
};
}
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index 63625996a3a98be237fab81f12b0bb6749bf8f40..bfc52156f2a649e628d9a8a89b6d155a824bfb66 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -72,16 +72,18 @@ class BoxModel
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ historySize = GET_PROP_VALUE(TypeTag, PTAG(TimeDiscHistorySize)),
dim = GridView::dimension
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
@@ -93,6 +95,7 @@ class BoxModel
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
@@ -106,9 +109,7 @@ public:
* \brief The constructor.
*/
BoxModel()
- {
- enableHints_ = GET_PARAM(TypeTag, bool, EnableHints);
- }
+ { }
~BoxModel()
{ delete jacAsm_; }
@@ -123,9 +124,11 @@ public:
{
problemPtr_ = &prob;
+ updateBoundaryTypes_();
+
int nDofs = asImp_().numDofs();
- uCur_.resize(nDofs);
- uPrev_.resize(nDofs);
+ for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
+ solution_[historyIdx].resize(nDofs);
boxVolume_.resize(nDofs);
localJacobian_.init(problem_());
@@ -135,10 +138,10 @@ public:
asImp_().applyInitialSolution_();
// resize the hint vectors
- if (enableHints_) {
+ if (enableHints_()) {
int nVerts = gridView_().size(dim);
- curHints_.resize(nVerts);
- prevHints_.resize(nVerts);
+ for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
+ hints_[historyIdx].resize(nVerts);
hintsUsable_.resize(nVerts);
std::fill(hintsUsable_.begin(),
hintsUsable_.end(),
@@ -147,74 +150,67 @@ public:
// also set the solution of the "previous" time step to the
// initial solution.
- uPrev_ = uCur_;
+ solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
}
- void setHints(const Element &elem,
- ElementVolumeVariables &prevVolVars,
- ElementVolumeVariables &curVolVars) const
+ void setAllHints(ElementVariables &elemVars) const
{
- if (!enableHints_)
+ if (!enableHints_())
return;
- int n = elem.template count<dim>();
- prevVolVars.resize(n);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int globalIdx = vertexMapper().map(elem, i, dim);
-
- if (!hintsUsable_[globalIdx]) {
- curVolVars[i].setHint(NULL);
- prevVolVars[i].setHint(NULL);
- }
- else {
- curVolVars[i].setHint(&curHints_[globalIdx]);
- prevVolVars[i].setHint(&prevHints_[globalIdx]);
- }
- }
- };
+ for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
+ setHints(elemVars, historyIdx);
+ }
- void setHints(const Element &elem,
- ElementVolumeVariables &curVolVars) const
+ void setHints(ElementVariables &elemVars, int historyIdx) const
{
- if (!enableHints_)
+ if (!enableHints_())
return;
- int n = elem.template count<dim>();
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int globalIdx = vertexMapper().map(elem, i, dim);
-
- if (!hintsUsable_[globalIdx])
- curVolVars[i].setHint(NULL);
- else
- curVolVars[i].setHint(&curHints_[globalIdx]);
+ int numScv = elemVars.numScv();
+ for (int scvIdx = 0; scvIdx < numScv; ++scvIdx) {
+ setHintsScv(elemVars, historyIdx, scvIdx);
}
+ }
+
+ void setHintsScv(ElementVariables &elemVars, int historyIdx, int scvIdx) const
+ {
+ int globalIdx = vertexMapper().map(elemVars.element(), scvIdx, dim);
+
+ if (!hintsUsable_[globalIdx])
+ elemVars.volVars(scvIdx, historyIdx).setHint(NULL);
+ else
+ elemVars.volVars(scvIdx, historyIdx)
+ .setHint(&hints_[historyIdx][globalIdx]);
};
- void updatePrevHints()
+ void updateAllHints(const ElementVariables &elemVars) const
{
- if (!enableHints_)
+ if (!enableHints_())
return;
- prevHints_ = curHints_;
+ for (int historyIdx = 0; historyIdx < historySize; ++historyIdx) {
+ updateHints(elemVars, historyIdx);
+ }
};
- void updateCurHints(const Element &elem,
- const ElementVolumeVariables &ev) const
+ void updateHints(const ElementVariables &elemVars, int historyIdx) const
{
- if (!enableHints_)
+ if (!enableHints_())
return;
-
- for (int i = 0; i < ev.size(); ++i) {
- int globalIdx = vertexMapper().map(elem, i, dim);
- curHints_[globalIdx] = ev[i];
- if (!hintsUsable_[globalIdx])
- prevHints_[globalIdx] = ev[i];
+
+ for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
+ int globalIdx = vertexMapper().map(elemVars.element(), scvIdx, dim);
+ hints_[historyIdx][globalIdx] = elemVars.volVars(scvIdx, historyIdx);
hintsUsable_[globalIdx] = true;
}
};
+ void shiftHints()
+ {
+ for (int historyIdx = 0; historyIdx < historySize - 1; ++ historyIdx)
+ hints_[historyIdx + 1] = hints_[historyIdx];
+ };
/*!
* \brief Compute the global residual for an arbitrary solution
@@ -226,10 +222,10 @@ public:
Scalar globalResidual(SolutionVector &dest,
const SolutionVector &u)
{
- SolutionVector tmp(curSol());
- curSol() = u;
+ SolutionVector tmp(solution(/*historyIdx=*/0));
+ solution(/*historyIdx=*/0) = u;
Scalar res = globalResidual(dest);
- curSol() = tmp;
+ solution(/*historyIdx=*/0) = tmp;
return res;
}
@@ -243,14 +239,16 @@ public:
{
dest = 0;
+ ElementVariables elemVars(this->problem_());
ElementIterator elemIt = gridView_().template begin<0>();
const ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().eval(*elemIt);
+ elemVars.updateAll(*elemIt);
+ localResidual().eval(elemVars);
- for (int i = 0; i < elemIt->template count<dim>(); ++i) {
- int globalI = vertexMapper().map(*elemIt, i, dim);
- dest[globalI] += localResidual().residual(i);
+ for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
+ int globalI = vertexMapper().map(*elemIt, scvIdx, dim);
+ dest[globalI] += localResidual().residual(scvIdx);
}
};
@@ -279,11 +277,14 @@ public:
void globalStorage(PrimaryVariables &dest)
{
dest = 0;
-
+
+ ElementVariables elemVars(this->problem_());
ElementIterator elemIt = gridView_().template begin<0>();
const ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
- localResidual().evalStorage(*elemIt);
+ elemVars.updateFVElemGeom(*elemIt);
+ elemVars.updateScvVars(/*historyIdx=*/0);
+ localResidual().evalStorage(elemVars, /*historyIdx=*/0);
for (int i = 0; i < elemIt->template count<dim>(); ++i)
dest += localResidual().storageTerm()[i];
@@ -302,28 +303,16 @@ public:
{ return boxVolume_[globalIdx][0]; }
/*!
- * \brief Reference to the current solution as a block vector.
- */
- const SolutionVector &curSol() const
- { return uCur_; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- SolutionVector &curSol()
- { return uCur_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
+ * \brief Reference to the solution at a given history index as a block vector.
*/
- const SolutionVector &prevSol() const
- { return uPrev_; }
+ const SolutionVector &solution(int historyIdx) const
+ { return solution_[historyIdx]; }
/*!
- * \brief Reference to the previous solution as a block vector.
+ * \brief Reference to the solution at a given history index as a block vector.
*/
- SolutionVector &prevSol()
- { return uPrev_; }
+ SolutionVector &solution(int historyIdx)
+ { return solution_[historyIdx]; }
/*!
* \brief Returns the operator assembler for the global jacobian of
@@ -374,7 +363,8 @@ public:
*/
Scalar primaryVarWeight(int vertIdx, int pvIdx) const
{
- return 1.0/std::max(std::abs(this->prevSol()[vertIdx][pvIdx]), 1.0);
+ Scalar absPv = std::abs(this->solution(/*historyIdx=*/1)[vertIdx][pvIdx]);
+ return 1.0/std::max(absPv, 1.0);
}
/*!
@@ -419,8 +409,8 @@ public:
NewtonController &controller)
{
#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i)
- Valgrind::CheckDefined(curSol()[i]);
+ for (size_t i = 0; i < solution(/*historyIdx=*/0).size(); ++i)
+ Valgrind::CheckDefined(solution(/*historyIdx=*/0)[i]);
#endif // HAVE_VALGRIND
asImp_().updateBegin();
@@ -433,8 +423,8 @@ public:
asImp_().updateFailed();
#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i) {
- Valgrind::CheckDefined(curSol()[i]);
+ for (size_t i = 0; i < solution(/*historyIdx=*/0).size(); ++i) {
+ Valgrind::CheckDefined(solution(/*historyIdx=*/0)[i]);
}
#endif // HAVE_VALGRIND
@@ -448,7 +438,9 @@ public:
* which the actual model can overload.
*/
void updateBegin()
- { }
+ {
+ updateBoundaryTypes_();
+ }
/*!
@@ -469,8 +461,7 @@ public:
// Reset the current solution to the one of the
// previous time step so that we can start the next
// update at a physically meaningful solution.
- uCur_ = uPrev_;
- curHints_ = prevHints_;
+ hints_[/*historyIdx=*/0] = hints_[/*historyIdx=*/1];
jacAsm_->reassembleAll();
};
@@ -485,10 +476,10 @@ public:
void advanceTimeLevel()
{
// make the current solution the previous one.
- uPrev_ = uCur_;
- prevHints_ = curHints_;
+ solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
- updatePrevHints();
+ // shift the hints by one position in the history
+ shiftHints();
}
/*!
@@ -513,7 +504,7 @@ public:
void deserialize(Restarter &res)
{
res.template deserializeEntities<dim>(asImp_(), this->gridView_());
- prevSol() = curSol();
+ solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
}
/*!
@@ -538,7 +529,7 @@ public:
}
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << curSol()[vertIdx][eqIdx] << " ";
+ outstream << solution_[/*historyIdx=*/0][vertIdx][eqIdx] << " ";
}
};
@@ -560,7 +551,7 @@ public:
DUNE_THROW(Dune::IOError,
"Could not deserialize vertex "
<< vertIdx);
- instream >> curSol()[vertIdx][eqIdx];
+ instream >> solution_[/*historyIdx=*/0][vertIdx][eqIdx];
}
};
@@ -602,6 +593,19 @@ public:
jacAsm_->init(problem_());
}
+ /*!
+ * \brief Return a pointer to the BoundaryTypes for a given global
+ * vertex index or 0 if the vertex is not on the boundary.
+ */
+ const BoundaryTypes *boundaryTypes(int globalIdx) const
+ {
+ static BoundaryTypes dummy;
+ int bvertIdx = boundaryVertexIndex_[globalIdx];
+ if (bvertIdx < 0)
+ return &dummy;
+ return &boundaryTypes_[bvertIdx];
+ }
+
/*!
* \brief Update the weights of all primary variables within an
* element given the complete set of volume variables
@@ -609,8 +613,7 @@ public:
* \param element The DUNE codim 0 entity
* \param volVars All volume variables for the element
*/
- void updatePVWeights(const Element &element,
- const ElementVolumeVariables &volVars) const
+ void updatePVWeights(const ElementVariables &elemVars) const
{ };
/*!
@@ -641,10 +644,11 @@ public:
ScalarField* def[numEq];
ScalarField* delta[numEq];
ScalarField* x[numEq];
- for (int i = 0; i < numEq; ++i) {
- x[i] = writer.allocateManagedBuffer(numVertices);
- delta[i] = writer.allocateManagedBuffer(numVertices);
- def[i] = writer.allocateManagedBuffer(numVertices);
+ ScalarField* relError = writer.allocateManagedBuffer(numVertices);
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
+ x[pvIdx] = writer.allocateManagedBuffer(numVertices);
+ delta[pvIdx] = writer.allocateManagedBuffer(numVertices);
+ def[pvIdx] = writer.allocateManagedBuffer(numVertices);
}
VertexIterator vIt = this->gridView_().template begin<dim>();
@@ -652,17 +656,25 @@ public:
for (; vIt != vEndIt; ++ vIt)
{
int globalIdx = vertexMapper().map(*vIt);
- for (int i = 0; i < numEq; ++i) {
- (*x[i])[globalIdx] = u[globalIdx][i];
- (*delta[i])[globalIdx] = - deltaU[globalIdx][i];
- (*def[i])[globalIdx] = globalResid[globalIdx][i];
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
+ (*x[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
+ (*delta[pvIdx])[globalIdx] =
+ - deltaU[globalIdx][pvIdx];
+ (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
}
+
+ PrimaryVariables uOld(u[globalIdx]);
+ PrimaryVariables uNew(uOld - deltaU[globalIdx]);
+ (*relError)[globalIdx] = asImp_().relativeErrorVertex(globalIdx,
+ uOld,
+ uNew);
}
- for (int i = 0; i < numEq; ++i) {
- writer.attachVertexData(*x[i], (boost::format("x_%i")%i).str().c_str());
- writer.attachVertexData(*delta[i], (boost::format("delta_%i")%i).str().c_str());
- writer.attachVertexData(*def[i], (boost::format("defect_%i")%i).str().c_str());
+ writer.attachVertexData(*relError, "relative Error");
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
+ writer.attachVertexData(*x[pvIdx], (boost::format("x_%i")%pvIdx).str().c_str());
+ writer.attachVertexData(*delta[pvIdx], (boost::format("delta_%i")%pvIdx).str().c_str());
+ writer.attachVertexData(*def[pvIdx], (boost::format("defect_%i")%pvIdx).str().c_str());
}
asImp_().addOutputVtkFields(u, writer);
@@ -684,7 +696,7 @@ public:
void addOutputVtkFields(const SolutionVector &sol,
MultiWriter &writer)
{
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
// create the required scalar fields
unsigned numVertices = this->gridView_().size(dim);
@@ -716,6 +728,9 @@ public:
{ return problem_().gridView(); }
protected:
+ static bool enableHints_()
+ { return GET_PARAM(TypeTag, bool, EnableHints); }
+
/*!
* \brief A reference to the problem on which the model is applied.
*/
@@ -739,13 +754,83 @@ protected:
LocalResidual &localResidual_()
{ return localJacobian_.localResidual(); }
+ /*!
+ * \brief Updates the stuff which determines a vertex' or
+ * element's boundary type
+ */
+ void updateBoundaryTypes_()
+ {
+ // resize the vectors
+ boundaryVertexIndex_.resize(numDofs());
+ std::fill(boundaryVertexIndex_.begin(),
+ boundaryVertexIndex_.end(),
+ -1);
+
+ int numBoundaryVertices = 0;
+ ElementVariables elemVars(problem_());
+
+ // loop over all elements of the grid
+ ElementIterator elemIt = gridView_().template begin<0>();
+ const ElementIterator elemEndIt = gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ // do nothing if the element does not have boundary intersections
+ if (!elemIt->hasBoundaryIntersections())
+ continue;
+
+ // retrieve the reference element for the current element
+ const Element &elem = *elemIt;
+ Dune::GeometryType geoType = elem.geometry().type();
+ const ReferenceElement &refElem = ReferenceElements::general(geoType);
+
+ elemVars.updateFVElemGeom(*elemIt);
+
+ // loop over all intersections of the element
+ IntersectionIterator isIt = gridView_().ibegin(elem);
+ const IntersectionIterator &endIt = gridView_().iend(elem);
+ for (; isIt != endIt; ++isIt)
+ {
+ // do nothing if the face is _not_ on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ // loop over all vertices of the intersection
+ int faceIdx = isIt->indexInInside();
+ int numFaceVerts = refElem.size(faceIdx, 1, dim);
+ for (int faceVertIdx = 0;
+ faceVertIdx < numFaceVerts;
+ ++faceVertIdx)
+ {
+ // find the local element index of the face's
+ // vertex
+ int scvIdx = refElem.subEntity(/*entityIdx=*/faceIdx,
+ /*entityCodim=*/1,
+ /*subEntityIdx=*/faceVertIdx,
+ /*subEntityCodim=*/dim);
+ int globalIdx = vertexMapper().map(*elemIt, scvIdx, /*codim=*/dim);
+ if (boundaryVertexIndex_[globalIdx] >= 0)
+ continue; // vertex has already been visited
+
+ // add a BoundaryTypes object
+ if (boundaryTypes_.size() <= numBoundaryVertices)
+ boundaryTypes_.resize(numBoundaryVertices + 1);
+ BoundaryTypes &bTypes = boundaryTypes_[numBoundaryVertices];
+ boundaryVertexIndex_[globalIdx] = numBoundaryVertices;
+ ++numBoundaryVertices;
+
+ problem_().boundaryTypes(bTypes, elemVars, scvIdx);
+ } // loop over intersection's vertices
+ } // loop over intersections
+ } // loop over elements
+ };
+
/*!
* \brief Applies the initial solution for all vertices of the grid.
*/
void applyInitialSolution_()
{
// first set the whole domain to zero
- uCur_ = Scalar(0.0);
+ SolutionVector &uCur = solution(/*historyIdx=*/0);
+ uCur = Scalar(0.0);
boxVolume_ = Scalar(0.0);
FVElementGeometry fvElemGeom;
@@ -785,8 +870,8 @@ protected:
// will be the arithmetic mean.
initVal *= fvElemGeom.subContVol[scvIdx].volume;
boxVolume_[globalIdx] += fvElemGeom.subContVol[scvIdx].volume;
- uCur_[globalIdx] += initVal;
- Valgrind::CheckDefined(uCur_[globalIdx]);
+ uCur[globalIdx] += initVal;
+ Valgrind::CheckDefined(uCur[globalIdx]);
}
}
@@ -801,7 +886,7 @@ protected:
Dune::ForwardCommunication);
VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
- sumPVHandle(uCur_, vertexMapper());
+ sumPVHandle(uCur, vertexMapper());
gridView().communicate(sumPVHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
@@ -810,15 +895,14 @@ protected:
// divide all primary variables by the volume of their boxes
int n = gridView_().size(dim);
for (int i = 0; i < n; ++i) {
- uCur_[i] /= boxVolume(i);
+ uCur[i] /= boxVolume(i);
}
}
// the hint cache for the previous and the current volume
// variables
mutable std::vector<bool> hintsUsable_;
- mutable std::vector<VolumeVariables> curHints_;
- mutable std::vector<VolumeVariables> prevHints_;
+ mutable std::vector<VolumeVariables> hints_[historySize];
/*!
* \brief Returns whether messages should be printed
@@ -843,13 +927,13 @@ protected:
// cur is the current iterative solution, prev the converged
// solution of the previous time step
- SolutionVector uCur_;
- SolutionVector uPrev_;
+ SolutionVector solution_[historySize];
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
+ // all the index of the BoundaryTypes object for a vertex
+ std::vector<int> boundaryVertexIndex_;
+ std::vector<BoundaryTypes> boundaryTypes_;
-private:
- bool enableHints_;
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
};
}
diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index b8284c38f341b1f8a8fc5d49a4c94ea855524ceb..9a5fbaab9b1e5e503fcb0e0236ffcad75652af4c 100644
--- a/dumux/boxmodels/common/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -60,7 +60,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
@@ -139,6 +139,28 @@ public:
model().init(asImp_());
}
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param context The local context. Only the element() and the fvElemGeom() methods are valid at this point
+ * \param localIdx The local index of the entity neighboring the Dirichlet boundary.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ template <class Context>
+ DUNE_DEPRECATED
+ void boundaryTypes(BoundaryTypes &values,
+ const Context &context,
+ int localIdx) const
+ {
+ const auto vPtr = context.element().template subEntity<dim>(localIdx);
+
+ // forward it to the method which only takes the global coordinate
+ asImp_().boundaryTypes(values, *vPtr);
+ }
+
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
@@ -170,6 +192,27 @@ public:
"a boundaryTypes() method.");
}
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param context The local context. Only the element() and the fvElemGeom() methods are valid at this point
+ * \param localIdx The local index of the entity neighboring the Dirichlet boundary.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ template <class Context>
+ DUNE_DEPRECATED
+ void dirichlet(PrimaryVariables &values,
+ const Context &context,
+ int localIdx) const
+ {
+ const auto vPtr = context.element().template subEntity<dim>(localIdx);
+
+ // forward it to the method which only takes the global coordinate
+ asImp_().dirichlet(values, *vPtr);
+ }
/*!
* \brief Evaluate the boundary conditions for a dirichlet
@@ -180,6 +223,7 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
+ DUNE_DEPRECATED
void dirichlet(PrimaryVariables &values,
const Vertex &vertex) const
{
@@ -198,6 +242,7 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
+ DUNE_DEPRECATED
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
{
@@ -209,6 +254,25 @@ public:
"a dirichlet() method.");
}
+ template <class Context>
+ DUNE_DEPRECATED
+ // if you get an deprecated warning here, please use context
+ // objects to specify your problem!
+ void neumann(PrimaryVariables &priVars,
+ const Context &context,
+ const Intersection &is,
+ int localIdx,
+ int boundaryIndex) const
+ {
+ return asImp_().boxSDNeumann(priVars,
+ context.element(),
+ context.fvElemGeom(),
+ is,
+ localIdx,
+ boundaryIndex,
+ context);
+ };
+
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
@@ -234,7 +298,8 @@ public:
const Intersection &is,
int scvIdx,
int boundaryFaceIdx,
- const ElementVolumeVariables &elemVolVars) const
+ const ElementVariables &elemVars) const
+ DUNE_DEPRECATED
{
// forward it to the interface without the volume variables
asImp_().neumann(values,
@@ -265,6 +330,7 @@ public:
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
+ DUNE_DEPRECATED
{
// forward it to the interface with only the global position
asImp_().neumannAtPos(values, fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal);
@@ -282,6 +348,7 @@ public:
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
+ DUNE_DEPRECATED
{
// Throw an exception (there is no reasonable default value
// for Neumann conditions)
@@ -291,6 +358,21 @@ public:
"a neumannAtPos() method.");
}
+ template <class Context>
+ // if you get an deprecated warning here, please use context
+ // objects to specify your problem!
+ DUNE_DEPRECATED
+ void source(PrimaryVariables &priVars,
+ const Context &context,
+ int localIdx) const
+ {
+ return asImp_().boxSDSource(priVars,
+ context.element(),
+ context.fvElemGeom(),
+ localIdx,
+ context);
+ };
+
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
@@ -313,10 +395,11 @@ public:
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx,
- const ElementVolumeVariables &elemVolVars) const
+ const ElementVariables &elemVars) const
+ DUNE_DEPRECATED
{
// forward to solution independent, box specific interface
- asImp_().source(values, element, fvElemGeom, scvIdx);
+ asImp_().source(values, elemVars, scvIdx);
}
/*!
@@ -336,6 +419,7 @@ public:
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
+ DUNE_DEPRECATED
{
// forward to generic interface
asImp_().sourceAtPos(values, fvElemGeom.subContVol[scvIdx].global);
@@ -356,6 +440,7 @@ public:
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
+ DUNE_DEPRECATED
{
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
@@ -411,6 +496,28 @@ public:
* thought as pipes with a cross section of 1 m^2 and 2D problems
* are assumed to extend 1 m to the back.
*/
+ template <class Context>
+ // if you get an deprecated warning here, please use context
+ // objects to specify your problem!
+ DUNE_DEPRECATED
+ Scalar extrusionFactor(const Context &context,
+ int localIdx) const
+ {
+ return asImp_().boxExtrusionFactor(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+ };
+
+ /*!
+ * \brief Return how much the domain is extruded at a given sub-control volume.
+ *
+ * This means the factor by which a lower-dimensional (1D or 2D)
+ * entity needs to be expanded to get a full dimensional cell. The
+ * default is 1.0 which means that 1D problems are actually
+ * thought as pipes with a cross section of 1 m^2 and 2D problems
+ * are assumed to extend 1 m to the back.
+ */
+ DUNE_DEPRECATED
Scalar boxExtrusionFactor(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
@@ -768,7 +875,7 @@ public:
Scalar t = timeManager().time() + timeManager().timeStepSize();
createResultWriter_();
resultWriter_->beginWrite(t);
- model().addOutputVtkFields(model().curSol(), *resultWriter_);
+ model().addOutputVtkFields(model().solution(/*historyIdx=*/0), *resultWriter_);
asImp_().addOutputVtkFields();
resultWriter_->endWrite();
}
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index 466a281393a39b24f57f3278fa3ab2e2a0da051d..ae75bb159e69f9254a2047f5e42be85f6d10ea81 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -74,7 +74,7 @@ NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector
NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume
NEW_PROP_TAG(VolumeVariables); //!< The secondary variables within a sub-control volume
-NEW_PROP_TAG(ElementVolumeVariables); //!< The secondary variables of all sub-control volumes in an element
+NEW_PROP_TAG(ElementVariables); //!< The secondary variables of all sub-control volumes in an element
NEW_PROP_TAG(FluxVariables); //!< Data required to calculate a flux over a face
NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces (outflow)
@@ -121,12 +121,15 @@ NEW_PROP_TAG(EnableHints);
// mappers from local to global indices
-//! maper for vertices
+//! mapper for vertices
NEW_PROP_TAG(VertexMapper);
-//! maper for elements
+//! mapper for elements
NEW_PROP_TAG(ElementMapper);
-//! maper for degrees of freedom
+//! mapper for degrees of freedom
NEW_PROP_TAG(DofMapper);
+
+//! The history size required by the time discretization
+NEW_PROP_TAG(TimeDiscHistorySize);
}
}
diff --git a/dumux/boxmodels/common/boxpropertydefaults.hh b/dumux/boxmodels/common/boxpropertydefaults.hh
index f1517d599b050e9aa01eb02daa6e6af3fc25fcbd..c5253660718029c693d3857d5357434224733098 100644
--- a/dumux/boxmodels/common/boxpropertydefaults.hh
+++ b/dumux/boxmodels/common/boxpropertydefaults.hh
@@ -124,7 +124,7 @@ SET_TYPE_PROP(BoxModel, VolumeVariables, Dumux::BoxVolumeVariables<TypeTag>);
/*!
* \brief An array of secondary variable containers.
*/
-SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
+SET_TYPE_PROP(BoxModel, ElementVariables, Dumux::BoxElementVariables<TypeTag>);
/*!
* \brief Boundary types at a single degree of freedom.
@@ -185,6 +185,8 @@ SET_INT_PROP(BoxModel, LinearSolverMaxIterations, 250);
//! set number of equations of the mathematical model as default
SET_INT_PROP(BoxModel, LinearSolverBlockSize, GET_PROP_VALUE(TypeTag, PTAG(NumEq)));
+//! Set the history size of the time discretiuation to 2 (for implicit euler)
+SET_INT_PROP(BoxModel, TimeDiscHistorySize, 2);
} // namespace Properties
} // namespace Dumux
diff --git a/dumux/boxmodels/common/boxvolumevariables.hh b/dumux/boxmodels/common/boxvolumevariables.hh
index 626d7a9801129d1d7e1481966f329650742b52fe..cf85043fbc42f154f0837b096f54545d837b37b4 100644
--- a/dumux/boxmodels/common/boxvolumevariables.hh
+++ b/dumux/boxmodels/common/boxvolumevariables.hh
@@ -49,6 +49,7 @@ class BoxVolumeVariables
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVariables)) ElementVariables;
public:
// default constructor
@@ -118,15 +119,13 @@ public:
* phase state in the 2p2c model)
*/
void update(const PrimaryVariables &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
+ const ElementVariables &elemVars,
int scvIdx,
- bool isOldSol)
+ int historyIdx)
{
Valgrind::CheckDefined(priVars);
primaryVars_ = priVars;
- extrusionFactor_ = problem.boxExtrusionFactor(element, elemGeom, scvIdx);
+ extrusionFactor_ = elemVars.problem().extrusionFactor(elemVars, scvIdx);
}
/*!
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index 1bbf74c58004eb4382ac52ccd333ae44eac5bb05..a7d244a324b8585a0f3e273165420724cd8fec35 100644
--- a/dumux/common/timemanager.hh
+++ b/dumux/common/timemanager.hh
@@ -379,7 +379,7 @@ public:
if (verbose_) {
std::cout <<
boost::format("Time step %d done. Wall time:%.4g, time:%.4g, time step size:%.4g\n")
- %timeStepIndex()%timer.elapsed()%time()%dt;
+ %timeStepIndex()%double(timer.elapsed())%double(time())%double(dt);
}
// write restart file if mandated by the problem
diff --git a/dumux/linear/overlappingscalarproduct.hh b/dumux/linear/overlappingscalarproduct.hh
index 2f0360e938678d9022f6acc49d9333b37c7880fc..71bfaff5048ce2fb088c2219016e48d5802190fc 100644
--- a/dumux/linear/overlappingscalarproduct.hh
+++ b/dumux/linear/overlappingscalarproduct.hh
@@ -67,8 +67,8 @@ public:
double norm(const OverlappingBlockVector &x)
{
- double tmp = dot(x, x);
- return std::sqrt(tmp);
+ field_type tmp = std::sqrt(dot(x, x));
+ return tmp;
};
private:
diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh
index 28cf38036859772409466dfb3bda6bdfefc5868b..79e4cd8caf27a84114a62bbd492595c6e5969bd1 100644
--- a/dumux/material/spatialparameters/boxspatialparameters.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters.hh
@@ -72,20 +72,30 @@ public:
* \brief Function for defining the parameters needed by constitutive relationships (kr-Sw, pc-Sw, etc.).
*
* \return the material parameters object
- * \param element The element
*/
+ template <class Context>
+ DUNE_DEPRECATED
+ const MaterialLawParams& materialLawParams(const Context &context, int localIdx) const
+ {
+ return asImp_().materialLawParams(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+ }
+
+ DUNE_DEPRECATED
const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ const FVElementGeometry &fvElemGeom,
+ int localIdx) const
{
- return asImp_().materialLawParamsAtPos(element.geometry().center());
+ return asImp_().materialLawParams(fvElemGeom.subContVol[localIdx].global);
}
+
/*!
* \brief Function for defining the parameters needed by constitutive relationships (kr-Sw, pc-Sw, etc.).
*
- * \return the material parameters object
* \param globalPos The position of the center of the element
+ * \return the material parameters object
*/
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
diff --git a/dumux/material/spatialparameters/boxspatialparameters1p.hh b/dumux/material/spatialparameters/boxspatialparameters1p.hh
index a70a3a4ae56406a2fc0d5b1e0e99019af536efb1..fe7232fc781a5b104778f6707727c3c33ff0f45f 100644
--- a/dumux/material/spatialparameters/boxspatialparameters1p.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters1p.hh
@@ -66,6 +66,7 @@ class BoxSpatialParametersOneP
typedef typename GridView::ctype CoordScalar;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar,dimWorld> Vector;
public:
BoxSpatialParametersOneP(const GridView &gv)
@@ -74,53 +75,6 @@ public:
~BoxSpatialParametersOneP()
{}
- /*!
- * \brief Returns the factor by which the volume of a sub control
- * volume needs to be multiplied in order to get cubic
- * meters.
- *
- * \param element The current finite element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index sub-control volume face where the
- * factor ought to be calculated.
- *
- * By default that's just 1.0
- */
- Scalar extrusionFactorScv(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- DUNE_DEPRECATED // use the extrusion factor of the volume variables
- { return 1.0; }
-
- /*!
- * \brief Returns the factor by which the area of a sub control
- * volume face needs to be multiplied in order to get
- * square meters.
- *
- * \param element The current finite element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvfIdx The index sub-control volume face where the
- * factor ought to be calculated.
- *
- * By default it is the arithmetic mean of the extrusion factor of
- * the face's two sub-control volumes.
- */
- Scalar extrusionFactorScvf(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvfIdx) const
- DUNE_DEPRECATED // use the extrusion factor of the volume variables
- {
- return
- 0.5 *
- (asImp_().extrusionFactorScv(element,
- fvElemGeom,
- fvElemGeom.subContVolFace[scvfIdx].i)
- +
- asImp_().extrusionFactorScv(element,
- fvElemGeom,
- fvElemGeom.subContVolFace[scvfIdx].j));
- }
-
/*!
* \brief Averages the intrinsic permeability (Scalar).
* \param result averaged intrinsic permeability
@@ -156,15 +110,26 @@ public:
result[i][j] = harmonicMean(K1[i][j], K2[i][j]);
}
+ template <class Context>
+ const Tensor &intrinsicPermeability(const Context &context,
+ int localIdx) const
+ {
+ return intrinsicPermeability(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+
+ }
+
/*!
* \brief Function for defining the intrinsic (absolute) permeability.
*
* \return intrinsic (absolute) permeability
* \param element The element
*/
- const Tensor& intrinsicPermeability (const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ DUNE_DEPRECATED
+ const Tensor &intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
return asImp_().intrinsicPermeabilityAtPos(element.geometry().center());
}
@@ -175,6 +140,7 @@ public:
* \return intrinsic (absolute) permeability
* \param globalPos The position of the center of the element
*/
+ DUNE_DEPRECATED
const Tensor& intrinsicPermeabilityAtPos (const GlobalPosition& globalPos) const
{
DUNE_THROW(Dune::InvalidStateException,
@@ -182,15 +148,26 @@ public:
"a intrinsicPermeabilityAtPos() method.");
}
+ template <class Context>
+ Scalar porosity(const Context &context,
+ int localIdx) const
+ {
+ return porosity(context.element(),
+ context.fvElemGeom(),
+ localIdx);
+
+ }
+
/*!
* \brief Function for defining the porosity.
*
* \return porosity
* \param element The element
*/
+ DUNE_DEPRECATED
Scalar porosity(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
return asImp_().porosityAtPos(element.geometry().center());
}
@@ -201,6 +178,7 @@ public:
* \return porosity
* \param globalPos The position of the center of the element
*/
+ DUNE_DEPRECATED
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
DUNE_THROW(Dune::InvalidStateException,
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index a083ef605c76f3e23f4067e2378b859b1380a60d..2165c39f8b8f10ad9032585101681148b13e230e 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -116,7 +116,7 @@ public:
protected:
bool execute_(NewtonController &ctl)
{
- SolutionVector &uCurrentIter = model().curSol();
+ SolutionVector &uCurrentIter = model().solution(/*historyIdx=*/0);
SolutionVector uLastIter(uCurrentIter);
SolutionVector deltaU(uCurrentIter);
diff --git a/test/boxmodels/2p/lensproblem.hh b/test/boxmodels/2p/lensproblem.hh
index dab4dc78d7b4ae513a3526ccbcceb3453ff43b77..d8bc6f25412d514a453bb6b655d306444fda66ec 100644
--- a/test/boxmodels/2p/lensproblem.hh
+++ b/test/boxmodels/2p/lensproblem.hh
@@ -25,10 +25,11 @@
* \brief Soil contamination problem where DNAPL infiltrates a fully
* water saturated medium.
*/
-
#ifndef DUMUX_LENSPROBLEM_HH
#define DUMUX_LENSPROBLEM_HH
+//#include <dumux/common/quad.hh>
+
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
@@ -57,6 +58,8 @@ namespace Properties
{
NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP, LensSpatialParameters));
+//SET_TYPE_PROP(LensProblem, Scalar, quad);
+
// Set the grid type
#if 0// HAVE_UG
SET_TYPE_PROP(LensProblem, Grid, Dune::UGGrid<2>);
@@ -243,6 +246,7 @@ public:
*
* This problem assumes a uniform temperature of 10 degrees Celsius.
*/
+ using ParentType::temperature;
Scalar temperature() const
{ return temperature_; };
diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh
index d0a34f0f153adba60b1a2aa9e2c56d045399e6de..ad82f044eac3d4e8c74eb3815bcbba5b81b07787 100644
--- a/test/boxmodels/2p/lensspatialparameters.hh
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -135,27 +135,24 @@ public:
* \param scvIdx The index sub-control volume face where the
* intrinsic velocity ought to be calculated.
*/
- Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ template <class Context>
+ Scalar intrinsicPermeability(const Context &context, int localIdx) const
{
- const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ const GlobalPosition &globalPos = context.pos(localIdx);
if (isInLens_(globalPos))
return lensK_;
return outerK_;
}
- Scalar porosity(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ template <class Context>
+ Scalar porosity(const Context &context, int localIdx) const
{ return 0.4; }
// return the parameter object for the Brooks-Corey material law which depends on the position
- const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ template <class Context>
+ const MaterialLawParams& materialLawParams(const Context &context, int localIdx) const
{
- const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ const GlobalPosition &globalPos = context.pos(localIdx);
if (isInLens_(globalPos))
return lensMaterialParams_;