Commit 5bae6bf1 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[implicit] make vol vars cache optional

The volume variables caching can now be disabled in order for the
simulation to use less memory. The FVGeometry can be bound to an
element which then creates the geometries necessary for operations
on the element or the complete stencil. Several classes now need to
be friend with the model, otherwise the FVGeometry can not be bound
to the element. Maybe the bind function could return a geometry instead
of a reference, then access rights can be limited better. It would decrease
the efficiency though.
parent be79cf6d
......@@ -151,6 +151,9 @@ public:
*/
void assemble(const Element& element, JacobianMatrix& matrix)
{
this->model_().curVolVars_().bind(element);
this->model_().prevVolVars_().bindElement(element);
// set the current grid element and update the element's
// finite volume geometry
globalI_ = this->problem_().elementMapper().index(element);
......
......@@ -53,6 +53,8 @@ class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) };
enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache) };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
......@@ -215,7 +217,7 @@ protected:
// temporary vector to store the Dirichlet boundary fluxes
PrimaryVariables flux(0);
if (!constantBC)
if (!constantBC || !enableVolVarsCache)
{
// update corresponding boundary volume variables before flux calculation
const auto insideScvIdx = scvf.insideScvIdx();
......
......@@ -43,12 +43,17 @@ class FluxVariablesCacheVector
public:
template <typename T = TypeTag>
typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(const Problem& problem)
typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
{
fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
for (const auto& element : elements(problem.gridView()))
{
for (auto&& scvf : problem.model().fvGeometries(element).scvfs())
// bind the geometries and volume variables to the element (all the elements in stencil)
problem.model().fvGeometries_().bind(element);
problem.model().curVolVars_().bind(element);
const auto& fvGeometry = problem.model().fvGeometries(element);
for (const auto& scvf : fvGeometry.scvfs())
{
(*this)[scvf.index()].update(problem, element, scvf);
}
......@@ -56,7 +61,7 @@ public:
}
template <typename T = TypeTag>
typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(const Problem& problem)
typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
{}
template <typename T = TypeTag>
......
......@@ -91,6 +91,10 @@ public:
*/
void eval(const Element &element)
{
// make sure FVElementGeometry and volume variables are bound to the element
model_().curVolVars_().bind(element);
model_().prevVolVars_().bindElement(element);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
......@@ -110,6 +114,10 @@ public:
{
elemPtr_ = &element;
// make sure FVElementGeometry and volume variables are bound to the element
model_().curVolVars_().bindElement(element);
model_().prevVolVars_().bindElement(element);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
bcTypesPtr_ = &bcTypes;
......@@ -129,6 +137,10 @@ public:
{
elemPtr_ = &element;
// make sure FVElementGeometry and volume variables are bound to the element
model_().curVolVars_().bind(element);
model_().prevVolVars_().bindElement(element);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
......
......@@ -27,9 +27,11 @@
#include <dune/istl/bvector.hh>
#include "properties.hh"
#include "localresidual.hh"
#include <dumux/implicit/adaptive/gridadaptproperties.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/parallel/vertexhandles.hh>
#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
namespace Dumux
{
......@@ -46,9 +48,17 @@ template <class TypeTag> class BoxLocalResidual;
template<class TypeTag>
class ImplicitModel
{
friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
friend CCLocalResidual<TypeTag>;
friend BoxLocalResidual<TypeTag>;
friend ImplicitLocalResidual<TypeTag>;
friend PrimaryVariableSwitch<TypeTag>;
friend typename GET_PROP_TYPE(TypeTag, Problem);
friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
friend typename GET_PROP_TYPE(TypeTag, StencilsVector);
friend typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector);
friend typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector);
friend typename GET_PROP_TYPE(TypeTag, FluxVariablesCacheVector);
typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
......@@ -58,7 +68,8 @@ class ImplicitModel
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector) VolumeVariablesVector;
typedef typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector) CurrentVolumeVariablesVector;
typedef typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector) PreviousVolumeVariablesVector;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
......@@ -123,12 +134,12 @@ public:
// resize and update the volVars with the initial solution
curVolVarsVector_.update(problem_(), curSol());
// update the flux variables caches
fluxVarsCacheVector_.update(problem_());
// update stencils
stencilsVector_.update(problem_());
// update the flux variables caches
fluxVarsCacheVector_.update(problem_());
// initialize assembler and create matrix
localJacobian_.init(problem_());
jacAsm_ = std::make_shared<JacobianAssembler>();
......@@ -788,17 +799,19 @@ public:
{ return fvGeometries_->fvGeometry(eIdx); }
protected:
// only to be called by friends and deriving classes
VolumeVariablesVector& curVolVars_()
CurrentVolumeVariablesVector& curVolVars_()
{ return curVolVarsVector_; }
// only to be called by friends and deriving classes
VolumeVariables& curVolVars_(const SubControlVolume& scv)
{ return curVolVarsVector_[scv.index()]; }
VolumeVariables& curVolVars_(unsigned int scvIdx)
{ return curVolVarsVector_[scvIdx]; }
PreviousVolumeVariablesVector& prevVolVars_()
{ return prevVolVarsVector_; }
VolumeVariables& prevVolVars_(const SubControlVolume& scv)
{ return prevVolVarsVector_[scv.index()]; }
......@@ -959,8 +972,8 @@ protected:
SolutionVector uPrev_;
// the current and previous variables (primary and secondary variables)
VolumeVariablesVector curVolVarsVector_;
VolumeVariablesVector prevVolVarsVector_;
CurrentVolumeVariablesVector curVolVarsVector_;
PreviousVolumeVariablesVector prevVolVarsVector_;
// the flux variables cache vector vector
FluxVariablesCacheVector fluxVarsCacheVector_;
......
......@@ -71,8 +71,10 @@ NEW_PROP_TAG(PrimaryVariables); //!< A vector of primary variables within a sub-
NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector of the grid
NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume
NEW_PROP_TAG(EnableVolumeVariablesCache); //!< If disabled, the volume variables are not stored (reduces memory, but is slower)
NEW_PROP_TAG(VolumeVariables); //!< The secondary variables within a sub-control volume
NEW_PROP_TAG(VolumeVariablesVector); //!< The type for a container of volume variables
NEW_PROP_TAG(CurrentVolumeVariablesVector); //!< The type for a container of current volume variables
NEW_PROP_TAG(PreviousVolumeVariablesVector); //!< The type for a container of previous volume variables
NEW_PROP_TAG(FluxVariables); //!< Container storing the different types of flux variables
NEW_PROP_TAG(FluxVariablesCache); //!< Stores data associated with flux vars (if enabled)
NEW_PROP_TAG(FluxVariablesCacheVector); //!< The global vector of flux variable containers
......
......@@ -96,8 +96,11 @@ SET_TYPE_PROP(ImplicitBase, FVElementGeometry, Dumux::FVElementGeometry<TypeTag>
//! The volume variable class, to be overloaded by the model
SET_TYPE_PROP(ImplicitBase, VolumeVariables, ImplicitVolumeVariables<TypeTag>);
//! The global volume variables vector class
SET_TYPE_PROP(ImplicitBase, VolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag>);
//! The global current volume variables vector class
SET_TYPE_PROP(ImplicitBase, CurrentVolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache)>);
//! The global previous volume variables vector class
SET_TYPE_PROP(ImplicitBase, PreviousVolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache)>);
//! The global flux variables cache vector class
SET_TYPE_PROP(ImplicitBase, FluxVariablesCacheVector, Dumux::FluxVariablesCacheVector<TypeTag>);
......@@ -143,6 +146,9 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false);
//! disable partial reassembling by default
SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false);
//! We do not store the volume variables by default
SET_BOOL_PROP(ImplicitBase, EnableVolumeVariablesCache, false);
//! disable flux variables data caching by default
SET_BOOL_PROP(ImplicitBase, EnableFluxVariablesCache, false);
......
......@@ -32,21 +32,42 @@ namespace Dumux
* \ingroup ImplicitModel
* \brief Base class for the volume variables vector
*/
template<class TypeTag>
class VolumeVariablesVector : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables)>
template<class TypeTag, bool useOldSol, bool enableVolVarsCache>
class VolumeVariablesVector
{};
// specialization in case of storing the volume variables
template<class TypeTag, bool useOldSol>
class VolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true> : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables)>
{
friend VolumeVariablesVector<TypeTag, !useOldSol, true>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
VolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const VolumeVariablesVector<TypeTag, useOldSol, true>& other) = default;
VolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const VolumeVariablesVector<TypeTag, !useOldSol, true>& other)
{
// do the copy
numScvs_ = other.numScvs_;
numBoundaryScvs_ = other.numBoundaryScvs_;
volumeVariables_ = other.volumeVariables_;
boundaryVolumeVariables_ = other.boundaryVolumeVariables_;
// return the existing object
return *this;
};
public:
void update(const Problem& problem, const SolutionVector& sol)
void update(Problem& problem, const SolutionVector& sol)
{
numScvs_ = problem.model().fvGeometries().numScv();
if (!isBox)
......@@ -109,6 +130,12 @@ public:
return boundaryVolumeVariables_[scvIdx-numScvs_];
}
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
void bind(const Element& element) const {}
// function to prepare the vol vars within the element
void bindElement(const Element& element) const {}
private:
IndexType numScvs_;
IndexType numBoundaryScvs_;
......@@ -116,6 +143,255 @@ private:
std::vector<VolumeVariables> boundaryVolumeVariables_;
};
// Specialization when the current volume variables are not stored
template<class TypeTag>
class VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>
{
// prev vol vars have to be a friend class in order for the assignment operator to work
friend VolumeVariablesVector<TypeTag, true, false>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
VolumeVariablesVector& operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other) = default;
void operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>& other)
{
eIdxBound_ = -1;
}
VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
public:
void update(Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
}
// Binding of an element, prepares the volume variables within the element stencil
// called by the local jacobian to prepare element assembly
// specialization for cc models
template <typename T = TypeTag>
typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
bind(const Element& element)
{
const auto eIdx = problem_().elementMapper().index(element);
if (eIdx == eIdxBound_ && volVarIndices_.size() > 1)
return;
eIdxBound_ = eIdx;
const auto& fvGeometry = problem_().model().fvGeometries(eIdx);
// get stencil information
const auto& elementStencil = problem_().model().stencils(element).elementStencil();
const auto& neighborStencil = problem_().model().stencils(element).neighborStencil();
const auto numDofs = elementStencil.size();
// resize volume variables to the required size
volumeVariables_.resize(numDofs);
volVarIndices_.resize(numDofs);
int localIdx = 0;
// update the volume variables of the element at hand
const auto& scvI = problem_().model().fvGeometries().subControlVolume(eIdx);
const auto& solI = problem_().model().curSol()[eIdx];
volumeVariables_[localIdx].update(solI, problem_(), element, scvI);
volVarIndices_[localIdx] = scvI.index();
localIdx++;
// Update the volume variables of the neighboring elements and the boundary
for (auto globalJ : neighborStencil)
{
const auto& elementJ = problem_().model().fvGeometries().element(globalJ);
const auto& scvJ = problem_().model().fvGeometries().subControlVolume(globalJ);
const auto& solJ = problem_().model().curSol()[globalJ];
volumeVariables_[localIdx].update(solJ, problem_(), elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
localIdx++;
}
// Check for boundary volume variables
for (const auto& scvFace : fvGeometry.scvfs())
{
// if we are not on a boundary, skip the rest
if (!scvFace.boundary())
continue;
// When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries
const auto bcTypes = problem_().boundaryTypes(element, scvFace);
if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann()))
continue;
volumeVariables_.resize(localIdx+1);
volVarIndices_.resize(localIdx+1);
const auto dirichletPriVars = problem_().dirichlet(element, scvFace);
volumeVariables_[localIdx].update(dirichletPriVars, problem_(), element, scvI);
volVarIndices_[localIdx] = scvFace.outsideScvIdx();
localIdx++;
}
}
// Binding of an element, prepares only the volume variables of the element
// specialization for cc models
template <typename T = TypeTag>
typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
bindElement(const Element& element)
{
const auto eIdx = problem_().elementMapper().index(element);
if (eIdx == eIdxBound_ || std::find(volVarIndices_.begin(), volVarIndices_.end(), eIdx) != volVarIndices_.end())
return;
volumeVariables_.resize(1);
volVarIndices_.resize(1);
eIdxBound_ = eIdx;
// update the volume variables of the element at hand
const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
const auto& sol = problem_().model().curSol()[eIdx];
volumeVariables_[0].update(sol, problem_(), element, scv);
volVarIndices_[0] = scv.index();
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
return volumeVariables_[getLocalIdx_(scvIdx)];
}
VolumeVariables& operator [](IndexType scvIdx)
{
return volumeVariables_[getLocalIdx_(scvIdx)];
}
private:
void release_()
{
volumeVariables_.clear();
volVarIndices_.clear();
}
const int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
if (it != volVarIndices_.end())
return std::distance(volVarIndices_.begin(), it);
else
DUNE_THROW(Dune::InvalidStateException, "Could not find the current volume variables for volVarIdx = " << volVarIdx <<
", make sure to properly bind the volume variables to the element before using them");
}
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
IndexType eIdxBound_;
std::vector<IndexType> volVarIndices_;
std::vector<VolumeVariables> volumeVariables_;
};
// Specialization when the previous volume variables are not stored
template<class TypeTag>
class VolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>
{
// current vol vars have to be a friend class in order for the assignment operator to work
friend VolumeVariablesVector<TypeTag, false, false>;
friend ImplicitModel<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
VolumeVariablesVector& operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other)
{
release_();
problemPtr_ = other.problemPtr_;
eIdxBound_ = -1;
return *this;
}
VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
public:
void update(const Problem& problem, const SolutionVector& sol)
{
problemPtr_ = &problem;
}
// Binding of an element, prepares the volume variables of only the element
// specialization for cc models
template <typename T = TypeTag>
typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
bindElement(const Element& element)
{
const auto eIdx = problem_().elementMapper().index(element);
if (eIdx == eIdxBound_)
return;
volumeVariables_.resize(1);
volVarIndices_.resize(1);
eIdxBound_ = eIdx;
// update the volume variables of the element at hand
const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
const auto& sol = problem_().model().prevSol()[eIdx];
volumeVariables_[0].update(sol, problem_(), element, scv);
volVarIndices_[0] = scv.index();
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
return volumeVariables_[getLocalIdx_(scvIdx)];
}
VolumeVariables& operator [](IndexType scvIdx)
{
return volumeVariables_[getLocalIdx_(scvIdx)];
}
private:
void release_()
{
volumeVariables_.clear();
volVarIndices_.clear();
}
const int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
if (it != volVarIndices_.end())
return std::distance(volVarIndices_.begin(), it);
else
DUNE_THROW(Dune::InvalidStateException, "Could not find the previous volume variables for volVarIdx = " << volVarIdx <<
", make sure to properly bind the volume variables to the element before using them");
}
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
IndexType eIdxBound_;
std::vector<IndexType> volVarIndices_;
std::vector<VolumeVariables> volumeVariables_;
};
} // end namespace
#endif
......@@ -93,7 +93,10 @@ public:
int eIdx = this->elementMapper().index(element);
(*rank)[eIdx] = this->gridView_().comm().rank();
for (auto&& scv : this->fvGeometries(element).scvs())
this->curVolVars_().bindElement(element);
const auto& fvGeometry = this->fvGeometries(element);
for (const auto& scv : fvGeometry.scvs())
{
const auto& spatialParams = this->problem_().spatialParams();
auto dofIdxGlobal = scv.dofIndex();
......
......@@ -157,7 +157,10 @@ public:
int eIdx = this->elementMapper().index(element);
(*rank)[eIdx] = this->gridView_().comm().rank();
for (auto&& scv : this->fvGeometries(element).scvs())
this->curVolVars_().bindElement(element);
const auto& fvGeometry = this->fvGeometries(element);
for (const auto& scv : fvGeometry.scvs())
{
auto dofIdxGlobal = scv.dofIndex();
const auto& volVars = this->curVolVars(scv);
......
......@@ -138,30 +138,24 @@ public:
{
const auto& fvGeometry = this->fvGeometries(element);
for (auto&& scv : fvGeometry.scvs())
this->curVolVars_().bindElement(element);
const auto& fvGeometry = this->fvGeometries(element);
for (const auto& scv : fvGeometry.scvs())
{
auto eIdx = scv.elementIndex();
auto dofIdxGlobal = scv.dofIndex();
(*rank)[eIdx] = this->gridView_().comm().rank();
const auto& volVars = this->curVolVars(scv);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
{
(*saturation[phaseIdx])[dofIdxGlobal] = volVars.saturation(phaseIdx);
(*pressure[phaseIdx])[dofIdxGlobal] = volVars.pressure(phaseIdx);
(*density[phaseIdx])[dofIdxGlobal] = volVars.density(phaseIdx);
}
(*poro)[dofIdxGlobal] = volVars.porosity();
(*perm)[dofIdxGlobal] = volVars.permeability();
(*temperature)[dofIdxGlobal] = volVars.temperature();
// velocity output
// velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
// velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
// velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
(*saturation[phaseIdx])[dofIdxGlobal] = volVars.saturation(phaseIdx);
(*pressure[phaseIdx])[dofIdxGlobal] = volVars.pressure(phaseIdx);