Prérequis : Chapitre 3 — Méthode des volumes finis (semi-discrétisation §3.3, méthode des lignes, condition CFL), Chapitre 5 — Reconstructions d'ordre élevé (ordres spatiaux 2 à 5).
Au chapitre 3 (§3.3), la semi-discrétisation spatiale a transformé l'EDP d'Euler en un système d'équations différentielles ordinaires (EDO) :
L'opérateur
Cette approche, appelée méthode des lignes (method of lines), découple la discrétisation spatiale de la discrétisation temporelle. Il reste à choisir un intégrateur en temps pour faire avancer la solution de
Implémentation : euler1d/solver.py — fonctions _compute_rhs (calcul de $\Delta t \cdot \mathbf{L}(\mathbf{U})$) et run_simulation (boucle en temps).
La condition de Courant-Friedrichs-Lewy (CFL) impose une borne supérieure au pas de temps pour garantir la stabilité des schémas explicites :
où
Le dénominateur
Dans le solveur, le pas de temps est recalculé à chaque itération :
- Calcul de
$s_{\max} = \max_i(\lvert u_i \rvert + c_i)$ viamax_wave_speed(U, gamma). -
$\Delta t = \text{CFL} \cdot \Delta x / s_{\max}$ . - Ajustement final : si
$t + \Delta t > t_{\text{final}}$ , on réduit$\Delta t$ pour atteindre exactement$t_{\text{final}}$ .
Ce recalcul est essentiel car la vitesse maximale évolue au cours de la simulation (notamment lors du passage d'un choc).
La méthode la plus simple pour intégrer le système d'EDO est l'Euler explicite :
C'est une méthode à un seul étage, d'ordre 1 en temps.
L'Euler explicite est le choix par défaut pour :
- Les schémas spatiaux d'ordre 1 (Rusanov, HLL, HLLC, Roe, Godunov, Lax-Friedrichs) : la précision temporelle d'ordre 1 est cohérente avec l'ordre spatial.
- Lax-Wendroff : ce schéma centré inclut sa propre intégration temporelle via un prédicteur-correcteur de Richtmyer. L'intégrateur extérieur doit être RK1 pour éviter un double pas de temps. Le solveur force
RK1lorsquetime_integral_included = True.
La méthode de Heun (ou RK2 explicite) utilise deux étages :
C'est une méthode à deux étages, d'ordre 2 en temps.
RK2 est le choix par défaut pour les reconstructions d'ordre spatial 2 :
- MUSCL avec limiteur de pente (minmod, van Leer, superbee, MC)
- ENO2 (sélection de stencil)
L'ordre temporel 2 correspond à l'ordre spatial 2, assurant un schéma globalement d'ordre 2 en espace et en temps.
Les schémas TVD (Total Variation Diminishing) garantissent que la variation totale de la solution ne croît pas au cours du temps, empêchant ainsi l'apparition d'oscillations parasites. Cependant, cette propriété est établie pour l'Euler explicite (RK1) sous la condition CFL
Les méthodes SSP (Strong Stability Preserving) sont conçues pour répondre à ce problème : si l'Euler explicite préserve une propriété de stabilité (TVD, positivité, etc.) sous une certaine condition CFL, alors la méthode SSP la préserve aussi, éventuellement sous une condition CFL réduite.
La méthode SSP-RK3, introduite par Shu et Osher [Shu, Osher, 1988], s'écrit sous forme de combinaisons convexes d'étapes d'Euler explicite :
C'est une méthode à trois étages, d'ordre 3 en temps.
Chaque étage est une combinaison convexe (coefficients positifs sommant à 1) d'un pas d'Euler explicite. Si l'Euler explicite préserve une propriété sous CFL
SSP-RK3 est le choix par défaut pour WENO3 (ordres JS et Z). L'ordre temporel 3 correspond à l'ordre spatial 3, et la propriété SSP est particulièrement importante pour les reconstructions WENO qui reposent sur la préservation de la TVD.
La méthode de Runge-Kutta d'ordre 4 classique utilise quatre étages :
C'est une méthode à quatre étages, d'ordre 4 en temps.
- Pas SSP : RK4 classique n'est pas une méthode SSP. Il n'existe pas de méthode RK explicite SSP d'ordre 4 avec seulement 4 étages [Gottlieb, Shu, Tadmor, 2001]. En pratique, pour les schémas spatiaux d'ordre 4 (JST), la dissipation artificielle adaptative stabilise la solution sans nécessiter la propriété SSP.
- Large domaine de stabilité : la région de stabilité de RK4 dans le plan complexe est sensiblement plus grande que celles de RK1-RK3, ce qui peut autoriser des CFL légèrement plus élevés dans certaines configurations.
RK4 est le choix par défaut pour le schéma JST (Jameson-Schmidt-Turkel), dont la dissipation artificielle adaptative requiert un intégrateur d'ordre suffisant.
Pour les reconstructions spatiales d'ordre 5 (WENO5), un intégrateur temporel d'ordre au moins 5 est nécessaire pour que l'erreur temporelle ne domine pas l'erreur spatiale. La méthode de Dormand-Prince [Dormand, Prince, 1980] est un schéma de Runge-Kutta à 6 étages et d'ordre 5.
Les six étages sont définis par :
La solution avancée est :
On note que
Dormand-Prince est le choix par défaut pour WENO5 (variantes JS et Z), assurant que l'erreur temporelle d'ordre 5 est cohérente avec la précision spatiale d'ordre 5.
Remarque : la méthode de Dormand-Prince est aussi la base du solveur adaptatif
ode45de MATLAB etdopri5de SciPy. Dans notre solveur, le pas de temps est contrôlé par la condition CFL et non par un estimateur d'erreur embarqué (embedded method), donc seule la formule d'ordre 5 est utilisée.
L'ordre de convergence global d'un schéma de volumes finis est limité par le minimum de l'ordre spatial et de l'ordre temporel :
Utiliser un intégrateur RK1 (ordre 1 en temps) avec une reconstruction WENO5 (ordre 5 en espace) limiterait l'ordre global à 1 : toute la précision spatiale serait gaspillée. Inversement, utiliser RK5 avec un flux d'ordre 1 est un surcoût inutile (6 évaluations de flux par pas de temps au lieu d'une).
| Ordre spatial | Reconstruction | Intégrateur par défaut | Étages |
|---|---|---|---|
| 1 | Constante (Rusanov, HLL, HLLC, Roe, Godunov, Lax-Friedrichs) | RK1 (Euler explicite) | 1 |
| 2 | MUSCL, ENO2 | RK2 (Heun) | 2 |
| 3 | WENO3 (JS, Z) | RK3 (SSP-RK3) | 3 |
| 4 | JST | RK4 (classique) | 4 |
| 5 | WENO5 (JS, Z) | RK5 (Dormand-Prince) | 6 |
La propriété default_time_integrator de chaque schéma renvoie la méthode RK correspondant à son ordre spatial. Le solveur (run_simulation) applique la logique suivante :
- Si
time_integral_included = True(Lax-Wendroff) : forcer RK1 pour éviter un double pas de temps. - Si l'utilisateur spécifie un intégrateur (
time_integratornon nul) : utiliser celui-ci. - Sinon : utiliser
scheme.default_time_integrator.
L'utilisateur peut toujours forcer un intégrateur différent (par exemple RK3 avec MUSCL pour bénéficier de la propriété SSP), mais la correspondance par défaut assure la consistance espace-temps.
| Méthode | Étages | Ordre | SSP ? | Schémas spatiaux par défaut |
|---|---|---|---|---|
| RK1 (Euler explicite) | 1 | 1 | Oui (trivial) | Rusanov, HLL, HLLC, Roe, Godunov, Lax-Friedrichs, Lax-Wendroff |
| RK2 (Heun) | 2 | 2 | Non | MUSCL, ENO2 |
| RK3 (SSP-RK3, Shu-Osher) | 3 | 3 |
Oui ( |
WENO3 (JS, Z) |
| RK4 (classique) | 4 | 4 | Non | JST |
| RK5 (Dormand-Prince) | 6 | 5 | Non | WENO5 (JS, Z) |
Coût par pas de temps : chaque étage nécessite une évaluation complète de l'opérateur spatial
- [Shu, Osher, 1988] C.-W. Shu et S. Osher, Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes, Journal of Computational Physics, 77(2), pp. 439-471, 1988.
- [Dormand, Prince, 1980] J. R. Dormand et P. J. Prince, A Family of Embedded Runge-Kutta Formulae, Journal of Computational and Applied Mathematics, 6(1), pp. 19-26, 1980.
- [Toro, 2009] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3e édition, Springer, 2009, Chapitre 6.
- [Gottlieb, Shu, Tadmor, 2001] S. Gottlieb, C.-W. Shu et E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review, 43(1), pp. 89-112, 2001.