Classical and Quantum Ideal Gases

Exact Riemann solvers for classical and quantum Euler gases, 1-D Navier–Stokes–Fourier (NSF) solvers, and a fast polylogarithm kernel used to resolve the quantum equation of state.

This repository ports the MATLAB implementation found in this thesis to Python 3.11. The polylog function has been ported from the MATLAB implementation to C++ and the Toro exact Riemann solver has been extended to support Fermi–Dirac (FD), Bose–Einstein (BE), and Maxwell–Boltzmann (MB) statistics.

Requirements

• Python 3.11+
• pip (or uv)

Building from source additionally requires a C++17 compiler. See DEVELOPER_GUIDE.md.

Installation

pip install ideal-gases

For plotting (euler plot, interactive explorers):

pip install ideal-gases[plot]

After install, the euler command-line tool is available.

Interactive mode

Launch matplotlib widget explorers to build custom Riemann problems with sliders, statistic toggles (quantum), and Save/Reset controls. Y-axis limits autoscale automatically on each update.

euler interactive classical
euler interactive quantum

Seed the initial state from CLI flags or a JSON config (same fields as euler solve):

euler interactive classical --gamma 1.4 --t-end 0.5 --nx 101
euler interactive quantum --rho-l 2 --t-l 1.5 --n 3 --h 0.5
euler interactive classical --config case.json

Optional domain flags (--x-min, --x-max, --x0, --nx) default to an interactive Sod-tube layout (x in [-10, 10], discontinuity at x0=0, nx=1024). Use -f path.png to set the Save button target; nothing is written until you click Save.

Example usage

euler interactive quantum

Outputs a Sod shock tube problem resolved with a quantum Euler solver for all statistics. We deactivate the solutions of MB and BE to focus on the FD solution. Using the slider, we can vary the left and right states and the thermal scale parameter h and the number of degrees of freedom n of the gas.

In Fig. 5 of Hu and Jing (2010), a fictitious 2-d fermi gas degenerate regime is used to prove the accuracy of Kinetic Flux Vector Splitting schemes for quantum Euler equations. Using the interactive mode, we set n : 2 and set the left and right states ($\rho,u,\theta$). Using the h slider, we found that the degenerate gas is resolved approximately for h $\approx$ 3.71. As show in the following figure:

Command-line mode

Compute exact solution profiles, save plots to PNG, and write CSV/JSON files with the solution fields.

Classical Sod shock tube

euler solve classical \
  --rho-l 1 --u-l 0 --p-l 1 \
  --rho-r 0.125 --u-r 0 --p-r 0.1 \
  --t-end 0.25 --gamma 1.4 \
  --nx 101 -o sod.csv

Quantum Euler

euler solve quantum \
  --rho-l 1 --u-l 0 --t-l 1 \
  --rho-r 0.125 --u-r 0 --t-r 0.25 \
  --t-end 0.20 --n 2 --h 0.1 --statistic FD \
  -o euler_fd.csv

Write separate files for FD, MB, and BE with --all-statistics (e.g. euler_case7_FD.csv, euler_case7_MB.csv, euler_case7_BE.csv):

euler solve quantum ... --all-statistics -o euler_case7

Equilibrium inversions

Compute the fugacity from density and temperature:

euler fugacity --rho 1.0 --theta 1.0 --n 3 --h 1.0 --statistic FD

Recover fugacity, temperature and pressure from density and internal energy:

euler moments --rho 1.0 --e 1.5 --n 3 --h 1.0 --statistic FD

Use -o result.json to write JSON output instead of printing to stdout.

Built-in benchmarks

euler toro 1 -o toro_test1.csv
euler list --toro

euler quantum-example 7 --all-statistics -o euler_eg7
euler list --quantum

JSON config files

Define a problem in JSON and run it with euler run or pass --config to euler solve:

euler run --config case.json
euler solve classical --config case.json -o override.csv

Example case.json:

{
  "mode": "quantum",
  "left": {"rho": 1.0, "u": 0.0, "theta": 1.0},
  "right": {"rho": 0.125, "u": 0.0, "theta": 0.25},
  "t_end": 0.20,
  "n": 2.0,
  "h": 0.1,
  "statistic": "FD",
  "all_statistics": true,
  "format": "json",
  "output": "euler_case7",
  "domain": {"x_min": 0.0, "x_max": 1.0, "x0": 0.5, "nx": 101}
}

Use --format json (or a .json output path) for JSON instead of CSV. CLI flags override values from the config file.

Visualization

Save a classical Sod shock tube figure:

euler plot classical \
  --rho-l 1 --u-l 0 --p-l 1 \
  --rho-r 0.125 --u-r 0 --p-r 0.1 \
  --t-end 0.2 --gamma 1.4 --nx 101 \
  -f sod.png

Plot a single quantum statistic or compare FD/MB/BE:

euler plot quantum \
  --rho-l 1 --u-l 0 --t-l 1 \
  --rho-r 0.125 --u-r 0 --t-r 0.25 \
  --t-end 0.20 --n 2 --h 0.1 --statistic FD \
  -f qfd.png

euler plot quantum-example 7 --all-statistics -f eg7

With --all-statistics, -f eg7 writes eg7_panels.png (3×6 grid) and eg7_comparison.png (overlay). Use --layout panels|comparison|both to select one or both (default: both). Add --show for an interactive window, or -o to export CSV/JSON in the same run.

Example usage

In Filbet, Hu and Jing (2010), the authors use a Sod shock tube initial condition with a fictitious 2-d fermi and bose gas to prove the accuracy of their numerical scheme in classical and quantum hydronamic regimes. These are examples 7 and 8, respectively, in the CLI plot tool.

euler plot quantum-example 7 --all-statistics -f sod_2d_gas_classical --layout comparison --show

yields the following plot:

euler plot quantum-example 8 --all-statistics -f sod_2d_gas_quantum --layout comparison --show

yields the following plot:

Python module

Import ideal_gases to compute classical and quantum Euler and NSF solutions in your own scripts.

Classical Euler

import numpy as np
from ideal_gases import classical_euler

x = np.linspace(0.0, 1.0, 101)
result = classical_euler(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    gamma=1.4,
    x=x,
    x0=0.5,
)

Quantum Euler (FD / BE / MB)

Left and right states are given in terms of density rho, velocity u, and temperature theta (written t in the API). The solver converts these to effective pressures via the quantum EOS, then applies the Toro exact Riemann solver.

import numpy as np
from ideal_gases import quantum_euler

x = np.linspace(0.0, 1.0, 101)
result = quantum_euler(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    n=2.0,          # degrees of freedom; gamma = (n+2)/n
    h=0.1,          # thermal scale parameter
    statistic="FD", # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
)

This returns a RiemannResult object that contains the solution fields: x, rho, ux, p, e, z (fugacity), t (temperature), mach, entropy.

In the classical limit, MB statistics with h → 0 recover the ideal-gas behaviour (pressures p = rho * theta).

Classical NSF

1-D Navier–Stokes–Fourier for a monatomic ideal gas. Same Sod left/right states as the classical Euler example; dim in {1, 2, 3} sets γ = (dim+2)/dim (unlike Euler's free gamma). kn is the Knudsen number used by the Chapman–Enskog closure μ = kn ρ T.

import numpy as np
from ideal_gases import classical_nsf

x = np.linspace(0.0, 1.0, 101)
result = classical_nsf(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    dim=3,
    kn=0.01,
    x=x,
    x0=0.5,
)

This returns a ClassicalNSFResult object that contains the cell-centered fields: rho, u, t (temperature), p, q (heat flux).

Quantum NSF (FD / BE / MB)

1-D Navier–Stokes–Fourier with the quantum EOS. Same left/right states as the quantum Euler example; dim replaces Euler's n, and kn sets the Chapman–Enskog viscosity μ = kn p(z).

import numpy as np
from ideal_gases import quantum_nsf

x = np.linspace(0.0, 1.0, 101)
result = quantum_nsf(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    dim=2,
    h=0.1,
    kn=0.01,
    statistic="FD", # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
)

This returns a QuantumNSFResult object (NSFResult is an alias) that contains: rho, u, t, p, z (fugacity), q (heat flux).

Equilibrium inversions

Given density and temperature, recover the fugacity:

from ideal_gases import find_fugacity

z = find_fugacity(rho=1.0, T=1.0, dim=3, h=1.0, eta=-1)

Given density and internal energy, recover fugacity, temperature and pressure:

from ideal_gases import find_moments

z, T, p = find_moments(rho=1.0, e=1.5, dim=3, h=1.0, eta=-1)

The eta parameter selects the statistic: -1 Fermi, 0 classical (Maxwell-Boltzmann), +1 Bose.

Polylogarithm module

Quantum solvers (G, find_moments, quantum_euler) and polylog(n, z) use the unified C++ kernel: Fukushima minimax Fermi–Dirac / Bose–Einstein integrals for supported half-integer orders on z < 0 and 0 < z < 1, with Bhagat / integer analytic branches as fallback.

We can use the polylogarithm module on our scripts as follows:

import numpy as np
from ideal_gases import polylog

polylog(2, 0.5)                         # scalar
polylog(1.5, np.linspace(0.2, 0.9, 50)) # array

We can plot the polylogarithm function to verify the accuracy of the implementation for integer and half-integer orders as follows:

uv run scripts/plot_polylogarithms.py

yields the following plot:

Public API

from ideal_gases import (
    G,
    ClassicalNSFResult,
    QuantumNSFResult,
    RiemannResult,
    adiabatic_index,
    classical_euler,
    classical_nsf,
    equilibrium_moments,
    find_fugacity,
    find_moments,
    polylog,
    quantum_euler,
    quantum_nsf,
)

License

MIT License. See LICENSE for the full text.

Copyright (c) 2026 Manuel A. Diaz

For building from source, tests, linting, CI, and releases, see DEVELOPER_GUIDE.md.