| 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 ShallowWaterFlux | ||
| 10 | * \brief This file includes a function to construct the Riemann problem. | ||
| 11 | * | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH | ||
| 14 | #define DUMUX_FLUX_SHALLOW_WATER_RIEMANN_PROBLEM_HH | ||
| 15 | |||
| 16 | #include <dumux/common/parameters.hh> | ||
| 17 | #include <dumux/flux/shallowwater/fluxlimiterlet.hh> | ||
| 18 | #include <dumux/flux/shallowwater/exactriemann.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | namespace ShallowWater { | ||
| 22 | |||
| 23 | /*! | ||
| 24 | * \ingroup ShallowWaterFlux | ||
| 25 | * \brief Construct a Riemann problem and solve it. | ||
| 26 | * | ||
| 27 | * | ||
| 28 | * Riemann problem applies the hydrostatic reconstruction, uses the | ||
| 29 | * Riemann invariants to transform the two-dimensional problem to a | ||
| 30 | * one-dimensional problem, solves this new problem, and rotates | ||
| 31 | * the problem back. Further it applies a flux limiter for the water | ||
| 32 | * flux to handle drying elements. | ||
| 33 | * The correction of the bed slope source term leads to a | ||
| 34 | * non-symmetric flux term at the interface for the momentum equations. | ||
| 35 | * Since DuMuX computes the fluxes twice from each side this does not | ||
| 36 | * matter. | ||
| 37 | * | ||
| 38 | * So far this implements the exact Riemann solver (with reconstruction | ||
| 39 | * after Audusse). | ||
| 40 | * | ||
| 41 | * The computed water flux (localFlux[0]) is given in m^2/s, the | ||
| 42 | * momentum fluxes (localFlux[1], localFlux[2]) are given in m^3/s^2. | ||
| 43 | * Later this flux will be multiplied by the scvf.area() (given in m | ||
| 44 | * for a 2D problem) to get the flux over a face. | ||
| 45 | * | ||
| 46 | * \param waterDepthLeft water depth on the left side | ||
| 47 | * \param waterDepthRight water depth on the right side | ||
| 48 | * \param velocityXLeft veloctiyX on the left side | ||
| 49 | * \param velocityXRight velocityX on the right side | ||
| 50 | * \param velocityYLeft velocityY on the left side | ||
| 51 | * \param velocityYRight velocityY on the right side | ||
| 52 | * \param bedSurfaceLeft surface of the bed on the left side | ||
| 53 | * \param bedSurfaceRight surface of the bed on the right side | ||
| 54 | * \param gravity gravity constant | ||
| 55 | * \param nxy the normal vector | ||
| 56 | * | ||
| 57 | */ | ||
| 58 | template<class Scalar, class GlobalPosition> | ||
| 59 | 159362171 | std::array<Scalar,3> riemannProblem(const Scalar waterDepthLeft, | |
| 60 | const Scalar waterDepthRight, | ||
| 61 | Scalar velocityXLeft, | ||
| 62 | Scalar velocityXRight, | ||
| 63 | Scalar velocityYLeft, | ||
| 64 | Scalar velocityYRight, | ||
| 65 | const Scalar bedSurfaceLeft, | ||
| 66 | const Scalar bedSurfaceRight, | ||
| 67 | const Scalar gravity, | ||
| 68 | const GlobalPosition& nxy) | ||
| 69 | { | ||
| 70 | using std::max; | ||
| 71 | |||
| 72 | // hydrostatic reconstrucion after Audusse | ||
| 73 |
2/2✓ Branch 0 taken 33158762 times.
✓ Branch 1 taken 126203409 times.
|
159362171 | const Scalar dzl = max(0.0, bedSurfaceRight - bedSurfaceLeft); |
| 74 |
2/2✓ Branch 0 taken 145972315 times.
✓ Branch 1 taken 13389856 times.
|
159362171 | const Scalar waterDepthLeftReconstructed = max(0.0, waterDepthLeft - dzl); |
| 75 |
2/2✓ Branch 0 taken 33158855 times.
✓ Branch 1 taken 126203316 times.
|
159362171 | const Scalar dzr = max(0.0, bedSurfaceLeft - bedSurfaceRight); |
| 76 |
2/2✓ Branch 0 taken 145972315 times.
✓ Branch 1 taken 13389856 times.
|
159362171 | const Scalar waterDepthRightReconstructed = max(0.0, waterDepthRight - dzr); |
| 77 | |||
| 78 | // make rotation of the flux we compute an 1d flux | ||
| 79 | 159362171 | Scalar tempFlux = velocityXLeft; | |
| 80 | 318724342 | velocityXLeft = nxy[0] * tempFlux + nxy[1] * velocityYLeft; | |
| 81 | 318724342 | velocityYLeft = -nxy[1] * tempFlux + nxy[0] * velocityYLeft; | |
| 82 | |||
| 83 | 159362171 | tempFlux = velocityXRight; | |
| 84 | 318724342 | velocityXRight = nxy[0] * tempFlux + nxy[1] * velocityYRight; | |
| 85 | 318724342 | velocityYRight = -nxy[1] * tempFlux + nxy[0] * velocityYRight; | |
| 86 | |||
| 87 | 159362171 | auto riemannResult = ShallowWater::exactRiemann(waterDepthLeftReconstructed, | |
| 88 | waterDepthRightReconstructed, | ||
| 89 | velocityXLeft, | ||
| 90 | velocityXRight, | ||
| 91 | velocityYLeft, | ||
| 92 | velocityYRight, | ||
| 93 | gravity); | ||
| 94 | |||
| 95 | //redo rotation | ||
| 96 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
|
159362171 | tempFlux = riemannResult.flux[1]; |
| 97 |
8/8✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
✓ Branch 2 taken 16 times.
✓ Branch 3 taken 159362155 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 159362155 times.
✓ Branch 6 taken 16 times.
✓ Branch 7 taken 159362155 times.
|
637448684 | riemannResult.flux[1] = nxy[0] * tempFlux - nxy[1] * riemannResult.flux[2]; |
| 98 |
8/8✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
✓ Branch 2 taken 16 times.
✓ Branch 3 taken 159362155 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 159362155 times.
✓ Branch 6 taken 16 times.
✓ Branch 7 taken 159362155 times.
|
637448684 | riemannResult.flux[2] = nxy[1] * tempFlux + nxy[0] * riemannResult.flux[2]; |
| 99 | |||
| 100 | // Add reconstruction flux from Audusse reconstruction | ||
| 101 | 159362171 | const Scalar hgzl = 0.5 * (waterDepthLeftReconstructed + waterDepthLeft) * (waterDepthLeftReconstructed - waterDepthLeft); | |
| 102 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
|
159362171 | const Scalar hdxzl = gravity * nxy[0] * hgzl; |
| 103 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
|
159362171 | const Scalar hdyzl = gravity * nxy[1] * hgzl; |
| 104 | |||
| 105 | /*Right side is computed from the other side otherwise the | ||
| 106 | following "non-symetric" fluxes are needed: | ||
| 107 | |||
| 108 | Scalar hgzr = 0.5 * (waterDepthRightReconstructed + waterDepthRight) * (waterDepthRightReconstructed - waterDepthRight); | ||
| 109 | Scalar hdxzr = gravity * nxy[0] * hgzr; | ||
| 110 | Scalar hdyzrhdyzr = gravity * nxy[1] * hgzr; | ||
| 111 | */ | ||
| 112 | |||
| 113 | // compute the mobility of the flux with the fluxlimiter | ||
| 114 |
4/6✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
|
159362171 | static const Scalar upperWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.UpperWaterDepth", 1e-3); |
| 115 |
4/6✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
|
159362171 | static const Scalar lowerWaterDepthFluxLimiting = getParam<Scalar>("FluxLimiterLET.LowerWaterDepth", 1e-5); |
| 116 |
4/6✓ Branch 0 taken 16 times.
✓ Branch 1 taken 159362155 times.
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
|
159362171 | static const bool upwindWaterDepthFluxLimiting = getParam<bool>("FluxLimiterLET.UpwindFluxLimiting", false); |
| 117 | |||
| 118 | 159362171 | Scalar limitingDepth = (waterDepthLeftReconstructed + waterDepthRightReconstructed) * 0.5; | |
| 119 | |||
| 120 | //Using the upwind water depth from the flux direction can improve stability. | ||
| 121 |
2/2✓ Branch 0 taken 55220864 times.
✓ Branch 1 taken 104141307 times.
|
159362171 | if (upwindWaterDepthFluxLimiting) |
| 122 | { | ||
| 123 |
4/4✓ Branch 0 taken 27965820 times.
✓ Branch 1 taken 27255044 times.
✓ Branch 2 taken 27965820 times.
✓ Branch 3 taken 27255044 times.
|
110441728 | if (riemannResult.flux[0] < 0) |
| 124 | { | ||
| 125 | limitingDepth = waterDepthRightReconstructed; | ||
| 126 | }else | ||
| 127 | { | ||
| 128 | 27965820 | limitingDepth = waterDepthLeftReconstructed; | |
| 129 | } | ||
| 130 | } | ||
| 131 | |||
| 132 | 159362171 | const Scalar mobility = ShallowWater::fluxLimiterLET(limitingDepth, | |
| 133 | limitingDepth, | ||
| 134 | upperWaterDepthFluxLimiting, | ||
| 135 | lowerWaterDepthFluxLimiting); | ||
| 136 | std::array<Scalar, 3> localFlux; | ||
| 137 | 318724342 | localFlux[0] = riemannResult.flux[0] * mobility; | |
| 138 | 318724342 | localFlux[1] = (riemannResult.flux[1] - hdxzl); | |
| 139 | 318724342 | localFlux[2] = (riemannResult.flux[2] - hdyzl); | |
| 140 | |||
| 141 | 159362171 | return localFlux; | |
| 142 | } | ||
| 143 | |||
| 144 | } // end namespace ShallowWater | ||
| 145 | } // end namespace Dumux | ||
| 146 | |||
| 147 | #endif | ||
| 148 |