3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
staggeredupwindhelper.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
25#define DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH
26
27#include <array>
28#include <optional>
29
30#include <dumux/common/math.hh>
35
38#include "velocitygradients.hh"
39
40namespace Dumux {
41
46template<class TypeTag, int upwindSchemeOrder>
48{
50
51 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
52 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
53 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
54
55 using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
56 using UpwindScheme = typename GridFluxVariablesCache::UpwindScheme;
57 using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
58
59 using GridFaceVariables = typename GridVariables::GridFaceVariables;
60 using ElementFaceVariables = typename GridFaceVariables::LocalView;
61 using FaceVariables = typename GridFaceVariables::FaceVariables;
62
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;
73 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
75 using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes;
77
78 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
79 static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
80
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!");
90
91public:
92 StaggeredUpwindHelper(const Element& element,
93 const FVElementGeometry& fvGeometry,
94 const SubControlVolumeFace& scvf,
95 const ElementFaceVariables& elemFaceVars,
96 const ElementVolumeVariables& elemVolVars,
97 const UpwindScheme& upwindScheme)
98 : element_(element)
99 , fvGeometry_(fvGeometry)
100 , scvf_(scvf)
101 , elemFaceVars_(elemFaceVars)
102 , faceVars_(elemFaceVars[scvf])
103 , elemVolVars_(elemVolVars)
104 , upwindScheme_(upwindScheme)
105 {}
106
115 FacePrimaryVariables computeUpwindFrontalMomentum(const bool selfIsUpstream) const
116 {
117 const auto density = elemVolVars_[scvf_.insideScvIdx()].density();
118
119 // for higher order schemes do higher order upwind reconstruction
120 if constexpr (useHigherOrder)
121 {
122 // only second order is implemented so far
123 if (canDoFrontalSecondOrder_(selfIsUpstream))
124 {
125 const auto distances = getFrontalDistances_(selfIsUpstream);
126 const auto upwindMomenta = getFrontalSecondOrderUpwindMomenta_(density, selfIsUpstream);
127 return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach());
128 }
129 }
130
131 // otherwise apply first order upwind scheme
132 const auto upwindMomenta = getFrontalFirstOrderUpwindMomenta_(density, selfIsUpstream);
133 return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]);
134 }
135
145 FacePrimaryVariables computeUpwindLateralMomentum(const bool selfIsUpstream,
146 const SubControlVolumeFace& lateralFace,
147 const int localSubFaceIdx,
148 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
149 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const
150 {
151 // Check whether the own or the neighboring element is upstream.
152 // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
153 // of interest is located.
154 const Scalar insideDensity = elemVolVars_[lateralFace.insideScvIdx()].density();
155 const Scalar outsideDensity = elemVolVars_[lateralFace.outsideScvIdx()].density();
156
157 // for higher order schemes do higher order upwind reconstruction
158 if constexpr (useHigherOrder)
159 {
160 // only second order is implemented so far
161 if (canDoLateralSecondOrder_(selfIsUpstream, localSubFaceIdx))
162 {
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());
166 }
167 }
168
169 // otherwise apply first order upwind scheme
170 const auto upwindMomenta = getLateralFirstOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes);
171 return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]);
172 }
173
174private:
183 bool canDoFrontalSecondOrder_(bool selfIsUpstream) const
184 {
185 // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
186 return selfIsUpstream ? scvf_.hasForwardNeighbor(0) : scvf_.hasBackwardNeighbor(0);
187 }
188
192 std::array<Scalar, 3> getFrontalSecondOrderUpwindMomenta_(const Scalar density, bool selfIsUpstream) const
193 {
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;
197 if (selfIsUpstream)
198 return {momentumOpposite, momentumSelf, faceVars_.velocityForward(0)*density};
199 else
200 return {momentumSelf, momentumOpposite, faceVars_.velocityBackward(0)*density};
201 }
202
206 std::array<Scalar, 2> getFrontalFirstOrderUpwindMomenta_(const Scalar density, bool selfIsUpstream) const
207 {
208 const Scalar momentumSelf = faceVars_.velocitySelf() * density;
209 const Scalar momentumOpposite = faceVars_.velocityOpposite() * density;
210 if (selfIsUpstream)
211 return {momentumOpposite, momentumSelf};
212 else
213 return {momentumSelf, momentumOpposite};
214 }
215
221 std::array<Scalar, 3> getFrontalDistances_(const bool selfIsUpstream) const
222 {
223 static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled");
224 if (selfIsUpstream)
225 {
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]);
230 return distances;
231 }
232 else
233 {
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]);
238 return distances;
239 }
240 }
241
248 bool canDoLateralSecondOrder_(const bool selfIsUpstream, const int localSubFaceIdx) const
249 {
250 // If the self velocity is upstream, the downstream velocity can be assigned or retrieved
251 // from the boundary, even if there is no parallel neighbor.
252 // If the self velocity is downstream and there is no parallel neighbor I cannot use a second order approximation.
253 return selfIsUpstream || scvf_.hasParallelNeighbor(localSubFaceIdx, 0);
254 }
255
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
302 {
303 static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled");
304
305 // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
306 // thus we always use this value for the computation of the transported momentum.
307 if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0) || scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
308 {
309 if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
310 {
311 const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
312 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
313
314 return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
315 }
316 else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
317 {
318 const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
319 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
320
321 return {boundaryMomentum, boundaryMomentum, boundaryMomentum};
322 }
323 }
324
325 if (selfIsUpstream)
326 {
327 std::array<Scalar, 3> momenta;
328 momenta[1] = faceVars_.velocitySelf()*insideDensity;
329 momenta[0] = faceVars_.velocityParallel(localSubFaceIdx, 0)*insideDensity;
330
331 // The local index of the faces that is opposite to localSubFaceIdx
332 const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
333
334 // The "upstream-upstream" velocity is retrieved from the other parallel neighbor or from the boundary.
335 if (scvf_.hasParallelNeighbor(oppositeSubFaceIdx, 0))
336 momenta[2] = faceVars_.velocityParallel(oppositeSubFaceIdx, 0)*insideDensity;
337 else
338 momenta[2] = getParallelVelocityFromOppositeBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, oppositeSubFaceIdx)*insideDensity;
339 return momenta;
340 }
341 else
342 {
343 std::array<Scalar, 3> momenta;
344 momenta[0] = faceVars_.velocitySelf()*outsideDensity;
345 momenta[1] = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
346
347 // If there is another parallel neighbor I can assign the "upstream-upstream" velocity, otherwise I retrieve it from the boundary.
348 if (scvf_.hasParallelNeighbor(localSubFaceIdx, 1))
349 momenta[2] = faceVars_.velocityParallel(localSubFaceIdx, 1)*outsideDensity;
350 else
351 {
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;
358 }
359 return momenta;
360 }
361 }
362
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
372 {
373 // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
374 // thus we always use this value for the computation of the transported momentum.
375 if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx))
376 {
377 const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx);
378 const Scalar boundaryMomentum = boundaryVelocity*outsideDensity;
379 return {boundaryMomentum, boundaryMomentum};
380 }
381 else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0))
382 {
383 const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
384 const Scalar boundaryMomentum = boundaryVelocity*insideDensity;
385 return {boundaryMomentum, boundaryMomentum};
386 }
387
388 const Scalar momentumParallel = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity;
389 const Scalar momentumSelf = faceVars_.velocitySelf()*insideDensity;
390 if (selfIsUpstream)
391 return {momentumParallel, momentumSelf};
392 else
393 return {momentumSelf, momentumParallel};
394 }
395
400 std::array<Scalar, 3> getLateralDistances_(const int localSubFaceIdx, const bool selfIsUpstream) const
401 {
402 static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled");
403 if (selfIsUpstream)
404 {
405 // The local index of the faces that is opposite to localSubFaceIdx
406 const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
407
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];
413 else
414 distances[2] = scvf_.area() / 2.0;
415 return distances;
416 }
417 else
418 {
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];
423 return distances;
424 }
425 }
426
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
437 {
438 // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
439 // so the velocity at the boundary equal to that on the scvf.
440 const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
441 if (useZeroGradient)
442 return faceVars.velocitySelf();
443
444 const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
445 if (lateralFaceHasBJS)
446 return VelocityGradients::beaversJosephVelocityAtLateralScvf(elemVolVars_.gridVolVars().problem(), element, fvGeometry_, scvf, faceVars,
447 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
448
449 const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
450 if (lateralFaceHasDirichletVelocity)
451 {
452 // ________________
453 // ---------------o || frontal face of staggered half-control-volume
454 // | || # current scvf # ghostFace of interest, lies on boundary
455 // | || # x dof position
456 // | || x~~~~> vel.Self -- element boundaries
457 // | || # __ domain boundary
458 // | || # o position at which the boundary conditions will be evaluated
459 // ---------------#
460
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())];
465 }
466
467 // Neumann conditions are not well implemented
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());
470 }
471
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
480 {
481 // A ghost subface at the boundary is created, featuring the location of the sub face's center
482 const auto& lateralOppositeScvf = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localOppositeSubFaceIdx).localLateralFaceIdx);
483 GlobalPosition center = scvf.pairData(localOppositeSubFaceIdx).lateralStaggeredFaceCenter + lateralOppositeScvf.center();
484 center *= 0.5;
485
486 // lateral face # lateralFace currently being assembled
487 // --------######## || frontal face of staggered half-control-volume
488 // | || | current scvf % lateralOppositeBoundaryFace of interest, lies on boundary
489 // | || | x dof position
490 // | || x~~~~> vel.Self -- element boundaries
491 // | || | __ domain boundary
492 // | || | o position at which the boundary conditions will be evaluated
493 // --------%%%c%%%o
494 // ________________ c center of opposite boundary face
495
496
497 // Get the boundary types of the lateral opposite boundary face
498 const auto& problem = elemVolVars_.gridVolVars().problem();
499 const auto lateralOppositeBoundaryFace = makeStaggeredBoundaryFace(lateralOppositeScvf, center);
500 const auto lateralOppositeFaceBoundaryTypes = problem.boundaryTypes(element, lateralOppositeBoundaryFace);
501 return getParallelVelocityFromBoundary_(element, scvf, faceVars, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes, localOppositeSubFaceIdx);
502 }
503
523 const SubControlVolumeFace& boundaryScvf_(const int localSubFaceIdx) const
524 {
525 if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
526 {
527 return fvGeometry_.scvf(scvf_.outsideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
528 }
529 else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
530 {
531 /*
532 * ------------
533 * | xxxxxxxx o
534 * | xxxxxxxx o
535 * | xxxxxxxx o
536 * lllllllllll-bbbbbbbbbbb
537 * | yyyyyyyy p |
538 * | yyyyyyyy p |
539 * | yyyyyyyy p |
540 * -----------------------
541 *
542 * o: scvf_, l: lateralFace, p: parallelFace, b: returned scvf, x: scvf_ inside scv, y: lateralFace
543 * outside scv
544 */
545 const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
546 const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
547
548 const auto& localLateralIdx = scvf_.pairData(localSubFaceIdx).localLateralFaceIdx;
549 const auto& localLateralOppositeIdx = (localLateralIdx % 2) ? (localLateralIdx - 1) : (localLateralIdx + 1);
550
551 return fvGeometry_.scvf(parallelFace.outsideScvIdx(), localLateralOppositeIdx);
552 }
553 else
554 {
555 DUNE_THROW(Dune::InvalidStateException, "The function boundaryScvf_ should only be called when hasHalfParallelNeighbor or hasCornerParallelNeighbor.");
556 }
557 }
558
578 Element boundaryElement_(const int localSubFaceIdx) const
579 {
580 if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx))
581 {
582 return fvGeometry_.gridGeometry().element(scvf_.outsideScvIdx());
583 }
584 else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx))
585 {
586 const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx);
587 const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx());
588
589 return fvGeometry_.gridGeometry().element(parallelFace.outsideScvIdx());
590 }
591 else
592 {
593 DUNE_THROW(Dune::InvalidStateException, "When entering boundaryElement_ scvf_ should have either hasHalfParallelNeighbor or hasCornerParallelNeighbor true. Not the case here.");
594 }
595 }
596
616 bool dirichletParallelNeighbor_(const int localSubFaceIdx) const
617 {
618 const auto& problem = elemVolVars_.gridVolVars().problem();
619 const Element& boundaryElement = boundaryElement_(localSubFaceIdx);
620 const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx);
621
622 return problem.boundaryTypes(boundaryElement, boundaryScvf).isDirichlet(Indices::velocity(scvf_.directionIndex()));
623 }
624
643 Scalar getParallelVelocityFromCorner_(const int localSubFaceIdx) const
644 {
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);
649
650 return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf_.directionIndex())];
651 }
652
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_;
660};
661
662} // end namespace Dumux
663
664#endif
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
Definition: adapt.hh:29
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 > &currentScvfBoundaryTypes, 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 > &currentScvfBoundaryTypes, 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.