Commit b1acb894 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/freeflow-momentum-new-staggered' into 'master'

Feature/freeflow momentum new staggered

Closes #896

See merge request !2779
parents 75ccf3d1 b7a44b24
Pipeline #8482 passed with stages
in 0 seconds
This diff is collapsed.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Assembly
* \ingroup BoxDiscretization
* \brief Calculates the element-wise residual for the box scheme
*/
#ifndef DUMUX_FACECENTERED_LOCAL_RESIDUAL_HH
#define DUMUX_FACECENTERED_LOCAL_RESIDUAL_HH
#include <dune/geometry/type.hh>
#include <dune/istl/matrix.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
/*!
* \ingroup Assembly
* \ingroup BoxDiscretization
* \brief The element-wise residual for the box scheme
* \tparam TypeTag the TypeTag
*/
template<class TypeTag>
class FaceCenteredLocalResidual : public FVLocalResidual<TypeTag>
{
using ParentType = FVLocalResidual<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
using FVElementGeometry = typename GridGeometry::LocalView;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
using Extrusion = Extrusion_t<GridGeometry>;
public:
using ElementResidualVector = typename ParentType::ElementResidualVector;
using ParentType::ParentType;
//! evaluate flux residuals for one sub control volume face and add to residual
void evalFlux(ElementResidualVector& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
const auto flux = evalFlux(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
residual[insideScv.localDofIndex()] += flux;
}
//! evaluate flux residuals for one sub control volume face
NumEqVector evalFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
if (elemBcTypes.hasDirichlet())
return this->asImp().maybeHandleDirichletBoundary(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
if (elemBcTypes.hasNeumann())
return this->asImp().maybeHandleNeumannBoundary(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
return this->asImp().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache, elemBcTypes);
}
/*!
* \brief Calculate the source term of the equation
*
* \param problem The problem to solve
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
* \param fvGeometry The finite-volume geometry of the element
* \param elemVolVars The volume variables associated with the element stencil
* \param scv The sub-control volume over which we integrate the source term
* \note This is the default implementation for all models as sources are computed
* in the user interface of the problem
*
*/
NumEqVector computeSource(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{
NumEqVector source(0.0);
// add contributions from volume flux sources
source += problem.source(element, fvGeometry, elemVolVars, scv)[scv.dofAxis()];
// add contribution from possible point sources
source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv)[scv.dofAxis()];
return source;
}
using ParentType::evalStorage;
/*!
* \brief Compute the storage local residual, i.e. the deviation of the
* storage term from zero for instationary problems.
*
* \param residual The residual vector to fill
* \param problem The problem to solve
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
* \param fvGeometry The finite-volume geometry of the element
* \param prevElemVolVars The volume averaged variables for all
* sub-control volumes of the element at the previous time level
* \param curElemVolVars The volume averaged variables for all
* sub-control volumes of the element at the current time level
* \param scv The sub control volume the storage term is integrated over
*/
void evalStorage(ElementResidualVector& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
const SubControlVolume& scv) const
{
const auto& curVolVars = curElemVolVars[scv];
const auto& prevVolVars = prevElemVolVars[scv];
// mass balance within the element. this is the
// \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
// euler as time discretization.
//
//! Compute storage with the model specific storage residual
NumEqVector prevStorage = this->asImp().computeStorage(problem, scv, prevVolVars, true/*isPreviousStorage*/);
NumEqVector storage = this->asImp().computeStorage(problem, scv, curVolVars, false/*isPreviousStorage*/);
prevStorage *= prevVolVars.extrusionFactor();
storage *= curVolVars.extrusionFactor();
storage -= prevStorage;
storage *= Extrusion::volume(scv);
storage /= this->timeLoop().timeStepSize();
residual[scv.localDofIndex()] += storage;
}
};
} // end namespace Dumux
#endif
......@@ -37,6 +37,47 @@
#include "diffmethod.hh"
#include "boxlocalassembler.hh"
#include "cclocalassembler.hh"
#include "fclocalassembler.hh"
namespace Dumux::Detail {
template<DiscretizationMethod diffMethod>
struct LocalAssemblerChooser;
template<>
struct LocalAssemblerChooser<DiscretizationMethod::box>
{
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using type = BoxLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
template<>
struct LocalAssemblerChooser<DiscretizationMethod::ccmpfa>
{
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
template<>
struct LocalAssemblerChooser<DiscretizationMethod::cctpfa>
{
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
template<>
struct LocalAssemblerChooser<DiscretizationMethod::fcstaggered>
{
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using type = FaceCenteredLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using LocalAssemblerChooser_t = typename LocalAssemblerChooser<
GetPropType<TypeTag, Properties::GridGeometry>::discMethod
>::template type<TypeTag, Impl, diffMethod, isImplicit>;
} // end namespace Dumux::Detail
namespace Dumux {
......@@ -50,23 +91,22 @@ namespace Dumux {
template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
class FVAssembler
{
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>;
using GridView = typename GridGeo::GridView;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
using Element = typename GridView::template Codim<0>::Entity;
using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
static constexpr DiscretizationMethod discMethod = GetPropType<TypeTag, Properties::GridGeometry>::discMethod;
static constexpr bool isBox = discMethod == DiscretizationMethod::box;
static constexpr bool isBox = GridGeo::discMethod == DiscretizationMethod::box;
using ThisType = FVAssembler<TypeTag, diffMethod, isImplicit>;
using LocalAssembler = std::conditional_t<isBox, BoxLocalAssembler<TypeTag, ThisType, diffMethod, isImplicit>,
CCLocalAssembler<TypeTag, ThisType, diffMethod, isImplicit>>;
using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod, isImplicit>;
public:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridGeometry = GridGeo;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
......
......@@ -43,7 +43,7 @@ Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
pattern.resize(numDofs, numDofs);
// matrix pattern for implicit Jacobians
if (isImplicit)
if constexpr (isImplicit)
{
static constexpr int dim = std::decay_t<decltype(gridGeometry.gridView())>::dimension;
for (const auto& element : elements(gridGeometry.gridView()))
......@@ -204,6 +204,55 @@ template<bool isImplicit, class GridGeometry,
Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
{ return getFEJacobianPattern(gridGeometry.feBasis()); }
/*!
* \ingroup Assembly
* \brief Helper function to generate Jacobian pattern for the face-centered staggered method
*/
template<bool isImplicit, class GridGeometry,
typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethod::fcstaggered) ), int> = 0>
Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
{
// resize the jacobian and the residual
const auto numDofs = gridGeometry.numDofs();
Dune::MatrixIndexSet pattern(numDofs, numDofs);
const auto& connectivityMap = gridGeometry.connectivityMap();
auto fvGeometry = localView(gridGeometry);
// set the pattern
for (const auto& element : elements(gridGeometry.gridView()))
{
fvGeometry.bind(element);
for (const auto& scv : scvs(fvGeometry))
{
const auto globalI = scv.dofIndex();
pattern.add(globalI, globalI);
for (const auto& scvIdxJ : connectivityMap[scv.index()])
{
const auto globalJ = fvGeometry.scv(scvIdxJ).dofIndex();
pattern.add(globalI, globalJ);
if (gridGeometry.isPeriodic())
{
if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
{
const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
pattern.add(globalIP, globalI);
pattern.add(globalI, globalIP);
if (globalI > globalIP)
pattern.add(globalIP, globalJ);
}
}
}
}
}
return pattern;
}
} // namespace Dumux
#endif
......@@ -62,6 +62,15 @@ struct CheckOverlapSize<DiscretizationMethod::fem>
{ return feBasis.gridView().comm().size() <= 1 || feBasis.gridView().overlapSize(0) == 0; }
};
// fc staggered requires an overlap of exactly 1
template<>
struct CheckOverlapSize<DiscretizationMethod::fcstaggered>
{
template<class GridView>
static bool isValid(const GridView& gridView) noexcept
{ return gridView.comm().size() <= 1 || gridView.overlapSize(0) == 1; }
};
} // end namespace Dumux
#endif
......@@ -84,8 +84,7 @@ public:
const SubControlVolumeFace& lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
{
assert(scvf.isLateral());
const auto otherGlobalIdx = scvfIndices_()[GridGeometry::GeometryHelper::lateralOrthogonalScvfLocalIndex(scvf.localIndex())];
return gridGeometry().scvf(otherGlobalIdx);
return gridGeometry().scvf(gridGeometry().lateralOrthogonalScvf(scvf));
}
//! Return the frontal sub control volume face on a the boundary for a given sub control volume
......@@ -328,7 +327,7 @@ public:
const SubControlVolumeFace& lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
{
assert(scvf.isLateral());
const auto otherGlobalIdx = scvfIndices_()[GridGeometry::GeometryHelper::lateralOrthogonalScvfLocalIndex(scvf.localIndex())];
const auto otherGlobalIdx = gridGeometry().lateralOrthogonalScvf(scvf);
return scvfs_[findLocalIndex_(otherGlobalIdx, scvfIndices_())];
}
......
......@@ -138,7 +138,7 @@ public:
, intersectionMapper_(gridView)
{
// Check if the overlap size is what we expect
if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
if (!CheckOverlapSize<DiscretizationMethod::fcstaggered>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
......@@ -229,6 +229,10 @@ public:
const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
{ return periodicFaceMap_; }
//! get the global index of the orthogonal face sharing a common entity
GridIndexType lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
{ return lateralOrthogonalScvf_[scvf.index()]; }
private:
void update_()
......@@ -238,6 +242,7 @@ private:
scvfs_.clear();
scvIndicesOfElement_.clear();
scvfIndicesOfElement_.clear();
lateralOrthogonalScvf_.clear();
intersectionMapper_.update(this->gridView());
// determine size of containers
......@@ -246,6 +251,8 @@ private:
scvfIndicesOfElement_.resize(numElements);
hasBoundaryScvf_.resize(numElements, false);
scvfOfScvInfo_.resize(numElements);
// on frontal + maybe boundary per face and 2*(dim-1) laterals per frontal
lateralOrthogonalScvf_.resize(numElements*(2*dim*(2 + 2*(dim-1))));
outSideBoundaryVolVarIdx_ = 0;
numBoundaryScv_ = 0;
......@@ -269,6 +276,9 @@ private:
// keep track of frontal boundary scvfs
std::size_t numFrontalBoundaryScvfs = 0;
// a temporary map to store pairs of common entities
std::unordered_map<GridIndexType, Dune::ReservedVector<GridIndexType, 2>> commonEntityIdxToScvfsMap;
SmallLocalIndexType localScvIdx = 0;
for (const auto& intersection : intersections(this->gridView(), element))
{
......@@ -287,6 +297,12 @@ private:
continue;
scvfsIndexSet.push_back(scvfIdx);
const auto& lateralFacet = geometryHelper.facet(lateralFacetIndex, element);
const auto commonEntityIdx = geometryHelper.globalCommonEntityIndex(
geometryHelper.facet(localScvIdx, element), lateralFacet
);
commonEntityIdxToScvfsMap[commonEntityIdx].push_back(scvfIdx);
++scvfIdx;
}
}
......@@ -310,6 +326,18 @@ private:
const auto eIdx = this->elementMapper().index(element);
scvIndicesOfElement_[eIdx] = std::move(scvsIndexSet);
scvfIndicesOfElement_[eIdx] = std::move(scvfsIndexSet);
// create the bi-directional map
for (const auto& [_, scvfs] : commonEntityIdxToScvfsMap)
{
// TODO: this maybe be less than 2 if there is a processor boundary
// are we sure that the lateral orthogonal scvf is then never needed?
if (scvfs.size() == 2)
{
lateralOrthogonalScvf_[scvfs[0]] = scvfs[1];
lateralOrthogonalScvf_[scvfs[1]] = scvfs[0];
}
}
}
// reserve memory
......@@ -466,6 +494,9 @@ private:
}
connectivityMap_.update(*this);
// allow to free the non-needed memory
lateralOrthogonalScvf_.resize(scvfs_.size());
}
bool onDomainBoundary_(const typename GridView::Intersection& intersection) const
......@@ -495,6 +526,7 @@ private:
std::vector<bool> hasBoundaryScvf_;
std::vector<std::vector<SmallLocalIndexType>> scvfOfScvInfo_;
std::vector<GridIndexType> lateralOrthogonalScvf_;
std::vector<std::array<GridIndexType, numScvsPerElement>> scvIndicesOfElement_;
std::vector<std::vector<GridIndexType>> scvfIndicesOfElement_;
......@@ -556,7 +588,7 @@ public:
, intersectionMapper_(gridView)
{
// Check if the overlap size is what we expect
if (!CheckOverlapSize<DiscretizationMethod::staggered>::isValid(gridView))
if (!CheckOverlapSize<DiscretizationMethod::fcstaggered>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The staggered discretization method needs at least an overlap of 1 for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
......@@ -636,6 +668,10 @@ public:
const std::unordered_map<GridIndexType, GridIndexType>& periodicVertexMap() const
{ return periodicFaceMap_; }
//! get the global index of the orthogonal face sharing a common entity
GridIndexType lateralOrthogonalScvf(const SubControlVolumeFace& scvf) const
{ return lateralOrthogonalScvf_[scvf.index()]; }
private:
void update_()
......@@ -649,12 +685,15 @@ private:
numBoundaryScvf_ = 0;
hasBoundaryScvf_.clear();
scvfIndicesOfElement_.clear();
lateralOrthogonalScvf_.clear();
outsideVolVarIndices_.clear();
// determine size of containers
const auto numElements = this->gridView().size(0);
scvfIndicesOfElement_.resize(numElements);
hasBoundaryScvf_.resize(numElements, false);
// on frontal + maybe boundary per face and 2*(dim-1) laterals per frontal
lateralOrthogonalScvf_.resize(numElements*(2*dim*(2 + 2*(dim-1))));
GeometryHelper geometryHelper(this->gridView());
......@@ -675,8 +714,10 @@ private:
// keep track of frontal boundary scvfs
std::size_t numFrontalBoundaryScvfs = 0;
SmallLocalIndexType localScvIdx = 0;
// a temporary map to store pairs of common entities
std::unordered_map<GridIndexType, Dune::ReservedVector<GridIndexType, 2>> commonEntityIdxToScvfsMap;
SmallLocalIndexType localScvIdx = 0;
for (const auto& intersection : intersections(this->gridView(), element))
{
++numScvs_;
......@@ -697,6 +738,13 @@ private:
outsideVolVarIndices_[scvfIdx] = neighborVolVarIdx++;
scvfsIndexSet.push_back(scvfIdx);
const auto& lateralFacet = geometryHelper.facet(lateralFacetIndex, element);
const auto commonEntityIdx = geometryHelper.globalCommonEntityIndex(
geometryHelper.facet(localScvIdx, element), lateralFacet
);
commonEntityIdxToScvfsMap[commonEntityIdx].push_back(scvfIdx);
++scvfIdx;
}
}
......@@ -762,6 +810,18 @@ private:
// Save the scv indices belonging to this element to build up fv element geometries fast
scvfIndicesOfElement_[eIdx] = std::move(scvfsIndexSet);
// create the bi-directional map
for (const auto& [_, scvfs] : commonEntityIdxToScvfsMap)
{
// TODO: this maybe be less than 2 if there is a processor boundary
// are we sure that the lateral orthogonal scvf is then never needed?
if (scvfs.size() == 2)
{
lateralOrthogonalScvf_[scvfs[0]] = scvfs[1];
lateralOrthogonalScvf_[scvfs[1]] = scvfs[0];
}
}
}
connectivityMap_.update(*this);
......@@ -793,6 +853,7 @@ private:
std::size_t numBoundaryScvf_;
std::vector<bool> hasBoundaryScvf_;
std::vector<GridIndexType> lateralOrthogonalScvf_;
std::vector<std::vector<GridIndexType>> scvfIndicesOfElement_;
// a map for periodic boundary vertices
......
......@@ -79,72 +79,6 @@ public:
return table[ownLocalFaceIndex];
}
//! Return the local index of a lateral orthogonal scvf
static constexpr int lateralOrthogonalScvfLocalIndex(const SmallLocalIndexType ownLocalScvfIndex)
{
if constexpr(GridView::Grid::dimension == 1)
{
assert(false && "There are no lateral scvfs in 1D");
return -1;
}
if constexpr (GridView::Grid::dimension == 2)
{
switch (ownLocalScvfIndex)
{
case 1: return 7;
case 7: return 1;
case 2: return 10;
case 10: return 2;
case 4: return 8;
case 8: return 4;
case 5: return 11;
case 11: return 5;
default:
{
assert(false && "No lateral orthogonal scvf found");
return -1;
}
}
}
else
{
switch (ownLocalScvfIndex)
{
case 1: return 11;
case 11: return 1;
case 2: return 16;
case 16: return 2;
case 3: return 21;
case 21: return 3;
case 4: return 26;
case 26: return 4;
case 6: return 12;
case 12: return 6;
case 7: return 17;
case 17: return 7;
case 8: return 22;
case 22: return 8;
case 9: return 27;
case 27: return 9;