24#ifndef DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
25#define DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
46template<
class TypeTag,
int upwindSchemeOrder>
51 using GridVolumeVariables =
typename GridVariables::GridVolumeVariables;
52 using ElementVolumeVariables =
typename GridVolumeVariables::LocalView;
53 using VolumeVariables =
typename GridVolumeVariables::VolumeVariables;
55 using GridFluxVariablesCache =
typename GridVariables::GridFluxVariablesCache;
57 using FluxVariablesCache =
typename GridFluxVariablesCache::FluxVariablesCache;
59 using GridFaceVariables =
typename GridVariables::GridFaceVariables;
60 using ElementFaceVariables =
typename GridFaceVariables::LocalView;
61 using FaceVariables =
typename GridFaceVariables::FaceVariables;
65 using FVElementGeometry =
typename GridGeometry::LocalView;
66 using GridView =
typename GridGeometry::GridView;
67 using Element =
typename GridView::template Codim<0>::Entity;
70 using Indices =
typename ModelTraits::Indices;
71 using SubControlVolume =
typename FVElementGeometry::SubControlVolume;
72 using SubControlVolumeFace =
typename FVElementGeometry::SubControlVolumeFace;
78 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
79 static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
81 static_assert(upwindSchemeOrder <= 2,
"Not implemented: Order higher than 2!");
82 static_assert(upwindSchemeOrder <= 1 || GridFluxVariablesCache::cachingEnabled,
83 "Higher order upwind method requires caching enabled for the GridFluxVariablesCache!");
84 static_assert(upwindSchemeOrder <= 1 || GridGeometry::cachingEnabled,
85 "Higher order upwind method requires caching enabled for the GridGeometry!");
86 static_assert(upwindSchemeOrder <= 1 || GridFaceVariables::cachingEnabled,
87 "Higher order upwind method requires caching enabled for the GridFaceVariables!");
88 static_assert(upwindSchemeOrder <= 1 || GridVolumeVariables::cachingEnabled,
89 "Higher order upwind method requires caching enabled for the GridGeometry!");
93 const FVElementGeometry& fvGeometry,
94 const SubControlVolumeFace& scvf,
95 const ElementFaceVariables& elemFaceVars,
96 const ElementVolumeVariables& elemVolVars,
97 const UpwindScheme& upwindScheme)
99 , fvGeometry_(fvGeometry)
101 , elemFaceVars_(elemFaceVars)
102 , faceVars_(elemFaceVars[scvf])
103 , elemVolVars_(elemVolVars)
104 , upwindScheme_(upwindScheme)
117 const auto density = elemVolVars_[scvf_.insideScvIdx()].density();
120 if constexpr (useHigherOrder)
123 if (canDoFrontalSecondOrder_(selfIsUpstream))
125 const auto distances = getFrontalDistances_(selfIsUpstream);
126 const auto upwindMomenta = getFrontalSecondOrderUpwindMomenta_(
density, selfIsUpstream);
127 return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach());
132 const auto upwindMomenta = getFrontalFirstOrderUpwindMomenta_(
density, selfIsUpstream);
133 return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]);
146 const SubControlVolumeFace& lateralFace,
147 const int localSubFaceIdx,
148 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
149 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
const
154 const Scalar insideDensity = elemVolVars_[lateralFace.insideScvIdx()].density();
155 const Scalar outsideDensity = elemVolVars_[lateralFace.outsideScvIdx()].density();
158 if constexpr (useHigherOrder)
161 if (canDoLateralSecondOrder_(selfIsUpstream, localSubFaceIdx))
163 const auto upwindMomenta = getLateralSecondOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes);
164 const auto distances = getLateralDistances_(localSubFaceIdx, selfIsUpstream);
165 return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach());
170 const auto upwindMomenta = getLateralFirstOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes);
171 return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]);
183 bool canDoFrontalSecondOrder_(
bool selfIsUpstream)
const
186 return selfIsUpstream ? scvf_.hasForwardNeighbor(0) : scvf_.hasBackwardNeighbor(0);
192 std::array<Scalar, 3> getFrontalSecondOrderUpwindMomenta_(
const Scalar
density,
bool selfIsUpstream)
const
194 static_assert(useHigherOrder,
"Should only be reached if higher order methods are enabled");
195 const Scalar momentumSelf = faceVars_.velocitySelf() *
density;
196 const Scalar momentumOpposite = faceVars_.velocityOpposite() *
density;
198 return {momentumOpposite, momentumSelf, faceVars_.velocityForward(0)*
density};
200 return {momentumSelf, momentumOpposite, faceVars_.velocityBackward(0)*
density};
206 std::array<Scalar, 2> getFrontalFirstOrderUpwindMomenta_(
const Scalar
density,
bool selfIsUpstream)
const
208 const Scalar momentumSelf = faceVars_.velocitySelf() *
density;
209 const Scalar momentumOpposite = faceVars_.velocityOpposite() *
density;
211 return {momentumOpposite, momentumSelf};
213 return {momentumSelf, momentumOpposite};
221 std::array<Scalar, 3> getFrontalDistances_(
const bool selfIsUpstream)
const
223 static_assert(useHigherOrder,
"Should only be reached if higher order methods are enabled");
226 std::array<Scalar, 3> distances;
227 distances[0] = scvf_.selfToOppositeDistance();
228 distances[1] = scvf_.axisData().inAxisForwardDistances[0];
229 distances[2] = 0.5 * (scvf_.axisData().selfToOppositeDistance + scvf_.axisData().inAxisBackwardDistances[0]);
234 std::array<Scalar, 3> distances;
235 distances[0] = scvf_.selfToOppositeDistance();
236 distances[1] = scvf_.axisData().inAxisBackwardDistances[0];
237 distances[2] = 0.5 * (scvf_.axisData().selfToOppositeDistance + scvf_.axisData().inAxisForwardDistances[0]);
248 bool canDoLateralSecondOrder_(
const bool selfIsUpstream,
const int localSubFaceIdx)
const
253 return selfIsUpstream || scvf_.hasParallelNeighbor(localSubFaceIdx, 0);
296 std::array<Scalar, 3> getLateralSecondOrderUpwindMomenta_(
const Scalar insideDensity,
297 const Scalar outsideDensity,
298 const bool selfIsUpstream,
299 const int localSubFaceIdx,
300 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
301 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
const
303 static_assert(useHigherOrder,
"Should only be reached if higher order methods are enabled");
307 if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0) || scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
309 if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
311 const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
312 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
314 return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
316 else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
318 const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
319 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
321 return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
327 std::array<Scalar, 3> momenta;
328 momenta[1] = faceVars_.velocitySelf()*insideDensity;
329 momenta[0] = faceVars_.velocityParallel(localSubFaceIdx, 0)*insideDensity;
332 const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
335 if (scvf_.hasParallelNeighbor(oppositeSubFaceIdx, 0))
336 momenta[2] = faceVars_.velocityParallel(oppositeSubFaceIdx, 0)*insideDensity;
338 momenta[2] = getParallelVelocityFromOppositeBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, oppositeSubFaceIdx)*insideDensity;
343 std::array<Scalar, 3> momenta;
344 momenta[0] = faceVars_.velocitySelf()*outsideDensity;
345 momenta[1] = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
348 if (scvf_.hasParallelNeighbor(localSubFaceIdx, 1))
349 momenta[2] = faceVars_.velocityParallel(localSubFaceIdx, 1)*outsideDensity;
352 const auto& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
353 const auto& elementParallel = fvGeometry_.gridGeometry().element(lateralFace.outsideScvIdx());
354 const auto& firstParallelScvf = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
355 const auto& problem = elemVolVars_.gridVolVars().problem();
356 const auto& boundaryTypes = problem.boundaryTypes(elementParallel, firstParallelScvf);
357 momenta[2] = getParallelVelocityFromOppositeBoundary_(elementParallel, firstParallelScvf, elemFaceVars_[firstParallelScvf], boundaryTypes, localSubFaceIdx)*outsideDensity;
366 std::array<Scalar, 2> getLateralFirstOrderUpwindMomenta_(
const Scalar insideDensity,
367 const Scalar outsideDensity,
368 const bool selfIsUpstream,
369 const int localSubFaceIdx,
370 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
371 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
const
375 if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
377 const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
378 const Scalar boundaryMomentum = boundaryVelocity*outsideDensity;
379 return {boundaryMomentum, boundaryMomentum};
381 else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
383 const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
384 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
385 return {boundaryMomentum, boundaryMomentum};
388 const Scalar momentumParallel = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
389 const Scalar momentumSelf = faceVars_.velocitySelf()*insideDensity;
391 return {momentumParallel, momentumSelf};
393 return {momentumSelf, momentumParallel};
400 std::array<Scalar, 3> getLateralDistances_(
const int localSubFaceIdx,
const bool selfIsUpstream)
const
402 static_assert(useHigherOrder,
"Should only be reached if higher order methods are enabled");
406 const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
408 std::array<Scalar, 3> distances;
409 distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0);
410 distances[1] = scvf_.parallelDofsDistance(oppositeSubFaceIdx, 0);
411 if (scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
412 distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0];
414 distances[2] = scvf_.area() / 2.0;
419 std::array<Scalar, 3> distances;
420 distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0);
421 distances[1] = scvf_.parallelDofsDistance(localSubFaceIdx, 1);
422 distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0];
431 Scalar getParallelVelocityFromBoundary_(
const Element& element,
432 const SubControlVolumeFace& scvf,
433 const FaceVariables& faceVars,
434 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
435 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
436 const int localSubFaceIdx)
const
440 const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
442 return faceVars.velocitySelf();
444 const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
445 if (lateralFaceHasBJS)
447 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
449 const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
450 if (lateralFaceHasDirichletVelocity)
461 const auto& lateralFace = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
462 const auto ghostFace =
makeStaggeredBoundaryFace(lateralFace, scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
463 const auto& problem = elemVolVars_.gridVolVars().problem();
464 return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
468 DUNE_THROW(Dune::InvalidStateException,
"Something went wrong with the boundary conditions for the momentum equations at global position "
469 << fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx).center());
475 Scalar getParallelVelocityFromOppositeBoundary_(
const Element& element,
476 const SubControlVolumeFace& scvf,
477 const FaceVariables& faceVars,
478 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
479 const int localOppositeSubFaceIdx)
const
482 const auto& lateralOppositeScvf = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localOppositeSubFaceIdx).localLateralFaceIdx);
483 GlobalPosition center = scvf.pairData(localOppositeSubFaceIdx).lateralStaggeredFaceCenter + lateralOppositeScvf.center();
498 const auto& problem = elemVolVars_.gridVolVars().problem();
500 const auto lateralOppositeFaceBoundaryTypes = problem.boundaryTypes(element, lateralOppositeBoundaryFace);
501 return getParallelVelocityFromBoundary_(element, scvf, faceVars, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes, localOppositeSubFaceIdx);
523 const SubControlVolumeFace& boundaryScvf_(
const int localSubFaceIdx)
const
525 if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
527 return fvGeometry_.scvf(scvf_.outsideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
529 else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
545 const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
546 const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
548 const auto& localLateralIdx = scvf_.pairData(localSubFaceIdx).localLateralFaceIdx;
549 const auto& localLateralOppositeIdx = (localLateralIdx % 2) ? (localLateralIdx - 1) : (localLateralIdx + 1);
551 return fvGeometry_.scvf(parallelFace.outsideScvIdx(), localLateralOppositeIdx);
555 DUNE_THROW(Dune::InvalidStateException,
"The function boundaryScvf_ should only be called when hasHalfParallelNeighbor or hasCornerParallelNeighbor.");
578 Element boundaryElement_(
const int localSubFaceIdx)
const
580 if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
582 return fvGeometry_.gridGeometry().element(scvf_.outsideScvIdx());
584 else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
586 const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
587 const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
589 return fvGeometry_.gridGeometry().element(parallelFace.outsideScvIdx());
593 DUNE_THROW(Dune::InvalidStateException,
"When entering boundaryElement_ scvf_ should have either hasHalfParallelNeighbor or hasCornerParallelNeighbor true. Not the case here.");
616 bool dirichletParallelNeighbor_(
const int localSubFaceIdx)
const
618 const auto& problem = elemVolVars_.gridVolVars().problem();
619 const Element& boundaryElement = boundaryElement_(localSubFaceIdx);
620 const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx);
622 return problem.boundaryTypes(boundaryElement, boundaryScvf).isDirichlet(Indices::velocity(scvf_.directionIndex()));
643 Scalar getParallelVelocityFromCorner_(
const int localSubFaceIdx)
const
645 const auto& problem = elemVolVars_.gridVolVars().problem();
646 const Element& boundaryElement = boundaryElement_(localSubFaceIdx);
647 const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx);
648 const auto ghostFace =
makeStaggeredBoundaryFace(boundaryScvf, scvf_.pairData(localSubFaceIdx).lateralStaggeredFaceCenter);
650 return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf_.directionIndex())];
653 const Element& element_;
654 const FVElementGeometry& fvGeometry_;
655 const SubControlVolumeFace& scvf_;
656 const ElementFaceVariables& elemFaceVars_;
657 const FaceVariables& faceVars_;
658 const ElementVolumeVariables& elemVolVars_;
659 const UpwindScheme& upwindScheme_;
Some exceptions thrown in DuMux
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
This file contains different higher order methods for approximating the velocity.
SubControlVolumeFace makeStaggeredBoundaryFace(const SubControlVolumeFace &scvf, const typename SubControlVolumeFace::GlobalPosition &newCenter)
Helper function to turn a given cell scvface into a fake boundary face.
Definition: discretization/staggered/freeflow/subcontrolvolumeface.hh:104
UpwindSchemeImpl< GridGeometry, GridGeometry::discMethod > UpwindScheme
The upwind scheme used for the advective fluxes. This depends on the chosen discretization method.
Definition: flux/upwindscheme.hh:42
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
typename Detail::template ProblemTraits< Problem, GridGeometry::discMethod >::BoundaryTypes BoundaryTypes
Definition: common/typetraits/problem.hh:46
The upwinding variables class for the Navier-Stokes model using the staggered grid discretization.
Definition: staggeredupwindhelper.hh:48
FacePrimaryVariables computeUpwindFrontalMomentum(const bool selfIsUpstream) const
Returns the momentum in the frontal directon.
Definition: staggeredupwindhelper.hh:115
FacePrimaryVariables computeUpwindLateralMomentum(const bool selfIsUpstream, const SubControlVolumeFace &lateralFace, const int localSubFaceIdx, const std::optional< BoundaryTypes > ¤tScvfBoundaryTypes, const std::optional< BoundaryTypes > &lateralFaceBoundaryTypes) const
Returns the momentum in the lateral directon.
Definition: staggeredupwindhelper.hh:145
StaggeredUpwindHelper(const Element &element, const FVElementGeometry &fvGeometry, const SubControlVolumeFace &scvf, const ElementFaceVariables &elemFaceVars, const ElementVolumeVariables &elemVolVars, const UpwindScheme &upwindScheme)
Definition: staggeredupwindhelper.hh:92
Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered g...
Definition: velocitygradients.hh:39
static Scalar beaversJosephVelocityAtLateralScvf(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const SubControlVolumeFace &scvf, const FaceVariables &faceVars, const std::optional< BoundaryTypes > ¤tScvfBoundaryTypes, const std::optional< BoundaryTypes > &lateralFaceBoundaryTypes, const std::size_t localSubFaceIdx)
Returns the Beavers-Jospeh slip velocity for a lateral scvf which lies on the boundary.
Definition: velocitygradients.hh:323
Declares all properties used in Dumux.
Type traits for problem classes.