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