Skip to content

GitLab

Commit fd7d7903 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[staggered][localresidual] Eval boundaries for cell center directly in flux calculation

parent dcf4df19
Hide whitespace changes
Inline Side-by-side
dumux/assembly/staggeredlocalresidual.hh
View file @ fd7d7903
......@@ -104,6 +104,8 @@ public:
{
if (!scvf.boundary())
residual += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
else
residual += asImp_().computeBoundaryFluxForCellCenter(problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
}
//! Evaluate the source terms for a cell center residual
......@@ -171,19 +173,6 @@ public:
residual += storage;
}
//! Evaluate the boundary conditions for a cell center residual
void evalBoundaryForCellCenter(CellCenterResidualValue& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
asImp_().evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, elemVolVars, elemFaceVars, bcTypes, elemFluxVarsCache);
}
//! for compatibility with FVLocalAssemblerBase
template<class... Args>
CellCenterResidualValue evalFluxAndSource(Args&&... args) const
......
dumux/freeflow/navierstokes/staggered/localresidual.hh
View file @ fd7d7903
......@@ -185,70 +185,63 @@ public:
* \brief Evaluate boundary conditions for a cell center dof
*/
template<class ElementBoundaryTypes>
void evalBoundaryForCellCenter_(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
CellCenterResidual computeBoundaryFluxForCellCenter(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
const auto bcTypes = problem.boundaryTypes(element, scvf);
CellCenterResidual result(0.0);
// no fluxes occur over symmetry boundaries
if (!bcTypes.isSymmetry())
{
const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
if (scvf.boundary())
{
const auto bcTypes = problem.boundaryTypes(element, scvf);
// treat Dirichlet and outflow BCs
FluxVariables fluxVars;
auto boundaryFlux = fluxVars.computeMassFlux(problem, element, fvGeometry, elemVolVars,
elemFaceVars, scvf, elemFluxVarsCache[scvf]);
// no fluxes occur over symmetry boundaries
if (!bcTypes.isSymmetry())
{
const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
// treat Dirichlet and outflow BCs
result = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
// treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
if(bcTypes.hasNeumann())
// treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
if(bcTypes.hasNeumann())
{
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
{
if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
{
boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
* extrusionFactor * scvf.area();
}
result[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
* extrusionFactor * scvf.area();
}
}
// account for wall functions, if used
incorporateWallFunction_(boundaryFlux, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars);
// add the flux over the boundary scvf to the residual
residual += boundaryFlux;
}
// account for wall functions, if used
incorporateWallFunction_(result, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars);
}
}
return result;
}
/*!
* \brief Evaluate boundary conditions for a face dof
* \brief Evaluate Dirichlet (fixed value) boundary conditions for a face dof
*/
template<class ElementBoundaryTypes>
void evalBoundaryForFace(FaceResidual& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
void evalDirichletBoundariesForFace(FaceResidual& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
if (scvf.boundary())
{
......@@ -262,17 +255,6 @@ public:
const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
residual = velocity - dirichletValue;
}
else if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
{
// the source term has already been accounted for, here we
// add a given Neumann flux for the face on the boundary itself ...
const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
residual += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
* extrusionFactor * scvf.area();
// ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
}
else if(bcTypes.isSymmetry())
{
// for symmetry boundary conditions, there is no flow accross the boundary and
......
dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
View file @ fd7d7903
......@@ -133,44 +133,6 @@ public:
return source;
}
protected:
/*!
* \brief Evaluate boundary conditions for a cell center dof
*/
template<class ElementBoundaryTypes>
void evalBoundaryForCellCenter_(CellCenterResidual& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
BaseLocalResidual::evalBoundaryForCellCenter_(residual, problem, element, fvGeometry,
elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
for (auto&& scvf : scvfs(fvGeometry))
{
unsigned int elementIdx = problem.fvGridGeometry().elementMapper().index(element);
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
// fixed value for the turbulent kinetic energy
if (insideVolVars.inNearWallRegion())
{
residual[Indices::turbulentKineticEnergyEqIdx - cellCenterOffset]
= insideVolVars.turbulentKineticEnergy() - problem.turbulentKineticEnergyWallFunction(elementIdx);
}
// fixed value for the dissipation
if (insideVolVars.inNearWallRegion() || insideVolVars.isMatchingPoint())
{
residual[Indices::dissipationEqIdx - cellCenterOffset]
= insideVolVars.dissipation() - problem.dissipationWallFunction(elementIdx);
}
}
}
};
} // end namespace Dumux
......
dumux/multidomain/subdomainstaggeredlocalassembler.hh
View file @ fd7d7903
......@@ -166,8 +166,6 @@ public:
if (!this->assembler().isStationaryProblem())
residual += evalLocalStorageResidualForCellCenter();
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);
......
test/freeflow/rans/problem.hh
View file @ fd7d7903
......@@ -297,6 +297,32 @@ public:
if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx)
return true;
return false;
}
#elif KEPSILON
/*!
* \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
{
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
// set a fixed turbulent kinetic energy and dissipation near the wall
if (this->inNearWallRegion(eIdx))
return pvIdx == Indices::turbulentKineticEnergyEqIdx || pvIdx == Indices::dissipationEqIdx;
// set a fixed dissipation at the matching point
if (this->isMatchingPoint(eIdx))
return pvIdx == Indices::dissipationEqIdx;
return false;
}
#endif
......@@ -338,10 +364,18 @@ public:
PrimaryVariables values(initialAtPos(globalPos));
#if KOMEGA
using std::pow;
unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element);
const auto wallDistance = ParentType::wallDistance_[elementIdx];
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
const auto wallDistance = ParentType::wallDistance_[eIdx];
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[eIdx]
/ (ParentType::betaOmega() * pow(wallDistance, 2));
#elif KEPSILON
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
// fixed value for the turbulent kinetic energy
values[Indices::turbulentKineticEnergyEqIdx] = this->turbulentKineticEnergyWallFunction(eIdx);
// fixed value for the dissipation
values[Indices::dissipationEqIdx] = this->dissipationWallFunction(eIdx);
#endif
return values;
}
......
test/freeflow/ransnc/problem.hh
View file @ fd7d7903
......@@ -323,6 +323,32 @@ public:
if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx)
return true;
return false;
}
#elif KEPSILON
/*!
* \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
{
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
// set a fixed turbulent kinetic energy and dissipation near the wall
if (this->inNearWallRegion(eIdx))
return pvIdx == Indices::turbulentKineticEnergyEqIdx || pvIdx == Indices::dissipationEqIdx;
// set a fixed dissipation at the matching point
if (this->isMatchingPoint(eIdx))
return pvIdx == Indices::dissipationEqIdx;
return false;
}
#endif
......@@ -371,10 +397,18 @@ public:
PrimaryVariables values(initialAtPos(globalPos));
#if KOMEGA
using std::pow;
unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element);
const auto wallDistance = ParentType::wallDistance_[elementIdx];
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
unsigned int eIdx = this->fvGridGeometry().elementMapper().index(element);
const auto wallDistance = ParentType::wallDistance_[eIdx];
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[eIdx]
/ (ParentType::betaOmega() * pow(wallDistance, 2));
#elif KEPSILON
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
// fixed value for the turbulent kinetic energy
values[Indices::turbulentKineticEnergyEqIdx] = this->turbulentKineticEnergyWallFunction(eIdx);
// fixed value for the dissipation
values[Indices::dissipationEqIdx] = this->dissipationWallFunction(eIdx);
#endif
return values;
}
......
Supports Markdown
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