Coverage for pyxfem/fluidmechanics/exact_riemann_solver.py: 0%
248 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1"""
2@author: chevaugeon-n
31D Exact Riemann Solver for the Euler Equations.
5This module provides a robust, high-performance implementation of the exact
6Riemann solver for the one-dimensional Euler equations of gas dynamics,
7based on the methodologies described in E.F. Toro's textbook.
9It is designed for High-Performance Computing (HPC) applications and features:
10 - A `PerfectGaz` thermodynamic closure model.
11 - Hardware-accelerated JIT-compiled critical paths using Numba (`@njit` and `@cfunc`),
12 allowing zero-overhead integration into broader finite volume schemes.
13 - An `ExactRiemannSolver` that computes the exact *star region* state using a
14 Newton-Raphson iterative scheme with robust initial guessing (PVRS, Two-Shock,
15 Two-Rarefaction).
16 - Full wave sampling across rarefaction fans, contact discontinuities, and shock waves.
17 - Built-in classical Riemann test cases (e.g., Sod shock tube).
19Dependencies:
20 - numpy
21 - numba (optional but strongly recommended for performance)
23Example:
24 from pyxfem.core.riemann import ExactRiemannSolver
25 solver = ExactRiemannSolver(rho0=1.0, u0=0.0, p0=1.0,
26 ... rho1=0.125, u1=0.0, p1=0.1)
27 x, rho, u, p = solver.solution_sample(x0=0.5, t=0.2, l=1.0, n=100)
29References:
30 Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics:
31 A Practical Introduction. Springer Science & Business Media.
32"""
34from typing import Tuple, Optional, Any
35import numpy as np
36from pyxfem.core.config import USE_NUMBA, FloatArray
38if USE_NUMBA :
39 import numba
40 @numba.njit(cache=True, fastmath=True)
41 def pg_soundspeed_numba(p: float, rho: float, gamma: float) -> float:
42 """
43 Numba-compiled function to compute the sound speed of a perfect gas.
45 Args:
46 p (float): Pressure.
47 rho (float): Density.
48 gamma (float): Heat capacity ratio (adiabatic index).
50 Returns:
51 float: Sound speed.
52 """
53 return np.sqrt(gamma*p/rho)
55 @numba.njit(cache=True, fastmath=True)
56 def pg_pressure_numba(e: float, rho: float, gamma: float) -> float:
57 """
58 Numba-compiled function to compute the pressure of a perfect gas.
60 Args:
61 e (float): Specific internal energy.
62 rho (float): Density.
63 gamma (float): Heat capacity ratio (adiabatic index).
65 Returns:
66 float: Pressure.
67 """
68 return (gamma-1.)*e*rho
70 nbfloat = numba.types.float64
71 c_sig = nbfloat(nbfloat, nbfloat, nbfloat)
73 @numba.cfunc(c_sig, nopython=True, cache=True, fastmath=True)
74 def pg_pressure_c(e: float, rho: float, gamma: float) -> float:
75 """
76 C-callable Numba function to compute the pressure of a perfect gas.
78 Args:
79 e (float): Specific internal energy.
80 rho (float): Density.
81 gamma (float): Heat capacity ratio (adiabatic index).
83 Returns:
84 float: Pressure.
85 """
86 return (gamma - 1.0) * e * rho
88 @numba.cfunc(c_sig, nopython=True, cache=True, fastmath=True)
89 def pg_soundspeed_c(p: float, rho: float, gamma: float) -> float:
90 """
91 C-callable Numba function to compute the sound speed of a perfect gas.
93 Args:
94 p (float): Pressure.
95 rho (float): Density.
96 gamma (float): Heat capacity ratio (adiabatic index).
98 Returns:
99 float: Sound speed.
100 """
101 return np.sqrt(gamma * p / rho)
103class PerfectGaz:
104 """
105 Thermodynamic model for a perfect gas.
107 Attributes:
108 gamma (float): Heat capacity ratio (adiabatic index).
109 """
110 def __init__(self, gamma: float = 1.4):
111 """
112 Initializes the PerfectGaz model.
114 Args:
115 gamma (float, optional): Heat capacity ratio. Defaults to 1.4.
116 """
117 self.gamma: float = gamma
119 def soundspeed(self, p: FloatArray, rho: FloatArray, out: Optional[FloatArray] = None)\
120 -> FloatArray:
121 """
122 Computes the sound speed as an array operation.
124 Args:
125 p (FloatVars): Pressure array or float
126 rho (FloatVars): Density array or float
127 out (Optional[FloatArray], optional): Output array to store results. Defaults to None.
129 Raises:
130 ValueError: If pressure or density lead to an invalid (NaN) sound speed calculation.
132 Returns:
133 FloatVars: Sound speed array (or float).
134 """
135 gamma = self.gamma
136 if out is not None:
137 cc = np.divide(p, rho, out=out)
138 cc *= gamma
139 cc = np.sqrt(cc, out=out)
140 else:
141 cc = np.sqrt(gamma * (p / rho))
142 if np.any(np.isnan(cc)):
143 print("Error : Invalid pressure or density in cells :")
144 ibad = np.argwhere(np.isnan(cc))
145 print(ibad)
146 print("pressure :", np.asarray(p)[ibad])
147 print("density :", np.asarray(rho)[ibad])
148 raise ValueError("invalid pressure or density when computing sound speed")
149 return cc
151 def soundspeed_scalar(self, p:float, rho: float) -> float:
152 """
153 Computes the sound speed as an array operation.
155 Args:
156 p (FloatVars): Pressure array or float
157 rho (FloatVars): Density array or float
158 out (Optional[FloatArray], optional): Output array to store results. Defaults to None.
160 Raises:
161 ValueError: If pressure or density lead to an invalid (NaN) sound speed calculation.
163 Returns:
164 FloatVars: Sound speed array (or float).
165 """
166 gamma = self.gamma
167 cc = np.sqrt(gamma * (p / rho))
168 # Vérification robuste (np.any fonctionne sur scalaires et tableaux)
169 if np.any(np.isnan(cc)):
170 print("Error : Invalid pressure or density in cells :")
171 print(f"pressure : {p}")
172 print(f"density : {rho}")
173 raise ValueError("invalid pressure or density when computing sound speed")
174 return cc
176 def internal_specific_energy(self, p: FloatArray, rho: FloatArray) -> FloatArray:
177 """
178 Returns the specific internal energy per unit of mass.
180 Args:
181 p (FloatArray): Pressure array.
182 rho (FloatArray): Density array.
184 Returns:
185 FloatArray: Specific internal energy array.
186 """
187 return p/rho/(self.gamma-1.)
189 def pressure(self, e: FloatArray, rho: FloatArray, out: Optional[FloatArray] = None)\
190 -> FloatArray:
191 """
192 Returns pressure as a function of internal energy and density.
194 Args:
195 e (FloatArray): Specific internal energy array.
196 rho (FloatArray): Density array.
197 out (Optional[FloatArray], optional): Output array to store results. Defaults to None.
199 Returns:
200 FloatArray: Pressure array.
201 """
202 p = np.multiply(e,rho, out = out)
203 p *=(self.gamma-1)
204 return p
206 if USE_NUMBA :
207 def get_numba_c_funcs(self) -> Tuple[Any, Any, float]:
208 """
209 Returns ctypes pointers for the Numba C functions and the gamma parameter.
211 Returns:
212 Tuple[ctypes._CFuncPtr, ctypes._CFuncPtr, float]:
213 Tuple containing ctypes pointers for pressure and sound speed functions, and gamma.
214 """
215 return pg_pressure_c.ctypes, pg_soundspeed_c.ctypes, self.gamma
217 # def make_numba_pressure_func(self):
218 # gamma = self.gamma
219 # @numba.njit(cache =True)
220 # def numba_pressure_func(e, rho):
221 # return (gamma-1.)*e*rho
222 # return numba_pressure_func
224 # def make_numba_soundspeed_func(self):
225 # gamma = self.gamma
226 # @numba.njit(cache =True)
227 # def numba_soundspeed_func(p, rho):
228 # return np.sqrt(gamma*p/rho)
229 # return numba_soundspeed_func
231 # def get_numba_funcs(self):
232 # # On retourne les pointeurs de fonction ET le tuple de paramètres
233 # return pg_pressure_numba, pg_soundspeed_numba, self.gamma
237class ExactRiemannSolver():
238 """ Implementation of the exact Riemann solver for the 1D Euler equations,
239 following E.F. Toro's textbook.
241 Attributes:
242 gamma (float): Heat capacity ratio.
243 gs (FloatArray): Array of constant gamma-dependent thermodynamic factors.
244 rho0 (float): Density of the left state.
245 rho1 (float): Density of the right state.
246 u0 (float): Velocity of the left state.
247 u1 (float): Velocity of the right state.
248 p0 (float): Pressure of the left state.
249 p1 (float): Pressure of the right state.
250 gaz (PerfectGaz): Perfect gas model instance.
251 c0 (float): Sound speed of the left state.
252 c1 (float): Sound speed of the right state.
253 pm (float): Resolved pressure in the star region.
254 um (float): Resolved velocity in the star region.
255 """
256 def __init__(self, gamma: float = 1.4, rho0: float = 1.0, rho1: float = 0.125,
257 u0: float = 0.0, u1: float = 0.0, p0: float = 1.0, p1: float = 0.1):
258 """
259 Initializes the exact Riemann solver and computes the star state (pm, um).
261 Args:
262 gamma (float, optional): Heat capacity ratio. Defaults to 1.4.
263 rho0 (float, optional): Left state density. Defaults to 1.0.
264 rho1 (float, optional): Right state density. Defaults to 0.125.
265 u0 (float, optional): Left state velocity. Defaults to 0.0.
266 u1 (float, optional): Right state velocity. Defaults to 0.0.
267 p0 (float, optional): Left state pressure. Defaults to 1.0.
268 p1 (float, optional): Right state pressure. Defaults to 0.1.
270 Raises:
271 ValueError: If initial data generates a vacuum.
272 """
273 self.gamma = gamma
274 g = self.gamma
275 # 0 1 2 3 4
276 self.gs = np.array([g, (g-1)/2./g, (g+1)/2./g, 2*g/(g-1), 2./(g-1), \
277 # 5 6 7 8
278 2./(g+1), (g-1)/(g+1), (g-1)/2, g-1])
279 self.rho0 = rho0
280 self.rho1 = rho1
281 self.u0 = u0
282 self.u1 = u1
283 self.p0 = p0
284 self.p1 = p1
285 self.gaz = PerfectGaz(gamma)
286 self.c0 = self.gaz.soundspeed_scalar(p0, rho0)
287 self.c1 = self.gaz.soundspeed_scalar(p1, rho1)
288 if 2./(gamma-1.)*(self.c0+self.c1) < (u1-u0) :
289 print("Generation of vaccum with initial data !")
290 raise ValueError()
291 self.pm, self.um = self.star_pu()
293 def guess_p(self) -> float:
294 """
295 Provides an initial guess for the pressure in the star region using
296 approximate Riemann solvers (PVRS, Two-Rarefaction, Two-Shock).
298 Returns:
299 float: Initial guess for pressure.
300 """
301 p0, p1 = self.p0, self.p1
302 u0, u1 = self.u0, self.u1
303 rho0, rho1 = self.rho0, self.rho1
304 c0, c1 = self.c0, self.c1
305 cup = 0.25*(rho0+rho1)*(c0+c1)
306 tmp = 0.5*(p0+p1) + 0.5*(u0-u1)*cup
307 ppv = max(tmp, 0.)
308 pmin = min(p0,p1)
309 pmax = max(p0,p1)
310 qmax = pmax/pmin
311 quser = 2.
312 if qmax <= quser and pmin <= ppv and ppv<= pmax :
313 # PVRS Riemann Solver
314 return ppv
315 if ppv <= pmin :
316 # two rarefaction rieman Solver
317 g1, g3, g4, g7 = self.gs[[1,3,4,7]]
318 pq = (p0/p1)**g1
319 um = (pq*u0/c0 + u1/c1 + g4*(pq-1.))/(pq/c0 + 1./c1)
320 pt0 = 1. + g7*(u0-um)/c0
321 pt1 = 1. + g7*(um-u1)/c1
322 return 0.5*(p0*pt0**g3 + p1*pt1**g3)
323 # two shock estimate
324 g5, g6 = self.gs[[5,6]]
325 ge0 = np.sqrt( (g5/rho0)/(g6*p0 + ppv) )
326 ge1 = np.sqrt( (g5/rho1)/(g6*p1 + ppv) )
327 return (ge0*p0+ge1*p1 - (u1 -u0))/(ge0+ge1)
329 def prefun(self, p: float, dk: float, pk: float, ck: float) -> Tuple[float, float]:
330 """
331 Computes the Riemann function and its derivative for a given state.
333 Args:
334 p (float): Current guessed pressure.
335 dk (float): State density (left or right).
336 pk (float): State pressure (left or right).
337 ck (float): State sound speed (left or right).
339 Returns:
340 Tuple[float, float]: Riemann function value (f) and its derivative (fd).
341 """
342 g1,g2,g4,g5,g6 = self.gs[[1,2,4, 5, 6]]
343 if p <= pk :
344 # rarefaction
345 rp = p/pk
346 f = g4*ck*(rp**g1 - 1.)
347 fd = (1./(dk*ck))*rp**(-g2)
348 else :
349 # shockwave
350 ak = g5/dk
351 bk = g6*pk
352 qrt = np.sqrt(ak/(bk + p))
353 f = (p-pk)*qrt
354 fd = (1.- 0.5*(p-pk)/(bk+p))*qrt
355 return f, fd
357 def star_pu(self) -> Tuple[float, float]:
358 """
359 Solves for pressure and velocity in the star region using Newton-Raphson iteration.
361 Raises:
362 ValueError: If the Newton scheme fails to converge within nitermax.
364 Returns:
365 Tuple[float, float]: Pressure (p) and velocity (u) in the star region.
366 """
367 p0, p1 = self.p0, self.p1
368 u0, u1 = self.u0, self.u1
369 rho0, rho1 = self.rho0, self.rho1
370 c0, c1 = self.c0, self.c1
371 pold = self.guess_p()
372 du = u1 - u0
373 nitermax = 100
374 tol = 1.e-6
375 for i in range(1, nitermax+1):
376 f0, f0d = self.prefun( pold, rho0, p0, c0)
377 f1, f1d = self.prefun( pold, rho1, p1, c1)
378 p = pold - (f0 + f1 + du)/( f0d+f1d)
379 change = 2.*np.abs( (p - pold)/(p+pold))
380 if change <= tol :
381 break
382 #print(i, p)
383 pold = p
384 if i >= nitermax:
385 print("divergence of newton scheme !")
386 raise ValueError()
387 u = 0.5*(u0+u1 + f1 -f0)
388 return p, u
390 def sample(self, pm: float, um: float, s: FloatArray)\
391 -> Tuple[FloatArray, FloatArray, FloatArray]:
392 """
393 Samples the complete Riemann solution across a range of wave speeds.
395 Args:
396 pm (float): Resolved star pressure.
397 um (float): Resolved star velocity.
398 s (FloatArray): Array of wave speeds (x/t) to sample.
400 Returns:
401 Tuple[FloatArray, FloatArray, FloatArray]:
402 Sampled density (rho), velocity (u),
403 and pressure (p) arrays corresponding to input speeds.
404 """
405 s = np.asarray(s)
406 p0, p1 = self.p0, self.p1
407 u0, u1 = self.u0, self.u1
408 rho0, rho1 = self.rho0, self.rho1
409 c0, c1 = self.c0, self.c1
410 g, g1, g2, g3, g4, g5, g6, g7, _ = self.gs[:]
411 rho = np.zeros(s.shape)
412 u = np.zeros(s.shape)
413 p = np.zeros(s.shape)
414 s_le_um = s<=um
415 #indices = np.where(s<=um)
416 # point left to the contact discontinuity
417 if pm <= p0:
418 # in left rarefaction
419 sh0 = u0 - c0
420 #lefts tate
421 s_le_sh0 = s <= sh0
422 indices = np.where(np.logical_and(s_le_um, s_le_sh0))
423 rho[indices] = rho0
424 u[indices] = u0
425 p[indices] = p0
426 cm0 = c0*(pm/p0)**g1
427 st0 = um - cm0
428 s_le_st0 = s <= st0
429 s_le_um_and_s_gt_sh0 = np.logical_and(s_le_um, np.logical_not(s_le_sh0))
430 indices = np.where(np.logical_and( s_le_um_and_s_gt_sh0, np.logical_not(s_le_st0)))
431 #point in star left
432 rho[indices] = rho0*(pm/p0)**(1./g)
433 u[indices] = um
434 p[indices] = pm
435 indices = np.where(np.logical_and(s_le_um_and_s_gt_sh0, s_le_st0))
436 # point inside left fan
437 c = g5*(c0+ g7*(u0 -s[indices]))
438 u[indices]= g5*(c0+ g7*u0 +s[indices])
439 rho[indices] = rho0*(c/c0)**g4
440 p[indices]= p0*(c/c0)**g3
441 #left shock
442 pm0 = pm/p0
443 s0 = u0 - c0*np.sqrt(g2*pm0+g1)
444 else : #pm > p0
445 pm0 = pm/p0
446 s0 = u0 - c0*np.sqrt(g2*pm0+g1)
447 indice = np.where(np.logical_and(s_le_um, s<=s0))
448 # point is left data state
449 rho[indice], u[indice], p[indice] = rho0, u0, p0
450 indice = np.where(np.logical_and(s_le_um, s>s0))
451 # point is star left state
452 rho[indice], u[indice], p[indice] = rho0*(pm0 + g6)/(pm0*g6 + 1.), um, pm
454 #indices where s>um
455 #sampig point is to the right of the contact discontinuity
456 if pm > p1 :
457 #right shock
458 pm1 = pm/p1
459 s1 = u1 + c1*np.sqrt(g2*pm1 + g1)
460 indice = np.where(np.logical_and(s > um, s >= s1))
461 rho[indice], u[indice], p[indice] = rho1, u1, p1
462 indice = np.where(np.logical_and(s > um, s < s1))
463 # sampled point is star right state
464 rho[indice], u[indice], p[indice] = rho1*(pm1 + g6)/(pm1*g6 + 1.), um, pm
465 else : #pm<= p1
466 #right rarefaction
467 sh1 = u1 +c1
468 indice = np.where( np.logical_and(s > um, s >= sh1))
469 # right state !
470 rho[indice], u[indice], p[indice] = rho1, u1, p1
471 cm1 = c1*(pm/p1)**g1
472 st1 = um + cm1
473 s_gt_um_and_s_lt_sh1 = np.logical_and(s > um, s < sh1)
474 indice = np.where(np.logical_and(s_gt_um_and_s_lt_sh1, s<= st1))
475 #point in star right
476 rho[indice], u[indice], p[indice]= rho1*(pm/p1)**(1./g), um, pm
477 indice = np.where(np.logical_and(s_gt_um_and_s_lt_sh1, s > st1))
478 # point inside right fan
479 c = g5*(c1 - g7*(u1 -s[indice]))
480 rho[indice] = rho1*(c/c1)**g4
481 u[indice] = g5*(-c1 + g7*u1 +s[indice])
482 p[indice] = p1*(c/c1)**g3
483 return rho, u, p
485 def solution_sample(self, x0: float, t: float, l: float, n: int) \
486 -> Tuple[FloatArray, FloatArray, FloatArray, FloatArray]:
487 """
488 Samples the complete solution profile over physical space at a specific time.
490 Args:
491 x0 (float): Initial position of the discontinuity.
492 t (float): Evaluation time.
493 l (float): Domain length (from 0 to l).
494 n (int): Number of spatial intervals.
496 Returns:
497 Tuple[FloatArray, FloatArray, FloatArray, FloatArray]:
498 Spatial grid (x), Density (rho), Velocity (u), and Pressure (p) arrays.
499 """
500 x = np.linspace(0,l,n+1)
501 if t != 0:
502 s = (x - x0)/t
503 rho = np.zeros(s.shape)
504 u = np.zeros(s.shape)
505 p = np.zeros(s.shape)
506 rho, u, p = self.sample(self.pm, self.um, s)
507 else:
508 p0, p1 = self.p0, self.p1
509 u0, u1 = self.u0, self.u1
510 rho0, rho1 = self.rho0, self.rho1
511 rho = np.where(x <= x0, rho0, rho1)
512 u = np.where(x <= x0, u0, u1)
513 p = np.where(x <= x0, p0, p1)
514 return x, rho, u, p
516def init_shock_tube(pbname: str = 'sod') \
517 -> Tuple[Tuple[float, float, float], Tuple[float, float, float], float]:
518 """
519 Returns initial conditions for classical 1D shock tube tests.
521 Args:
522 pbname (str, optional): Name of the problem ('sod', 'toro1', 'toro2', 'toro3').
523 Defaults to 'sod'.
525 Raises:
526 ValueError: If an unknown problem name is provided.
528 Returns:
529 Tuple[Tuple[float, float, float], Tuple[float, float, float], float]:
530 (rho_L, u_L, p_L), (rho_R, u_R, p_R), and end time (tend).
531 """
532 if pbname == "sod":
533 rhol, rhor = 1., 0.125
534 ul, ur = 0. , 0.
535 pl, pr = 1., 0.1
536 tend = 0.2
537 elif pbname == "toro1":
538 rhol, rhor = 1., 0.125
539 ul, ur = 0.75, 0.
540 pl, pr = 1., 0.1
541 tend = 0.2
542 elif pbname =="toro2":
543 rhol, rhor = 1., 1.
544 ul, ur = -2. , 2.
545 pl, pr = 0.4, 0.4
546 tend = 0.3
547 elif pbname =="toro3":
548 rhol, rhor = 1., 1.
549 ul, ur = 0. , 0.
550 pl, pr = 100., 0.1
551 tend = 0.03
552 else :
553 print(f"pb {pbname} unknown")
554 raise ValueError()
555 # elif pb =="toro4":
556 # rhol, rhor = 1., 1.
557 # ul, ur = 0. , 0.
558 # pl, pr = 1000., 0.01
559 return (rhol, ul, pl), (rhor, ur, pr), tend
561def test_Riemann_Sod(plot = False) -> None: #pylint: disable=invalid-name
562 """
563 Tests the exact Riemann solver on the standard Sod shock tube problem and plots results.
564 """
565 import matplotlib.pyplot as plt #pylint: disable=import-outside-toplevel
566 # revert sod :
567 rs = ExactRiemannSolver(rho0 = 0.125, p0 = 0.1, rho1 =1., p1 = 1.)
568 # sod :
569 rs = ExactRiemannSolver()
570 x, rho, u, p = rs.solution_sample(0.5, 0.1, 1., 50)
571 if plot :
572 plt.plot(x, rho)
573 plt.ylabel('rho')
574 plt.figure()
575 plt.plot(x, u )
576 plt.ylabel('u')
577 plt.figure()
578 plt.plot(x, p )
579 plt.ylabel('p')
580 plt.show()
582if __name__ == "__main__":
583 test_Riemann_Sod(plot = True)