3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
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 * 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 2 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_FLUX_SHALLOW_WATER_EXACT_RIEMANN_HH
25#define DUMUX_FLUX_SHALLOW_WATER_EXACT_RIEMANN_HH
26
27#include <array>
28#include <cmath>
29
30namespace Dumux {
31namespace ShallowWater{
32
33template<typename Scalar>
34struct RiemannSolution {
35 std::array<Scalar,3> flux;
36 Scalar waterDepth;
37 Scalar velocityX;
38 Scalar velocityY;
39};
40
41
60template<class Scalar>
61RiemannSolution<Scalar> exactRiemann(const Scalar dl,
62 const Scalar dr,
63 const Scalar ul,
64 const Scalar ur,
65 const Scalar vl,
66 const Scalar vr,
67 const Scalar grav,
68 const Scalar s = 0.0)
69{
70 RiemannSolution<Scalar> sol;
71 sol.waterDepth = 0.0; // d (Toro)
72 sol.velocityX = 0.0; // u (Toro)
73 sol.velocityY = 0.0; // v (Toro)
74
75 constexpr int maxsteps = 200; //maximum steps for the Newton method
76 constexpr Scalar tol = 1.0E-12; //tolerance for the Newton method
77
78 using std::sqrt;
79 using std::max;
80 using std::min;
81 using std::abs;
82
83 //================== Exact Riemann solver
84 const auto cl = sqrt(grav * max(dl, 0.0)); //celerity sqrt(g*h) left side
85 const auto cr = sqrt(grav * max(dr, 0.0)); //celerity sqrt(g*h) right side
86 const auto dcrit = (ur - ul) - 2.0*(cl + cr); //critical water depth
87
88 // dry case
89 if ((dl <= 0) || (dr <= 0) || (dcrit >= 0))
90 {
91 //left side is dry
92 if(dl <= 0)
93 {
94 const auto shr = ur + cr; //wave speed at head right
95
96 if (s >= shr)
97 {
98 sol.waterDepth = dr;
99 sol.velocityX = ur;
100 sol.velocityY = vr;
101
102 }
103
104 else
105 {
106
107 const auto str = ur - 2.0*cr; //wave speed at tail right
108
109 if(s >= str){
110 sol.velocityX = (ur - 2.0 *cr + 2.0*s)/3.0;
111 const auto c = (-ur + 2.0*cr +s )/3.0;
112 sol.waterDepth = c*c/grav;
113 sol.velocityY = vr;
114
115 }
116
117 else
118 {
119 sol.waterDepth = dl;
120 sol.velocityX = ul;
121 sol.velocityY = vl;
122 }
123 }
124 }
125
126 else
127 {
128 //right side is dry
129 if(dr <= 0)
130 {
131 const auto shl = ul-cl; //wave speed at head left
132
133 if(s <= shl)
134 {
135 sol.waterDepth = dl;
136 sol.velocityX = ul;
137 sol.velocityY = vl;
138 }
139
140 else
141 {
142 const auto stl = ul + 2.0 *cl; //wave speed at tail left
143
144 if ( s <= stl)
145 {
146 sol.velocityX = (ul + 2.0 *cl + 2.0*s)/3.0;
147 const auto c = (ul + 2.0*cl - s)/3.0;
148 sol.waterDepth = c*c/grav;
149 sol.velocityY = vl;
150 }
151
152 else
153 {
154 sol.waterDepth = dr;
155 sol.velocityX = ur;
156 sol.velocityY = vr;
157 }
158 }
159 }
160
161 //middle state is dry
162 else
163 {
164 const auto shl = ul - cl; //wave speed at head left
165 const auto shr = ur + cr; //wave speed at head right
166 const auto ssl = ul + 2.0 * cl;
167 const auto ssr = ur - 2.0 * cr;
168
169 if(s <= shl)
170 {
171 sol.waterDepth = dl;
172 sol.velocityX = ul;
173 sol.velocityY = vl;
174 }
175
176 if((s > shl) && (s <= ssl))
177 {
178 sol.velocityX = (ul + 2.0 * cl + 2.0*s)/3.0;
179 const auto c = (ul + 2.0*cl - s)/3.0;
180 sol.waterDepth = c*c/grav;
181 sol.velocityY = vl;
182 }
183
184 if((s > ssl) && (s <= ssr))
185 {
186 sol.waterDepth = 0.0;
187 sol.velocityX = 0.0;
188 sol.velocityY = 0.0;
189 }
190
191 if ((s > ssr ) && (s <= shr))
192 {
193 sol.velocityX = (ur - 2.0*cr + 2.0*s)/3.0;
194 const auto c = (-ur + 2.0*cr + s)/3.0;
195 sol.waterDepth = c*c/grav;
196 sol.velocityY = vr;
197 }
198
199 if(s > shr)
200 {
201 sol.waterDepth = dr;
202 sol.velocityX = ur;
203 sol.velocityY = vr;
204 }
205 }
206 }
207 }
208
209 // wet case
210 else
211 {
212 // Get a starting value for the Newton-Raphson method
213 const auto tmp = 0.5 * (cl + cr) - 0.25 * (ul - ur);
214 auto ds = (1.0/grav) * tmp * tmp; // water depth
215
216 //default is Two-Rarefaction as starting value, if ds > dmin we use two-shock
217 const auto dmin = min(dl, dr); // minimum depth
218
219 if (ds > dmin)
220 {
221 const auto gel = sqrt(0.5*grav * (ds+dl)/(ds*dl));
222 const auto ger = sqrt(0.5*grav * (ds+dr)/(ds*dr));
223 ds = (gel * dl + ger * dr - (ur-ul))/(gel + ger);
224 }
225
226 //Start Newton-Raphson loop ds = hstar
227 ds = max(ds, tol);
228 auto d0 = ds;
229
230 Scalar fl = 0.0, fr = 0.0, fld = 0.0, frd = 0.0;
231
232 for (int i=0; i <= maxsteps; ++i)
233 {
234 //Compute fl,fr,fld,frd
235
236 //first left side
237 if (ds <= dl)
238 {
239 const auto c = sqrt(grav * ds);
240 fl = 2.0 * (c-cl);
241 fld = grav/c;
242 }
243
244 else
245 {
246 const auto ges = sqrt(0.5 * grav * (ds+dl)/(ds*dl));
247 fl = (ds - dl) * ges;
248 fld = ges - 0.25 * grav * (ds - dl)/(ges * ds * ds);
249 }
250
251 //second right side
252 if (ds <= dr)
253 {
254 const auto c = sqrt(grav * ds);
255 fr = 2.0 * (c-cr);
256 frd = grav/c;
257 }
258
259 else
260 {
261 const auto ges = sqrt(0.5 * grav * (ds+dr)/(ds*dr));
262 fr = (ds - dr) * ges;
263 frd = ges - 0.25 * grav * (ds - dr)/(ges * ds * ds);
264 }
265
266 ds -= (fl + fr + ur - ul)/(fld + frd);
267 const auto cha = abs(ds - d0)/(0.5*(ds + d0));
268
269 if (cha <= tol)
270 {
271 break;
272 }
273
274 if (ds < 0.0)
275 {
276 ds = tol;
277 }
278
279 d0 = ds;
280 }
281
282 //compute ustar (us)
283 if (ds <= dl)
284 {
285 const auto c = sqrt(grav * ds);
286 fl = 2.0 * (c-cl);
287 }
288
289 else
290 {
291 const auto ges = sqrt(0.5*grav*(ds + dl)/(ds*dl));
292 fl = (ds - dl) * ges;
293 }
294
295 if (ds <= dr)
296 {
297 const auto c = sqrt(grav * ds);
298 fr = 2.0 * (c-cr);
299 }
300
301 else
302 {
303 const auto ges = sqrt(0.5*grav*(ds + dr)/(ds*dr));
304 fr = (ds - dr) * ges;
305 }
306
307 //exact Riemann solver sstar = ustar
308 const auto us = 0.5*(ul + ur) + 0.5*(fr - fl); // u_star
309 const auto cs = sqrt(grav*ds); // c_star
310
311 /*********************** computation of u and d *******************/
312 //left wave
313 if (s <= us)
314 {
315 sol.velocityY = vl;
316 //left shock
317 if (ds >= dl)
318 {
319 const auto ql = sqrt((ds + dl)*ds/(2.0*dl*dl)); // flux left
320 const auto sl = ul - cl*ql; // shock left
321
322 //left side of shock
323 if(s <= sl)
324 {
325 sol.waterDepth = dl;
326 sol.velocityX = ul;
327 }
328
329 //right side of shock
330 else
331 {
332 sol.waterDepth = ds;
333 sol.velocityX = us;
334 }
335 }
336
337 //left rarefaction
338 else
339 {
340 const auto shl = ul-cl; //wave speed at head left
341
342 //right side of rarefaction
343 if (s <= shl)
344 {
345 sol.waterDepth = dl;
346 sol.velocityX = ul;
347 }
348
349 else
350 {
351 const auto stl = us -cs; //wave speed at tail left
352
353 //inside the rarefaction
354 if(s <= stl)
355 {
356 sol.velocityX = (ul + 2.0*cl + 2.0 * s)/3.0;
357 const auto c = (ul + 2.0* cl - s)/3.0;
358 sol.waterDepth = c*c/grav;
359 }
360
361 //inside the star region
362 else
363 {
364 sol.waterDepth = ds;
365 sol.velocityX = us;
366 }
367 }
368 }
369 }
370
371 //right wave
372 else
373 {
374 sol.velocityY = vr;
375
376 //right shock
377 if(ds >= dr)
378 {
379 const auto qr = sqrt((ds + dr)*ds /(2.0*dr*dr)); // flux right
380 const auto sr = ur + cr*qr; // shock right
381
382 //right side of shock
383 if(s >= sr)
384 {
385 sol.waterDepth = dr;
386 sol.velocityX = ur;
387 }
388
389 //left side of shock
390 else
391 {
392 sol.waterDepth = ds;
393 sol.velocityX = us;
394 }
395
396 //right rarefaction
397 }
398
399 else
400 {
401 const auto shr = ur + cr; //wave speed at head right
402
403 //right side of Rarefaction
404 if (s >= shr)
405 {
406 sol.waterDepth = dr;
407 sol.velocityX = ur;
408 }
409
410 else
411 {
412 const auto str = us + cs; //wave speed at tail right
413
414 //inside the rarefaction
415 if(s>=str)
416 {
417 sol.velocityX = (ur - 2.0*cr + 2.0*s)/3.0;
418 const auto c = (-ur + 2.0*cr + s)/3.0;
419 sol.waterDepth = c*c/grav;
420
421 }
422
423 //inside the star region
424 else
425 {
426 sol.waterDepth = ds;
427 sol.velocityX = us;
428 }
429 }
430 }
431 }
432 }
433 //============================== Flux computation ===================
434
435 sol.flux[0] = sol.waterDepth * sol.velocityX;
436 sol.flux[1] = sol.waterDepth * sol.velocityX * sol.velocityX + 0.5 * grav * sol.waterDepth * sol.waterDepth;
437 sol.flux[2] = sol.waterDepth * sol.velocityX * sol.velocityY;
438
439 return sol;
440}
441
442} // end namespace ShallowWater
443} // end namespace Dumux
444
445#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 Shallow water equations.
Definition: exactriemann.hh:61
Dumux
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Dumux::ShallowWater::RiemannSolution
Definition: exactriemann.hh:34
Dumux::ShallowWater::RiemannSolution::waterDepth
Scalar waterDepth
Definition: exactriemann.hh:36
Dumux::ShallowWater::RiemannSolution::flux
std::array< Scalar, 3 > flux
Definition: exactriemann.hh:35
Dumux::ShallowWater::RiemannSolution::velocityY
Scalar velocityY
Definition: exactriemann.hh:38
Dumux::ShallowWater::RiemannSolution::velocityX
Scalar velocityX
Definition: exactriemann.hh:37
Generated by doxygen 1.9.3