Séminaire Calcul Scientifique et Modélisation

We discuss a new swarm-based gradient descent (SBGD) method for non-convex optimization. The swarm consists of agents, each is identified with position $x$ and mass $m$. There are three key aspects to the SBGD dynamics: (i) persistent transition of mass from agents at high to lower ground; (ii) a random marching direction, aligned with the steepest gradient descent; and (iii) a time stepping protocol which decreases with $m$.

The interplay between positions and masses leads to dynamic distinction between `heavier leaders’ near local minima, and `lighter explorers’ which explore for improved position with large(r) time steps. Convergence analysis and numerical simulations demonstrate the effectiveness of SBGD method as a global optimizer. 

The coupling of coastal wave models, such as Boussinesq-type (BT) and Saint-Venant (SV) equations, has been explored since the 1990s. Despite numerous models and coupling examples, the literature exhibits significant disagreement regarding induced artifacts and methods for their analysis. This work aims to elucidate these issues, proposing explanations and a method for evaluating and comparing coupling techniques. We ground our explanation in the mathematical properties of each model's Cauchy and half-line problems, highlighting the sensitivity of these models to numerical artifacts. Additionally, we demonstrate how one-way models provide insights into expected physical effects, unexpected artifacts, and errors relative to 3D models. We demonstrate this analysis with linearized models, where we establish the well-posedness of a popular coupling, characterize analytically the "coupling error" in terms of wave reflections, and prove its asymptotic behavior in shallow water. We will discuss how these insights can be applied to other linear/nonlinear models, providing a foundation for the evaluation and comparison of new coupled coastal wave models.

Simuler numériquement de manière précise l'évolution des interfaces séparant différents milieux est un enjeu crucial dans de nombreuses applications (multi-fluides, fluide-structure, etc). La méthode MOF (moment-of-fluid), extension de la méthode VOF (volume-of-fluid), utilise une reconstruction affine des interfaces par cellule basée sur les fractions volumiques et les centroïdes de chaque phase. Cette reconstruction d'interface est solution d'un problème de minimisation sous contrainte de volume. Ce problème est résolu dans la littérature par des calculs géométriques sur des polyèdres qui ont un coût important en 3D. On propose dans cet exposé une nouvelle approche du calcul de la fonction objectif et de ses dérivées de manière complètement analytique dans le cas de cellules hexaédriques rectangulaires et tétraédriques en 3D. Les résultats numériques montrent un gain important en temps de calcul.

We present a new model for heat transfer in compressible fluid flows. The model is derived from Hamilton’s principle of stationary action in Eulerian coordinates, in a setting where the entropy conservation is recovered as an Euler–Lagrange equation. A sufficient criterion for the hyperbolicity of the model is formulated. The governing equations are asymptotically consistent with the Euler equations for compressible heat conducting fluids, provided the addition of suitable relaxation terms. A study of the Rankine–Hugoniot conditions and Clausius–Duhem inequality is performed for a specific choice of the equation of state. In particular, this reveals that contact discontinuities cannot exist while expansion waves and compression fans are possible solutions to the governing equations. Evidence of these properties is provided on a set of numerical test cases.

This talk presents a family of algebraically constrained finite element schemes for hyperbolic conservation laws. The validity of generalized discrete maximum principles is enforced using monolithic convex limiting (MCL), a new flux correction procedure based on representation of spatial semi-discretizations in terms of admissible intermediate states. Semi-discrete entropy stability is enforced using a limiter-based fix. Time integration is performed using explicit or implicit Runge-Kutta methods, which can also be equipped with property-preserving flux limiters. In MCL schemes for nonlinear systems, problem-dependent inequality constraints are imposed on scalar functions of conserved variables to ensure physical and numerical admissibility of approximate solutions. After explaining the design philosophy behind our flux-corrected finite element approximations and showing some numerical examples, we turn to the analysis of consistency and convergence. For the Euler equations of gas dynamics, we prove weak convergence to a dissipative weak solution. The convergence analysis to be presented in this talk is joint work with Maria Lukáčová-Medvid’ová and Philipp Öffner.

Networks of hyperbolic PDEs arise in different applications, e.g. modeling water- or gas-networks or road traffic. In the first part of this talk we discuss modeling aspects of coupling conditions for hyperbolic PDEs.

Starting from an kinetic description we derive coupling conditions for the associated macroscopic equations. For this process a detailed description of the boundary layer is important. In the second part appropriate numerical methods are considered.

Different high order approaches are compared and applications to district heating or water networks are discussed.

Multidimensional simulations of magnetohydrodynamic phenomena occurring in stellar interiors are essential for understanding how stars evolve and die. The highly subsonic flow regimes found in the regions deep inside stars pose severe challenges to conventional methods of computational MHD, such as the popular "high-resolution shock-capturing'' schemes. After giving a brief overview of work on astrophysical simulations (including also supernova explosions and common-envelope evolution) in our group at Heidelberg, we summarize the challenges and present suitable numerical solvers optimized for magnetized, low-Mach-number stellar flows, implemented in our Seven-League Hydro code. We show how the choice of the numerical method can drastically affect both the performance of the code and its accuracy in real astrophysical simulations.

We study non-conservative hyperbolic systems of balance laws and are interested in development of well-balanced (WB) numerical methods for such systems. One of the ways to enforce the balance between the flux terms and source and non-conservative product terms is to rewrite the studied system in a quasi-conservative form by incorporating the latter terms into the modified global flux. The resulting system can be quite easily solved by Riemann-problem-solver-free central-upwind (CU) schemes. This approach, however, does not allow to accurately treat non-conservative products. We therefore apply a path-conservative (PC) integration technique and develop a very robust and accurate path-conservative central-upwind schemes (PCCU) based on flux globalization. I will demonstrate the performance of the WB PCCU schemes on a wide variety of examples.

This talk explores two advanced numerical methods for solving compressible two-phase flows modelled using the conservative Symmetric Hyperbolic Thermodynamically Compatible (SHTC) model proposed by Romenski et al. I first address the weak hyperbolicity of the original model in multidimensional cases by restoring strong hyperbolicity through two distinct approaches: the explicit symmetrization of the system and the hyperbolic Generalized Lagrangian Multiplier (GLM) curl-cleaning approach. Then, I will present two numerical methods to solve the proposed problem: a high-order ADER Discontinuous Galerkin (ADER-DG) scheme with an a posteriori sub-cell finite volume limiter and an exactly curl-free finite volume scheme to handle the curl involution in the relative velocity field. The latter method uses a staggered grid discretization and defines a proper compatible gradient and a curl operator to achieve a curl-free discrete solution. Extensive numerical test cases in one and multiple dimensions validate both methods' accuracy and stability. 

A Coulter counter is an impedance measurement system widely used in blood analyzers to count and size red blood cells, thus providing information about the most numerous cells of the body. In Coulter counters, cells flow through a detection zone where an electric field is imposed, which is disturbed when a cell passes through. The number of these impedance signals yield the red blood cell count, while the cell volume is supposed to be proportional to the amplitude of the signals. However, in real systems, the red blood cells trajectories in the system does not allow to verify the assumptions necessary to provide an accurate volume measurement. For a few years, IMAG has been developing the YALES2BIO solver for the prediction of red blood cell dynamics under flow. In this presentation, I will describe the fluid-structure problem and the numerical method used, then share how numerical simulation has been used to understand the signals in industrial Coulter counters and to improve the measurements of red blood cell volumes rendered by such systems. In addition, I will discuss how the mechanical properties of RBCs impact the measurements. This work has been performed during the PhD theses of Pierre Taraconat and Pierre Pottier (Horiba Medical & IMAG).

This talk is concerned with representation issues associated with the numerical solution of a unified mathematical model of continuum mechanics, due to Godunov, Peshkov, and Romenski, which can describe ideal fluids, viscous fluids and elastoplastic solids as special cases of a general continuum. The different regimes are characterized solely by the choice of material parameters and the resulting PDE system is of hyperbolic nature, with clearly defined finite wave speeds, in contrast to the standard formulation of viscous fluids via the Navier–Stokes equations. The description of such a general continuum hinges on the evolution of a matrix-valued field called distortion, which is a generalization of the inverse deformation gradient in solid mechanics. In the fluid regime, this quantity can no longer be recovered as a gradient of displacements and encodes very rich information, in particular due to the different orientations that ideal fluid parcels can be found in. The fine features of the distortion field can be challenging (or outright impossible) to resolve with standard well-tested Finite Volume methods. Degenerate situations are routinely encountered where unphysical states are generated simply as a result of taking a convex combination of two data points. We show how changing to an alternative representation of the same distortion field, obtained via polar decomposition, can be used to solve such discretization issues. Instead of the original PDE system, one can instead evolve the rotational and stretch components of the distortion matrix separately, which allows the description of the rotational components through a quaternion-valued partial differential equation. We discuss the peculiarities of quaternion PDEs and some of the discretization strategies that they enable. We present numerical examples of high-Reynolds number simulations which could not be carried out with the previous formulation of the model.

In this talk, we propose a novel collocation-based Model Order Reduction (cMOR) strategy for solving parametric advection-diffusion PDEs on moving Chimera grids. Unlike traditional projection-based MOR, cMOR solves the High-Dimensional Model on a small subset of collocation points and extends the solution to the entire domain using a global reduced basis. By leveraging the ADER method on unsteady Chimera meshes, cMOR addresses the computational challenges posed by convective dominated problems, particularly the Kolmogorov N-width barrier. Our results demonstrate the efficiency of cMOR in reducing the computational cost while maintaining accuracy.

I will present a work in collaboration with M. Bessemoulin-Chatard and T. Rey, in which we consider a non-linear kinetic model describing the reaction of two-species. It can be considered as a simplified version of models describing the generation and recombination of electron-hole pairs in semiconductors. I will introduce a finite volume discretization of this model for which we can prove an exponential decay towards the steady state using discrete hypocoercivity methods. After presenting the ideas of the proof in the continuous framework, I will highlight the main difficulties induced by the discretization process. The properties of the method will then be illustrated by several numerical examples.

Dans ce travail, nous proposons une analyse mathématique et numérique de problèmes intervenant en électrophysiologie cardiaque, notamment les modèles bidomaine et tridomaine.

La modélisation mathématique du tissu cardiaque nécessite des géométries complexes multi-échelles pour tenir compte de la taille du cœur et des processus biologiques au niveau microscopique. À l'échelle microscopique, le tissu cardiaque est soumis à des phénomènes particulièrement complexes, et il est très difficile de comprendre et de prédire son comportement à l'échelle macroscopique.

En se basant sur la loi d'Ohm de la conduction électrique et sur la conservation de la charge électrique, nous obtenons le modèle microscopique qui donne une description détaillée de l'activité électrique dans les cellules responsables de la contraction du muscle cardiaque. Ce modèle est de type elliptique, couplé à un système d'EDO non linéaire.

Grâce à des techniques d'homogénéisation et de développements asymptotiques à partir du modèle microscopique, nous obtenons un modèle modèle macroscopique. Ce dernier, de type parabolique, permet à son tour de décrire la propagation des ondes électriques dans le cœur tout entier. Nous apportons en plus, pour ces deux formulations, des résultats d'existence et d'unicité de la solution.

Dans une seconde partie, nous traitons le cas particulier du modèle monodomaine. Ce dernier étant une simplification du modèle bidomaine.

Ce modèle est représenté par un système d'équations aux dérivées partielles de type réaction diffusion non linéaire couplé à des EDO dont la résolution numérique est très coûteuse. Nous utilisons le schéma d'Euler implicite pour la discrétisation de la variable temporelle et la méthode des éléments finis pour la discrétisation en espace. Le système non linéaire obtenu est résolu grâce à la méthode de Newton.

Ensuite, en utilisant la méthodologie de la reconstruction des flux équilibrés dans l’espace $\textbf{H}(div, \Omega)$, nous établissons une estimation d’erreur a posteriori entre la solution exacte et la solution approchée à chaque pas du solveur de Newton. Cette estimation permet de certifier la solution approchée à chaque pas du solveur de linéarisation et permet également de distinguer les différentes composantes de l’erreur de la simulation numérique à savoir l’erreur de discrétisation par éléments finis et l’erreur de linéarisation issue de la méthode de Newton.

I will consider the simulation of slender structures immersed in a three-dimensional (3D) flow. By exploiting the special geometric configuration of the slender structures, this problem can be modeled by mixed-dimensional coupled equations (3D for the fluid and 1D for the solid). Several challenges must be faced when dealing with this type of problems. From a mathematical point of view, these include defining well posed trace operators of codimension two. On the computational standpoint, the non-standard mathematical formulation makes it difficult to ensure the accuracy of the solutions obtained with the mixed-dimensional discrete formulation as compared to a fully resolved one. I define the continuous formulation using the Navier-Stokes equations for the fluid and a Timoshenko beam model for the structure. I complement these models with a mixed-dimensional version of the fluid-structure interface conditions, based on the projection of kinematic coupling conditions on a finite-dimensional Fourier space. One of the fundamental advantages of this approach is that it enables the approximation of the problem within the framework of the Finite Element Method (FEM). I establish the energy stability of the discrete formulation and provide extensive numerical evidence of the accuracy of the mixed-dimensional model, notably with respect to a fully resolved (ALE based) model.

Many physical or economical applications rely on Monte-Carlo (MC) codes to solve deterministic partial differential equations (PDEs). This is the case for example for (non-exhaustive list) neutronics and photonics. The Monte-Carlo resolution implies the sampling of the physical variables: x the position, v the velocity and t the time. The simulations are costly but the MC resolution is competitive with respect to other methods due to the high dimensional (7d) deterministic problem. The numerical parameter controling the accuracy is N_MC, the number of particles. The larger N_MC, the more accurate the results. The convergence rate obeys the central limit theorem: it is O(1/sqrt(N_MC)).

Obviously, propagating uncertainties (for sensitivity analysis etc.) with respect to different parameters X \in R^d is of great interest in every of the aforementioned applications (uncertain cross-sections etc.). In fact, in our physical applications, we would like to be able to perform systematic uncertainty propagations. As a consequence, we often face a 7+d dimensional problem. Non-intrusive methods are usually applied (use of black box codes). But it demands a high number N of evaluations. In our MC resolution context, each one of them is costly. One accurate run can take several hours on hundreds of processors.

When applying any non-intrusive method to propagate uncertainties through the linear Boltzmann equation solved with an MC code, basically, the physical space (x,t,v) and the uncertain space (X) are both explored thanks to two different MC experimental designs. The first one has particles to explore the space of physical variables (x, t, v), the second one has N runs for the space of the uncertain variable X. In this non-intrusive context, the two MC samplings are tensorised in the sense we process N_MC x N = 1e9 -- 1e15 particles for an overall O()1/sqrt(N_MC) error. An uncertainty propagation study is consequently costly. The main idea of the present work comes from the fact that MC experimental designs should allow avoiding the tensorisation of the N_MC particles and N runs [1,2,3,4]. For this, we sample the whole space relative to (x,t,v,X) within the same MC design. This implies sampling the uncertain parameters X within the code, hence the intrusiveness of the approach. In practice in [1], fast convergence rates have been observed with respect to the polynomial Chaos truncation order P: the method is efficient for the linear [2], nonlinear [4] Boltzman equation and keff computations [3]. The aim of the talk is to present the details of the uncertain MC solver.

[1] G. Poëtte. A gPC-intrusive Monte-Carlo scheme for the resolution of the uncertain linear Boltzmann equation. Journal of Computational Physics, 385:135 – 162, 2019

[2] G. Poëtte. Spectral convergence of the generalized Polynomial Chaos reduced model obtained from the uncertain linear Boltzmann equation. Preprint submitted to Mathematics and Computers in Simulation , 2019.

[3] G. Poëtte and E. Brun. Efficient uncertain k eff computations with the Monte Carlo resolution of generalised Polynomial Chaos Based reduced models. Preprint,November 2020.

[4] G. Poëtte. Efficient uncertainty propagation for photonics: combining Implicit Semi-analog Monte Carlo (ISMC) and Monte Carlo generalised Polynomial Chaos (MC-gPC). Preprint, 2020.

La motilité cellulaire est un phénomène impliqué dans de nombreux processus biologiques tels que la propagation des cancers, la réponse immunitaire, la cicatrisation ou le développement embryonnaire. Après avoir présenté le contexte biologique, je présenterai un modèle à frontière libre en dimension 2 modélisant ce phénomène. Je présenterai des résultats sur l'existence et la stabilité d'états stationnaires. Enfin je présenterai un schéma numérique aux éléments finis permettant de réaliser des simulations numériques mettant en avant l'influence du noyau sur la motilité cellulaire.

Cette étude s'intéresse aux phénomènes géophysiques, et plus particulièrement aux courants de densité pyroclastiques, des mélanges complexes composés de pyroclastes, de fragments rocheux et d'air. Ces phénomènes destructeurs, capables de parcourir de grandes distances et d'impacter des zones urbanisées, se distinguent par leur capacité à se propager même sur des terrains à faible pente. La fluidisation et la dilatation de ces matériaux granulaires denses semblent jouer un rôle clé dans ces dynamiques. Des approches de modélisation ont ainsi été développées pour approfondir leur compréhension.

Un modèle de mélange fluide-solide, adapté pour intégrer les propriétés spécifiques du gaz interstitiel, a été utilisé. La compressibilité du gaz permet de reformuler l'équation de conservation de la masse de la phase gazeuse en une équation dépendant de la pression. Pour décrire la dynamique de la phase solide, l'équation de quantité de mouvement de cette phase est complétée par des lois constitutives basées sur une rhéologie seuil et une fonction de dilatance. La divergence du champ de vitesse, qui reflète la capacité de l'écoulement granulaire à s'étendre ou à se comprimer, dépend ainsi de la fraction volumique, la pression, le taux de déformation et le nombre inertiel. Ce cadre théorique fournit une description réaliste et robuste des écoulements granulaires non-isochore fluidisés et constitue une base solide pour des études numériques.

We characterised all non-zero vector-fields S ∈ [W^{−1,p}(Ω)]^n , 1 < p < ∞, n ≥ 3, whose potential, φ, linked to S by the elliptic problem ∇·(M ∇φ)= ∇·S, attains a constant value on each of the finitely many connected components of Rn \Ω, where M a symmetric positive definite matrix. Our characterisation states that such S posses a Stokes decomposition and when such S are extended by zero to R^n their Stokes decomposition vanishes identically outside Ω. We also showed that given S ∈ [W^{−1,p} (Ω)]^n there is a unique S_{nm} ∈ [W^{−1,p}(Ω)]^n of minimum norm among all vector-fields that generate the same potential as S on R^n\Ω modulo constants. We showed that when Ω admits the Gauss divergence theorem there is a unique h∗ ∈ W 2,q (Ω) such that S_{nm} = ⟨S, ∇h∗⟩∆q ∇h∗ where q = (p−1)/p and ∆q is the vector q−Laplacian hence each vector-field S ∈ [W^{−1,p}(Ω)]^n can be written as S = ∆v + ∇ψ − ⟨S, ∇h∗⟩∆q ∇h∗ for unique v ∈ [W 1,p (Ω)]n and ψ ∈ Lp (Ω). Finally, we showed that when Ω is Lipschitz, under certain circumstances it is possible to determine S_{nm} from φ on ∂Ω.

This work is devoted to the numerical simulation of low Mach number flows, modeled by the compressible Euler system. Commonly used solvers for discretizing this model are Godunov-type schemes. These schemes exhibit poor performance at low Mach number in terms of efficiency and accuracy.

Indeed, regarding the accuracy problem observed with quadrangular grids, it arises from the fact that the discrete solution fails to converge to the incompressible solution as the Mach number tends to zero. To overcome this accuracy problem, many fixes have been developed and consist in modifying the numerical diffusion of the original scheme. These corrections improve the accuracy of compressible schemes as the Mach number goes to zero. Unfortunately they introduce other problems, such as the appearance of numerical oscillations (checkerboard modes on a Cartesian grid) in the numerical solution, or the damping of acoustic waves as the Mach number goes to zero. Efficiency is also compromised as these schemes are stable under a more restrictive CFL condition compared to the original scheme.

In this talk, we propose to study the phenomenon of oscillations that plagues some of the fixes proposed in the literature. We focus on Roe-type fixes, in particular those that reduce the numerical diffusion on the jump of the normal velocity. The asymptotic analysis of these schemes leads to a discretization of a wave system in which the pressure gradient is centered. To better understand the phenomenon, we focus on the linear wave system. We then show that this fix is not TVD, unlike the Godunov scheme, which explains the appearance of numerical oscillations in the unsteady solution. Next, we study the long-time behavior of the numerical solution. It turns out that spurious stationary oscillations appear on the velocity field, preventing mesh convergence of the numerical solution. Moreover, the dimension of the space of spurious elements is huge, making it difficult to consider any form of filtering for these modes. We conclude this work by noting that the modification of the numerical diffusion makes it difficult to develop numerical methods that are both stable and accurate for low Mach number flows.

Je montrerai tout d'abord comment à partir de considérations énergétiques et du principe de moindre action on peut utiliser une fonction indicatrice de phase (fonction couleur) pour modéliser l'effet des forces capillaires sur la dynamique d'un film mince. Je montrerai ensuite comment construire un solveur HLLC pour discrétiser le système d'edp obtenu dans le cas où on peut négliger les effets liés à la pression capillaire (donc liés à la courbure de l'interface film-air). Je terminerai par l'état d'avancement de nos travaux dans le cas où toutes les forces capillaires sont prises en compte dans le modèle. Des résultats numériques permettront d'illustrer la présentation et de mettre en évidence ce qui marche mais aussi ce qui ne marche pas encore ... Il s'agit d'un travail commun réalisé avec B. Delacroix, G. Blanchard, M. Bouyges et C. Laurent.

Aortic stenosis (AS) is a heart disease characterized by the narrowing of the aortic valve, restricting the blood flow from the left ventricle. Without timely intervention, severe AS can lead to life-threatening complications. Aortic valve replacement (AVR) is the primary treatment of AS, replacing the diseased valve with a mechanical or biological prosthesis to restore normal blood flow and improve the heart function.

AVR is associated with ventricular remodeling which includes a reduction in cardiac mass and an increase in extracellular volume (ECV). Both aspects of the remodeling process can be studied by magnetic resonance imaging (MRI), with ECV being recognized as a marker of cardiac fibrosis and a predictor of mortality in AS patients.

We developed personalized computational models of the human ventricles using MRI data of 12 AS patients before and 3 months after AVR. Using these models, we provide new insights into the remodeling process that follows AVR. We demonstrate that an ECV increase does not have to be associated with changes in cardiac fibrosis, and that the remodeling process does not affect conduction velocity, an important parameter in cardiac electrophysiology.

In this talk, I review recent efforts on the development of registration methods for parametric model order reduction (MOR), with emphasis on advection-dominated flows. In computer vision and pattern recognition, registration refers to the process of finding a parametric transformation that aligns two datasets; in model order reduction, registration methods seek a parametric bijection that tracks coherent structures (e.g., shocks, shear layers) of the solution field. We integrate registration in the offline/online model reduction framework to tackle problems with parameter-dependent discontinuities. Our approach combines registration with three additional building blocks: (i) an hyper-reduced least-squares Petrov-Galerkin (LSPG) reduced-order model, to estimate the mapped solution; (ii) a parametric mesh adaptation procedure to build a parsimonious yet accurate representation of the solution field; and (iii) a multi-fidelity strategy to reduce offline training costs. We present numerical results for several two-dimensional inviscid compressible flows, to show the potential of the method.

On the mesoscopic level, motion of individual particles can be modeled by a kinetic transport equation for the population density f(t,x,v) as a function of time t, space x and velocity v \in V. A relaxation term on the right hand side accounts for scattering due to self-induced velocity changes and typically involves a parameter K(x,v,v') encoding the probability of changing from velocity v' to v at location x:

\partial_t f(t,x,v) + v \cdot abla f(t,x,v) = \int K(x,v,v') f(t,x,v') - K(x,v',v)f(t,x,v) dv'

This hyperbolic model is widely used to model bacterial motion, called chemotaxis.

We study the inverse parameter reconstruction problem whose aim is to recover the scattering parameter $K$ and that has to be solved when fitting the model to a real situation. We restrict ourselves to macroscopic, i.e. velocity averaged data $\rho = \int f dv$ as a basis of our reconstruction. This introduces additional difficulties, which can be overcome by the use of short time interior domain data. In this way, we can establish theoretical existence and uniqueness of the reconstruction, study its macroscopic limiting behavior and numerically conduct the inversion under suitable data generating experimental designs.

This work based on a collaboration with Kathrin Hellmuth (Würzburg, Germany), Qin Li (Madison, Wisc., USA) and Min Tang (Shanghai, China).

The efficient numerical solution of compressible Euler equations of gas dynamics poses several computational challenges. Severe time restrictions are required by standard explicit time discretization techniques, in particular for flow regimes characterized by low Mach number values. We analyze schemes based on a general Implicit-Explicit (IMEX) Runge-Kutta (RK) time discretization for the compressible Euler equations of gas dynamics. We prove the asymptotic-preserving (AP) property of the numerical scheme in the low Mach number limit for both a single asymptotic length scale and two length scales. The analysis is carried out considering a general equation of state (EOS) and therefore it is not restricted to the ideal gas law as done for standard AP schemes. We couple implicitly the energy equation to the momentum one, while treating the continuity equation in an explicit fashion.

We also present an alternative strategy, which consists of employing a semi-implicit temporal integrator based on IMEX-RK methods (SI-IMEX-RK). The stiff dependence is carefully analyzed, so as to avoid the solution of a nonlinear equation for the pressure also for equations of state (EOS) of non-ideal gases.

The asymptotic-preserving (AP) and the asymptotically-accurate (AA) properties of the two approaches are assessed on a number of classical benchmarks for ideal gases and on their extension to non-ideal gases. Particular attention is devoted to the robustness of the numerical methods with respect to the boundary conditions and to non well-prepared initial conditions.

This work is based on a collaboration with Luca Bonaventura (Politecnico di Milano), Sebastiano Boscarino and Giovanni Russo (Università di Catania).