Learning the Exact Flux: Neural Riemann Solvers with Hard Constraints

Abstract

Godunov-type methods, which obtain numerical fluxes through local Riemann problems at cell interfaces, are among the most fundamental and widely used numerical methods in computational fluid dynamics. Exact Riemann solvers faithfully solve the underlying equations, but can be computationally expensive due to the iterative root-finding procedures they often require. Consequently, most practical computations rely on classical approximate Riemann solvers, such as Rusanov and Roe, which trade accuracy for computational speed. Neural networks have recently shown promise as an alternative for approximating exact Riemann solvers, but most existing approaches are data-driven or impose weak constraints. This may result in problems with maintaining balanced states, symmetry breaking, and conservation errors when integrated into a Godunov-type scheme. To address these issues, we propose a hard-constrained neural Riemann solver (HCNRS) and enforce five constraints: positivity, consistency, mirror symmetry, Galilean invariance, and scaling invariance. Numerical experiments are carried out for the shallow water and ideal-gas Euler equations on standard benchmark problems. In the absence of hard constraints, violations of the well-balanced property, mass conservation, and symmetry are observed. Notably, in the Euler implosion problem, the exact Riemann solver with MUSCL-Hancock captures the jet structure well, whereas the Rusanov flux is too diffusive and smears it out. HCNRS accurately reproduces the solution obtained by the exact Riemann solver. In contrast, an unconstrained neural formulation lacks mirror symmetry, which makes the solution depend on the choice of flux normal direction. As a result, the jet is either shifted or lost, along with diagonal symmetry.

Figures (6)

• Figure 1: Typical solutions of the Riemann problem. (a) Shallow water equations: a left-going shock and a right-going rarefaction wave. (b) Euler equations: a left-going rarefaction, a contact discontinuity, and a right-going shock. $h$ and $u$ denote water depth and velocity for the shallow water equations, while $\rho$, $u$, and $p$ denote density, velocity, and pressure for the Euler equations. Subscripts $L$ and $R$ denote the left and right initial states, respectively, and the subscript $*$ indicates intermediate (star) states. $*L$ and $*R$ denote left and right star-region states. For the Euler equations, the velocity and pressure are uniform across the contact discontinuity, while the density may differ.
• Figure 2: Still water test at $t = 0.1\,\mathrm{s}$. Top: initial condition showing bathymetry and water surface elevation. Middle: deviation in water depth $\Delta h = h - h_{\mathrm{ref}}$, where $h_{\mathrm{ref}}$ denotes the equilibrium depth. Bottom: deviation in velocity $\Delta u = u - u_{\mathrm{ref}}$, with $u_{\mathrm{ref}} = 0$.
• Figure 3: Absolute value of the mass flux returned by UCNRS for states with mirrored velocities, with $h\in[5,10]$, $u_L\in[-1,1]$, and $u_R=-u_L$. The nonzero flux indicates that UCNRS does not enforce zero mass flux at the wall.
• Figure 4: Two-dimensional dam break test at $t = 5$. Top row: spatial error in water depth, $\Delta h = h_{\mathrm{Exact RS}} - h$, for the Rusanov, UCNRS, and HCNRS. Bottom-left: midline transect of the solution, $h(x, y = 0)$. Bottom-middle: corresponding transect errors $\Delta h(x,0)$ for all methods. Bottom-right: transect error for HCNRS only.
• Figure 5: Euler implosion at $T=2.5$. Pressure is shown as a color map, overlaid with 31 density contours ($0.35$--$1.1$, step $0.025$) and velocity vectors. Formatted for comparison with Fig. 4.11 of Liska and Wendroff. Top row: Exact RS, HCNRS, Rusanov. Bottom row: UCNRS, UCNRS(rotate), HCNRS(rotate). Here "rotate" denotes the solution obtained from a $90^\circ$ counterclockwise rotation of the initial condition.