Commit 4aff4c15 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[staggered] Improve handling of fixed Dirichlet cell

* allow setting a fixed value anywhere in the domain
* do not use boundaryType but a separate function to check for fixed cell
* set 1.0 on the diagonal in the system matrix and 0.0 for the off-diagonal entries
  to improve numerical accuracy
parent ae23491c
......@@ -82,6 +82,19 @@ public:
: ParentType(fvGridGeometry, paramGroup)
{ }
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{ return false; }
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume (-face).
......
......@@ -55,35 +55,10 @@ public:
ParentType::resetEq(eqIdx);
boundaryInfo_[eqIdx].visited = false;
boundaryInfo_[eqIdx].isDirichletCell = false;
boundaryInfo_[eqIdx].isSymmetry = false;
boundaryInfo_[eqIdx].isBJS = false;
}
/*!
* \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
* This is a provisional alternative to setting the Dirichlet value on the boundary directly.
*
* \param eqIdx The index of the equation which should used to set
* the Dirichlet condition
*/
void setDirichletCell(int eqIdx)
{
resetEq(eqIdx);
boundaryInfo_[eqIdx].visited = true;
boundaryInfo_[eqIdx].isDirichletCell = true;
}
/*!
* \brief Returns true if an equation is used to specify a
* Dirichlet condition.
*
* \param eqIdx The index of the equation
*/
bool isDirichletCell(unsigned eqIdx) const
{ return boundaryInfo_[eqIdx].isDirichletCell; }
/*!
* \brief Sets a symmetry boundary condition for all equations
*/
......@@ -158,7 +133,6 @@ public:
protected:
struct StaggeredFreeFlowBoundaryInfo {
bool visited;
bool isDirichletCell;
bool isSymmetry;
bool isBJS;
};
......
......@@ -97,37 +97,6 @@ public:
return storage;
}
/*!
* \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
* This is a provisional alternative to setting the Dirichlet value on the boundary directly.
*/
template<class ElementVolumeVariables, class BoundaryTypes>
void setFixedCell(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
const SubControlVolume& insideScv,
const ElementVolumeVariables& elemVolVars,
const BoundaryTypes& bcTypes) const
{
ParentType::setFixedCell(residual, problem, element, insideScv, elemVolVars, bcTypes);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// get equation index
const auto eqIdx = Indices::conti0EqIdx + compIdx;
// set a fixed mole fraction for cells
if(eqIdx != Indices::conti0EqIdx && bcTypes.isDirichletCell(eqIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(compIdx) : insideVolVars.massFraction(compIdx);
residual[eqIdx - cellCenterOffset] = massOrMoleFraction - problem.dirichlet(element, insideScv)[eqIdx];
}
}
}
};
}
......
......@@ -185,26 +185,6 @@ public:
return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
}
/*!
* \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
* This is a provisional alternative to setting the Dirichlet value on the boundary directly.
*/
template<class BoundaryTypes>
void setFixedCell(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
const SubControlVolume& insideScv,
const ElementVolumeVariables& elemVolVars,
const BoundaryTypes& bcTypes) const
{
// set a fixed pressure for cells adjacent to a wall
if(bcTypes.isDirichletCell(Indices::pressureIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
residual[Indices::pressureIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichlet(element, insideScv)[Indices::pressureIdx];
}
}
protected:
/*!
......@@ -258,10 +238,6 @@ protected:
// add the flux over the boundary scvf to the residual
residual += boundaryFlux;
}
// if specified, set a fixed value at the center of a cell at the boundary
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
asImp_().setFixedCell(residual, problem, element, scv, elemVolVars, bcTypes);
}
}
}
......
......@@ -128,25 +128,6 @@ public:
return source;
}
/*!
* \brief Sets a fixed Dirichlet value for a cell (such as for dissipation) at the boundary.
*/
template<class BoundaryTypes>
void setFixedCell(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
const SubControlVolume& insideScv,
const ElementVolumeVariables& elemVolVars,
const BoundaryTypes& bcTypes) const
{
// set a fixed dissipation for cells adjacent to a wall
if(bcTypes.isDirichletCell(Indices::dissipationEqIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
residual[dissipationEqIdx] = insideVolVars.dissipation() - problem.dirichlet(element, insideScv)[Indices::dissipationEqIdx];
}
}
};
}
......
......@@ -28,7 +28,6 @@
#include <dune/common/reservedvector.hh>
#include <dune/grid/common/gridenums.hh> // for GhostEntity
#include <dune/istl/matrixindexset.hh>
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
......@@ -57,18 +56,12 @@ class SubDomainStaggeredLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag
using ParentType = FVLocalAssemblerBase<TypeTag, Assembler,Implementation, isImplicit>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using ElementResidualVector = typename LocalResidual::ElementResidualVector;
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using SolutionVector = typename Assembler::SolutionVector;
using SubSolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using Scalar = typename GridVariables::Scalar;
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
......@@ -77,7 +70,6 @@ class SubDomainStaggeredLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag
using FVGridGeometry = typename GridVariables::GridGeometry;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
......@@ -91,7 +83,8 @@ public:
static constexpr auto cellCenterId = typename Dune::index_constant<0>();
static constexpr auto faceId = typename Dune::index_constant<1>();
static constexpr auto faceOffset = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
static constexpr auto numEqCellCenter = CellCenterResidualValue::dimension;
static constexpr auto faceOffset = numEqCellCenter;
using ParentType::ParentType;
......@@ -170,7 +163,17 @@ public:
if (!this->assembler().isStationaryProblem())
residual += evalLocalStorageResidualForCellCenter();
this->localResidual().evalBoundaryForCellCenter(residual, this->problem(), this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache());
this->localResidual().evalBoundaryForCellCenter(residual, problem(), this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache());
// handle cells with a fixed Dirichlet value
const auto cellCenterGlobalI = problem().fvGridGeometry().elementMapper().index(this->element());
const auto& scvI = this->fvGeometry().scv(cellCenterGlobalI);
for (int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx)
{
static constexpr auto offset = numEq - numEqCellCenter;
if (this->problem().isDirichletCell(this->element(), this->fvGeometry(), scvI, pvIdx + offset))
residual[pvIdx] = this->curSol()[cellCenterId][cellCenterGlobalI][pvIdx] - this->problem().dirichlet(this->element(), scvI)[pvIdx + offset];
}
return residual;
}
......@@ -243,7 +246,7 @@ public:
if (!this->assembler().isStationaryProblem())
residual += evalLocalStorageResidualForFace(scvf);
this->localResidual().evalBoundaryForFace(residual, this->problem(), this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache(), scvf);
this->localResidual().evalBoundaryForFace(residual, problem(), this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache(), scvf);
return residual;
}
......@@ -321,7 +324,7 @@ private:
template<class SubSol>
void assembleResidualImpl_(Dune::index_constant<1>, SubSol& res)
{
for(auto&& scvf : scvfs(this->fvGeometry()))
for (auto&& scvf : scvfs(this->fvGeometry()))
res[scvf.dofIndex()] += this->asImp_().assembleFaceResidualImpl(scvf);
}
......@@ -342,6 +345,9 @@ private:
{
this->asImp_().assembleJacobianCellCenterCoupling(domainJ, jacRow[domainJ], residual, gridVariablesI);
});
// handle cells with a fixed Dirichlet value
incorporateDirichletCells_(jacRow);
}
//! Assembles the residuals and derivatives for the face dofs.
......@@ -363,6 +369,34 @@ private:
});
}
//! If specified in the problem, a fixed Dirichlet value can be assigned to cell centered unknows such as pressure
template<class JacobianMatrixRow>
void incorporateDirichletCells_(JacobianMatrixRow& jacRow)
{
const auto cellCenterGlobalI = problem().fvGridGeometry().elementMapper().index(this->element());
// overwrite the partial derivative with zero in case a fixed Dirichlet BC is used
static constexpr auto offset = numEq - numEqCellCenter;
for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
if (problem().isDirichletCell(this->element(), this->fvGeometry(), this->fvGeometry().scv(cellCenterGlobalI), eqIdx + offset))
{
using namespace Dune::Hybrid;
forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
{
auto& ccRowI = jacRow[i][cellCenterGlobalI];
for (auto col = ccRowI.begin(); col != ccRowI.end(); ++col)
{
ccRowI[col.index()][eqIdx] = 0.0;
// set the diagonal entry to 1.0
if ((i == domainId) && (col.index() == cellCenterGlobalI))
ccRowI[col.index()][eqIdx][eqIdx] = 1.0;
}
});
}
}
}
ElementFaceVariables curElemFaceVars_;
ElementFaceVariables prevElemFaceVars_;
CouplingManager& couplingManager_; //!< the coupling manager
......@@ -449,7 +483,6 @@ class SubDomainStaggeredLocalAssembler<id, TypeTag, Assembler, DiffMethod::numer
using FaceVariables = typename ElementFaceVariables::FaceVariables;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
......@@ -502,9 +535,6 @@ public:
const auto& fvGridGeometry = this->problem().fvGridGeometry();
const auto& curSol = this->curSol()[domainI];
// get the vecor of the acutal element residuals
// const auto origResiduals = this->evalLocalResidual();
const auto cellCenterGlobalI = fvGridGeometry.elementMapper().index(element);
const auto origResidual = this->evalLocalResidualForCellCenter();
......@@ -916,9 +946,7 @@ public:
template<class JacobianMatrixDiagBlock, class GridVariables>
void evalAdditionalDerivatives(const std::vector<std::size_t>& additionalDofDependencies,
JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
{
}
{ }
/*!
* \brief Updates the current global Jacobian matrix with the
......@@ -930,7 +958,7 @@ public:
const int globalI,
const int globalJ,
const int pvIdx,
const CCOrFacePrimaryVariables &partialDeriv)
const CCOrFacePrimaryVariables& partialDeriv)
{
for (int eqIdx = 0; eqIdx < partialDeriv.size(); eqIdx++)
{
......
......@@ -171,10 +171,27 @@ public:
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
bool isDirichletCell(const Element& element,
const typename FVGridGeometry::LocalView& fvGeometry,
const typename FVGridGeometry::SubControlVolume& scv,
int pvIdx) const
{
// set a fixed pressure in all cells at the boundary
values.setDirichletCell(Indices::pressureIdx);
for (const auto& scvf : scvfs(fvGeometry))
if (scvf.boundary())
return true;
return values;
return false;
}
/*!
......
......@@ -91,6 +91,7 @@ public:
using CellArray = std::array<unsigned int, dimWorld>;
const CellArray numCells = getParam<CellArray>("Grid.Cells");
cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0];
cellSizeY_ = this->fvGridGeometry().bBoxMax()[1] / numCells[1];
}
/*!
......@@ -143,13 +144,26 @@ public:
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
// set a fixed pressure in one cell
if (isLowerLeftCell_(globalPos))
values.setDirichletCell(Indices::pressureIdx);
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
template<class Element, class FVElementGeometry, class SubControlVolume>
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{
// set a fixed pressure in one cell
return (isLowerLeftCell_(scv.center()) && pvIdx == Indices::pressureIdx);
}
/*!
* \brief Return dirichlet boundary values at a given position
*
......@@ -189,12 +203,13 @@ private:
bool isLowerLeftCell_(const GlobalPosition& globalPos) const
{
return globalPos[0] < (0.5*cellSizeX_ + eps_) && globalPos[1] < eps_;
return globalPos[0] < (0.5*cellSizeX_ + eps_) && globalPos[1] < (0.5*cellSizeY_ + eps_);
}
Scalar eps_;
Scalar lidVelocity_;
Scalar cellSizeX_;
Scalar cellSizeY_;
};
} //end namespace
......
......@@ -174,11 +174,28 @@ public:
// set Dirichlet values for the velocity and pressure everywhere
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setDirichletCell(Indices::pressureIdx);
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
bool isDirichletCell(const Element& element,
const typename FVGridGeometry::LocalView& fvGeometry,
const typename FVGridGeometry::SubControlVolume& scv,
int pvIdx) const
{
bool onBoundary = false;
for (const auto& scvf : scvfs(fvGeometry))
onBoundary = std::max(onBoundary, scvf.boundary());
return onBoundary;
}
/*!
* \brief Return dirichlet boundary values at a given position
*
......
......@@ -100,7 +100,7 @@ public:
using CellArray = std::array<unsigned int, dimWorld>;
const auto numCells = getParam<CellArray>("Grid.Cells");
cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0];
cellSizeX_ = (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]) / numCells[0];
createAnalyticalSolution_();
}
......@@ -172,13 +172,27 @@ public:
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
// set a fixed pressure in one cell
if (isLowerLeftCell_(globalPos))
values.setDirichletCell(Indices::pressureIdx);
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
template<class Element, class FVElementGeometry, class SubControlVolume>
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{
// set fixed pressure in all cells at the left boundary
auto isAtLeftBoundary = [&](const auto& globalPos) { return globalPos[0] < (this->fvGridGeometry().bBoxMin()[0] + 0.6*cellSizeX_ + eps_); };
return (isAtLeftBoundary(scv.center()) && pvIdx == Indices::pressureIdx);
}
/*!
* \brief Return dirichlet boundary values at a given position
*
......@@ -293,11 +307,6 @@ private:
}
}
bool isLowerLeftCell_(const GlobalPosition& globalPos) const
{
return globalPos[0] < (this->fvGridGeometry().bBoxMin()[0] + 0.5*cellSizeX_ + eps_);
}
Scalar eps_;
Scalar cellSizeX_;
......
......@@ -26,6 +26,7 @@
#ifndef DUMUX_DONEA_TEST_PROBLEM_HH
#define DUMUX_DONEA_TEST_PROBLEM_HH
#include <dune/common/fmatrix.hh>
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/components/constant.hh>
......@@ -90,9 +91,10 @@ class NavierStokesAnalyticProblem : public NavierStokesProblem<TypeTag>
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
static constexpr auto dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
using DimMatrix = Dune::FieldVector<Dune::FieldVector<Scalar, dimWorld>, dimWorld>;
using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using DimVector = GlobalPosition;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
public:
NavierStokesAnalyticProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
......@@ -199,13 +201,32 @@ public:
{
BoundaryTypes values;
// set Dirichlet values for the velocity and pressure everywhere
values.setDirichletCell(Indices::conti0EqIdx);
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::momentumXBalanceIdx);
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
bool isDirichletCell(const Element& element,
const typename FVGridGeometry::LocalView& fvGeometry,
const typename FVGridGeometry::SubControlVolume& scv,
int pvIdx) const
{
// set a fixed pressure in all cells at the boundary
for (const auto& scvf : scvfs(fvGeometry))
if (scvf.boundary())
return true;
return false;
}
/*!
* \brief Return dirichlet boundary values at a given position
*
......
......@@ -103,10 +103,6 @@ public:
useWholeLength_ = getParam<bool>("Problem.UseWholeLength");
FluidSystem::init();
deltaRho_.resize(this->fvGridGeometry().numCellCenterDofs());
using CellArray = std::array<unsigned int, dimWorld>;
const CellArray numCells = getParam<CellArray>("Grid.Cells");
cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0];
}
/*!
......@@ -159,9 +155,6 @@ public:
values.setNeumann(Indices::conti0EqIdx);
values.setNeumann(transportEqIdx);
if(isLowerLeftCell_(globalPos))
values.setDirichletCell(Indices::pressureIdx);
if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_)
{
if(useWholeLength_)
......@@ -174,6 +167,23 @@ public:
return values;
}
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
*/
template<class FVElementGeometry, class SubControlVolume>
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{