diff --git a/dumux/boxmodels/2p/2pfluxvariables.hh b/dumux/boxmodels/2p/2pfluxvariables.hh
index 082d62298db5175332c4ea8370ac2e9e309516cd..4ea7096c07d4c285f68fd0dc19ac8e1bcdabc179 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
enum {
dim = GridView::dimension,
@@ -77,90 +77,100 @@ 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()
- {}
-
-#warning Docme
- void update(const ElementVariables &elemVars, int scvfIdx)
+ TwoPFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int faceIdx,
+ const ElementVolumeVariables &elemDat)
+ : fvElemGeom_(elemGeom)
{
- insideScvIdx_ = elemVars.fvElemGeom().subContVolFace[scvfIdx].i;
- outsideScvIdx_ = elemVars.fvElemGeom().subContVolFace[scvfIdx].j;
+ scvfIdx_ = faceIdx;
+
+ for (int phase = 0; phase < numPhases; ++phase) {
+ potentialGrad_[phase] = Scalar(0);
+ }
- calculateGradients_(elemVars, scvfIdx);
- calculateNormalFlux_(elemVars, scvfIdx);
+ calculateGradients_(problem, element, elemDat);
+ calculateK_(problem, element, elemDat);
};
- /*!
- * \brief Return the extrusion factor of the SCVF.
+public:
+ /*
+ * \brief Return the intrinsic permeability.
*/
- Scalar extrusionFactor() const
- { return 1.0; }
+ const Tensor &intrinsicPermeability() const
+ { return K_; }
/*!
- * \brief Return a phase's pressure potential gradient.
+ * \brief Return the 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 phaseIdx The index of the fluid phase for which the downstream
- * direction is requested.
+ * \param normalFlux The normal flux i.e. the given intrinsic permeability
+ * times the pressure potential gradient and SCV face normal.
*/
- int downstreamIdx(int phaseIdx) const
- { return (normalFlux_[phaseIdx] > 0)?outsideScvIdx_:insideScvIdx_; }
+ int downstreamIdx(Scalar normalFlux) const
+ { return (normalFlux >= 0)?face().j:face().i; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase as a function of the normal flux.
*
- * \param phaseIdx The index of the fluid phase for which the upstream
- * direction is requested.
+ * \param normalFlux The normal flux i.e. the given intrinsic permeability
+ * times the pressure potential gradient and SCV face normal.
*/
- int upstreamIdx(int phaseIdx) const
- { return (normalFlux_[phaseIdx] > 0)?insideScvIdx_:outsideScvIdx_; }
+ 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_]; }
protected:
- 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];
+ 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)
+ {
// calculate gradients
- for (int scvIdx = 0;
- scvIdx < elemVars.numScv();
- scvIdx ++) // loop over adjacent vertices
+ for (int idx = 0;
+ idx < fvElemGeom_.numVertices;
+ idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
- const Vector &feGrad = scvf.grad[scvIdx];
+ const Vector &feGrad = face().grad[idx];
// compute sum of pressure gradients for each phase
for (int phase = 0; phase < numPhases; phase++)
{
// the pressure gradient
Vector tmp(feGrad);
- tmp *= elemVars.volVars(scvIdx, /*historyIdx=*/0).pressure(phase);
+ tmp *= elemVolVars[idx].pressure(phase);
potentialGrad_[phase] += tmp;
}
}
@@ -172,22 +182,18 @@ protected:
{
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
- Vector g(elemVars.problem().boxGravity(elemVars.element(),
- elemVars.fvElemGeom(),
- insideScvIdx_));
- g += elemVars.problem().boxGravity(elemVars.element(),
- elemVars.fvElemGeom(),
- outsideScvIdx_);
+ Vector g(problem.boxGravity(element, fvElemGeom_, face().i));
+ g += problem.boxGravity(element, fvElemGeom_, face().j);
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 = 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 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 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)
@@ -206,52 +212,20 @@ protected:
}
}
- void calculateNormalFlux_(const ElementVariables &elemVars,
- int scvfIdx)
+ void calculateK_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
{
- const SpatialParameters &spatialParams = elemVars.problem().spatialParameters();
-
+ const SpatialParameters &spatialParams = problem.spatialParameters();
// calculate the intrinsic permeability
- 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;
- }
+ spatialParams.meanK(K_,
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().i),
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().j));
}
-
- // 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 8c45c76ba8bc1a42b0893da20ed72b049a0c8ba4..0b554f26fdb2beb559e25dff62c57e86428bf911 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
@@ -107,26 +107,23 @@ 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,
- const ElementVariables &elemVars,
- int scvIdx,
- int historyIdx) const
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
{
- // retrieve the volume variables for the SCV at the specified
- // point in time
- const VolumeVariables &volVars = elemVars.volVars(scvIdx, historyIdx);
+ // 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];
// wetting phase mass
- result[contiWEqIdx] =
- volVars.density(wPhaseIdx)
- * volVars.porosity()
- * volVars.saturation(wPhaseIdx);
+ result[contiWEqIdx] = vertDat.density(wPhaseIdx) * vertDat.porosity()
+ * vertDat.saturation(wPhaseIdx);
// non-wetting phase mass
- result[contiNEqIdx] =
- volVars.density(nPhaseIdx)
- * volVars.porosity()
- * volVars.saturation(nPhaseIdx);
+ result[contiNEqIdx] = vertDat.density(nPhaseIdx) * vertDat.porosity()
+ * vertDat.saturation(nPhaseIdx);
}
/*!
@@ -136,13 +133,16 @@ 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,
- const ElementVariables &elemVars,
- int scvfIdx) const
+ void computeFlux(PrimaryVariables &flux, int faceIdx) const
{
+ FluxVariables vars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_());
flux = 0;
- asImp_()->computeAdvectiveFlux(flux, elemVars, scvfIdx);
- asImp_()->computeDiffusiveFlux(flux, elemVars, scvfIdx);
+ asImp_()->computeAdvectiveFlux(flux, vars);
+ asImp_()->computeDiffusiveFlux(flux, vars);
}
/*!
@@ -155,40 +155,31 @@ 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 ElementVariables &elemVars,
- int scvfIdx) const
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) 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)
{
- // 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
+ // calculate the flux in the normal direction of the
// current sub control volume face:
//
- // v_\alpha := - (K grad p) * n
+ // v = - (K grad p) * n
//
- // 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);
+ // (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));
// add advective flux of current component in current
// phase
@@ -213,11 +204,10 @@ public:
* non-isothermal two-phase models to calculate diffusive heat
* fluxes
*/
- void computeDiffusiveFlux(PrimaryVariables &flux,
- const ElementVariables &elemVars,
- int scvfIdx) const
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxData) const
{
- // no diffusive fluxes for immiscible models
+ // diffusive fluxes
+ flux += 0.0;
}
/*!
@@ -227,12 +217,14 @@ public:
* \param localVertexIdx The index of the SCV
*
*/
- void computeSource(PrimaryVariables &values,
- const ElementVariables &elemVars,
- int scvIdx) const
+ void computeSource(PrimaryVariables &q, int localVertexIdx) const
{
// retrieve the source term intrinsic to the problem
- elemVars.problem().source(values, elemVars, scvIdx);
+ this->problem_().boxSDSource(q,
+ this->elem_(),
+ this->fvElemGeom_(),
+ localVertexIdx,
+ this->curVolVars_());
}
diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 8c2126b9a3f41d834be180f831b9c798f75174f9..604d538d0d87294433bc5962c4fa453089356518 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
enum {
@@ -122,11 +122,8 @@ public:
*/
Scalar primaryVarWeight(int globalVertexIdx, int pvIdx) const
{
- if (Indices::pressureIdx == pvIdx) {
- Scalar absPv =
- std::abs(this->solution(/*historyIdx=*/1)[globalVertexIdx][pvIdx]);
- return std::min(10.0/absPv, 1.0);
- }
+ if (Indices::pressureIdx == pvIdx)
+ return std::min(10.0/this->prevSol()[globalVertexIdx][pvIdx], 1.0);
return 1;
}
@@ -176,7 +173,9 @@ public:
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
- ElementVariables elemVars(this->problem_());
+ FVElementGeometry fvElemGeom;
+ VolumeVariables volVars;
+ ElementVolumeVariables elemVolVars;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
@@ -189,14 +188,19 @@ public:
int idx = this->elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
- elemVars.updateFVElemGeom(*elemIt);
- elemVars.updateScvVars(/*historyIdx=*/0);
+ fvElemGeom.update(this->gridView_(), *elemIt);
- for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx)
+ int numVerts = elemIt->template count<dim> ();
+ for (int i = 0; i < numVerts; ++i)
{
- int globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
- const VolumeVariables &volVars = elemVars.volVars(scvIdx, /*historyIdx=*/0);
(*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
(*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
(*pC)[globalIdx] = volVars.capillaryPressure();
@@ -214,7 +218,6 @@ public:
}
};
-#if 0
if(velocityOutput)
{
// calculate vertex velocities
@@ -270,7 +273,6 @@ 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);
@@ -329,9 +331,7 @@ 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,13 +352,11 @@ 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 bb30c5d275166b2197d61e3679a6a4d0ef24a764..6981be5acc2f130a392677a9fa65c720929c1003 100644
--- a/dumux/boxmodels/2p/2pproblem.hh
+++ b/dumux/boxmodels/2p/2pproblem.hh
@@ -111,26 +111,6 @@ 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.
*
@@ -145,7 +125,6 @@ public:
Scalar boxTemperature(const Element &element,
const FVElementGeometry fvGeom,
int scvIdx) const
- DUNE_DEPRECATED // use context objects!
{ return asImp_().temperatureAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
@@ -157,7 +136,6 @@ 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(); }
/*!
@@ -168,30 +146,8 @@ 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$.
*
@@ -199,9 +155,8 @@ public:
* it just calls gravityAtPos().
*/
const Vector &boxGravity(const Element &element,
- const FVElementGeometry &fvGeom,
- int scvIdx) const
- DUNE_DEPRECATED // use context objects!
+ const FVElementGeometry &fvGeom,
+ int scvIdx) const
{ return asImp_().gravityAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index 52391ab80b82605d9c05af62eae5289ddd974dea..8612933a7e4ad1a204ccd5deb2d18fe7f2a60262 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -56,8 +56,6 @@ 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 {
@@ -91,26 +89,28 @@ public:
* \param isOldSol Evaluate function with solution of current or previous time step
*/
void update(const PrimaryVariables &priVars,
- const ElementVariables &elemVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
int scvIdx,
- int historyIdx)
+ bool isOldSol)
{
ParentType::update(priVars,
- elemVars,
+ problem,
+ element,
+ elemGeom,
scvIdx,
- historyIdx);
+ isOldSol);
asImp().updateTemperature_(priVars,
- elemVars,
+ element,
+ elemGeom,
scvIdx,
- historyIdx);
+ problem);
- const SpatialParameters &spatialParams =
- elemVars.problem().spatialParameters();
-
- // material law parameters
+ // material law parameters
const MaterialLawParams &materialParams =
- spatialParams.materialLawParams(elemVars, scvIdx);
+ problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
Scalar p[numPhases];
Scalar Sn;
@@ -147,7 +147,9 @@ public:
fluidState_);
// porosity
- porosity_ = spatialParams.porosity(elemVars, scvIdx);
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
}
/*!
@@ -216,12 +218,12 @@ public:
protected:
void updateTemperature_(const PrimaryVariables &priVars,
- const ElementVariables &elemVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
int scvIdx,
- int historyIdx)
+ const Problem &problem)
{
- temperature_ =
- elemVars.problem().temperature(elemVars, scvIdx);
+ temperature_ = problem.boxTemperature(element, elemGeom, scvIdx);
}
FluidState fluidState_;
diff --git a/dumux/boxmodels/common/boxassembler.hh b/dumux/boxmodels/common/boxassembler.hh
index ef31386e032977419c8bd9438fe49eddef70a7b5..045ad46f77fd5b9033efc95b492d44fb60eb403e 100644
--- a/dumux/boxmodels/common/boxassembler.hh
+++ b/dumux/boxmodels/common/boxassembler.hh
@@ -113,6 +113,9 @@ public:
BoxAssembler()
{
+ enableJacobianRecycling_ = GET_PARAM(TypeTag, bool, EnableJacobianRecycling);
+ enablePartialReassemble_ = GET_PARAM(TypeTag, bool, EnablePartialReassemble);
+
problemPtr_ = 0;
matrix_ = 0;
@@ -156,7 +159,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);
}
@@ -164,7 +167,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);
@@ -173,14 +176,14 @@ public:
}
/*!
- * \brief Assemble the global Jacobian of the residual and the residual for the current solution.
+ * \brief Assemble the local jacobian of the problem.
*
* 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_();
@@ -228,7 +231,7 @@ public:
*/
void setMatrixReuseable(bool yesno = true)
{
- if (enableJacobianRecycling_())
+ if (enableJacobianRecycling_)
reuseMatrix_ = yesno;
}
@@ -244,7 +247,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);
@@ -280,7 +283,7 @@ public:
void updateDiscrepancy(const SolutionVector &u,
const SolutionVector &uDelta)
{
- if (!enablePartialReassemble_())
+ if (!enablePartialReassemble_)
return;
// update the vector with the distances of the current
@@ -339,7 +342,7 @@ public:
*/
void computeColors(Scalar relTol)
{
- if (!enablePartialReassemble_())
+ if (!enablePartialReassemble_)
return;
ElementIterator elemIt = gridView_().template begin<0>();
@@ -483,7 +486,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);
@@ -497,7 +500,7 @@ public:
*/
int vertexColor(int globalVertIdx) const
{
- if (!enablePartialReassemble_())
+ if (!enablePartialReassemble_)
return Red; // reassemble unconditionally!
return vertexColor_[globalVertIdx];
}
@@ -509,7 +512,7 @@ public:
*/
int elementColor(const Element &element) const
{
- if (!enablePartialReassemble_())
+ if (!enablePartialReassemble_)
return Red; // reassemble unconditionally!
int globalIdx = elementMapper_().map(element);
@@ -523,7 +526,7 @@ public:
*/
int elementColor(int globalElementIdx) const
{
- if (!enablePartialReassemble_())
+ if (!enablePartialReassemble_)
return Red; // reassemble unconditionally!
return elementColor_[globalElementIdx];
}
@@ -540,27 +543,8 @@ 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(); }
+private:
// Construct the BCRS matrix for the global jacobian
void createMatrix_()
{
@@ -639,7 +623,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;
@@ -756,20 +740,21 @@ private:
residual_[globI] += model_().localJacobian().residual(i);
if (enableJacobianRecycling_) {
- // save storage term and its Jacobian in case we can
- // reuse the current linearization...
+ // save the flux term and the jacobian of the
+ // storage term in case we can reuse the current
+ // linearization...
storageJacobian_[globI] +=
- model_().localJacobian().jacobianStorage(i);
+ model_().localJacobian().storageJacobian(i);
storageTerm_[globI] +=
- model_().localJacobian().residualStorage(i);
+ model_().localJacobian().storageTerm(i);
}
// update the jacobian matrix
for (int j=0; j < numVertices; ++ j) {
int globJ = vertexMapper_().map(elem, j, dim);
(*matrix_)[globI][globJ] +=
- model_().localJacobian().jacobian(i, j);
+ model_().localJacobian().mat(i,j);
}
}
}
@@ -800,6 +785,22 @@ 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
@@ -826,6 +827,9 @@ 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 63ac9daf88d3e2406c314304cf51db6e79da28a6..8e692ab0d072b8e8e4d6446f7f3ac915241082a6 100644
--- a/dumux/boxmodels/common/boxelementvolumevariables.hh
+++ b/dumux/boxmodels/common/boxelementvolumevariables.hh
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2010-2011 by Andreas Lauser *
+ * Copyright (C) 2010 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_VARIABLES_HH
-#define DUMUX_BOX_ELEMENT_VARIABLES_HH
+#ifndef DUMUX_BOX_ELEMENT_VOLUME_VARIABLES_HH
+#define DUMUX_BOX_ELEMENT_VOLUME_VARIABLES_HH
#include "boxproperties.hh"
@@ -38,28 +38,16 @@ namespace Dumux
* volume variables object for each of the element's vertices
*/
template<class TypeTag>
-class BoxElementVariables
+class BoxElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- 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 std::vector<VolumeVariables> ParentType;
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;
@@ -69,107 +57,47 @@ class BoxElementVariables
enum { dim = GridView::dimension };
enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
- typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
-
public:
/*!
* \brief The constructor.
*/
- explicit BoxElementVariables(const Problem &problem)
- : gridView_(problem.gridView())
- {
- // remember the problem object
- problemPtr_ = &problem;
- modelPtr_ = &problem.model();
- }
+ BoxElementVolumeVariables()
+ { }
/*!
- * \brief Construct the volume variables of an element from scratch.
+ * \brief Construct the volume variables for all of vertices of an element.
*
* \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 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()
+ void update(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ bool oldSol)
{
- // resize the SCV and the SCVF arrays
- if (fvElemGeom_.numVertices > scvVars_.size()) {
- scvVars_.resize(fvElemGeom_.numVertices);
- boundaryTypes_.resize(fvElemGeom_.numVertices);
+ 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);
}
- 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.
@@ -184,216 +112,27 @@ public:
* \param elementSolVector The local solution for the element using PDELab ordering
*/
template<typename ElemSolVectorType>
- void updatePDELab(const Element &element,
- const ElemSolVectorType &elementSolVector)
+ void updatePDELab(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const ElemSolVectorType& elementSolVector)
{
- updateFVElemGeom(element);
- updateBoundaryTypes(element);
-
- // update the current time step's volume variables
- PrimaryVariables scvSol;
- for (int scvIdx = 0; scvIdx < numScv(); scvIdx++)
+ int n = element.template count<dim>();
+ this->resize(n);
+ for (int vertexIdx = 0; vertexIdx < n; vertexIdx++)
{
- // reorder the solution
- for (int eqnIdx = 0; eqnIdx < numEq; eqnIdx++)
- scvSol[eqnIdx] = elementSolVector[scvIdx + eqnIdx*numScv()];
+ 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);
- // 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 05fee5964a1c63a6dcf78e22c55fb9cbb0720c4b..9ccf1ffed81e285d28b74674d5e9d9fcff0eeab8 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -30,7 +30,6 @@
#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"
@@ -78,12 +77,17 @@ 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,
};
@@ -109,30 +113,23 @@ 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
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()
- {
- internalElemVars_ = 0;
- }
-
- ~BoxLocalJacobian()
{
- delete internalElemVars_;
+ numericDifferenceMethod_ = GET_PARAM(TypeTag, int, NumericDifferenceMethod);
+ Valgrind::SetUndefined(problemPtr_);
}
+
/*!
* \brief Initialize the local Jacobian object.
*
@@ -142,67 +139,140 @@ public:
* \param prob The problem which we want to simulate.
*/
void init(Problem &prob)
- {
+ {
problemPtr_ = &prob;
- modelPtr_ = &prob.model();
- internalElemVars_ = new ElementVariables(prob);
+ localResidual_.init(prob);
+
+ // assume quadrilinears as elements with most vertices
+ A_.setSize(2<<dim, 2<<dim);
+ storageJacobian_.resize(2<<dim);
}
/*!
* \brief Assemble an element's local Jacobian matrix of the
* defect.
*
- * This assembles the 'grad f(x^k)' and 'f(x^k)' part of the newton update
+ * \param element The DUNE Codim<0> entity which we look at.
*/
void assemble(const Element &element)
{
- internalElemVars_->updateAll(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_);
- 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(residual_, residualStorage_, elemVars);
+ localResidual().eval(elem_(),
+ fvElemGeom_,
+ prevVolVars_,
+ curVolVars_,
+ bcTypes_);
+ residual_ = localResidual().residual();
+ storageTerm_ = localResidual().storageTerm();
- // save all flux variables calculated using the unmodified
- // primary variables. This automatically makes these flux
- // variables the evaluation point.
- elemVars.saveScvfVars();
+ model_().updatePVWeights(elem_(), curVolVars_);
// calculate the local jacobian matrix
- int numScv = elemVars.numScv();
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
+ int numVertices = fvElemGeom_.numVertices;
+ ElementSolutionVector partialDeriv(numVertices);
+ PrimaryVariables storageDeriv;
+ for (int j = 0; j < numVertices; j++) {
for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
- asImp_().evalPartialDerivative_(elemVars,
- scvIdx,
+ asImp_().evalPartialDerivative_(partialDeriv,
+ storageDeriv,
+ j,
pvIdx);
- // update the local stiffness matrix with the current
- // partial derivatives
- updateLocalJacobian_(elemVars, scvIdx, pvIdx);
+ // update the local stiffness matrix with the current partial
+ // derivatives
+ updateLocalJacobian_(j,
+ pvIdx,
+ partialDeriv,
+ storageDeriv);
}
}
-
- // restore flux variables.
- //elemVars.restoreScvfVars(); // not necessary
}
+ /*!
+ * \brief Returns a reference to the object which calculates the
+ * local residual.
+ */
+ const LocalResidual &localResidual() const
+ { return localResidual_; }
+
+ /*!
+ * \brief Returns a reference to the object which calculates the
+ * local residual.
+ */
+ LocalResidual &localResidual()
+ { return localResidual_; }
+
+ /*!
+ * \brief Returns the Jacobian of the equations at vertex i to the
+ * primary variables at vertex j.
+ *
+ * \param i The local vertex (or sub-control volume) index on which
+ * the equations are defined
+ * \param j The local vertex (or sub-control volume) index which holds
+ * primary variables
+ */
+ const MatrixBlock &mat(int i, int j) const
+ { return A_[i][j]; }
+
+ /*!
+ * \brief Returns the Jacobian of the storage term at vertex i.
+ *
+ * \param i The local vertex (or sub-control volume) index
+ */
+ const MatrixBlock &storageJacobian(int i) const
+ { return storageJacobian_[i]; }
+
+ /*!
+ * \brief Returns the residual of the equations at vertex i.
+ *
+ * \param i The local vertex (or sub-control volume) index on which
+ * the equations are defined
+ */
+ const PrimaryVariables &residual(int i) const
+ { return residual_[i]; }
+
+ /*!
+ * \brief Returns the storage term for vertex i.
+ *
+ * \param i The local vertex (or sub-control volume) index on which
+ * the equations are defined
+ */
+ const PrimaryVariables &storageTerm(int i) const
+ { return storageTerm_[i]; }
+
/*!
* \brief Returns the epsilon value which is added and removed
* from the current solution.
@@ -211,8 +281,7 @@ public:
* which the local derivative ought to be calculated.
* \param pvIdx The index of the primary variable which gets varied
*/
- Scalar numericEpsilon(const ElementVariables &elemVars,
- int scvIdx,
+ Scalar numericEpsilon(int scvIdx,
int pvIdx) const
{
// define the base epsilon as the geometric mean of 1 and the
@@ -222,102 +291,71 @@ 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 = elemVars.volVars(scvIdx, /*historyIdx=*/0).primaryVar(pvIdx);
+ Scalar pv = this->curVolVars_[scvIdx].primaryVar(pvIdx);
return baseEps*(std::abs(pv) + 1);
}
- /*!
- * \brief Return reference to the local residual.
- */
- LocalResidual &localResidual()
- { return localResidual_; }
- const LocalResidual &localResidual() const
- { return localResidual_; }
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
/*!
- * \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
+ * \brief Returns a reference to the problem.
*/
- const MatrixBlock &jacobian(int domainScvIdx, int rangeScvIdx) const
- { return jacobian_[domainScvIdx][rangeScvIdx]; }
+ const Problem &problem_() const
+ {
+ Valgrind::CheckDefined(problemPtr_);
+ return *problemPtr_;
+ };
/*!
- * \brief Returns the local Jacobian matrix the storage term of a sub-control volume.
- *
- * \param scvIdx The local index of sub control volume
+ * \brief Returns a reference to the grid view.
*/
- const MatrixBlock &jacobianStorage(int scvIdx) const
- { return jacobianStorage_[scvIdx]; }
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
/*!
- * \brief Returns the local residual of a sub-control volume.
- *
- * \param scvIdx The local index of the sub control volume
+ * \brief Returns a reference to the element.
*/
- const VectorBlock &residual(int scvIdx) const
- { return residual_[scvIdx]; }
+ const Element &elem_() const
+ {
+ Valgrind::CheckDefined(elemPtr_);
+ return *elemPtr_;
+ };
/*!
- * \brief Returns the local storage term of a sub-control volume.
- *
- * \param scvIdx The local index of the sub control volume
+ * \brief Returns a reference to the model.
*/
- 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 *modelPtr_; }
+ { return problem_().model(); };
/*!
- * \brief Returns the numeric difference method which is applied.
+ * \brief Returns a reference to the jacobian assembler.
*/
- static int numericDifferenceMethod_()
- { return GET_PARAM(TypeTag, int, NumericDifferenceMethod); }
+ const JacobianAssembler &jacAsm_() const
+ { return model_().jacobianAssembler(); }
/*!
- * \brief Resize all internal attributes to the size of the
- * element.
+ * \brief Returns a reference to the vertex mapper.
*/
- 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);
- };
+ const VertexMapper &vertexMapper_() const
+ { return problem_().vertexMapper(); };
/*!
- * \brief Reset the all relevant internal attributes to 0
+ * \brief Reset the local jacobian matrix to 0
*/
- void reset_(const ElementVariables &elemVars)
+ void reset_()
{
- 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;
+ 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;
}
}
}
@@ -366,19 +404,21 @@ protected:
* for which the partial derivative ought to be
* calculated
*/
- void evalPartialDerivative_(ElementVariables &elemVars,
+ void evalPartialDerivative_(ElementSolutionVector &dest,
+ PrimaryVariables &destStorage,
int scvIdx,
int pvIdx)
{
- // save all quantities which depend on the specified primary
- // variable at the given sub control volume
- elemVars.saveScvVars(scvIdx);
+ int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
+
+ PrimaryVariables priVars(model_().curSol()[globalIdx]);
+ VolumeVariables origVolVars(curVolVars_[scvIdx]);
- PrimaryVariables priVars(elemVars.volVars(scvIdx, /*historyIdx=*/0).primaryVars());
- Scalar eps = asImp_().numericEpsilon(elemVars, scvIdx, pvIdx);
+ curVolVars_[scvIdx].setEvalPoint(&origVolVars);
+ Scalar eps = asImp_().numericEpsilon(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)
@@ -387,19 +427,32 @@ protected:
delta += eps;
// calculate the residual
- elemVars.updateScvVars(priVars, scvIdx, /*historyIdx=*/0);
- elemVars.updateAllScvfVars();
- localResidual_.eval(derivResidual_, derivStorage_, elemVars);
+ 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];
}
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)
- derivResidual_ = residual_;
- derivStorage_[scvIdx] = residualStorage_[scvIdx];
+ dest = residual_;
+ destStorage = storageTerm_[scvIdx];
}
- if (numericDifferenceMethod_() <= 0) {
+
+ if (numericDifferenceMethod_ <= 0) {
// we are not using forward differences, i.e. we don't
// need to calculate f(x - \epsilon)
@@ -407,35 +460,40 @@ protected:
priVars[pvIdx] -= delta + eps;
delta += eps;
- // 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];
+ // 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];
}
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)
- derivResidual_ -= residual_;
- derivStorage_[scvIdx] -= residualStorage_[scvIdx];
+ dest -= residual_;
+ destStorage -= storageTerm_[scvIdx];
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
- derivResidual_ /= delta;
- derivStorage_[scvIdx] /= delta;
+ dest /= delta;
+ destStorage /= delta;
- // restore the original state of the element's volume
- // variables
- elemVars.restoreScvVars(scvIdx);
+ // restore the original state of the element's volume variables
+ curVolVars_[scvIdx] = origVolVars;
-#ifndef NDEBUG
- for (unsigned i = 0; i < derivResidual_.size(); ++i)
- Valgrind::CheckDefined(derivResidual_[i]);
+#if HAVE_VALGRIND
+ for (unsigned i = 0; i < dest.size(); ++i)
+ Valgrind::CheckDefined(dest[i]);
#endif
}
@@ -444,44 +502,56 @@ protected:
* partial derivatives of all equations in regard to the
* primary variable 'pvIdx' at vertex 'scvIdx' .
*/
- void updateLocalJacobian_(const ElementVariables &elemVars,
- int scvIdx,
- int pvIdx)
+ void updateLocalJacobian_(int scvIdx,
+ int pvIdx,
+ const ElementSolutionVector &deriv,
+ const PrimaryVariables &storageDeriv)
{
// store the derivative of the storage term
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- jacobianStorage_[scvIdx][eqIdx][pvIdx] = derivStorage_[scvIdx][eqIdx];
+ storageJacobian_[scvIdx][eqIdx][pvIdx] = storageDeriv[eqIdx];
}
- int numScv = elemVars.numScv();
- for (int eqScvIdx = 0; eqScvIdx < numScv; eqScvIdx++)
+ for (int i = 0; i < fvElemGeom_.numVertices; i++)
{
+ if (jacAsm_().vertexColor(elem_(), i) == Green) {
+ // Green vertices are not to be changed!
+ continue;
+ }
+
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // 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]);
+ // 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]);
}
}
}
+ const Element *elemPtr_;
+ FVElementGeometry fvElemGeom_;
+
+ ElementBoundaryTypes bcTypes_;
+
+ // The problem we would like to solve
Problem *problemPtr_;
- Model *modelPtr_;
- ElementVariables *internalElemVars_;
+ // secondary variables at the previous and at the current time
+ // levels
+ ElementVolumeVariables prevVolVars_;
+ ElementVolumeVariables curVolVars_;
- LocalBlockMatrix jacobian_;
- LocalStorageMatrix jacobianStorage_;
+ LocalResidual localResidual_;
- LocalBlockVector residual_;
- LocalBlockVector residualStorage_;
+ LocalBlockMatrix A_;
+ std::vector<MatrixBlock> storageJacobian_;
- LocalBlockVector derivResidual_;
- LocalBlockVector derivStorage_;
+ ElementSolutionVector residual_;
+ ElementSolutionVector storageTerm_;
- LocalResidual localResidual_;
+ int numericDifferenceMethod_;
};
}
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index b3176627521054f4dd00e5e3a041a1efc27dda12..d22a826c8fa3dc6143c212761b96aac1e6918a3c 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/bvector.hh>
+#include <dune/istl/matrix.hh>
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
@@ -47,33 +47,48 @@ 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 { dim = GridView::dimension };
+
+ 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 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(ElementVariables)) ElementVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ 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(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
- enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
+ 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;
- 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()
@@ -82,217 +97,312 @@ 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
- * conservation equations from zero and store the results
- * internally.
+ * equations from zero.
*
- * The results can be requested afterwards using the residual()
- * and storageTerm() methods.
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
*/
- void eval(const ElementVariables &elemVars)
+ void eval(const Element &element)
{
- int numScv = elemVars.numScv();
- internalResidual_.resize(numScv);
- internalStorageTerm_.resize(numScv);
- eval(internalResidual_, internalStorageTerm_, 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);
+ }
/*!
- * \brief Return the result of the eval() call using internal
- * storage.
+ * \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
*/
- const LocalBlockVector &residual() const
- { return internalResidual_; }
- const VectorBlock &residual(int scvIdx) const
- { return internalResidual_[scvIdx]; }
+ 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);
+ }
/*!
- * \brief Return the storage term calculated using the last call
- * to eval() using internal storage.
+ * \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
*/
- const LocalBlockVector &storageTerm() const
- { return internalStorageTerm_; }
- const VectorBlock &storageTerm(int scvIdx) const
- { return internalStorageTerm_[scvIdx]; }
-
+ 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_();
+ }
+
/*!
* \brief Compute the local residual, i.e. the deviation of the
- * conservation equations from zero.
+ * equations from zero.
*
- * \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
+ * \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
*/
- void eval(LocalBlockVector &residual,
- LocalBlockVector &storageTerm,
- const ElementVariables &elemVars) const
+ void eval(const Element &element,
+ const FVElementGeometry &fvGeom,
+ const ElementVolumeVariables &prevVolVars,
+ const ElementVolumeVariables &curVolVars,
+ const ElementBoundaryTypes &bcTypes)
{
- residual = 0.0;
- storageTerm = 0.0;
+ //Valgrind::CheckDefined(fvGeom);
+ Valgrind::CheckDefined(prevVolVars);
+ Valgrind::CheckDefined(curVolVars);
- // evaluate the flux terms
- asImp_().evalFluxes(residual, elemVars);
+#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_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
- Valgrind::CheckDefined(residual[i]);
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
- // evaluate the storage and the source terms
- asImp_().evalVolumeTerms_(residual, storageTerm, elemVars);
+ asImp_().evalVolumeTerms_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
- Valgrind::CheckDefined(residual[i]);
-#endif // !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++) {
+ Valgrind::CheckDefined(residual_[i]);
+ }
+#endif // HAVE_VALGRIND
// evaluate the boundary conditions
- asImp_().evalBoundary_(residual, storageTerm, elemVars);
+ asImp_().evalBoundary_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < elemVars.fvElemGeom().numVertices; i++)
- Valgrind::CheckDefined(residual[i]);
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
}
/*!
- * \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.
+ * \brief Returns the local residual for a given sub-control
+ * volume of the element.
*/
- void evalStorage(const ElementVariables &elemVars,
- int historyIdx)
- {
- int numScv = elemVars.numScv();
- internalStorageTerm_.resize(numScv);
- evalStorage(internalStorageTerm_,
- elemVars,
- historyIdx);
- }
+ const ElementSolutionVector &residual() const
+ { return residual_; }
/*!
- * \brief Calculate the amount of all conservation quantities
- * stored in all element's sub-control volumes for a given
- * history index.
+ * \brief Returns the local residual for a given sub-control
+ * volume of the element.
*
- * This is used to figure out how much of each conservation
- * quantity is inside the element.
+ * \param scvIdx The local index of the sub-control volume
+ * (i.e. the element's local vertex index)
*/
- 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();
- }
- }
+ const PrimaryVariables &residual(int scvIdx) const
+ { return residual_[scvIdx]; }
/*!
- * \brief Add the flux term to a local residual.
+ * \brief Returns the storage term for all sub-control volumes of the
+ * element.
*/
- void evalFluxes(LocalBlockVector &residual,
- const ElementVariables &elemVars) const
- {
- 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 ElementSolutionVector &storageTerm() const
+ { return storageTerm_; }
- Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, /*context=*/elemVars, scvfIdx);
- flux *= elemVars.fluxVars(scvfIdx).extrusionFactor();
- Valgrind::CheckDefined(flux);
+protected:
+ Implementation &asImp_()
+ {
+ assert(static_cast<Implementation*>(this) != 0);
+ return *static_cast<Implementation*>(this);
+ }
- // 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;
- }
+ const Implementation &asImp_() const
+ {
+ assert(static_cast<const Implementation*>(this) != 0);
+ return *static_cast<const Implementation*>(this);
}
-protected:
/*!
- * \brief Evaluate the boundary conditions of an element.
+ * \brief Evaluate the boundary conditions
+ * of the current element.
*/
- void evalBoundary_(LocalBlockVector &residual,
- LocalBlockVector &storageTerm,
- const ElementVariables &elemVars) const
+ void evalBoundary_()
{
- if (!elemVars.onBoundary()) {
- return;
- }
-
- if (elemVars.hasNeumann())
- asImp_().evalNeumann_(residual, elemVars);
-
- if (elemVars.hasDirichlet())
- asImp_().evalDirichlet_(residual, storageTerm, elemVars);
+ 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_();
}
/*!
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
*/
- void evalDirichlet_(LocalBlockVector &residual,
- LocalBlockVector &storageTerm,
- const ElementVariables &elemVars) const
+ void evalDirichlet_()
{
PrimaryVariables tmp(0);
- for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
- const BoundaryTypes &bcTypes = elemVars.boundaryTypes(scvIdx);
- if (!bcTypes.hasDirichlet())
+ for (int i = 0; i < fvElemGeom_().numVertices; ++i) {
+ const BoundaryTypes &bcTypes = bcTypes_(i);
+ if (! bcTypes.hasDirichlet())
continue;
-
- // ask the problem for the dirichlet values
- asImp_().computeDirichlet_(tmp, elemVars, scvIdx);
- const PrimaryVariables &priVars =
- elemVars.volVars(scvIdx).primaryVars();
+ // ask the problem for the dirichlet values
+ const VertexPointer vPtr = elem_().template subEntity<dim>(i);
+ asImp_().problem_().dirichlet(tmp, *vPtr);
// 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(priVars);
+ Valgrind::CheckDefined(curPrimaryVar_(i, pvIdx));
Valgrind::CheckDefined(tmp[pvIdx]);
- residual[scvIdx][eqIdx] = priVars[pvIdx] - tmp[pvIdx];
- storageTerm[scvIdx][eqIdx] = 0.0;
+ this->residual_[i][eqIdx] =
+ curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
+
+ this->storageTerm_[i][eqIdx] = 0.0;
};
};
}
@@ -301,16 +411,13 @@ protected:
* \brief Add all Neumann boundary conditions to the local
* residual.
*/
- void evalNeumann_(LocalBlockVector &residual,
- const ElementVariables &elemVars) const
+ void evalNeumann_()
{
- const Element &elem = elemVars.element();
- Dune::GeometryType geoType = elem.geometry().type();
+ Dune::GeometryType geoType = elem_().geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
- const GridView &gridView = elemVars.gridView();
- IntersectionIterator isIt = gridView.ibegin(elem);
- const IntersectionIterator &endIt = gridView.iend(elem);
+ IntersectionIterator isIt = gridView_().ibegin(elem_());
+ const IntersectionIterator &endIt = gridView_().iend(elem_());
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -325,59 +432,112 @@ protected:
faceVertIdx < numFaceVerts;
++faceVertIdx)
{
- int scvIdx = refElem.subEntity(/*entityIdx=*/faceIdx,
- /*entityCodim=*/1,
- /*subEntityIdx=*/faceVertIdx,
- /*subEntityCodim=*/dim);
-
+ int elemVertIdx = refElem.subEntity(faceIdx,
+ 1,
+ faceVertIdx,
+ dim);
+
int boundaryFaceIdx =
- elemVars.fvElemGeom().boundaryFaceIndex(faceIdx, faceVertIdx);
+ fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
- evalNeumannSegment_(residual,
- elemVars,
- *isIt,
- scvIdx,
+ evalNeumannSegment_(isIt,
+ elemVertIdx,
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_(LocalBlockVector &residual,
- const ElementVariables &elemVars,
- const Intersection &is,
+ void evalNeumannSegment_(const IntersectionIterator &isIt,
int scvIdx,
- int boundaryFaceIdx) const
+ int boundaryFaceIdx)
{
// temporary vector to store the neumann boundary fluxes
- const BoundaryTypes &bcTypes = elemVars.boundaryTypes(scvIdx);
- PrimaryVariables values;
-
+ PrimaryVariables values(0.0);
+ const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
- elemVars.problem().neumann(values,
- elemVars,
- is,
- scvIdx,
- boundaryFaceIdx);
- Valgrind::CheckDefined(values);
-
+ problem_().boxSDNeumann(values,
+ elem_(),
+ fvElemGeom_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx,
+ curVolVars_());
+ if (!Valgrind::CheckDefined(values))
+ std::cerr << "undefined: " << values << "\n";
values *=
- elemVars.fvElemGeom().boundaryFace[boundaryFaceIdx].area
- * elemVars.volVars(scvIdx, /*historyIdx*/0).extrusionFactor();
+ fvElemGeom_().boundaryFace[boundaryFaceIdx].area
+ * curVolVars_(scvIdx).extrusionFactor();
Valgrind::CheckDefined(values);
- residual[scvIdx] += 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();
}
}
@@ -386,17 +546,15 @@ protected:
* to the local residual of all sub-control volumes of the
* current element.
*/
- void evalVolumeTerms_(LocalBlockVector &residual,
- LocalBlockVector &storageTerm,
- const ElementVariables &elemVars) const
+ void evalVolumeTerms_()
{
- PrimaryVariables tmp, tmp2;
-
// evaluate the volume terms (storage + source terms)
- for (int scvIdx=0; scvIdx < elemVars.numScv(); scvIdx++)
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
{
Scalar extrusionFactor =
- elemVars.volVars(scvIdx, /*historyIdx=*/0).extrusionFactor();
+ curVolVars_(i).extrusionFactor();
+
+ PrimaryVariables tmp(0);
// mass balance within the element. this is the
// $\frac{m}{\partial t}$ term if using implicit
@@ -404,51 +562,162 @@ protected:
//
// TODO (?): we might need a more explicit way for
// doing the time discretization...
- asImp_().computeStorage(tmp2,
- elemVars,
- scvIdx,
- /*historyIdx=*/1);
- asImp_().computeStorage(tmp,
- elemVars,
- scvIdx,
- /*historyIdx=*/0);
+ this->asImp_().computeStorage(storageTerm_[i], i, false);
+ this->asImp_().computeStorage(tmp, i, true);
- tmp -= tmp2;
- tmp *=
- elemVars.fvElemGeom().subContVol[scvIdx].volume
- / elemVars.problem().timeManager().timeStepSize()
+ storageTerm_[i] -= tmp;
+ storageTerm_[i] *=
+ fvElemGeom_().subContVol[i].volume
+ / problem_().timeManager().timeStepSize()
* extrusionFactor;
+ residual_[i] += storageTerm_[i];
- 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;
+ // subtract the source term from the local rate
+ this->asImp_().computeSource(tmp, i);
+ tmp *= fvElemGeom_().subContVol[i].volume * extrusionFactor;
+ residual_[i] -= tmp;
// make sure that only defined quantities were used
// to calculate the residual.
- Valgrind::CheckDefined(storageTerm[scvIdx]);
- Valgrind::CheckDefined(residual[scvIdx]);
+ Valgrind::CheckDefined(residual_[i]);
}
}
-private:
- Implementation &asImp_()
+ /*!
+ * \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
{
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
+ Valgrind::CheckDefined(elemPtr_);
+ return *elemPtr_;
}
- const Implementation &asImp_() const
+ /*!
+ * \brief Returns a reference to the current element's finite
+ * volume geometry.
+ */
+ const FVElementGeometry &fvElemGeom_() const
{
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
+ 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
+ {
+ Valgrind::CheckDefined(bcTypesPtr_);
+ return *bcTypesPtr_;
}
- LocalBlockVector internalResidual_;
- LocalBlockVector internalStorageTerm_;
+ /*!
+ * \brief Returns a reference to the boundary types of the i-th
+ * sub-control volume of the current element.
+ */
+ const BoundaryTypes &bcTypes_(int i) const
+ {
+ return bcTypes_()[i];
+ }
+
+protected:
+ ElementSolutionVector storageTerm_;
+ ElementSolutionVector residual_;
+
+ // The problem we would like to solve
+ Problem *problemPtr_;
+
+ const Element *elemPtr_;
+ const FVElementGeometry *fvElemGeomPtr_;
+
+ // current and previous secondary variables for the element
+ const ElementVolumeVariables *prevVolVarsPtr_;
+ const ElementVolumeVariables *curVolVarsPtr_;
+
+ const ElementBoundaryTypes *bcTypesPtr_;
};
}
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index bfc52156f2a649e628d9a8a89b6d155a824bfb66..63625996a3a98be237fab81f12b0bb6749bf8f40 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -72,18 +72,16 @@ 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
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;
@@ -95,7 +93,6 @@ 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;
@@ -109,7 +106,9 @@ public:
* \brief The constructor.
*/
BoxModel()
- { }
+ {
+ enableHints_ = GET_PARAM(TypeTag, bool, EnableHints);
+ }
~BoxModel()
{ delete jacAsm_; }
@@ -124,11 +123,9 @@ public:
{
problemPtr_ = &prob;
- updateBoundaryTypes_();
-
int nDofs = asImp_().numDofs();
- for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
- solution_[historyIdx].resize(nDofs);
+ uCur_.resize(nDofs);
+ uPrev_.resize(nDofs);
boxVolume_.resize(nDofs);
localJacobian_.init(problem_());
@@ -138,10 +135,10 @@ public:
asImp_().applyInitialSolution_();
// resize the hint vectors
- if (enableHints_()) {
+ if (enableHints_) {
int nVerts = gridView_().size(dim);
- for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
- hints_[historyIdx].resize(nVerts);
+ curHints_.resize(nVerts);
+ prevHints_.resize(nVerts);
hintsUsable_.resize(nVerts);
std::fill(hintsUsable_.begin(),
hintsUsable_.end(),
@@ -150,67 +147,74 @@ public:
// also set the solution of the "previous" time step to the
// initial solution.
- solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
+ uPrev_ = uCur_;
}
- void setAllHints(ElementVariables &elemVars) const
+ void setHints(const Element &elem,
+ ElementVolumeVariables &prevVolVars,
+ ElementVolumeVariables &curVolVars) const
{
- if (!enableHints_())
+ if (!enableHints_)
return;
- for (int historyIdx = 0; historyIdx < historySize; ++historyIdx)
- setHints(elemVars, historyIdx);
- }
+ 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]);
+ }
+ }
+ };
- void setHints(ElementVariables &elemVars, int historyIdx) const
+ void setHints(const Element &elem,
+ ElementVolumeVariables &curVolVars) const
{
- if (!enableHints_())
+ if (!enableHints_)
return;
- int numScv = elemVars.numScv();
- for (int scvIdx = 0; scvIdx < numScv; ++scvIdx) {
- setHintsScv(elemVars, historyIdx, scvIdx);
- }
- }
+ int n = elem.template count<dim>();
+ curVolVars.resize(n);
+ for (int i = 0; i < n; ++i) {
+ int globalIdx = vertexMapper().map(elem, i, dim);
- 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]);
+ if (!hintsUsable_[globalIdx])
+ curVolVars[i].setHint(NULL);
+ else
+ curVolVars[i].setHint(&curHints_[globalIdx]);
+ }
};
- void updateAllHints(const ElementVariables &elemVars) const
+ void updatePrevHints()
{
- if (!enableHints_())
+ if (!enableHints_)
return;
- for (int historyIdx = 0; historyIdx < historySize; ++historyIdx) {
- updateHints(elemVars, historyIdx);
- }
+ prevHints_ = curHints_;
};
- void updateHints(const ElementVariables &elemVars, int historyIdx) const
+ void updateCurHints(const Element &elem,
+ const ElementVolumeVariables &ev) const
{
- if (!enableHints_())
+ if (!enableHints_)
return;
-
- for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
- int globalIdx = vertexMapper().map(elemVars.element(), scvIdx, dim);
- hints_[historyIdx][globalIdx] = elemVars.volVars(scvIdx, historyIdx);
+
+ 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];
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
@@ -222,10 +226,10 @@ public:
Scalar globalResidual(SolutionVector &dest,
const SolutionVector &u)
{
- SolutionVector tmp(solution(/*historyIdx=*/0));
- solution(/*historyIdx=*/0) = u;
+ SolutionVector tmp(curSol());
+ curSol() = u;
Scalar res = globalResidual(dest);
- solution(/*historyIdx=*/0) = tmp;
+ curSol() = tmp;
return res;
}
@@ -239,16 +243,14 @@ public:
{
dest = 0;
- ElementVariables elemVars(this->problem_());
ElementIterator elemIt = gridView_().template begin<0>();
const ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
- elemVars.updateAll(*elemIt);
- localResidual().eval(elemVars);
+ localResidual().eval(*elemIt);
- for (int scvIdx = 0; scvIdx < elemVars.numScv(); ++scvIdx) {
- int globalI = vertexMapper().map(*elemIt, scvIdx, dim);
- dest[globalI] += localResidual().residual(scvIdx);
+ for (int i = 0; i < elemIt->template count<dim>(); ++i) {
+ int globalI = vertexMapper().map(*elemIt, i, dim);
+ dest[globalI] += localResidual().residual(i);
}
};
@@ -277,14 +279,11 @@ 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) {
- elemVars.updateFVElemGeom(*elemIt);
- elemVars.updateScvVars(/*historyIdx=*/0);
- localResidual().evalStorage(elemVars, /*historyIdx=*/0);
+ localResidual().evalStorage(*elemIt);
for (int i = 0; i < elemIt->template count<dim>(); ++i)
dest += localResidual().storageTerm()[i];
@@ -303,16 +302,28 @@ public:
{ return boxVolume_[globalIdx][0]; }
/*!
- * \brief Reference to the solution at a given history index as a block vector.
+ * \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.
*/
- const SolutionVector &solution(int historyIdx) const
- { return solution_[historyIdx]; }
+ const SolutionVector &prevSol() const
+ { return uPrev_; }
/*!
- * \brief Reference to the solution at a given history index as a block vector.
+ * \brief Reference to the previous solution as a block vector.
*/
- SolutionVector &solution(int historyIdx)
- { return solution_[historyIdx]; }
+ SolutionVector &prevSol()
+ { return uPrev_; }
/*!
* \brief Returns the operator assembler for the global jacobian of
@@ -363,8 +374,7 @@ public:
*/
Scalar primaryVarWeight(int vertIdx, int pvIdx) const
{
- Scalar absPv = std::abs(this->solution(/*historyIdx=*/1)[vertIdx][pvIdx]);
- return 1.0/std::max(absPv, 1.0);
+ return 1.0/std::max(std::abs(this->prevSol()[vertIdx][pvIdx]), 1.0);
}
/*!
@@ -409,8 +419,8 @@ public:
NewtonController &controller)
{
#if HAVE_VALGRIND
- for (size_t i = 0; i < solution(/*historyIdx=*/0).size(); ++i)
- Valgrind::CheckDefined(solution(/*historyIdx=*/0)[i]);
+ for (size_t i = 0; i < curSol().size(); ++i)
+ Valgrind::CheckDefined(curSol()[i]);
#endif // HAVE_VALGRIND
asImp_().updateBegin();
@@ -423,8 +433,8 @@ public:
asImp_().updateFailed();
#if HAVE_VALGRIND
- for (size_t i = 0; i < solution(/*historyIdx=*/0).size(); ++i) {
- Valgrind::CheckDefined(solution(/*historyIdx=*/0)[i]);
+ for (size_t i = 0; i < curSol().size(); ++i) {
+ Valgrind::CheckDefined(curSol()[i]);
}
#endif // HAVE_VALGRIND
@@ -438,9 +448,7 @@ public:
* which the actual model can overload.
*/
void updateBegin()
- {
- updateBoundaryTypes_();
- }
+ { }
/*!
@@ -461,7 +469,8 @@ 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.
- hints_[/*historyIdx=*/0] = hints_[/*historyIdx=*/1];
+ uCur_ = uPrev_;
+ curHints_ = prevHints_;
jacAsm_->reassembleAll();
};
@@ -476,10 +485,10 @@ public:
void advanceTimeLevel()
{
// make the current solution the previous one.
- solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
+ uPrev_ = uCur_;
+ prevHints_ = curHints_;
- // shift the hints by one position in the history
- shiftHints();
+ updatePrevHints();
}
/*!
@@ -504,7 +513,7 @@ public:
void deserialize(Restarter &res)
{
res.template deserializeEntities<dim>(asImp_(), this->gridView_());
- solution_[/*historyIdx=*/1] = solution_[/*historyIdx=*/0];
+ prevSol() = curSol();
}
/*!
@@ -529,7 +538,7 @@ public:
}
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << solution_[/*historyIdx=*/0][vertIdx][eqIdx] << " ";
+ outstream << curSol()[vertIdx][eqIdx] << " ";
}
};
@@ -551,7 +560,7 @@ public:
DUNE_THROW(Dune::IOError,
"Could not deserialize vertex "
<< vertIdx);
- instream >> solution_[/*historyIdx=*/0][vertIdx][eqIdx];
+ instream >> curSol()[vertIdx][eqIdx];
}
};
@@ -593,19 +602,6 @@ 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
@@ -613,7 +609,8 @@ public:
* \param element The DUNE codim 0 entity
* \param volVars All volume variables for the element
*/
- void updatePVWeights(const ElementVariables &elemVars) const
+ void updatePVWeights(const Element &element,
+ const ElementVolumeVariables &volVars) const
{ };
/*!
@@ -644,11 +641,10 @@ public:
ScalarField* def[numEq];
ScalarField* delta[numEq];
ScalarField* x[numEq];
- 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);
+ for (int i = 0; i < numEq; ++i) {
+ x[i] = writer.allocateManagedBuffer(numVertices);
+ delta[i] = writer.allocateManagedBuffer(numVertices);
+ def[i] = writer.allocateManagedBuffer(numVertices);
}
VertexIterator vIt = this->gridView_().template begin<dim>();
@@ -656,25 +652,17 @@ public:
for (; vIt != vEndIt; ++ vIt)
{
int globalIdx = vertexMapper().map(*vIt);
- 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];
+ 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];
}
-
- PrimaryVariables uOld(u[globalIdx]);
- PrimaryVariables uNew(uOld - deltaU[globalIdx]);
- (*relError)[globalIdx] = asImp_().relativeErrorVertex(globalIdx,
- uOld,
- uNew);
}
- 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());
+ 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());
}
asImp_().addOutputVtkFields(u, writer);
@@ -696,7 +684,7 @@ public:
void addOutputVtkFields(const SolutionVector &sol,
MultiWriter &writer)
{
- typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
// create the required scalar fields
unsigned numVertices = this->gridView_().size(dim);
@@ -728,9 +716,6 @@ 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.
*/
@@ -754,83 +739,13 @@ 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
- SolutionVector &uCur = solution(/*historyIdx=*/0);
- uCur = Scalar(0.0);
+ uCur_ = Scalar(0.0);
boxVolume_ = Scalar(0.0);
FVElementGeometry fvElemGeom;
@@ -870,8 +785,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]);
}
}
@@ -886,7 +801,7 @@ protected:
Dune::ForwardCommunication);
VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
- sumPVHandle(uCur, vertexMapper());
+ sumPVHandle(uCur_, vertexMapper());
gridView().communicate(sumPVHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
@@ -895,14 +810,15 @@ 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> hints_[historySize];
+ mutable std::vector<VolumeVariables> curHints_;
+ mutable std::vector<VolumeVariables> prevHints_;
/*!
* \brief Returns whether messages should be printed
@@ -927,13 +843,13 @@ protected:
// cur is the current iterative solution, prev the converged
// solution of the previous time step
- SolutionVector solution_[historySize];
-
- // all the index of the BoundaryTypes object for a vertex
- std::vector<int> boundaryVertexIndex_;
- std::vector<BoundaryTypes> boundaryTypes_;
+ SolutionVector uCur_;
+ SolutionVector uPrev_;
Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
+
+private:
+ bool enableHints_;
};
}
diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index 9a5fbaab9b1e5e503fcb0e0236ffcad75652af4c..b8284c38f341b1f8a8fc5d49a4c94ea855524ceb 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(ElementVariables)) ElementVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
@@ -139,28 +139,6 @@ 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.
@@ -192,27 +170,6 @@ 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
@@ -223,7 +180,6 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
- DUNE_DEPRECATED
void dirichlet(PrimaryVariables &values,
const Vertex &vertex) const
{
@@ -242,7 +198,6 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
- DUNE_DEPRECATED
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
{
@@ -254,25 +209,6 @@ 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.
@@ -298,8 +234,7 @@ public:
const Intersection &is,
int scvIdx,
int boundaryFaceIdx,
- const ElementVariables &elemVars) const
- DUNE_DEPRECATED
+ const ElementVolumeVariables &elemVolVars) const
{
// forward it to the interface without the volume variables
asImp_().neumann(values,
@@ -330,7 +265,6 @@ 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);
@@ -348,7 +282,6 @@ public:
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
- DUNE_DEPRECATED
{
// Throw an exception (there is no reasonable default value
// for Neumann conditions)
@@ -358,21 +291,6 @@ 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.
@@ -395,11 +313,10 @@ public:
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx,
- const ElementVariables &elemVars) const
- DUNE_DEPRECATED
+ const ElementVolumeVariables &elemVolVars) const
{
// forward to solution independent, box specific interface
- asImp_().source(values, elemVars, scvIdx);
+ asImp_().source(values, element, fvElemGeom, scvIdx);
}
/*!
@@ -419,7 +336,6 @@ public:
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
- DUNE_DEPRECATED
{
// forward to generic interface
asImp_().sourceAtPos(values, fvElemGeom.subContVol[scvIdx].global);
@@ -440,7 +356,6 @@ public:
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &pos) const
- DUNE_DEPRECATED
{
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
@@ -496,28 +411,6 @@ 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
@@ -875,7 +768,7 @@ public:
Scalar t = timeManager().time() + timeManager().timeStepSize();
createResultWriter_();
resultWriter_->beginWrite(t);
- model().addOutputVtkFields(model().solution(/*historyIdx=*/0), *resultWriter_);
+ model().addOutputVtkFields(model().curSol(), *resultWriter_);
asImp_().addOutputVtkFields();
resultWriter_->endWrite();
}
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index ae75bb159e69f9254a2047f5e42be85f6d10ea81..466a281393a39b24f57f3278fa3ab2e2a0da051d 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(ElementVariables); //!< The secondary variables of all sub-control volumes in an element
+NEW_PROP_TAG(ElementVolumeVariables); //!< 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,15 +121,12 @@ NEW_PROP_TAG(EnableHints);
// mappers from local to global indices
-//! mapper for vertices
+//! maper for vertices
NEW_PROP_TAG(VertexMapper);
-//! mapper for elements
+//! maper for elements
NEW_PROP_TAG(ElementMapper);
-//! mapper for degrees of freedom
+//! maper 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 c5253660718029c693d3857d5357434224733098..f1517d599b050e9aa01eb02daa6e6af3fc25fcbd 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, ElementVariables, Dumux::BoxElementVariables<TypeTag>);
+SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
/*!
* \brief Boundary types at a single degree of freedom.
@@ -185,8 +185,6 @@ 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 cf85043fbc42f154f0837b096f54545d837b37b4..626d7a9801129d1d7e1481966f329650742b52fe 100644
--- a/dumux/boxmodels/common/boxvolumevariables.hh
+++ b/dumux/boxmodels/common/boxvolumevariables.hh
@@ -49,7 +49,6 @@ 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
@@ -119,13 +118,15 @@ public:
* phase state in the 2p2c model)
*/
void update(const PrimaryVariables &priVars,
- const ElementVariables &elemVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
int scvIdx,
- int historyIdx)
+ bool isOldSol)
{
Valgrind::CheckDefined(priVars);
primaryVars_ = priVars;
- extrusionFactor_ = elemVars.problem().extrusionFactor(elemVars, scvIdx);
+ extrusionFactor_ = problem.boxExtrusionFactor(element, elemGeom, scvIdx);
}
/*!
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index a7d244a324b8585a0f3e273165420724cd8fec35..1bbf74c58004eb4382ac52ccd333ae44eac5bb05 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()%double(timer.elapsed())%double(time())%double(dt);
+ %timeStepIndex()%timer.elapsed()%time()%dt;
}
// write restart file if mandated by the problem
diff --git a/dumux/linear/overlappingscalarproduct.hh b/dumux/linear/overlappingscalarproduct.hh
index 71bfaff5048ce2fb088c2219016e48d5802190fc..2f0360e938678d9022f6acc49d9333b37c7880fc 100644
--- a/dumux/linear/overlappingscalarproduct.hh
+++ b/dumux/linear/overlappingscalarproduct.hh
@@ -67,8 +67,8 @@ public:
double norm(const OverlappingBlockVector &x)
{
- field_type tmp = std::sqrt(dot(x, x));
- return tmp;
+ double tmp = dot(x, x);
+ return std::sqrt(tmp);
};
private:
diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh
index 79e4cd8caf27a84114a62bbd492595c6e5969bd1..28cf38036859772409466dfb3bda6bdfefc5868b 100644
--- a/dumux/material/spatialparameters/boxspatialparameters.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters.hh
@@ -72,30 +72,20 @@ 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 localIdx) const
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
- return asImp_().materialLawParams(fvElemGeom.subContVol[localIdx].global);
+ return asImp_().materialLawParamsAtPos(element.geometry().center());
}
-
/*!
* \brief Function for defining the parameters needed by constitutive relationships (kr-Sw, pc-Sw, etc.).
*
- * \param globalPos The position of the center of the element
* \return the material parameters object
+ * \param globalPos The position of the center of the element
*/
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
diff --git a/dumux/material/spatialparameters/boxspatialparameters1p.hh b/dumux/material/spatialparameters/boxspatialparameters1p.hh
index fe7232fc781a5b104778f6707727c3c33ff0f45f..a70a3a4ae56406a2fc0d5b1e0e99019af536efb1 100644
--- a/dumux/material/spatialparameters/boxspatialparameters1p.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters1p.hh
@@ -66,7 +66,6 @@ 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)
@@ -75,6 +74,53 @@ 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
@@ -110,26 +156,15 @@ 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
*/
- DUNE_DEPRECATED
- const Tensor &intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ const Tensor& intrinsicPermeability (const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
return asImp_().intrinsicPermeabilityAtPos(element.geometry().center());
}
@@ -140,7 +175,6 @@ 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,
@@ -148,26 +182,15 @@ 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());
}
@@ -178,7 +201,6 @@ 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 2165c39f8b8f10ad9032585101681148b13e230e..a083ef605c76f3e23f4067e2378b859b1380a60d 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -116,7 +116,7 @@ public:
protected:
bool execute_(NewtonController &ctl)
{
- SolutionVector &uCurrentIter = model().solution(/*historyIdx=*/0);
+ SolutionVector &uCurrentIter = model().curSol();
SolutionVector uLastIter(uCurrentIter);
SolutionVector deltaU(uCurrentIter);
diff --git a/test/boxmodels/2p/lensproblem.hh b/test/boxmodels/2p/lensproblem.hh
index d8bc6f25412d514a453bb6b655d306444fda66ec..dab4dc78d7b4ae513a3526ccbcceb3453ff43b77 100644
--- a/test/boxmodels/2p/lensproblem.hh
+++ b/test/boxmodels/2p/lensproblem.hh
@@ -25,11 +25,10 @@
* \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
@@ -58,8 +57,6 @@ 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>);
@@ -246,7 +243,6 @@ 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 ad82f044eac3d4e8c74eb3815bcbba5b81b07787..d0a34f0f153adba60b1a2aa9e2c56d045399e6de 100644
--- a/test/boxmodels/2p/lensspatialparameters.hh
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -135,24 +135,27 @@ public:
* \param scvIdx The index sub-control volume face where the
* intrinsic velocity ought to be calculated.
*/
- template <class Context>
- Scalar intrinsicPermeability(const Context &context, int localIdx) const
+ Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
- const GlobalPosition &globalPos = context.pos(localIdx);
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
if (isInLens_(globalPos))
return lensK_;
return outerK_;
}
- template <class Context>
- Scalar porosity(const Context &context, int localIdx) const
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{ return 0.4; }
// return the parameter object for the Brooks-Corey material law which depends on the position
- template <class Context>
- const MaterialLawParams& materialLawParams(const Context &context, int localIdx) const
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
- const GlobalPosition &globalPos = context.pos(localIdx);
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
if (isInLens_(globalPos))
return lensMaterialParams_;