| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup NavierStokesModel | ||
| 10 | * \copydoc Dumux::StaggeredVelocityReconstruction | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYRECONSTRUCTION_HH | ||
| 13 | #define DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYRECONSTRUCTION_HH | ||
| 14 | |||
| 15 | #include <numeric> | ||
| 16 | #include <algorithm> | ||
| 17 | #include <dune/common/reservedvector.hh> | ||
| 18 | #include <dumux/discretization/method.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup NavierStokesModel | ||
| 24 | * \brief Helper class for reconstructing the velocity. | ||
| 25 | */ | ||
| 26 | struct StaggeredVelocityReconstruction | ||
| 27 | { | ||
| 28 | //! Return the velocity vector at the center of the primal grid. | ||
| 29 | template<class VelocityHelper, class FVElementGeometry> | ||
| 30 | 922347 | static auto cellCenterVelocity(const VelocityHelper& getFaceVelocity, | |
| 31 | const FVElementGeometry& fvGeometry) | ||
| 32 | { | ||
| 33 | static_assert(FVElementGeometry::GridGeometry::discMethod == DiscretizationMethods::cctpfa); | ||
| 34 | using VelocityVector = typename FVElementGeometry::GridGeometry::GlobalCoordinate; | ||
| 35 | 922347 | VelocityVector result(0.0); | |
| 36 | |||
| 37 | 3693024 | const auto directionIndex = [&](const auto& vector) | |
| 38 | { | ||
| 39 |
14/36✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 128 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 128 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3728 times.
✓ Branch 17 taken 7456 times.
✓ Branch 18 taken 3728 times.
✓ Branch 19 taken 7456 times.
✓ Branch 20 taken 1844548 times.
✓ Branch 21 taken 1844676 times.
✓ Branch 22 taken 1844548 times.
✓ Branch 23 taken 1844548 times.
✓ Branch 24 taken 1844548 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 1844548 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 3692824 times.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
|
25861864 | return std::find_if(vector.begin(), vector.end(), [](const auto& x) { return std::abs(x) > 1e-8; } ) - vector.begin(); |
| 40 | }; | ||
| 41 | |||
| 42 |
4/4✓ Branch 0 taken 3692988 times.
✓ Branch 1 taken 922347 times.
✓ Branch 2 taken 3692988 times.
✓ Branch 3 taken 922347 times.
|
9230670 | for (const auto& scvf : scvfs(fvGeometry)) |
| 43 | { | ||
| 44 | 7385976 | const auto dirIdx = directionIndex(scvf.unitOuterNormal()); | |
| 45 | 3693096 | result[dirIdx] += 0.5*getFaceVelocity(fvGeometry, scvf)[dirIdx]; | |
| 46 | } | ||
| 47 | 922347 | return result; | |
| 48 | } | ||
| 49 | |||
| 50 | //! Return the velocity vector at dof position given an scv | ||
| 51 | template<class SubControlVolume, class FVElementGeometry, class VelocityHelper> | ||
| 52 | 90380 | static auto faceVelocity(const SubControlVolume& scv, | |
| 53 | const FVElementGeometry& fvGeometry, | ||
| 54 | const VelocityHelper& getNormalVelocityDofValue) | ||
| 55 | { | ||
| 56 | 90380 | int axis = scv.dofAxis(); | |
| 57 | 90380 | const int dim = FVElementGeometry::GridGeometry::GridView::dimension; | |
| 58 | using Scalar = typename SubControlVolume::Traits::Scalar; | ||
| 59 | using GlobalPosition = typename FVElementGeometry::GridGeometry::GlobalCoordinate; | ||
| 60 | using VelocityVector = typename FVElementGeometry::GridGeometry::GlobalCoordinate; | ||
| 61 | 90380 | VelocityVector faceVelocityVector(0.0); | |
| 62 | |||
| 63 | // per dimension, we have at max two velocities from which we'll compute an average | ||
| 64 |
4/4✓ Branch 0 taken 90380 times.
✓ Branch 1 taken 90380 times.
✓ Branch 2 taken 90380 times.
✓ Branch 3 taken 90380 times.
|
271140 | std::array<Dune::ReservedVector<Scalar, 2>, dim> normalVelocities; |
| 65 |
4/6✓ Branch 0 taken 576 times.
✓ Branch 1 taken 89804 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 576 times.
✓ Branch 5 taken 576 times.
✗ Branch 6 not taken.
|
90380 | if (scv.boundary() && !fvGeometry.gridGeometry().dofOnPeriodicBoundary(scv.dofIndex())) |
| 66 | { | ||
| 67 | // iterate through the inner lateral velocities, | ||
| 68 |
4/4✓ Branch 1 taken 2304 times.
✓ Branch 2 taken 576 times.
✓ Branch 3 taken 2304 times.
✓ Branch 4 taken 576 times.
|
2880 | for (const auto& scvf : scvfs(fvGeometry, scv)) |
| 69 | { | ||
| 70 |
2/2✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 1152 times.
|
2304 | if (scvf.isFrontal()) |
| 71 | continue; | ||
| 72 | |||
| 73 | // at a lateral velocity, find the inner and outer normal velocities | ||
| 74 | 1152 | const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf); | |
| 75 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1152 times.
|
2304 | const auto& lateralScv = fvGeometry.scv(orthogonalScvf.insideScvIdx()); |
| 76 | 1152 | auto lateralAxis = lateralScv.dofAxis(); | |
| 77 | 2304 | normalVelocities[lateralAxis].push_back( getNormalVelocityDofValue(lateralScv) ) ; | |
| 78 | } | ||
| 79 | } | ||
| 80 | else | ||
| 81 | { | ||
| 82 | // Find the location of interpolation, if periodic, from both sides | ||
| 83 | 89804 | const GlobalPosition selfPosition = scv.dofPosition(); | |
| 84 | 89804 | GlobalPosition outsideSelfPosition = selfPosition; | |
| 85 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 89804 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 89792 times.
|
89804 | if (fvGeometry.gridGeometry().dofOnPeriodicBoundary(scv.dofIndex())) |
| 86 | 12 | outsideSelfPosition = fvGeometry.outsidePeriodicScv(scv).dofPosition(); | |
| 87 | |||
| 88 | // iterate through the lateral velocities to get to the normal locations | ||
| 89 |
4/4✓ Branch 1 taken 269412 times.
✓ Branch 2 taken 89804 times.
✓ Branch 3 taken 269412 times.
✓ Branch 4 taken 89804 times.
|
359216 | for (const auto& scvf : scvfs(fvGeometry, scv)) |
| 90 | { | ||
| 91 |
2/2✓ Branch 0 taken 179608 times.
✓ Branch 1 taken 89804 times.
|
269412 | if (scvf.isFrontal()) |
| 92 | continue; | ||
| 93 | |||
| 94 | // at a lateral velocity, find the inner and outer normal velocities | ||
| 95 | 179608 | const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf); | |
| 96 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 179608 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 179608 times.
|
359216 | const auto& insideLateralScv = fvGeometry.scv(orthogonalScvf.insideScvIdx()); |
| 97 | 359216 | const auto& outsideLateralScv = fvGeometry.scv(orthogonalScvf.outsideScvIdx()); | |
| 98 | 179608 | const auto& lateralAxis = insideLateralScv.dofAxis(); | |
| 99 | |||
| 100 | // Find the inside normal velocities | ||
| 101 | 179608 | const auto& insideNormalVelocity = getNormalVelocityDofValue(insideLateralScv); | |
| 102 | 359216 | const auto& insideNormalPosition = insideLateralScv.dofPosition()[axis]; | |
| 103 | |||
| 104 | // Find the outside normal velocities | ||
| 105 | 179608 | const auto& outsideNormalVelocity = getNormalVelocityDofValue(outsideLateralScv); | |
| 106 | 359216 | const auto& outsideNormalPosition = outsideLateralScv.dofPosition()[axis]; | |
| 107 | |||
| 108 | // Linear interpolation at the face plane and add to normal velocity collection | ||
| 109 | 359216 | const auto& innerDistance = std::abs(insideNormalPosition - selfPosition[axis]); | |
| 110 | 359216 | const auto& totalDistance = std::abs(outsideNormalPosition - outsideSelfPosition[axis]) + innerDistance; | |
| 111 | 179608 | const auto& velDiff = outsideNormalVelocity - insideNormalVelocity; | |
| 112 | |||
| 113 | 538824 | normalVelocities[lateralAxis].push_back(insideNormalVelocity + (velDiff * innerDistance / totalDistance)); | |
| 114 | } | ||
| 115 | } | ||
| 116 | |||
| 117 |
2/2✓ Branch 0 taken 180760 times.
✓ Branch 1 taken 90380 times.
|
271140 | for (int i = 0; i < faceVelocityVector.size(); i++) |
| 118 | { | ||
| 119 |
2/2✓ Branch 0 taken 90380 times.
✓ Branch 1 taken 90380 times.
|
180760 | if (i == axis) |
| 120 | 90452 | faceVelocityVector[i] = getNormalVelocityDofValue(scv); | |
| 121 | else | ||
| 122 | 723040 | faceVelocityVector[i] = std::accumulate(normalVelocities[i].begin(), normalVelocities[i].end(), 0.0) / normalVelocities[i].size(); | |
| 123 | } | ||
| 124 | |||
| 125 | 90380 | return faceVelocityVector; | |
| 126 | } | ||
| 127 | |||
| 128 | }; | ||
| 129 | |||
| 130 | } // end namespace Dumux | ||
| 131 | |||
| 132 | #endif // DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYRECONSTRUCTION_HH | ||
| 133 |