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

[FreeFlow] Revise boundary conditions and indices

* index of the pressure has to equal the one of the main component
* advection law does not need outflow information
* fix some reference solutions
parent 1211f9f3
......@@ -165,26 +165,10 @@ public:
if (!scvf.boundary())
continue;
const auto bcTypes = problem.boundaryTypes(element, scvf);
CellCenterPrimaryVariables boundaryPriVars(0.0);
for(int eqIdx = 0; eqIdx < CellCenterPrimaryVariables::dimension; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
boundaryPriVars[eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[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);
auto boundaryPriVars = GVV::Traits::getBoundaryPriVars(problem, sol, element, scvf);
auto elemSol = elementSolution<FVElementGeometry>(std::move(boundaryPriVars));
volumeVariables_[localIdx].update(elemSol,
problem,
......
......@@ -49,9 +49,11 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::staggered >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
......@@ -81,6 +83,7 @@ public:
using Cache = FluxVariablesCaching::EmptyDiffusionCache;
static CellCenterPrimaryVariables flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
......@@ -102,7 +105,7 @@ public:
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
const auto bcTypes = problem.boundaryTypes(element, scvf);
if(bcTypes.isOutflow(eqIdx))
return flux;
}
......
......@@ -44,6 +44,7 @@ template <class TypeTag>
class FouriersLawImplementation<TypeTag, DiscretizationMethod::staggered >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
......@@ -63,13 +64,17 @@ public:
using Cache = FluxVariablesCaching::EmptyDiffusionCache;
//! calculate the diffusive energy fluxes
static Scalar flux(const Element& element,
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
{
Scalar flux(0.0);
if(scvf.boundary() && problem.boundaryTypes(element, scvf).isOutflow(Indices::energyBalanceIdx))
return flux;
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
......
......@@ -38,9 +38,75 @@ struct StaggeredGridDefaultGridVolumeVariablesTraits
{
using Problem = P;
using VolumeVariables = VV;
using PrimaryVariables = typename VV::PrimaryVariables;
template<class GridVolumeVariables, bool cachingEnabled>
using LocalView = StaggeredElementVolumeVariables<GridVolumeVariables, cachingEnabled>;
//! Returns the primary variales used for the boundary volVars and checks for admissable
//! combinations for boundary conditions.
template<class Problem, class SolutionVector, class Element, class SubControlVolumeFace>
static PrimaryVariables getBoundaryPriVars(const Problem& problem,
const SolutionVector& sol,
const Element& element,
const SubControlVolumeFace& scvf)
{
using CellCenterPrimaryVariables = typename SolutionVector::value_type;
using Indices = typename VolumeVariables::Indices;
static constexpr auto dim = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
static constexpr auto offset = dim;
const auto bcTypes = problem.boundaryTypes(element, scvf);
PrimaryVariables boundaryPriVars(0.0);
// make sure to not use outflow BC for momentum balance
for(int i = 0; i < dim; ++i)
{
if(bcTypes.isOutflow(Indices::velocity(i)))
DUNE_THROW(Dune::InvalidStateException, "Outflow condition cannot be used for velocity. Set only a Dirichlet value for pressure instead.");
}
if(bcTypes.isOutflow(Indices::pressureIdx))
DUNE_THROW(Dune::InvalidStateException, "Outflow condition cannot be used for pressure. Set only a Dirichlet value for velocity instead.");
// Determine the pressure value at a boundary with a Dirichlet condition for velocity.
// This just takes the value of the adjacent inner cell.
if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
{
if(bcTypes.isDirichlet(Indices::pressureIdx))
DUNE_THROW(Dune::InvalidStateException, "A Dirichlet condition for velocity must not be combined with a Dirichlet condition for pressure");
else
boundaryPriVars[Indices::pressureIdx] = sol[scvf.insideScvIdx()][Indices::pressureIdx - offset];
// TODO: pressure could be extrapolated to the boundary
}
// Determine the pressure value for a boundary with a Dirichlet condition for pressure.
// Takes a value specified in the problem.
if(bcTypes.isDirichlet(Indices::pressureIdx))
{
if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
DUNE_THROW(Dune::InvalidStateException, "A Dirichlet condition for velocity must not be combined with a Dirichlet condition for pressure");
else
boundaryPriVars[Indices::pressureIdx] = problem.dirichlet(element, scvf)[Indices::pressureIdx];
}
// Return for isothermal single-phase systems ...
if(CellCenterPrimaryVariables::dimension == 1)
return boundaryPriVars;
// ... or handle values for components, temperature, etc.
for(int eqIdx = offset; eqIdx < PrimaryVariables::dimension; ++eqIdx)
{
if(eqIdx == Indices::pressureIdx)
continue;
if(bcTypes.isDirichlet(eqIdx))
boundaryPriVars[eqIdx] = problem.dirichlet(element, scvf)[eqIdx];
else if(bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry() || bcTypes.isNeumann(eqIdx))
boundaryPriVars[eqIdx] = sol[scvf.insideScvIdx()][eqIdx - offset];
}
return boundaryPriVars;
}
};
/*!
......@@ -81,7 +147,6 @@ public:
template<class FVGridGeometry, class SolutionVector>
void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
{
using CellCenterPrimaryVariables = typename SolutionVector::value_type;
auto numScv = fvGridGeometry.numScv();
auto numBoundaryScvf = fvGridGeometry.numBoundaryScvf();
......@@ -108,25 +173,9 @@ public:
if (!scvf.boundary())
continue;
const auto bcTypes = problem().boundaryTypes(element, scvf);
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
PrimaryVariables boundaryPriVars(0.0);
constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension;
for(int eqIdx = offset; eqIdx < PrimaryVariables::dimension; ++eqIdx)
{
if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
boundaryPriVars[eqIdx] = problem().dirichlet(element, scvf)[eqIdx];
else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
boundaryPriVars[eqIdx] = sol[scvf.insideScvIdx()][eqIdx - offset];
//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.");
}
auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(boundaryPriVars));
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
auto boundaryPriVars = Traits::getBoundaryPriVars(problem(), sol, element, scvf);
const auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(boundaryPriVars));
volumeVariables_[scvf.outsideScvIdx()].update(elemSol, problem(), element, insideScv);
}
}
......
......@@ -47,10 +47,12 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::staggered >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
......@@ -85,6 +87,7 @@ public:
using CacheFiller = FluxVariablesCaching::EmptyCacheFiller;
static CellCenterPrimaryVariables flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
......@@ -119,7 +122,7 @@ public:
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
const auto bcTypes = problem.boundaryTypes(element, scvf);
if(bcTypes.isOutflow(eqIdx) && eqIdx != pressureIdx)
return componentFlux;
}
......@@ -183,7 +186,7 @@ public:
for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
{
componentFlux[compIdx] = reducedFlux[compIdx];
componentFlux[numComponents-1] -=reducedFlux[compIdx];
componentFlux[numComponents-1] -= reducedFlux[compIdx];
}
return componentFlux ;
......
......@@ -44,6 +44,9 @@ public:
//! The index of the main component
static constexpr int mainCompIdx = fluidSystemPhaseIdx;
//! Index of the pressure has to equal the one of the main component
static constexpr int pressureIdx = FreeflowIndices::conti0EqIdx + mainCompIdx;
//! The index of the component whose mass balance will be replaced by the total one
static constexpr int replaceCompEqIdx = theReplaceCompEqIdx;
};
......
......@@ -64,28 +64,17 @@ public:
*/
template<class ElementVolumeVariables, class ElementFaceVariables, class FluxVariablesCache>
CellCenterPrimaryVariables computeMassFlux(const Problem& problem,
const Element &element,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache)
{
CellCenterPrimaryVariables flux(0.0);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// get equation index
const auto eqIdx = Indices::conti0EqIdx + compIdx;
bool isOutflow = false;
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypes(element, scvf);
if(bcTypes.isOutflow(eqIdx))
isOutflow = true;
}
auto upwindTerm = [compIdx](const auto& volVars)
{
const auto density = useMoles ? volVars.molarDensity() : volVars.density();
......@@ -93,10 +82,10 @@ public:
return density * fraction;
};
flux[compIdx] = ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
flux[compIdx] = ParentType::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm);
}
flux += MolecularDiffusionType::flux(problem, fvGeometry, elemVolVars, scvf);
flux += MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
// in case one balance is substituted by the total mass balance
if (Indices::replaceCompEqIdx < numComponents)
......
......@@ -120,11 +120,9 @@ public:
if(compIdx == Indices::mainCompIdx)
continue;
int offset = Indices::mainCompIdx != 0 ? 1 : 0;
// temporary add 1.0 to remove spurious differences in mole fractions
// which are below the numerical accuracy
Scalar moleOrMassFraction = elemSol[0][Indices::conti0EqIdx+compIdx+offset] + 1.0;
Scalar moleOrMassFraction = elemSol[0][Indices::conti0EqIdx+compIdx] + 1.0;
moleOrMassFraction = moleOrMassFraction - 1.0;
if(useMoles)
fluidState.setMoleFraction(fluidSystemPhaseIdx, compIdx, moleOrMassFraction);
......
......@@ -91,21 +91,19 @@ public:
* \param elemFaceVars The face variables
* \param scvf The sub control volume face
* \param upwindTerm The uwind term (i.e. the advectively transported quantity)
* \param isOutflow Determines if we are on an outflow boundary
*/
template<class UpwindTerm>
static Scalar advectiveFluxForCellCenter(const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
UpwindTerm upwindTerm,
bool isOutflow = false)
UpwindTerm upwindTerm)
{
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const bool insideIsUpstream = scvf.directionSign() == sign(velocity);
static const Scalar upWindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;
......@@ -143,21 +141,15 @@ public:
// The advectively transported quantity (i.e density for a single-phase system).
auto upwindTerm = [](const auto& volVars) { return volVars.density(); };
// Check if we are on an outflow boundary.
const bool isOutflow = scvf.boundary()
? problem.boundaryTypes(element, scvf).isOutflow(Indices::conti0EqIdx)
: false;
// Call the generic flux function.
CellCenterPrimaryVariables result(0.0);
result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm);
return result;
}
/*!
* \brief Returns the momentum flux over all staggered faces.
*
*/
FacePrimaryVariables computeMomentumFlux(const Problem& problem,
const Element& element,
......@@ -252,16 +244,11 @@ public:
// (pointing in opposite direction of the scvf's one).
frontalFlux += pressure * -1.0 * scvf.directionSign();
// Treat outflow conditions.
// Handle inflow or outflow conditions.
// Treat the staggered half-volume adjacent to the boundary as if it was on the opposite side of the boundary.
// The respective face's outer normal vector will point in the same direction as the scvf's one.
if(scvf.boundary())
{
if(problem.boundaryTypes(element, scvf).isOutflow(Indices::velocity(scvf.directionIndex())))
{
// Treat the staggered half-volume adjacent to the boundary as if it was on the opposite side of the boundary.
// The respective face's outer normal vector will point in the same direction as the scvf's one.
frontalFlux += outflowBoundaryFlux_(problem, element, scvf, elemVolVars, elemFaceVars);
}
}
frontalFlux += inflowOutflowBoundaryFlux_(problem, element, scvf, elemVolVars, elemFaceVars);
// Account for the staggered face's area. For rectangular elements, this equals the area of the scvf
// our velocity dof of interest lives on.
......@@ -468,36 +455,27 @@ private:
if (!enableUnsymmetrizedVelocityGradient)
{
// For the normal gradient, get the velocities perpendicular to the velocity at the current scvf.
// The inner one is located at staggered face within the own element,
// the outer one at the respective staggered face of the element on the other side of the
// current scvf.
const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
// Lambda to conveniently get the outer normal velocity for faces that are on the boundary
// and therefore have no neighbor. Calls the problem to retrieve a fixed value set on the boundary.
auto getNormalVelocityFromBoundary = [&]()
// If we are at a boundary, a gradient of zero is implictly assumed for all velocities,
// thus no further calculations are required.
if (!scvf.boundary())
{
if (problem.boundaryTypesAtPos(scvf.center()).isOutflow(Indices::velocity(normalFace.directionIndex())))
return faceVars.velocityNormalInside(localSubFaceIdx);
const auto ghostFace = makeNormalGhostFace_(scvf, localSubFaceIdx);
return problem.dirichlet(element, ghostFace)[Indices::velocity(normalFace.directionIndex())];
};
const Scalar outerNormalVelocity = scvf.hasFrontalNeighbor(localSubFaceIdx)
? faceVars.velocityNormalOutside(localSubFaceIdx)
: getNormalVelocityFromBoundary();
// Calculate the velocity gradient in positive coordinate direction.
const Scalar normalDeltaV = scvf.normalInPosCoordDir()
? (outerNormalVelocity - innerNormalVelocity)
: (innerNormalVelocity - outerNormalVelocity);
const Scalar normalGradient = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
// Account for the orientation of the staggered normal face's outer normal vector.
normalDiffusiveFlux -= muAvg * normalGradient * normalFace.directionSign();
// For the normal gradient, get the velocities perpendicular to the velocity at the current scvf.
// The inner one is located at staggered face within the own element,
// the outer one at the respective staggered face of the element on the other side of the
// current scvf.
const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
const Scalar outerNormalVelocity = faceVars.velocityNormalOutside(localSubFaceIdx);
// Calculate the velocity gradient in positive coordinate direction.
const Scalar normalDeltaV = scvf.normalInPosCoordDir()
? (outerNormalVelocity - innerNormalVelocity)
: (innerNormalVelocity - outerNormalVelocity);
const Scalar normalGradient = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
// Account for the orientation of the staggered normal face's outer normal vector.
normalDiffusiveFlux -= muAvg * normalGradient * normalFace.directionSign();
}
}
// For the parallel derivative, get the velocities at the current (own) scvf
......@@ -528,7 +506,7 @@ private:
}
/*!
* \brief Returns the momentum flux over an outflow boundary face.
* \brief Returns the momentum flux over an inflow or outflow boundary face.
*
* \verbatim
* scvf //
......@@ -543,15 +521,15 @@ private:
* // boundary
* \endverbatim
*/
FacePrimaryVariables outflowBoundaryFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars)
FacePrimaryVariables inflowOutflowBoundaryFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars)
{
FacePrimaryVariables outflow(0.0);
FacePrimaryVariables inOrOutflow(0.0);
// Advective momentum outflow.
// Advective momentum flux.
if(enableInertiaTerms)
{
const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
......@@ -560,7 +538,7 @@ private:
const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
insideVolVars : outsideVolVars;
outflow += velocitySelf * velocitySelf * upVolVars.density();
inOrOutflow += velocitySelf * velocitySelf * upVolVars.density();
}
// Apply a pressure at the boundary.
......@@ -568,10 +546,10 @@ private:
? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
problem.initial(scvf)[Indices::pressureIdx])
: problem.dirichlet(element, scvf)[Indices::pressureIdx];
outflow += boundaryPressure;
inOrOutflow += boundaryPressure;
// Account for the orientation of the face at the boundary,
return outflow * scvf.directionSign();
return inOrOutflow * scvf.directionSign();
}
private:
......@@ -582,12 +560,6 @@ private:
return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()}, ownScvf.dofIndex(), ownScvf.index());
};
//! helper function to conveniently create a ghost face which is outside the domain, normal to the scvf of interest
SubControlVolumeFace makeNormalGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
{
return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualOuterNormalFaceDofPos);
};
//! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
{
......@@ -601,7 +573,6 @@ private:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
}
};
} // end namespace
......
......@@ -195,10 +195,10 @@ public:
const BoundaryTypes& bcTypes) const
{
// set a fixed pressure for cells adjacent to a wall
if(bcTypes.isDirichletCell(Indices::conti0EqIdx))
if(bcTypes.isDirichletCell(Indices::pressureIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
residual[Indices::conti0EqIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichlet(element, insideScv)[Indices::pressureIdx];
residual[Indices::pressureIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichlet(element, insideScv)[Indices::pressureIdx];
}
}
......@@ -273,18 +273,16 @@ protected:
// handle the actual boundary conditions:
const auto bcTypes = problem.boundaryTypes(element, scvf);
// set a fixed value for the velocity for Dirichlet boundary conditions
if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
{
// set a fixed value for the velocity for Dirichlet boundary conditions
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
residual = velocity - dirichletValue;
}
// handle Neumann boundary conditions
if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
else if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
{
// set the given Neumann flux for the face on the boundary itself
// set 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();
......@@ -293,24 +291,21 @@ protected:
FluxVariables fluxVars;
residual += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
}
// For symmetry boundary conditions, there is no flow accross the boundary and
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
if(bcTypes.isSymmetry())
else if(bcTypes.isSymmetry())
{
// For symmetry boundary conditions, there is no flow accross the boundary and
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar fixedValue = 0.0;
residual = velocity - fixedValue;
}
// outflow condition for the momentum balance equation
if(bcTypes.isOutflow(Indices::velocity(scvf.directionIndex())))
else if(bcTypes.isDirichlet(Indices::pressureIdx))
{
if(bcTypes.isDirichlet(Indices::conti0EqIdx))
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
else
DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center() << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
// If none of the above conditions apply, we are at an "fixed pressure" boundary for which the velocity needs to be assembled
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
}
else
DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations.");
}
}
......