Exact Riemann Solver for the 1-D Euler Equations

A MATLAB implementation of the exact Riemann solver for the one-dimensional Euler equations of gas dynamics. The code finds the intermediate pressure (p*) using Newton–Raphson iteration, then samples the full wave structure (shocks, rarefactions, and contact discontinuity) to produce the exact solution at any given time.

The default test case is the classic Sod shock-tube problem.

Background

The Riemann problem asks: given two constant gas states separated by a membrane at x = 0, what happens when the membrane is removed? The solution consists of three waves — a left wave (shock or rarefaction), a contact discontinuity, and a right wave (shock or rarefaction) — that separate four constant-state regions.

The key unknown is the pressure p* in the intermediate region. It satisfies

f_L(p*) + f_R(p*) + u_R − u_L = 0

where f_L and f_R take different forms depending on whether the corresponding wave is a shock or a rarefaction. This code solves the equation above with Newton–Raphson and then reconstructs the full spatial profile of density, velocity, pressure, and specific internal energy.

Reference: Toro, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd edition, Springer.

Files

File	Description
s_start_main.m	Main driver script — sets initial conditions (Sod problem by default), calls the solver, and saves a figure.
f_euler_riemann_root_finder.m	Newton–Raphson iteration to find p*.
f_fofp.m	Evaluates the pressure function f(p) = f_L + f_R + u_R − u_L. Contains local sub-functions for the left and right contributions, each branching on shock vs. rarefaction.
f_fofpprime.m	Evaluates f'(p), the derivative needed by Newton–Raphson, again branching on shock vs. rarefaction for each side.
f_plot_riemann.m	Samples the exact solution on a spatial grid at time T, handling all four regions and the fan inside any rarefaction wave.
f_plot_solution.m	Produces a 2×2 subplot figure of density, velocity, pressure, and specific internal energy.

Quick Start

1. Clone or download this repository.
2. Open MATLAB and navigate to the project folder.
3. Run:

   s_start_main

   This will print nothing to the command window (the iteration is silent) and save a figure figure1_sod_solution.png showing the four solution profiles at T = 0.25.

Changing the Problem

Edit the left and right primitive-variable vectors in s_start_main.m:

WL = [rho_L, u_L, p_L];
WR = [rho_R, u_R, p_R];

You can also change the domain limits (LX0, LX1), number of grid points (N), and solution time (T). The initial pressure guess pguess is set to the arithmetic mean of p_L and p_R, which works well for most moderate-strength problems.

License

This project is licensed under the GNU General Public License v3.0 — see the LICENSE file for details.