diff --git a/dumux/implicit/cellcentered/elementvolumevariables.hh b/dumux/implicit/cellcentered/elementvolumevariables.hh
index dfa3ae46af0d1c5a494e19ad34832b689d3a13e2..c09efacc7cfa010884b16827877f1629683f2bdc 100644
--- a/dumux/implicit/cellcentered/elementvolumevariables.hh
+++ b/dumux/implicit/cellcentered/elementvolumevariables.hh
@@ -34,12 +34,13 @@ namespace Dumux
* volume variables object for each of the element's vertices
*/
template<class TypeTag>
-class CCElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) >
+class CCElementVolumeVariables
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -55,35 +56,21 @@ public:
* \param fvGeometry The finite volume geometry of the element
* \param oldSol Tells whether the model's previous or current solution should be used.
*/
- void update(const Problem &problem,
+ void update(Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
bool oldSol)
{
- const SolutionVector &globalSol =
- oldSol?
- problem.model().prevSol():
- problem.model().curSol();
+ oldSol_ = oldSol;
+ problem_ = &problem;
int numNeighbors = fvGeometry.numNeighbors;
- this->resize(numNeighbors);
+ stencil_.resize(numNeighbors);
for (int i = 0; i < numNeighbors; i++)
{
const Element& neighbor = fvGeometry.neighbors[i];
-
- const PrimaryVariables &solI
- = globalSol[problem.elementMapper().index(neighbor)];
-
- FVElementGeometry neighborFVGeom;
- neighborFVGeom.updateInner(neighbor);
-
- (*this)[i].update(solI,
- problem,
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- oldSol);
+ stencil_[i] = problem.elementMapper().index(neighbor);
}
// only treat boundary if current solution is evaluated
@@ -96,40 +83,69 @@ public:
|| elemBCTypes.hasNeumann()
|| elemBCTypes.hasOutflow())
{
- this->resize(numNeighbors + element.subEntities(1));
+ ghostVolVars_ = std::vector<VolumeVariables>(element.subEntities(1));
+ isDirichlet_ = std::vector<bool>(element.subEntities(1), false);
// add volume variables for the boundary faces
for (const auto& intersection : intersections(problem.gridView(), element))
{
- if (intersection.boundary())
+ if (!intersection.boundary())
+ continue;
+
+ BoundaryTypes bcTypes;
+ problem.boundaryTypes(bcTypes, intersection);
+
+ int fIdx = intersection.indexInInside();
+ if (bcTypes.hasDirichlet())
{
- BoundaryTypes bcTypes;
- problem.boundaryTypes(bcTypes, intersection);
-
- int fIdx = intersection.indexInInside();
- int indexInVariables = numNeighbors + fIdx;
-
- if (bcTypes.hasDirichlet())
- {
- PrimaryVariables dirichletValues;
- problem.dirichlet(dirichletValues, intersection);
-
- (*this)[indexInVariables].update(dirichletValues,
- problem,
- element,
- fvGeometry,
- /*scvIdx=*/0,
- oldSol);
- }
- else
- {
- (*this)[indexInVariables] = (*this)[0];
- }
+ PrimaryVariables dirichletValues;
+ problem.dirichlet(dirichletValues, intersection);
+
+ ghostVolVars_[fIdx].update(dirichletValues,
+ problem,
+ element,
+ fvGeometry,
+ /*scvIdx=*/0,
+ oldSol);
+ isDirichlet_[fIdx] = true;
}
}
}
}
}
+
+ const VolumeVariables& operator [](const unsigned int i) const
+ {
+ if (i >= stencil_.size())
+ {
+ assert(!oldSol_); // old boundary volVars don't exist
+ if (isDirichlet_[i - stencil_.size()]) // ghost volVars only exist for Dirichlet boundaries
+ return ghostVolVars_[i - stencil_.size()];
+ else
+ return problem_->model().volVars(stencil_[0]);
+ }
+ else
+ {
+ if(!oldSol_)
+ return problem_->model().volVars(stencil_[i]);
+ else
+ return problem_->model().prevVolVars(stencil_[i]);
+ }
+ }
+
+ VolumeVariables& operator [](const unsigned int i)
+ {
+ assert(!oldSol_); // old volVars should never be modified
+ assert(i < stencil_.size()); // boundary volVars should not be modified
+ return problem_->model().volVars(stencil_[i]);
+ }
+
+private:
+ bool oldSol_ = false;
+ std::vector<unsigned int> stencil_;
+ std::vector<VolumeVariables> ghostVolVars_;
+ std::vector<bool> isDirichlet_;
+ Problem* problem_;
};
} // namespace Dumux
diff --git a/dumux/implicit/localjacobian.hh b/dumux/implicit/localjacobian.hh
index 4ba53f954b9386f13cfb224e592c6182307612d6..b5d27dc017023ed930f8f7f8daf77d236c9f5f30 100644
--- a/dumux/implicit/localjacobian.hh
+++ b/dumux/implicit/localjacobian.hh
@@ -23,6 +23,7 @@
#ifndef DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
#define DUMUX_IMPLICIT_LOCAL_JACOBIAN_HH
+#include <dune/istl/io.hh>
#include <dune/istl/matrix.hh>
#include <dumux/common/math.hh>
@@ -145,32 +146,31 @@ public:
// set the current grid element and update the element's
// finite volume geometry
elemPtr_ = &element;
- fvElemGeom_.update(gridView_(), element);
reset_();
- bcTypes_.update(problem_(), element_(), fvElemGeom_);
+ bcTypes_.update(problem_(), element_(), fvElemGeom_());
// set the hints for the volume variables
- model_().setHints(element, prevVolVars_, curVolVars_);
+ // model_().setHints(element, prevVolVars_, curVolVars_);
// update the secondary variables for the element at the last
// and the current time levels
prevVolVars_.update(problem_(),
element_(),
- fvElemGeom_,
+ fvElemGeom_(),
true /* isOldSol? */);
curVolVars_.update(problem_(),
element_(),
- fvElemGeom_,
+ fvElemGeom_(),
false /* isOldSol? */);
// update the hints of the model
- model_().updateCurHints(element, curVolVars_);
+ // model_().updateCurHints(element, curVolVars_);
// calculate the local residual
localResidual().eval(element_(),
- fvElemGeom_,
+ fvElemGeom_(),
prevVolVars_,
curVolVars_,
bcTypes_);
@@ -183,7 +183,7 @@ public:
int numRows, numCols;
if (isBox)
{
- numRows = numCols = fvElemGeom_.numScv;
+ numRows = numCols = fvElemGeom_().numScv;
// resize for hanging nodes or lower dimensional grids
if(numRows > 1<<dim || numCols > 1<<dim)
A_.setSize(numRows, numCols);
@@ -191,7 +191,7 @@ public:
else
{
numRows = 1;
- numCols = fvElemGeom_.numNeighbors;
+ numCols = fvElemGeom_().numNeighbors;
// resize for hanging nodes or lower dimensional grids
if(numCols > 2*dim + 1)
A_.setSize(numRows, numCols);
@@ -309,6 +309,15 @@ protected:
return *problemPtr_;
}
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ Problem &problem_()
+ {
+ Valgrind::CheckDefined(problemPtr_);
+ return *problemPtr_;
+ }
+
/*!
* \brief Returns a reference to the grid view.
*/
@@ -413,15 +422,17 @@ protected:
}
else
{
- neighbor = fvElemGeom_.neighbors[col];
+ neighbor = fvElemGeom_().neighbors[col];
neighborFVGeom.updateInner(neighbor);
dofIdxGlobal = problemPtr_->elementMapper().index(neighbor);
}
- PrimaryVariables priVars(model_().curSol()[dofIdxGlobal]);
- VolumeVariables origVolVars(curVolVars_[col]);
+ auto priVars = model_().curSol()[dofIdxGlobal];
+ // std::cout << "Copy old volVars: " << std::endl;
+ auto origVolVars = curVolVars_[col];
+ // std::cout << "...end copy old volVars." << std::endl;
- curVolVars_[col].setEvalPoint(&origVolVars);
+ // curVolVars_[col].setEvalPoint(&origVolVars);
Scalar eps = asImp_().numericEpsilon(col, pvIdx);
Scalar delta = 0;
@@ -434,11 +445,12 @@ protected:
delta += eps;
// calculate the residual
+ // std::cout << "Compute forward-deflected residual: " << std::endl;
if (isBox)
curVolVars_[col].update(priVars,
problem_(),
element_(),
- fvElemGeom_,
+ fvElemGeom_(),
col,
false);
else
@@ -450,10 +462,11 @@ protected:
false);
localResidual().eval(element_(),
- fvElemGeom_,
+ fvElemGeom_(),
prevVolVars_,
curVolVars_,
bcTypes_);
+ // std::cout << "...end compute forward-deflected residual." << std::endl;
// store the residual and the storage term
partialDeriv = localResidual().residual();
@@ -479,11 +492,12 @@ protected:
delta += eps;
// calculate residual again
+ // std::cout << "Compute backward-deflected residual: " << std::endl;
if (isBox)
curVolVars_[col].update(priVars,
problem_(),
element_(),
- fvElemGeom_,
+ fvElemGeom_(),
col,
false);
else
@@ -495,10 +509,11 @@ protected:
false);
localResidual().eval(element_(),
- fvElemGeom_,
+ fvElemGeom_(),
prevVolVars_,
curVolVars_,
bcTypes_);
+ // std::cout << "...end compute backward-deflected residual." << std::endl;
partialDeriv -= localResidual().residual();
if (isBox || col == 0)
storageDeriv -= localResidual().storageTerm()[col];
@@ -518,7 +533,9 @@ protected:
storageDeriv /= delta;
// restore the original state of the element's volume variables
+ // std::cout << "Restore to orgininal volVars: " << std::endl;
curVolVars_[col] = origVolVars;
+ // std::cout << "...end Restore to orgininal volVars." << std::endl << std::endl;
#if HAVE_VALGRIND
for (unsigned i = 0; i < partialDeriv.size(); ++i)
@@ -544,7 +561,7 @@ protected:
}
}
- for (int i = 0; i < fvElemGeom_.numScv; i++)
+ for (int i = 0; i < fvElemGeom_().numScv; i++)
{
// Green vertices are not to be changed!
if (!isBox || jacAsm_().vertexColor(element_(), i) != Green) {
@@ -558,11 +575,19 @@ protected:
}
}
}
+
+ // std::cout << std::endl << std::endl;
+ // Dune::printmatrix(std::cout, A_, "THE LOCAL JACOBIAN MATRIX", "");
+ // std::cout << std::endl << std::endl;
+
}
- const Element *elemPtr_;
- FVElementGeometry fvElemGeom_;
+ const FVElementGeometry& fvElemGeom_() const
+ {
+ return model_().fvGeometries(problem_().elementMapper().index(element_()));
+ }
+ const Element *elemPtr_;
ElementBoundaryTypes bcTypes_;
// The problem we would like to solve
diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh
index f48ca6355b1b4600f224a8bd43b80789c9f5cc70..98d9f3fe991ac6d6815476a298a3e7cc7b29aee4 100644
--- a/dumux/implicit/localresidual.hh
+++ b/dumux/implicit/localresidual.hh
@@ -85,14 +85,11 @@ public:
*/
void eval(const Element &element)
{
- FVElementGeometry fvGeometry;
-
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
+ fvElemGeomPtr_ = &model_().fvGeometries(model_().elementMapper().index(element));
ElementVolumeVariables volVarsPrev, volVarsCur;
// update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
+ // model_().setHints(element, volVarsPrev, volVarsCur);
volVarsPrev.update(problem_(),
element,
@@ -121,10 +118,7 @@ public:
void evalStorage(const Element &element)
{
elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
+ fvElemGeomPtr_ = &model_().fvGeometries(model_().elementMapper().index(element));
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
@@ -136,7 +130,7 @@ public:
ElementVolumeVariables volVars;
// update the hints
- model_().setHints(element, volVars);
+ // model_().setHints(element, volVars);
// calculate volume current variables
volVars.update(problem_(), element, fvGeometry_(), false);
@@ -157,10 +151,7 @@ public:
const ElementVolumeVariables &curVolVars)
{
elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
+ fvElemGeomPtr_ = &model_().fvGeometries(model_().elementMapper().index(element));
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
@@ -424,6 +415,12 @@ protected:
const Problem &problem_() const
{ return *problemPtr_; }
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ Problem &problem_()
+ { return *problemPtr_; }
+
/*!
* \brief Returns a reference to the model.
*/
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index 40b2d0cec6d2da4667bc65d77c0171a0e4789c3e..f1de5b9228972818167c9bced1c2ca24c0b2b2f6 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -84,7 +84,8 @@ public:
ImplicitModel()
: problemPtr_(0)
{
- enableHints_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableHints);
+ // enableHints_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableHints);
+ // std::cout << enableHints_ << std::endl;
}
/*!
@@ -101,7 +102,8 @@ public:
int numDofs = asImp_().numDofs();
uCur_.resize(numDofs);
- uPrev_.resize(numDofs);
+ volVars_.resize(numDofs);
+ fvGeometries_.resize(numDofs);
if (isBox)
boxVolume_.resize(numDofs);
@@ -111,89 +113,102 @@ public:
asImp_().applyInitialSolution_();
- // resize the hint vectors
- if (isBox && enableHints_) {
- int numVertices = gridView_().size(dim);
- curHints_.resize(numVertices);
- prevHints_.resize(numVertices);
- hintsUsable_.resize(numVertices);
- std::fill(hintsUsable_.begin(),
- hintsUsable_.end(),
- false);
+ // update the volVars with the initial solution
+ for (const auto& element : elements(gridView_()))
+ {
+ int dofIdxGlobal = elementMapper().index(element);
+ volVars_[dofIdxGlobal].update(uCur_[dofIdxGlobal],
+ this->problem_(),
+ element,
+ fvGeometries_[dofIdxGlobal],
+ 0,
+ false);
}
+ // // resize the hint vectors
+ // if (isBox && enableHints_) {
+ // int numVertices = gridView_().size(dim);
+ // curHints_.resize(numVertices);
+ // prevHints_.resize(numVertices);
+ // hintsUsable_.resize(numVertices);
+ // std::fill(hintsUsable_.begin(),
+ // hintsUsable_.end(),
+ // false);
+ // }
+
// also set the solution of the "previous" time step to the
// initial solution.
uPrev_ = uCur_;
+ prevVolVars_ = volVars_;
}
- void setHints(const Element &element,
- ElementVolumeVariables &prevVolVars,
- ElementVolumeVariables &curVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- int n = element.subEntities(dim);
- prevVolVars.resize(n);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i)
- {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
-
- if (!hintsUsable_[vIdxGlobal]) {
- curVolVars[i].setHint(NULL);
- prevVolVars[i].setHint(NULL);
- }
- else {
- curVolVars[i].setHint(&curHints_[vIdxGlobal]);
- prevVolVars[i].setHint(&prevHints_[vIdxGlobal]);
- }
- }
- }
-
- void setHints(const Element &element,
- ElementVolumeVariables &curVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- int n = element.subEntities(dim);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i)
- {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
-
- if (!hintsUsable_[vIdxGlobal])
- curVolVars[i].setHint(NULL);
- else
- curVolVars[i].setHint(&curHints_[vIdxGlobal]);
- }
- }
-
- void updatePrevHints()
- {
- if (!isBox || !enableHints_)
- return;
-
- prevHints_ = curHints_;
- }
-
- void updateCurHints(const Element &element,
- const ElementVolumeVariables &elemVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- for (unsigned int i = 0; i < elemVolVars.size(); ++i)
- {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
- curHints_[vIdxGlobal] = elemVolVars[i];
- if (!hintsUsable_[vIdxGlobal])
- prevHints_[vIdxGlobal] = elemVolVars[i];
- hintsUsable_[vIdxGlobal] = true;
- }
- }
+ // void setHints(const Element &element,
+ // ElementVolumeVariables &prevVolVars,
+ // ElementVolumeVariables &curVolVars) const
+ // {
+ // if (!isBox || !enableHints_)
+ // return;
+
+ // int n = element.subEntities(dim);
+ // prevVolVars_.resize(n);
+ // curVolVars_.resize(n);
+ // for (int i = 0; i < n; ++i)
+ // {
+ // int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
+
+ // if (!hintsUsable_[vIdxGlobal]) {
+ // curVolVars[i].setHint(NULL);
+ // prevVolVars[i].setHint(NULL);
+ // }
+ // else {
+ // curVolVars[i].setHint(&curHints_[vIdxGlobal]);
+ // prevVolVars[i].setHint(&prevHints_[vIdxGlobal]);
+ // }
+ // }
+ // }
+
+ // void setHints(const Element &element,
+ // ElementVolumeVariables &curVolVars) const
+ // {
+ // if (!isBox || !enableHints_)
+ // return;
+
+ // int n = element.subEntities(dim);
+ // curVolVars.resize(n);
+ // for (int i = 0; i < n; ++i)
+ // {
+ // int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
+
+ // if (!hintsUsable_[vIdxGlobal])
+ // curVolVars[i].setHint(NULL);
+ // else
+ // curVolVars[i].setHint(&curHints_[vIdxGlobal]);
+ // }
+ // }
+
+ // void updatePrevHints()
+ // {
+ // if (!isBox || !enableHints_)
+ // return;
+
+ // prevHints_ = curHints_;
+ // }
+
+ // void updateCurHints(const Element &element,
+ // const ElementVolumeVariables &elemVolVars) const
+ // {
+ // if (!isBox || !enableHints_)
+ // return;
+
+ // for (unsigned int i = 0; i < elemVolVars.size(); ++i)
+ // {
+ // int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
+ // curHints_[vIdxGlobal] = elemVolVars[i];
+ // if (!hintsUsable_[vIdxGlobal])
+ // prevHints_[vIdxGlobal] = elemVolVars[i];
+ // hintsUsable_[vIdxGlobal] = true;
+ // }
+ // }
/*!
@@ -456,6 +471,7 @@ public:
if(GET_PROP_VALUE(TypeTag, AdaptiveGrid) && problem_().gridAdapt().wasAdapted())
{
uPrev_ = uCur_;
+ prevVolVars_ = volVars_;
updateBoundaryIndices_();
@@ -466,17 +482,19 @@ public:
jacAsm_->init(problem_());
- // resize the hint vectors
- if (isBox && enableHints_) {
- int numVertices = gridView_().size(dim);
- curHints_.resize(numVertices);
- prevHints_.resize(numVertices);
- hintsUsable_.resize(numVertices);
- std::fill(hintsUsable_.begin(),
- hintsUsable_.end(),
- false);
- }
+ // // resize the hint vectors
+ // if (isBox && enableHints_) {
+ // int numVertices = gridView_().size(dim);
+ // curHints_.resize(numVertices);
+ // prevHints_.resize(numVertices);
+ // hintsUsable_.resize(numVertices);
+ // std::fill(hintsUsable_.begin(),
+ // hintsUsable_.end(),
+ // false);
+ // }
+
}
+
}
@@ -488,6 +506,21 @@ public:
void updateSuccessful()
{ }
+ void updateVolVars()
+ {
+ // update the volVars with the computed solution
+ for (const auto& element : Dune::elements(gridView_()))
+ {
+ int dofIdxGlobal = dofMapper().subIndex(element, 0, dofCodim);
+ this->volVars(dofIdxGlobal).update(curSol()[dofIdxGlobal],
+ this->problem_(),
+ element,
+ fvGeometries_[dofIdxGlobal],
+ 0,
+ false);
+ }
+ }
+
/*!
* \brief Called by the update() method if it was
* unsuccessful. This is primarily a hook which the actual
@@ -499,8 +532,9 @@ public:
// previous time step so that we can start the next
// update at a physically meaningful solution.
uCur_ = uPrev_;
- if (isBox)
- curHints_ = prevHints_;
+ volVars_ = prevVolVars_;
+ // if (isBox)
+ // curHints_ = prevHints_;
jacAsm_->reassembleAll();
}
@@ -516,10 +550,12 @@ public:
{
// make the current solution the previous one.
uPrev_ = uCur_;
- if (isBox)
- prevHints_ = curHints_;
+ prevVolVars_ = volVars_;
+
+ // if (isBox)
+ // prevHints_ = curHints_;
- updatePrevHints();
+ // updatePrevHints();
}
/*!
@@ -846,6 +882,25 @@ public:
VolumeVariables::completeFluidState(priVars, problem, element,
fvGeometry, scvIdx, fluidState);
}
+
+ const VolumeVariables& volVars(unsigned int dofIdxGlobal) const
+ { return volVars_[dofIdxGlobal]; }
+
+ VolumeVariables& volVars(unsigned int dofIdxGlobal)
+ { return volVars_[dofIdxGlobal]; }
+
+ const VolumeVariables& prevVolVars(unsigned int dofIdxGlobal) const
+ { return prevVolVars_[dofIdxGlobal]; }
+
+ VolumeVariables& prevVolVars(unsigned int dofIdxGlobal)
+ { return prevVolVars_[dofIdxGlobal]; }
+
+ const FVElementGeometry& fvGeometries(unsigned int dofIdxGlobal) const
+ { return fvGeometries_[dofIdxGlobal]; }
+
+ FVElementGeometry& fvGeometries(unsigned int dofIdxGlobal)
+ { return fvGeometries_[dofIdxGlobal]; }
+
protected:
/*!
* \brief A reference to the problem on which the model is applied.
@@ -882,14 +937,13 @@ protected:
uCur_ = Scalar(0.0);
boxVolume_ = Scalar(0.0);
- FVElementGeometry fvGeometry;
-
// iterate through leaf grid and evaluate initial
// condition at the center of each sub control volume
for (const auto& element : elements(gridView_())) {
// deal with the current element
- fvGeometry.update(gridView_(), element);
-
+ int eIdx = elementMapper().index(element);
+ fvGeometries_[eIdx].update(gridView_(), element);
+ const auto& fvGeometry = fvGeometries_[eIdx];
// loop over all element vertices, i.e. sub control volumes
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; scvIdx++)
{
@@ -991,9 +1045,9 @@ protected:
// the hint cache for the previous and the current volume
// variables
- mutable std::vector<bool> hintsUsable_;
- mutable std::vector<VolumeVariables> curHints_;
- mutable std::vector<VolumeVariables> prevHints_;
+ // mutable std::vector<bool> hintsUsable_;
+ // mutable std::vector<VolumeVariables> curHints_;
+ // mutable std::vector<VolumeVariables> prevHints_;
// the problem we want to solve. defines the constitutive
// relations, matxerial laws, etc.
@@ -1008,6 +1062,10 @@ protected:
// the set of all indices of vertices on the boundary
std::vector<bool> boundaryIndices_;
+ std::vector<VolumeVariables> volVars_;
+ std::vector<VolumeVariables> prevVolVars_;
+ std::vector<FVElementGeometry> fvGeometries_;
+
// cur is the current iterative solution, prev the converged
// solution of the previous time step
SolutionVector uCur_;
@@ -1027,7 +1085,7 @@ private:
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
- bool enableHints_;
+ // bool enableHints_;
};
} // end namespace Dumux
diff --git a/dumux/implicit/volumevariables.hh b/dumux/implicit/volumevariables.hh
index 99789c2f7689071fe47bc0c37d343cc4b95650f7..d16e258963238472dd8004b6c1203019297f37d0 100644
--- a/dumux/implicit/volumevariables.hh
+++ b/dumux/implicit/volumevariables.hh
@@ -85,20 +85,20 @@ public:
Valgrind::CheckDefined(evalPoint_);
}
- /*!
- * \brief Returns the evaluation point used by the local jacobian.
- *
- * The evaluation point is only used by semi-smooth models.
- */
- const Implementation &evalPoint() const
- { return (evalPoint_ == 0)?asImp_():*evalPoint_; }
-
- /*!
- * \brief Set the volume variables which should be used as initial
- * conditions for complex calculations.
- */
- void setHint(const Implementation *hint)
- {}
+ // /*!
+ // * \brief Returns the evaluation point used by the local jacobian.
+ // *
+ // * The evaluation point is only used by semi-smooth models.
+ // */
+ // const Implementation &evalPoint() const
+ // { return (evalPoint_ == 0)?asImp_():*evalPoint_; }
+
+ // /*!
+ // * \brief Set the volume variables which should be used as initial
+ // * conditions for complex calculations.
+ // */
+ // void setHint(const Implementation *hint)
+ // {}
/*!
* \brief Update all quantities for a given control volume
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 6f43f1b0fd2e456615b91efcbb107172158930d7..e4186aabea3005a0e0fe4c9e67383370ae892809 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -471,6 +471,8 @@ public:
}
}
+ this->model_().updateVolVars();
+
}
/*!
diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh
index 5ea13e75476ec0f22cba24a769cf7ff2e334f09d..2af79e8a4697a0feadbcbb03df6e551a494ff1c4 100644
--- a/dumux/porousmediumflow/2p/implicit/model.hh
+++ b/dumux/porousmediumflow/2p/implicit/model.hh
@@ -155,19 +155,33 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
- int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
- (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
- (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
- (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
- (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
- (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
- (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
- (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
- (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
- (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
- (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
- (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+ int eIdx = this->elementMapper().index(element);
+
+ (*rank)[eIdx] = this->gridView_().comm().rank();
+
+ const auto& fvGeometry = this->fvGeometries(eIdx);
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
+ element,
+ fvGeometry,
+ false /* oldSol? */);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+
+ (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+ (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+ (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
+ (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+ (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+ (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
+ (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
+ (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+ (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+ (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+ (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+ }
}
// velocity output