Skip to content
Snippets Groups Projects
Commit 85c6691e authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[staggered] Allow switching off globalVolVarCaching

parent 6126af8b
No related branches found
No related tags found
3 merge requests!617[WIP] Next,!576Feature/port staggered ff to next next,!571Cleanup/next
Hide whitespace changes
Inline Side-by-side
dumux/discretization/staggered/elementvolumevariables.hh
+ 71
− 27
View file @ 85c6691e
......@@ -93,8 +93,11 @@ class StaggeredElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -103,6 +106,12 @@ class StaggeredElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
public:
//! Constructor
......@@ -115,30 +124,35 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
auto eIdx = globalVolVars().problem_().elementMapper().index(element);
clear();
// stencil information
const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil();
const auto numDofs = neighborStencil.size() + 1;
const auto& problem = globalVolVars().problem();
const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
const auto globalI = fvGridGeometry.elementMapper().index(element);
const auto map = fvGridGeometry.connectivityMap();
const auto& connectivityMapI = map(cellCenterIdx, cellCenterIdx, globalI);
const auto numDofs = connectivityMapI.size();
auto&& scvI = fvGeometry.scv(globalI);
// resize local containers to the required size (for internal elements)
volumeVariables_.resize(numDofs);
volVarIndices_.resize(numDofs);
int localIdx = 0;
// update the volume variables of the element at hand
auto&& scvI = fvGeometry.scv(eIdx);
volumeVariables_[localIdx].update(sol[eIdx], globalVolVars().problem_(), element, scvI);
volVarIndices_[localIdx] = scvI.index();
++localIdx;
// Update the volume variables of the neighboring elements
for (auto globalJ : neighborStencil)
// Update the volume variables of the element at hand and the neighboring elements
for (auto globalJ : connectivityMapI)
{
const auto& elementJ = fvGeometry.fvGridGeometry().element(globalJ);
const auto& elementJ = fvGridGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
volumeVariables_[localIdx].update(sol[globalJ], globalVolVars().problem_(), elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
PrimaryVariables priVars(0.0);
priVars[cellCenterIdx] = sol[cellCenterIdx][globalJ];
ElementSolution elemSol{std::move(priVars)};
volumeVariables_[localIdx].update(elemSol,
problem,
elementJ,
scvJ);
volVarIndices_[localIdx] = scvJ.dofIndex();
++localIdx;
}
......@@ -149,18 +163,33 @@ public:
if (!scvf.boundary())
continue;
// check if boundary is a pure dirichlet boundary
const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf);
const auto bcTypes = problem.boundaryTypes(element, scvf);
PrimaryVariables boundaryPriVars(0.0);
volumeVariables_.resize(localIdx+1);
volVarIndices_.resize(localIdx+1);
volumeVariables_[localIdx].update(dirichletPriVars, globalVolVars().problem_(), element, scvI);
volVarIndices_[localIdx] = scvf.outsideScvIdx();
++localIdx;
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
//TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
// could be made more general by allowing a non-zero-gradient, provided in problem file
else
if(eqIdx == Indices::pressureIdx)
DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
}
volumeVariables_.resize(localIdx+1);
volVarIndices_.resize(localIdx+1);
ElementSolution elemSol{std::move(boundaryPriVars)};
volumeVariables_[localIdx].update(elemSol,
problem,
element,
scvI);
volVarIndices_[localIdx] = scvf.outsideScvIdx();
++localIdx;
}
}
......@@ -170,13 +199,21 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
auto eIdx = globalVolVars().problem_().elementMapper().index(element);
clear();
const auto eIdx = fvGeometry.fvGridGeometry().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element
auto&& scv = fvGeometry.scv(eIdx);
volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
PrimaryVariables priVars(0.0);
priVars[cellCenterIdx] = sol[cellCenterIdx][eIdx];
ElementSolution elemSol{std::move(priVars)};
volumeVariables_[0].update(elemSol,
globalVolVars().problem(),
element,
scv);
volVarIndices_[0] = scv.dofIndex();
}
......@@ -196,6 +233,13 @@ public:
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
//! Clear all local storage
void clear()
{
volVarIndices_.clear();
volumeVariables_.clear();
}
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
......
dumux/discretization/staggered/globalvolumevariables.hh
+ 8
− 13
View file @ 85c6691e
......@@ -23,9 +23,6 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
#include <dumux/common/basicproperties.hh>
#include <dumux/discretization/staggered/elementvolumevariables.hh>
namespace Dumux
{
......@@ -41,9 +38,6 @@ class StaggeredGlobalVolumeVariables
template<class TypeTag>
class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
{
// The local class needs to access and change volVars
friend StaggeredElementVolumeVariables<TypeTag, true>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
......@@ -154,15 +148,15 @@ private:
template<class TypeTag>
class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
// local class needs access to the problem
friend StaggeredElementVolumeVariables<TypeTag, false>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
public:
void update(Problem& problem, const SolutionVector& sol)
{ problemPtr_ = &problem; }
StaggeredGlobalVolumeVariables(const Problem& problem) : problemPtr_(&problem) {}
void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol) {}
/*!
* \brief Return a local restriction of this global object
......@@ -172,11 +166,12 @@ public:
friend inline ElementVolumeVariables localView(const StaggeredGlobalVolumeVariables& global)
{ return ElementVolumeVariables(global); }
private:
Problem& problem_() const
const Problem& problem() const
{ return *problemPtr_;}
Problem* problemPtr_;
private:
const Problem* problemPtr_;
};
} // end namespace
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment