| 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 Function to compute the Riemann flux at the interface | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_FLUX_SHALLOW_WATER_EXACT_RIEMANN_HH | ||
| 13 | #define DUMUX_FLUX_SHALLOW_WATER_EXACT_RIEMANN_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <cmath> | ||
| 17 | |||
| 18 | namespace Dumux { | ||
| 19 | namespace ShallowWater{ | ||
| 20 | |||
| 21 | template<typename Scalar> | ||
| 22 | struct RiemannSolution { | ||
| 23 | std::array<Scalar,3> flux; | ||
| 24 | Scalar waterDepth; | ||
| 25 | Scalar velocityX; | ||
| 26 | Scalar velocityY; | ||
| 27 | }; | ||
| 28 | |||
| 29 | |||
| 30 | /*! | ||
| 31 | * \ingroup ShallowWaterFlux | ||
| 32 | * \brief Exact Riemann solver for the shallow water equations. | ||
| 33 | * | ||
| 34 | * The flux of the 2D shallow water equations must be rotated | ||
| 35 | * to a 1D problem before the Riemann solver can be applied. | ||
| 36 | * The computed water flux is given in m^2/s, the momentum | ||
| 37 | * fluxes are given in m^3/s^2. | ||
| 38 | * | ||
| 39 | * This Riemann solver is described in the book | ||
| 40 | * "Shock-capturing methods for free-surface shallow flows" | ||
| 41 | * from Toro, 2001. We keep the notation for the variables | ||
| 42 | * after Toro. | ||
| 43 | * | ||
| 44 | * \param dl water depth on the left side | ||
| 45 | * \param dr water depth on the right side | ||
| 46 | * \param ul veloctiyX on the left side | ||
| 47 | * \param ur velocityX on the right side | ||
| 48 | * \param vl velocityY on the left side | ||
| 49 | * \param vr velocityY on the right side | ||
| 50 | * \param grav gravity constant | ||
| 51 | * \param s sample point (default = 0 since x = 0 for flux computation) | ||
| 52 | */ | ||
| 53 | template<class Scalar> | ||
| 54 | 160377221 | RiemannSolution<Scalar> exactRiemann(const Scalar dl, | |
| 55 | const Scalar dr, | ||
| 56 | const Scalar ul, | ||
| 57 | const Scalar ur, | ||
| 58 | const Scalar vl, | ||
| 59 | const Scalar vr, | ||
| 60 | const Scalar grav, | ||
| 61 | const Scalar s = 0.0) | ||
| 62 | { | ||
| 63 | RiemannSolution<Scalar> sol; | ||
| 64 | 160377221 | sol.waterDepth = 0.0; // d (Toro) | |
| 65 | 160377221 | sol.velocityX = 0.0; // u (Toro) | |
| 66 | 160377221 | sol.velocityY = 0.0; // v (Toro) | |
| 67 | |||
| 68 | 160377221 | constexpr int maxsteps = 200; //maximum steps for the Newton method | |
| 69 | 160377221 | constexpr Scalar tol = 1.0E-12; //tolerance for the Newton method | |
| 70 | |||
| 71 | using std::sqrt; | ||
| 72 | using std::max; | ||
| 73 | using std::min; | ||
| 74 | using std::abs; | ||
| 75 | |||
| 76 | //================== Exact Riemann solver | ||
| 77 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 160377221 times.
|
160377221 | const auto cl = sqrt(grav * max(dl, 0.0)); //celerity sqrt(g*h) left side |
| 78 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 160377221 times.
|
160377221 | const auto cr = sqrt(grav * max(dr, 0.0)); //celerity sqrt(g*h) right side |
| 79 | 160377221 | const auto dcrit = (ur - ul) - 2.0*(cl + cr); //critical water depth | |
| 80 | |||
| 81 | // dry case | ||
| 82 |
6/6✓ Branch 0 taken 146987365 times.
✓ Branch 1 taken 13389856 times.
✓ Branch 2 taken 133597509 times.
✓ Branch 3 taken 13389856 times.
✓ Branch 4 taken 512 times.
✓ Branch 5 taken 133596997 times.
|
160377221 | if ((dl <= 0) || (dr <= 0) || (dcrit >= 0)) |
| 83 | { | ||
| 84 | //left side is dry | ||
| 85 |
2/2✓ Branch 0 taken 13389856 times.
✓ Branch 1 taken 13390368 times.
|
26780224 | if(dl <= 0) |
| 86 | { | ||
| 87 | 13389856 | const auto shr = ur + cr; //wave speed at head right | |
| 88 | |||
| 89 |
2/2✓ Branch 0 taken 32768 times.
✓ Branch 1 taken 13357088 times.
|
13389856 | if (s >= shr) |
| 90 | { | ||
| 91 | 32768 | sol.waterDepth = dr; | |
| 92 | 32768 | sol.velocityX = ur; | |
| 93 | 32768 | sol.velocityY = vr; | |
| 94 | |||
| 95 | } | ||
| 96 | |||
| 97 | else | ||
| 98 | { | ||
| 99 | |||
| 100 | 13357088 | const auto str = ur - 2.0*cr; //wave speed at tail right | |
| 101 | |||
| 102 |
2/2✓ Branch 0 taken 13350048 times.
✓ Branch 1 taken 7040 times.
|
13357088 | if(s >= str){ |
| 103 | 13350048 | sol.velocityX = (ur - 2.0 *cr + 2.0*s)/3.0; | |
| 104 | 13350048 | const auto c = (-ur + 2.0*cr +s )/3.0; | |
| 105 | 13350048 | sol.waterDepth = c*c/grav; | |
| 106 | 13350048 | sol.velocityY = vr; | |
| 107 | |||
| 108 | } | ||
| 109 | |||
| 110 | else | ||
| 111 | { | ||
| 112 | 7040 | sol.waterDepth = dl; | |
| 113 | 7040 | sol.velocityX = ul; | |
| 114 | 7040 | sol.velocityY = vl; | |
| 115 | } | ||
| 116 | } | ||
| 117 | } | ||
| 118 | |||
| 119 | else | ||
| 120 | { | ||
| 121 | //right side is dry | ||
| 122 |
2/2✓ Branch 0 taken 13389856 times.
✓ Branch 1 taken 512 times.
|
13390368 | if(dr <= 0) |
| 123 | { | ||
| 124 | 13389856 | const auto shl = ul-cl; //wave speed at head left | |
| 125 | |||
| 126 |
2/2✓ Branch 0 taken 32768 times.
✓ Branch 1 taken 13357088 times.
|
13389856 | if(s <= shl) |
| 127 | { | ||
| 128 | 32768 | sol.waterDepth = dl; | |
| 129 | 32768 | sol.velocityX = ul; | |
| 130 | 32768 | sol.velocityY = vl; | |
| 131 | } | ||
| 132 | |||
| 133 | else | ||
| 134 | { | ||
| 135 | 13357088 | const auto stl = ul + 2.0 *cl; //wave speed at tail left | |
| 136 | |||
| 137 |
2/2✓ Branch 0 taken 13350048 times.
✓ Branch 1 taken 7040 times.
|
13357088 | if ( s <= stl) |
| 138 | { | ||
| 139 | 13350048 | sol.velocityX = (ul + 2.0 *cl + 2.0*s)/3.0; | |
| 140 | 13350048 | const auto c = (ul + 2.0*cl - s)/3.0; | |
| 141 | 13350048 | sol.waterDepth = c*c/grav; | |
| 142 | 13350048 | sol.velocityY = vl; | |
| 143 | } | ||
| 144 | |||
| 145 | else | ||
| 146 | { | ||
| 147 | 7040 | sol.waterDepth = dr; | |
| 148 | 7040 | sol.velocityX = ur; | |
| 149 | 7040 | sol.velocityY = vr; | |
| 150 | } | ||
| 151 | } | ||
| 152 | } | ||
| 153 | |||
| 154 | //middle state is dry | ||
| 155 | else | ||
| 156 | { | ||
| 157 | 512 | const auto shl = ul - cl; //wave speed at head left | |
| 158 | 512 | const auto shr = ur + cr; //wave speed at head right | |
| 159 | 512 | const auto ssl = ul + 2.0 * cl; // wave speed of a wet/dry front for a left dry state | |
| 160 | 512 | const auto ssr = ur - 2.0 * cr; // wave speed of a wet/dry front for a right dry state | |
| 161 | |||
| 162 |
2/2✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
|
512 | if(s <= shl) |
| 163 | { | ||
| 164 | 256 | sol.waterDepth = dl; | |
| 165 | 256 | sol.velocityX = ul; | |
| 166 | 256 | sol.velocityY = vl; | |
| 167 | } | ||
| 168 | |||
| 169 |
3/4✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 256 times.
|
512 | if((s > shl) && (s <= ssl)) |
| 170 | { | ||
| 171 | ✗ | sol.velocityX = (ul + 2.0 * cl + 2.0*s)/3.0; | |
| 172 | ✗ | const auto c = (ul + 2.0*cl - s)/3.0; | |
| 173 | ✗ | sol.waterDepth = c*c/grav; | |
| 174 | ✗ | sol.velocityY = vl; | |
| 175 | } | ||
| 176 | |||
| 177 |
3/4✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 256 times.
|
512 | if((s > ssl) && (s <= ssr)) |
| 178 | { | ||
| 179 | ✗ | sol.waterDepth = 0.0; | |
| 180 | ✗ | sol.velocityX = 0.0; | |
| 181 | ✗ | sol.velocityY = 0.0; | |
| 182 | } | ||
| 183 | |||
| 184 |
3/4✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 256 times.
|
512 | if ((s > ssr ) && (s <= shr)) |
| 185 | { | ||
| 186 | ✗ | sol.velocityX = (ur - 2.0*cr + 2.0*s)/3.0; | |
| 187 | ✗ | const auto c = (-ur + 2.0*cr + s)/3.0; | |
| 188 | ✗ | sol.waterDepth = c*c/grav; | |
| 189 | ✗ | sol.velocityY = vr; | |
| 190 | } | ||
| 191 | |||
| 192 |
2/2✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
|
512 | if(s > shr) |
| 193 | { | ||
| 194 | 256 | sol.waterDepth = dr; | |
| 195 | 256 | sol.velocityX = ur; | |
| 196 | 256 | sol.velocityY = vr; | |
| 197 | } | ||
| 198 | } | ||
| 199 | } | ||
| 200 | } | ||
| 201 | |||
| 202 | // wet case | ||
| 203 | else | ||
| 204 | { | ||
| 205 | // Get a starting value for the Newton-Raphson method | ||
| 206 | 133596997 | const auto tmp = 0.5 * (cl + cr) - 0.25 * (ul - ur); | |
| 207 | 133596997 | auto ds = (1.0/grav) * tmp * tmp; // water depth | |
| 208 | |||
| 209 | //default is Two-Rarefaction as starting value, if ds > dmin we use two-shock | ||
| 210 |
2/2✓ Branch 0 taken 57336037 times.
✓ Branch 1 taken 76260960 times.
|
133596997 | const auto dmin = min(dl, dr); // minimum depth |
| 211 | |||
| 212 |
2/2✓ Branch 0 taken 90031483 times.
✓ Branch 1 taken 43565514 times.
|
133596997 | if (ds > dmin) |
| 213 | { | ||
| 214 | // function g in (Toro, 2001) | ||
| 215 | 90031483 | const auto gel = sqrt(0.5*grav * (ds+dl)/(ds*dl)); | |
| 216 | 90031483 | const auto ger = sqrt(0.5*grav * (ds+dr)/(ds*dr)); | |
| 217 | 90031483 | ds = (gel * dl + ger * dr - (ur-ul))/(gel + ger); | |
| 218 | } | ||
| 219 | |||
| 220 | //Start Newton-Raphson loop ds = hstar | ||
| 221 |
2/2✓ Branch 0 taken 4928 times.
✓ Branch 1 taken 133592069 times.
|
133596997 | ds = max(ds, tol); |
| 222 | 133596997 | auto d0 = ds; // initial guess water depth | |
| 223 | |||
| 224 | // helper functions (left and right) and their derivatives | ||
| 225 | 133596997 | Scalar fl = 0.0, fr = 0.0, fld = 0.0, frd = 0.0; | |
| 226 | |||
| 227 |
1/2✓ Branch 0 taken 280974922 times.
✗ Branch 1 not taken.
|
280974922 | for (int i=0; i <= maxsteps; ++i) |
| 228 | { | ||
| 229 | //Compute fl,fr,fld,frd | ||
| 230 | |||
| 231 | //first left side | ||
| 232 |
2/2✓ Branch 0 taken 151874743 times.
✓ Branch 1 taken 129100179 times.
|
280974922 | if (ds <= dl) // wave is rarefaction wave (or depression) |
| 233 | { | ||
| 234 | 151874743 | const auto c = sqrt(grav * ds); | |
| 235 | 151874743 | fl = 2.0 * (c-cl); | |
| 236 | 151874743 | fld = grav/c; | |
| 237 | } | ||
| 238 | |||
| 239 | else // wave is shock wave (or bore) | ||
| 240 | { | ||
| 241 | 129100179 | const auto ges = sqrt(0.5 * grav * (ds+dl)/(ds*dl)); // function g in (Toro, 2001) | |
| 242 | 129100179 | fl = (ds - dl) * ges; | |
| 243 | 129100179 | fld = ges - 0.25 * grav * (ds - dl)/(ges * ds * ds); | |
| 244 | } | ||
| 245 | |||
| 246 | //second right side | ||
| 247 |
2/2✓ Branch 0 taken 147853276 times.
✓ Branch 1 taken 133121646 times.
|
280974922 | if (ds <= dr) |
| 248 | { | ||
| 249 | 147853276 | const auto c = sqrt(grav * ds); | |
| 250 | 147853276 | fr = 2.0 * (c-cr); | |
| 251 | 147853276 | frd = grav/c; | |
| 252 | } | ||
| 253 | |||
| 254 | else | ||
| 255 | { | ||
| 256 | 133121646 | const auto ges = sqrt(0.5 * grav * (ds+dr)/(ds*dr)); | |
| 257 | 133121646 | fr = (ds - dr) * ges; | |
| 258 | 133121646 | frd = ges - 0.25 * grav * (ds - dr)/(ges * ds * ds); | |
| 259 | } | ||
| 260 | |||
| 261 | 280974922 | ds -= (fl + fr + ur - ul)/(fld + frd); | |
| 262 |
2/2✓ Branch 0 taken 147377925 times.
✓ Branch 1 taken 133596997 times.
|
280974922 | const auto cha = abs(ds - d0)/(0.5*(ds + d0)); |
| 263 | |||
| 264 |
2/2✓ Branch 0 taken 147377925 times.
✓ Branch 1 taken 133596997 times.
|
280974922 | if (cha <= tol) |
| 265 | { | ||
| 266 | break; | ||
| 267 | } | ||
| 268 | |||
| 269 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 147377925 times.
|
147377925 | if (ds < 0.0) |
| 270 | { | ||
| 271 | ✗ | ds = tol; | |
| 272 | } | ||
| 273 | |||
| 274 | 147377925 | d0 = ds; | |
| 275 | } | ||
| 276 | |||
| 277 | //compute ustar (us) | ||
| 278 |
2/2✓ Branch 0 taken 69546710 times.
✓ Branch 1 taken 64050287 times.
|
133596997 | if (ds <= dl) |
| 279 | { | ||
| 280 | 69546710 | const auto c = sqrt(grav * ds); | |
| 281 | 69546710 | fl = 2.0 * (c-cl); | |
| 282 | } | ||
| 283 | |||
| 284 | else | ||
| 285 | { | ||
| 286 | 64050287 | const auto ges = sqrt(0.5*grav*(ds + dl)/(ds*dl)); | |
| 287 | 64050287 | fl = (ds - dl) * ges; | |
| 288 | } | ||
| 289 | |||
| 290 |
2/2✓ Branch 0 taken 68550745 times.
✓ Branch 1 taken 65046252 times.
|
133596997 | if (ds <= dr) |
| 291 | { | ||
| 292 | 68550745 | const auto c = sqrt(grav * ds); | |
| 293 | 68550745 | fr = 2.0 * (c-cr); | |
| 294 | } | ||
| 295 | |||
| 296 | else | ||
| 297 | { | ||
| 298 | 65046252 | const auto ges = sqrt(0.5*grav*(ds + dr)/(ds*dr)); | |
| 299 | 65046252 | fr = (ds - dr) * ges; | |
| 300 | } | ||
| 301 | //exact Riemann solver sstar = ustar | ||
| 302 | 133596997 | const auto us = 0.5*(ul + ur) + 0.5*(fr - fl); // velocity in the star region (u_star) | |
| 303 | 133596997 | const auto cs = sqrt(grav*ds); // wave speed in the star region | |
| 304 | |||
| 305 | |||
| 306 | /*********************** computation of u and d *******************/ | ||
| 307 | //left wave | ||
| 308 |
2/2✓ Branch 0 taken 68350026 times.
✓ Branch 1 taken 65246971 times.
|
133596997 | if (s <= us) |
| 309 | { | ||
| 310 | 68350026 | sol.velocityY = vl; | |
| 311 | //left shock | ||
| 312 |
2/2✓ Branch 0 taken 25333497 times.
✓ Branch 1 taken 43016529 times.
|
68350026 | if (ds >= dl) |
| 313 | { | ||
| 314 | 25333497 | const auto ql = sqrt((ds + dl)*ds/(2.0*dl*dl)); // flux left | |
| 315 | 25333497 | const auto sl = ul - cl*ql; // shock left | |
| 316 | |||
| 317 | //left side of shock | ||
| 318 |
2/2✓ Branch 0 taken 11111 times.
✓ Branch 1 taken 25322386 times.
|
25333497 | if(s <= sl) |
| 319 | { | ||
| 320 | 11111 | sol.waterDepth = dl; | |
| 321 | 11111 | sol.velocityX = ul; | |
| 322 | } | ||
| 323 | |||
| 324 | //right side of shock | ||
| 325 | else | ||
| 326 | { | ||
| 327 | 25322386 | sol.waterDepth = ds; | |
| 328 | 25322386 | sol.velocityX = us; | |
| 329 | } | ||
| 330 | } | ||
| 331 | |||
| 332 | //left rarefaction | ||
| 333 | else | ||
| 334 | { | ||
| 335 | 43016529 | const auto shl = ul-cl; //wave speed at head left | |
| 336 | |||
| 337 | //right side of rarefaction | ||
| 338 |
2/2✓ Branch 0 taken 689192 times.
✓ Branch 1 taken 42327337 times.
|
43016529 | if (s <= shl) |
| 339 | { | ||
| 340 | 689192 | sol.waterDepth = dl; | |
| 341 | 689192 | sol.velocityX = ul; | |
| 342 | } | ||
| 343 | |||
| 344 | else | ||
| 345 | { | ||
| 346 | 42327337 | const auto stl = us -cs; //wave speed at tail left | |
| 347 | |||
| 348 | //inside the rarefaction | ||
| 349 |
2/2✓ Branch 0 taken 313416 times.
✓ Branch 1 taken 42013921 times.
|
42327337 | if(s <= stl) |
| 350 | { | ||
| 351 | 313416 | sol.velocityX = (ul + 2.0*cl + 2.0 * s)/3.0; | |
| 352 | 313416 | const auto c = (ul + 2.0* cl - s)/3.0; | |
| 353 | 313416 | sol.waterDepth = c*c/grav; | |
| 354 | } | ||
| 355 | |||
| 356 | //inside the star region | ||
| 357 | else | ||
| 358 | { | ||
| 359 | 42013921 | sol.waterDepth = ds; | |
| 360 | 42013921 | sol.velocityX = us; | |
| 361 | } | ||
| 362 | } | ||
| 363 | } | ||
| 364 | } | ||
| 365 | |||
| 366 | //right wave | ||
| 367 | else | ||
| 368 | { | ||
| 369 | 65246971 | sol.velocityY = vr; | |
| 370 | |||
| 371 | //right shock | ||
| 372 |
2/2✓ Branch 0 taken 24401920 times.
✓ Branch 1 taken 40845051 times.
|
65246971 | if(ds >= dr) |
| 373 | { | ||
| 374 | 24401920 | const auto qr = sqrt((ds + dr)*ds /(2.0*dr*dr)); // flux right | |
| 375 | 24401920 | const auto sr = ur + cr*qr; // shock right | |
| 376 | |||
| 377 | //right side of shock | ||
| 378 |
2/2✓ Branch 0 taken 365624 times.
✓ Branch 1 taken 24036296 times.
|
24401920 | if(s >= sr) |
| 379 | { | ||
| 380 | 365624 | sol.waterDepth = dr; | |
| 381 | 365624 | sol.velocityX = ur; | |
| 382 | } | ||
| 383 | |||
| 384 | //left side of shock | ||
| 385 | else | ||
| 386 | { | ||
| 387 | 24036296 | sol.waterDepth = ds; | |
| 388 | 24036296 | sol.velocityX = us; | |
| 389 | } | ||
| 390 | |||
| 391 | //right rarefaction | ||
| 392 | } | ||
| 393 | |||
| 394 | else | ||
| 395 | { | ||
| 396 | 40845051 | const auto shr = ur + cr; //wave speed at head right | |
| 397 | |||
| 398 | //right side of Rarefaction | ||
| 399 |
2/2✓ Branch 0 taken 344199 times.
✓ Branch 1 taken 40500852 times.
|
40845051 | if (s >= shr) |
| 400 | { | ||
| 401 | 344199 | sol.waterDepth = dr; | |
| 402 | 344199 | sol.velocityX = ur; | |
| 403 | } | ||
| 404 | |||
| 405 | else | ||
| 406 | { | ||
| 407 | 40500852 | const auto str = us + cs; //wave speed at tail right | |
| 408 | |||
| 409 | //inside the rarefaction | ||
| 410 |
2/2✓ Branch 0 taken 190896 times.
✓ Branch 1 taken 40309956 times.
|
40500852 | if(s>=str) |
| 411 | { | ||
| 412 | 190896 | sol.velocityX = (ur - 2.0*cr + 2.0*s)/3.0; | |
| 413 | 190896 | const auto c = (-ur + 2.0*cr + s)/3.0; | |
| 414 | 190896 | sol.waterDepth = c*c/grav; | |
| 415 | |||
| 416 | } | ||
| 417 | |||
| 418 | //inside the star region | ||
| 419 | else | ||
| 420 | { | ||
| 421 | 40309956 | sol.waterDepth = ds; | |
| 422 | 40309956 | sol.velocityX = us; | |
| 423 | } | ||
| 424 | } | ||
| 425 | } | ||
| 426 | } | ||
| 427 | } | ||
| 428 | //============================== Flux computation =================== | ||
| 429 | |||
| 430 | 320754442 | sol.flux[0] = sol.waterDepth * sol.velocityX; | |
| 431 | 320754442 | sol.flux[1] = sol.waterDepth * sol.velocityX * sol.velocityX + 0.5 * grav * sol.waterDepth * sol.waterDepth; | |
| 432 | 320754442 | sol.flux[2] = sol.waterDepth * sol.velocityX * sol.velocityY; | |
| 433 | |||
| 434 | 160377221 | return sol; | |
| 435 | } | ||
| 436 | |||
| 437 | } // end namespace ShallowWater | ||
| 438 | } // end namespace Dumux | ||
| 439 | |||
| 440 | #endif | ||
| 441 |