version 3.10-dev
exactriemann.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// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
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
18namespace Dumux {
19namespace ShallowWater{
20
21template<typename Scalar>
22struct RiemannSolution {
23 std::array<Scalar,3> flux;
24 Scalar waterDepth;
25 Scalar velocityX;
26 Scalar velocityY;
27};
28
29
53template<class Scalar>
54RiemannSolution<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 sol.waterDepth = 0.0; // d (Toro)
65 sol.velocityX = 0.0; // u (Toro)
66 sol.velocityY = 0.0; // v (Toro)
67
68 constexpr int maxsteps = 200; //maximum steps for the Newton method
69 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 const auto cl = sqrt(grav * max(dl, 0.0)); //celerity sqrt(g*h) left side
78 const auto cr = sqrt(grav * max(dr, 0.0)); //celerity sqrt(g*h) right side
79 const auto dcrit = (ur - ul) - 2.0*(cl + cr); //critical water depth
80
81 // dry case
82 if ((dl <= 0) || (dr <= 0) || (dcrit >= 0))
83 {
84 //left side is dry
85 if(dl <= 0)
86 {
87 const auto shr = ur + cr; //wave speed at head right
88
89 if (s >= shr)
90 {
91 sol.waterDepth = dr;
92 sol.velocityX = ur;
93 sol.velocityY = vr;
94
95 }
96
97 else
98 {
99
100 const auto str = ur - 2.0*cr; //wave speed at tail right
101
102 if(s >= str){
103 sol.velocityX = (ur - 2.0 *cr + 2.0*s)/3.0;
104 const auto c = (-ur + 2.0*cr +s )/3.0;
105 sol.waterDepth = c*c/grav;
106 sol.velocityY = vr;
107
108 }
109
110 else
111 {
112 sol.waterDepth = dl;
113 sol.velocityX = ul;
114 sol.velocityY = vl;
115 }
116 }
117 }
118
119 else
120 {
121 //right side is dry
122 if(dr <= 0)
123 {
124 const auto shl = ul-cl; //wave speed at head left
125
126 if(s <= shl)
127 {
128 sol.waterDepth = dl;
129 sol.velocityX = ul;
130 sol.velocityY = vl;
131 }
132
133 else
134 {
135 const auto stl = ul + 2.0 *cl; //wave speed at tail left
136
137 if ( s <= stl)
138 {
139 sol.velocityX = (ul + 2.0 *cl + 2.0*s)/3.0;
140 const auto c = (ul + 2.0*cl - s)/3.0;
141 sol.waterDepth = c*c/grav;
142 sol.velocityY = vl;
143 }
144
145 else
146 {
147 sol.waterDepth = dr;
148 sol.velocityX = ur;
149 sol.velocityY = vr;
150 }
151 }
152 }
153
154 //middle state is dry
155 else
156 {
157 const auto shl = ul - cl; //wave speed at head left
158 const auto shr = ur + cr; //wave speed at head right
159 const auto ssl = ul + 2.0 * cl; // wave speed of a wet/dry front for a left dry state
160 const auto ssr = ur - 2.0 * cr; // wave speed of a wet/dry front for a right dry state
161
162 if(s <= shl)
163 {
164 sol.waterDepth = dl;
165 sol.velocityX = ul;
166 sol.velocityY = vl;
167 }
168
169 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 if((s > ssl) && (s <= ssr))
178 {
179 sol.waterDepth = 0.0;
180 sol.velocityX = 0.0;
181 sol.velocityY = 0.0;
182 }
183
184 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 if(s > shr)
193 {
194 sol.waterDepth = dr;
195 sol.velocityX = ur;
196 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 const auto tmp = 0.5 * (cl + cr) - 0.25 * (ul - ur);
207 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 const auto dmin = min(dl, dr); // minimum depth
211
212 if (ds > dmin)
213 {
214 // function g in (Toro, 2001)
215 const auto gel = sqrt(0.5*grav * (ds+dl)/(ds*dl));
216 const auto ger = sqrt(0.5*grav * (ds+dr)/(ds*dr));
217 ds = (gel * dl + ger * dr - (ur-ul))/(gel + ger);
218 }
219
220 //Start Newton-Raphson loop ds = hstar
221 ds = max(ds, tol);
222 auto d0 = ds; // initial guess water depth
223
224 // helper functions (left and right) and their derivatives
225 Scalar fl = 0.0, fr = 0.0, fld = 0.0, frd = 0.0;
226
227 for (int i=0; i <= maxsteps; ++i)
228 {
229 //Compute fl,fr,fld,frd
230
231 //first left side
232 if (ds <= dl) // wave is rarefaction wave (or depression)
233 {
234 const auto c = sqrt(grav * ds);
235 fl = 2.0 * (c-cl);
236 fld = grav/c;
237 }
238
239 else // wave is shock wave (or bore)
240 {
241 const auto ges = sqrt(0.5 * grav * (ds+dl)/(ds*dl)); // function g in (Toro, 2001)
242 fl = (ds - dl) * ges;
243 fld = ges - 0.25 * grav * (ds - dl)/(ges * ds * ds);
244 }
245
246 //second right side
247 if (ds <= dr)
248 {
249 const auto c = sqrt(grav * ds);
250 fr = 2.0 * (c-cr);
251 frd = grav/c;
252 }
253
254 else
255 {
256 const auto ges = sqrt(0.5 * grav * (ds+dr)/(ds*dr));
257 fr = (ds - dr) * ges;
258 frd = ges - 0.25 * grav * (ds - dr)/(ges * ds * ds);
259 }
260
261 ds -= (fl + fr + ur - ul)/(fld + frd);
262 const auto cha = abs(ds - d0)/(0.5*(ds + d0));
263
264 if (cha <= tol)
265 {
266 break;
267 }
268
269 if (ds < 0.0)
270 {
271 ds = tol;
272 }
273
274 d0 = ds;
275 }
276
277 //compute ustar (us)
278 if (ds <= dl)
279 {
280 const auto c = sqrt(grav * ds);
281 fl = 2.0 * (c-cl);
282 }
283
284 else
285 {
286 const auto ges = sqrt(0.5*grav*(ds + dl)/(ds*dl));
287 fl = (ds - dl) * ges;
288 }
289
290 if (ds <= dr)
291 {
292 const auto c = sqrt(grav * ds);
293 fr = 2.0 * (c-cr);
294 }
295
296 else
297 {
298 const auto ges = sqrt(0.5*grav*(ds + dr)/(ds*dr));
299 fr = (ds - dr) * ges;
300 }
301 //exact Riemann solver sstar = ustar
302 const auto us = 0.5*(ul + ur) + 0.5*(fr - fl); // velocity in the star region (u_star)
303 const auto cs = sqrt(grav*ds); // wave speed in the star region
304
305
306 /*********************** computation of u and d *******************/
307 //left wave
308 if (s <= us)
309 {
310 sol.velocityY = vl;
311 //left shock
312 if (ds >= dl)
313 {
314 const auto ql = sqrt((ds + dl)*ds/(2.0*dl*dl)); // flux left
315 const auto sl = ul - cl*ql; // shock left
316
317 //left side of shock
318 if(s <= sl)
319 {
320 sol.waterDepth = dl;
321 sol.velocityX = ul;
322 }
323
324 //right side of shock
325 else
326 {
327 sol.waterDepth = ds;
328 sol.velocityX = us;
329 }
330 }
331
332 //left rarefaction
333 else
334 {
335 const auto shl = ul-cl; //wave speed at head left
336
337 //right side of rarefaction
338 if (s <= shl)
339 {
340 sol.waterDepth = dl;
341 sol.velocityX = ul;
342 }
343
344 else
345 {
346 const auto stl = us -cs; //wave speed at tail left
347
348 //inside the rarefaction
349 if(s <= stl)
350 {
351 sol.velocityX = (ul + 2.0*cl + 2.0 * s)/3.0;
352 const auto c = (ul + 2.0* cl - s)/3.0;
353 sol.waterDepth = c*c/grav;
354 }
355
356 //inside the star region
357 else
358 {
359 sol.waterDepth = ds;
360 sol.velocityX = us;
361 }
362 }
363 }
364 }
365
366 //right wave
367 else
368 {
369 sol.velocityY = vr;
370
371 //right shock
372 if(ds >= dr)
373 {
374 const auto qr = sqrt((ds + dr)*ds /(2.0*dr*dr)); // flux right
375 const auto sr = ur + cr*qr; // shock right
376
377 //right side of shock
378 if(s >= sr)
379 {
380 sol.waterDepth = dr;
381 sol.velocityX = ur;
382 }
383
384 //left side of shock
385 else
386 {
387 sol.waterDepth = ds;
388 sol.velocityX = us;
389 }
390
391 //right rarefaction
392 }
393
394 else
395 {
396 const auto shr = ur + cr; //wave speed at head right
397
398 //right side of Rarefaction
399 if (s >= shr)
400 {
401 sol.waterDepth = dr;
402 sol.velocityX = ur;
403 }
404
405 else
406 {
407 const auto str = us + cs; //wave speed at tail right
408
409 //inside the rarefaction
410 if(s>=str)
411 {
412 sol.velocityX = (ur - 2.0*cr + 2.0*s)/3.0;
413 const auto c = (-ur + 2.0*cr + s)/3.0;
414 sol.waterDepth = c*c/grav;
415
416 }
417
418 //inside the star region
419 else
420 {
421 sol.waterDepth = ds;
422 sol.velocityX = us;
423 }
424 }
425 }
426 }
427 }
428 //============================== Flux computation ===================
429
430 sol.flux[0] = sol.waterDepth * sol.velocityX;
431 sol.flux[1] = sol.waterDepth * sol.velocityX * sol.velocityX + 0.5 * grav * sol.waterDepth * sol.waterDepth;
432 sol.flux[2] = sol.waterDepth * sol.velocityX * sol.velocityY;
433
434 return sol;
435}
436
437} // end namespace ShallowWater
438} // end namespace Dumux
439
440#endif
Dumux::ShallowWater::exactRiemann
RiemannSolution< Scalar > exactRiemann(const Scalar dl, const Scalar dr, const Scalar ul, const Scalar ur, const Scalar vl, const Scalar vr, const Scalar grav, const Scalar s=0.0)
Exact Riemann solver for the shallow water equations.
Definition: exactriemann.hh:54
Dumux
Definition: adapt.hh:17
Dumux::ShallowWater::RiemannSolution
Definition: exactriemann.hh:22
Dumux::ShallowWater::RiemannSolution::waterDepth
Scalar waterDepth
Definition: exactriemann.hh:24
Dumux::ShallowWater::RiemannSolution::flux
std::array< Scalar, 3 > flux
Definition: exactriemann.hh:23
Dumux::ShallowWater::RiemannSolution::velocityY
Scalar velocityY
Definition: exactriemann.hh:26
Dumux::ShallowWater::RiemannSolution::velocityX
Scalar velocityX
Definition: exactriemann.hh:25