Skip to content
GitLab
Commit 00794bdc authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[mpfa] complete makeover and l-method introduction

parent a35aac59
Expand all Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/darcyslaw.hh
View file @ 00794bdc
......@@ -33,6 +33,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/cellcentered/mpfa/methods.hh>
namespace Dumux
......@@ -53,6 +54,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
......@@ -64,6 +66,8 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
static const MpfaMethods method = GET_PROP_VALUE(TypeTag, MpfaMethod);
public:
static Scalar flux(const Problem& problem,
......@@ -79,19 +83,10 @@ public:
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions(phaseIdx);
const auto& tij = fluxVarsCache.advectionTij(phaseIdx);
// interface density as arithmetic mean of the two neighbors (when gravity is on)
// interface density as arithmetic mean of the neighbors (when gravity is on)
Scalar rho;
if (gravity)
{
if (!scvf.boundary())
{
rho = elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
rho += elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
rho /= 2.0;
}
else
rho = elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
}
rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
// calculate Tij*pj
Scalar flux(0.0);
......@@ -110,7 +105,6 @@ public:
h -= rho*(g*x);
}
flux += tij[localIdx++]*h;
}
......@@ -130,6 +124,27 @@ public:
else
return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices();
}
private:
static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx)
{
// use arithmetic mean of the densities around the scvf
if (!scvf.boundary())
{
const auto& outsideScvIndices = scvf.outsideScvIndices();
Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
for (auto outsideIdx : outsideScvIndices)
rho += elemVolVars[outsideIdx].density(phaseIdx);
rho /= outsideScvIndices.size();
return rho;
}
else
return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
}
};
} // end namespace
......
dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
View file @ 00794bdc
......@@ -44,6 +44,9 @@ class CCMpfaElementFluxVariablesCache;
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, true>
{
// the local jacobian needs to be able to update the cache during assembly
friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
......@@ -68,11 +71,6 @@ public:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
// Specialization for the global caching being enabled - do nothing here
void update(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
// Specialization for the global caching being enabled - do nothing here
void bindScvf(const Element& element,
const FVElementGeometry& fvGeometry,
......@@ -89,6 +87,11 @@ public:
private:
const GlobalFluxVariablesCache* globalFluxVarsCachePtr_;
// Specialization for the global caching being enabled - do nothing here
void update(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
};
/*!
......@@ -98,6 +101,9 @@ private:
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, false>
{
// the local jacobian needs to be able to update the cache during assembly
friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
......@@ -119,6 +125,7 @@ public:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
// TODO
DUNE_THROW(Dune::InvalidStateException, "Does local flux var cache binding make sense in general?");
}
......@@ -129,52 +136,53 @@ public:
const ElementVolumeVariables& elemVolVars)
{
clear();
const auto& problem = globalFluxVarsCache().problem_();
const auto& globalFvGeometry = problem.model().globalFvGeometry();
const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
const auto& assemblyMap = problem.model().localJacobian().assemblyMap();
const auto globalI = problem.elementMapper().index(element);
// // first prepare all the global indices
// reserve initial guess of memory (won't be enough though - several scvfs per neighbor will be required)
globalScvfIndices_.reserve(fvGeometry.numScvf() + assemblyMap[globalI].size());
// first add all the indices inside the element
for (auto&& scvf : scvfs(fvGeometry))
{
const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
globalScvfIndices_.push_back(scvf.index());
// for the indices in the neighbors, use assembly map of the local jacobian
for (unsigned int j = 0; j < neighborStencil.size(); ++j)
for (auto fluxVarIdx : assemblyMap[globalI][j])
globalScvfIndices_.push_back(fluxVarIdx);
if (boundary)
{
const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
globalScvfIndices_.push_back(scvfIdx);
}
else
{
const auto& ivSeed = globalFvGeometry.interactionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
globalScvfIndices_.push_back(scvfIdx);
}
}
// make global indices unique
std::sort(globalScvfIndices_.begin(), globalScvfIndices_.end());
globalScvfIndices_.erase(std::unique(globalScvfIndices_.begin(), globalScvfIndices_.end()), globalScvfIndices_.end());
// prepare all the caches of the scvfs inside the corresponding interaction volume using helper class
// prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class
fluxVarsCache_.resize(globalScvfIndices_.size());
for (auto&& scvf : scvfs(fvGeometry))
{
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this);
}
}
// This function updates the transmissibilities after the solution has been deflected during jacobian assembly
void update(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
for (auto&& scvf : scvfs(fvGeometry))
(*this)[scvf].setOutdated();
for (auto&& scvf : scvfs(fvGeometry))
// prepare the caches in the remaining neighbors
unsigned int j = 0;
for (auto globalJ : neighborStencil)
{
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
for (auto fluxVarIdx : assemblyMap[globalI][j])
{
const auto& scvf = fvGeometry.scvf(fluxVarIdx);
if (!(*this)[scvf].isUpdated())
{
auto elementJ = globalFvGeometry.element(globalJ);
FluxVariablesCacheFiller::fillFluxVarCache(problem, elementJ, fvGeometry, elemVolVars, scvf, *this);
}
}
// increment counter
j++;
}
}
......@@ -183,6 +191,7 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
{
// TODO
DUNE_THROW(Dune::InvalidStateException, "Does binding of one scvf make sense in general?");
}
......@@ -212,11 +221,25 @@ private:
globalScvfIndices_.clear();
}
// This function updates the transmissibilities after the solution has been deflected during jacobian assembly
void update(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
for (auto&& scvf : scvfs(fvGeometry))
(*this)[scvf].setUpdateStatus(false);
for (auto&& scvf : scvfs(fvGeometry))
{
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
}
}
// get index of scvf in the local container
int getLocalScvfIdx_(const int scvfIdx) const
{
auto it = std::lower_bound(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
assert(it != globalScvfIndices_.end() && "Could not find the flux vars cache for scvfIdx");
assert(globalScvfIndices_[std::distance(globalScvfIndices_.begin(), it)] == scvfIdx && "Could not find the flux vars cache for scvfIdx");
return std::distance(globalScvfIndices_.begin(), it);
}
......
dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
View file @ 00794bdc
......@@ -91,6 +91,7 @@ class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
......@@ -114,7 +115,7 @@ public:
const SolutionVector& sol)
{
const auto& problem = globalVolVars().problem_();
auto eIdx = problem.elementMapper().index(element);
const auto& globalFvGeometry = problem.model().globalFvGeometry();
// stencil information
const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
......@@ -126,6 +127,7 @@ public:
int localIdx = 0;
// update the volume variables of the element at hand
auto eIdx = problem.elementMapper().index(element);
auto&& scvI = fvGeometry.scv(eIdx);
volumeVariables_[localIdx].update(sol[eIdx], problem, element, scvI);
volVarIndices_[localIdx] = scvI.index();
......@@ -134,86 +136,92 @@ public:
// Update the volume variables of the neighboring elements
for (auto globalJ : neighborStencil)
{
const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ);
const auto& elementJ = globalFvGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
volumeVariables_[localIdx].update(sol[globalJ], problem, elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
++localIdx;
}
std::size_t neighborScvfEstimate = neighborStencil.size()*dim*(2*dim-2);
std::vector<IndexType> finishedBoundaries;
finishedBoundaries.reserve(neighborScvfEstimate);
// if the element is connected to a boundary, additionally prepare BCs
bool boundary = element.hasBoundaryIntersections();
if (!boundary)
for (auto&& scvf : scvfs(fvGeometry))
if (globalFvGeometry.scvfTouchesBoundary(scvf))
{
boundary = true;
break;
}
// if the element is connected to a boundary, additionally treat the BC
if (element.hasBoundaryIntersections())
if (boundary)
{
// reserve enough space for the boundary volume variables
volumeVariables_.reserve(numDofs+neighborScvfEstimate);
volVarIndices_.reserve(numDofs+neighborScvfEstimate);
for (auto&& scvf : scvfs(fvGeometry))
// reserve memory (12 is the case of 3d hexahedron with one face on boundary)
std::size_t estimate = 12;
std::vector<IndexType> finishedBoundaries;
finishedBoundaries.reserve(estimate);
volumeVariables_.reserve(numDofs+estimate);
volVarIndices_.reserve(numDofs+estimate);
// treat the BCs inside the element
if (element.hasBoundaryIntersections())
{
// if we are not on a boundary, skip to the next scvf
if (!scvf.boundary())
continue;
// boundary volume variables
VolumeVariables dirichletVolVars;
const auto dirichletPriVars = problem.dirichlet(element, scvf);
dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
// boundary vol var index
auto bVolVarIdx = scvf.outsideScvIdx();
volVarIndices_.push_back(bVolVarIdx);
finishedBoundaries.push_back(bVolVarIdx);
for (auto&& scvf : scvfs(fvGeometry))
{
// if we are not on a boundary, skip to the next scvf
if (!scvf.boundary())
continue;
// boundary volume variables
VolumeVariables dirichletVolVars;
const auto dirichletPriVars = problem.dirichlet(element, scvf);
dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
// boundary vol var index
auto bVolVarIdx = scvf.outsideScvIdx();
volVarIndices_.push_back(bVolVarIdx);
finishedBoundaries.push_back(bVolVarIdx);
}
}
}
// Update boundary volume variables in the neighbors
const auto& globalFvGeometry = problem.model().globalFvGeometry();
for (auto&& scvf : scvfs(fvGeometry))
{
// reserve enough space for the boundary volume variables
volumeVariables_.reserve(numDofs+neighborScvfEstimate);
volVarIndices_.reserve(numDofs+neighborScvfEstimate);
// skip the rest if the scvf does not touch a boundary
const bool boundary = globalFvGeometry.scvfTouchesBoundary(scvf);
if (!boundary)
continue;
// loop over all the scvfs in the interaction region
const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
// Update boundary volume variables in the neighbors
for (auto&& scvf : scvfs(fvGeometry))
{
auto&& ivScvf = fvGeometry.scvf(scvfIdx);
if (!ivScvf.boundary() || contains(ivScvf.outsideScvIdx(), finishedBoundaries))
// skip the rest if the scvf does not touch a boundary
if (!globalFvGeometry.scvfTouchesBoundary(scvf))
continue;
// that means we are on a not yet handled boundary scvf
auto insideScvIdx = ivScvf.insideScvIdx();
auto&& ivScv = fvGeometry.scv(insideScvIdx);
auto ivElement = globalFvGeometry.element(insideScvIdx);
// boundary volume variables
VolumeVariables dirichletVolVars;
const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf);
dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
// boundary vol var index
auto bVolVarIdx = ivScvf.outsideScvIdx();
volVarIndices_.push_back(bVolVarIdx);
finishedBoundaries.push_back(bVolVarIdx);
// loop over all the scvfs in the interaction region
const auto& ivSeed = globalFvGeometry.boundaryInteractionVolumeSeed(scvf);
for (auto scvfIdx : ivSeed.globalScvfIndices())
{
auto&& ivScvf = fvGeometry.scvf(scvfIdx);
if (!ivScvf.boundary() || MpfaHelper::contains(finishedBoundaries, ivScvf.outsideScvIdx()))
continue;
// that means we are on a not yet handled boundary scvf
auto insideScvIdx = ivScvf.insideScvIdx();
auto&& ivScv = fvGeometry.scv(insideScvIdx);
auto ivElement = globalFvGeometry.element(insideScvIdx);
// boundary volume variables
VolumeVariables dirichletVolVars;
const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf);
dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
// boundary vol var index
auto bVolVarIdx = ivScvf.outsideScvIdx();
volVarIndices_.push_back(bVolVarIdx);
finishedBoundaries.push_back(bVolVarIdx);
}
}
}
// free unused memory
volumeVariables_.shrink_to_fit();
volVarIndices_.shrink_to_fit();
// free unused memory
volumeVariables_.shrink_to_fit();
volVarIndices_.shrink_to_fit();
}
}
// Binding of an element, prepares only the volume variables of the element
......@@ -258,10 +266,6 @@ private:
return std::distance(volVarIndices_.begin(), it);
}
const bool contains(const IndexType idx,
const std::vector<IndexType>& indices) const
{ return std::find(indices.begin(), indices.end(), idx) != indices.end(); }
std::vector<IndexType> volVarIndices_;
std::vector<VolumeVariables> volumeVariables_;
};
......
dumux/discretization/cellcentered/mpfa/facetypes.hh
View file @ 00794bdc
......@@ -29,9 +29,7 @@ namespace Dumux
{
interior,
neumann,
dirichlet,
interiorNeumann,
interiorDirichlet
dirichlet
};
} // end namespace
......
dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
View file @ 00794bdc
......@@ -108,7 +108,7 @@ public:
{
if (!updateNeumannOnly)
cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
cache.setUpdated();
cache.setUpdateStatus(true);
}
// update flux variable caches of the other scvfs of the interaction volume
......@@ -124,7 +124,7 @@ public:
{
if (!updateNeumannOnly)
cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
cacheJ.setUpdated();
cacheJ.setUpdateStatus(true);
}
}
}
......@@ -140,7 +140,7 @@ public:
const auto scvfIdx = scvf.index();
auto& cache = fluxVarsCache[scvfIdx];
cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
cache.setUpdated();
cache.setUpdateStatus(true);
// update flux variable caches of the other scvfs of the interaction volume
for (const auto scvfIdxJ : iv.globalScvfs())
......@@ -151,7 +151,7 @@ public:
const auto elementJ = problem.model().globalFvGeometry().element(scvfJ.insideScvIdx());
auto& cacheJ = fluxVarsCache[scvfIdxJ];
cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
cacheJ.setUpdated();
cacheJ.setUpdateStatus(true);
}
}
}
......@@ -168,7 +168,7 @@ public:
{
if (!solDependentParams)
{
fluxVarsCache[scvf.index()].setUpdated();
fluxVarsCache[scvf.index()].setUpdateStatus(true);
// if we do not use tpfa on the boundary, we have to update the neumann fluxes
if (!useTpfaBoundary)
......
dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
View file @ 00794bdc
......@@ -30,8 +30,6 @@
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include "mpfageometryhelper.hh"
namespace Dumux
{
/*!
......@@ -148,6 +146,7 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
......@@ -157,10 +156,10 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
using CoordScalar = typename GridView::ctype;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
using MpfaGeometryHelper = Dumux::MpfaGeometryHelper<GridView, dim>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
//! Constructor
......@@ -226,36 +225,34 @@ public:
{
bindElement(element);
// get some references for convenience
const auto& problem = globalFvGeometry().problem_();
const auto& assemblyMap = problem.model().localJacobian().assemblyMap();
const auto& neighborStencil = problem.model().stencils(element).neighborStencil();
const auto globalI = problem.elementMapper().index(element);
// reserve memory
auto numNeighbors = globalFvGeometry().problem_().model().stencils(element).neighborStencil().size();
auto numNeighborScvf = numNeighbors*dim*(2*dim-2);
const auto numNeighbors = neighborStencil.size();
const auto estimate = numNeighbors*dim*(4*dim-4); // estimate holds for quadrilaterals in 2D/3D, overestimates else
neighborScvs_.reserve(numNeighbors);
neighborScvIndices_.reserve(numNeighbors);
neighborScvfIndices_.reserve(numNeighborScvf);
neighborScvfs_.reserve(numNeighborScvf);
// build scvfs connected to the interaction regions around the element vertices
std::vector<IndexType> finishedVertices;
finishedVertices.reserve(element.geometry().corners());
const auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
for (auto&& scvf : scvfs_)
neighborScvfIndices_.reserve(estimate);
neighborScvfs_.reserve(estimate);
// build scvfs in neighbors
// use the assembly map to determine which faces are necessary
unsigned int j = 0;
for (auto globalJ : neighborStencil)
{
// skip the rest if scvf has been handled already
if (contains(scvf.vertexIndex(), finishedVertices))
continue;
// get data on the neighbor
auto elementJ = globalFvGeometry().element(globalJ);
const auto& scvfIndices = assemblyMap[globalI][j];
if (globalFvGeometry().scvfTouchesBoundary(scvf))
{
const auto& ivSeed = globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
handleInteractionRegion(ivSeed);
}
else
{
const auto& ivSeed = globalFvGeometry().interactionVolumeSeed(scvf);
handleInteractionRegion(ivSeed);
}
// make neighbor geometries
makeNeighborGeometries(elementJ, globalJ, scvfIndices);