Add a 1D macroscopic traffic-flow solver (`Traffic1D`) via CESE

[tigercosmos]

Summary

Add a 1D macroscopic traffic-flow solver (Traffic1D) to modmesh, built on the CESE machinery in cpp/modmesh/onedim/. The Lighthill-Whitham-Richards (LWR) model is a scalar nonlinear hyperbolic conservation law, so Traffic1D is the first scalar solver in onedim/ -- a single-variable (NVAR=1, state u = rho) analog of the existing Euler1DCore.

The onedim/ module currently solves only the 1D Euler equations (Euler1DCore); a legacy Burgers solver exists under the deprecated spacetime/ module but is not the basis for this work. The deliverable targets the canonical congestion patterns of a first-order model: queue formation at a red light (a shock) and queue discharge at a green light (a rarefaction).

Model

The LWR model conserves vehicle density rho(x, t) along a road coordinate x:

d(rho)/dt + d(q)/dx = 0,    q(rho) = rho * V(rho)

With the Greenshields closure V(rho) = v_max * (1 - rho/rho_max), the flux is quadratic and concave:

q(rho)   = v_max * rho * (1 - rho/rho_max)     (quadratic, concave)
q'(rho)  = v_max * (1 - 2*rho/rho_max)         (kinematic wave speed)
q''(rho) = -2*v_max/rho_max < 0                (concavity)

Key derived quantities, used directly for validation:

critical density  rho_c = rho_max / 2          (q'(rho_c) = 0, standing wave)
road capacity     q_max = v_max * rho_max / 4   (peak of the fundamental diagram)

Waves travel downstream for rho < rho_c and upstream for rho > rho_c, with the wave speed q'(rho) spanning [-v_max, +v_max].

Relation to inviscid Burgers

This is the same class as inviscid Burgers -- a scalar law with a quadratic flux -- but with the opposite convexity: Burgers' flux u^2/2 is convex, whereas the Greenshields flux is concave. Consequently the Riemann wave types are reversed, so Burgers' shock/rarefaction logic must not be reused verbatim:

Riemann data	Burgers (convex)	LWR/Greenshields (concave)
rho_l < rho_r	rarefaction	shock
rho_l > rho_r	shock	rarefaction

Equivalently, the characteristic speed c = q'(rho) itself satisfies Burgers' equation, but c is a decreasing function of rho, which is what flips the role of the left and right states.

Physics-to-PDE map

Traffic phenomenon	PDE feature
Red light: queue tail	shock (rho_l < rho_r); s = [q]/[rho] < 0, upstream
Green light: discharge	rarefaction fan (rho_l > rho_r)
Vehicle conservation	exact flux balance

CESE captures shocks without an explicit artificial-viscosity term; the necessary numerical dissipation is supplied by the a-scheme's gradient-amplitude weighting (the ALPHA parameter in march_half_so1_alpha). Note that the scheme becomes very dissipative at small CFL, so dt should not be over-conservative or the shock will smear.

Scope and approach

Euler1DCore (cpp/modmesh/onedim/Euler1DCore.hpp) is the structural template, but it is hard-wired for Euler (NVAR=3, m_gamma, the Euler flux/Jacobian, and 3-component boundary updates). So Traffic1D is best done as a specialized scalar copy/adaptation (NVAR=1, state u = rho), not a "swap the flux" change and not an up-front refactor to a generic base.

The marching structure (update_cfl, march_half_so0, march_half_so1_alpha, treat_boundary_*, march_alpha) is reused as a pattern; the new physics is the scalar flux:

// Greenshields flux for the LWR model (scalar; NVAR = 1, u[0] = rho)
double flux(double rho)  const { return m_vmax * rho * (1.0 - rho / m_rhomax); }
double dflux(double rho) const { return m_vmax * (1.0 - 2.0 * rho / m_rhomax); }
// PDE: rho_t = -q_x  =>  ut[0] = -f_x ;  ft[0] = dflux(rho) * ut[0]
// CFL: match Euler1DCore::update_cfl, which uses the HALF step:
//   cfl = |dflux(rho)| * (dt/2) / dx,  dx = min(dxpos, dxneg)

Parameters v_max (free-flow speed) and rho_max (jam density) replace m_gamma; rho is in physical units (vehicles per length). Examples may normalize to rho_max = v_max = 1, giving q = rho*(1-rho), q' = 1-2*rho, rho_c = 0.5, and q_max = 0.25.

Two inherited constraints must be kept. The CESE staggered mesh requires odd ncoord (Euler1DCore::initialize_data rejects even sizes), and the per-half-step CFL must use the half time increment to match Euler1DCore::update_cfl (and the legacy spacetime Burgers core, which uses field().hdt() likewise).

Files and integration points

cpp/modmesh/onedim/Traffic1DCore.hpp / .cpp   # the scalar CESE core

• cpp/modmesh/onedim/CMakeLists.txt -- add the new sources.
• cpp/modmesh/onedim/pymod/wrap_onedim.cpp -- add a WrapTraffic1DCore (the wrapper file is centralized; no new wrap file).
• modmesh/onedim/__init__.py -- export a traffic1d module with a Traffic1DSolver (mirroring euler1d.Euler1DSolver).
• (optional, follow-up) modmesh/pilot/_traffic1d.py + registration in modmesh/pilot/_gui.py, on OneDimBaseApp, to animate rho, V, q.

v1 scope: initial / boundary conditions

• Initial conditions: Riemann data (density jumps) for the canonical red-light stop (rho_l < rho_r -> shock) and green-light discharge (rho_l > rho_r -> rarefaction).
• Boundaries: closed / extrapolation only, matching the current solver.
• Defer real inflow/outflow flux boundaries: admissible LWR boundary data depends on wave direction (supply/demand), which needs its own design.

Validation

• Shock speed vs Rankine-Hugoniot s = [q]/[rho] for the red-light problem (expect s < 0, i.e. the jam front propagates upstream).
• Rarefaction-fan shape for the green-light problem.
• Exact vehicle-count conservation (integral of rho over the domain).
• Recover the Greenshields fundamental diagram q(rho), including capacity q_max = v_max*rho_max/4 at rho_c = rho_max/2.
• Monitor positivity 0 <= rho <= rho_max as a test: the exact entropy solution obeys a maximum principle, but the CESE a-scheme is not provably monotone, so overshoots near the shock are a real failure mode to watch -- hence a test, not an assumed invariant.
• Reuse tests/test_onedim_euler.py as a structural template; add tests/test_onedim_traffic.py.

Non-goals / follow-ups

• Real inflow/outflow (supply-demand) boundary conditions.
• Spontaneous stop-and-go / phantom jams and string-stability studies: a first-order LWR model only transports the imposed initial condition and cannot nucleate jams on its own, so these require a second-order model (Aw-Rascle-Zhang / Payne-Whitham).
• Multi-lane coupling.

Plan

1. Traffic1DCore (LWR + Greenshields) + Python Traffic1DSolver + Riemann and conservation tests.
2. Optional follow-up: pilot app for live visualization.

References

• R. J. LeVeque, D. I. Ketcheson, K. T. Mandli, Riemann Problems and Jupyter Solutions, "Traffic flow: the Lighthill-Whitham-Richards model" -- concave LWR flux, decreasing characteristic speed, red-light shock moving upstream. http://www.clawpack.org/riemann_book/html/Traffic_flow.html
• M. Treiber and A. Kesting, Traffic Flow Dynamics, Lecture 06 "The LWR Model" (TU Dresden) -- Greenshields fundamental diagram, wave speed c = Q'(rho), shock equation c_12 = (Q_2-Q_1)/(rho_2-rho_1), and the traffic-light shock/rarefaction picture. https://mtreiber.de/Vkmod_Skript/Lecture06_Macro_LWR.pdf
• X.-Y. Wang, A Summary of the Space-Time Conservation Element and Solution Element (CESE) Method, NASA/TM-2015-218743, 2015 -- nondissipative baseline with controlled dissipation, the alpha = 0,1,2 gradient limiter, no Riemann solver, CFL < 1. https://ntrs.nasa.gov/api/citations/20160001345/downloads/20160001345.pdf
• B. D. Greenshields, "A study of traffic capacity," Highway Research Board Proceedings, 1935 -- the linear speed-density closure.
• M. J. Lighthill and G. B. Whitham, "On kinematic waves II: A theory of traffic flow on long crowded roads," Proc. R. Soc. Lond. A, 1955; P. I. Richards, "Shock waves on the highway," Oper. Res., 1956 -- the LWR model.
• S. C. Chang, "The method of space-time conservation element and solution element," J. Comput. Phys. 119(2), 1995 -- the original CESE scheme and a-scheme dissipation.
• A. Aw and M. Rascle, SIAM J. Appl. Math., 2000; H. M. Zhang, Transp. Res. B, 2002 -- second-order (ARZ) models that, unlike first-order LWR, reproduce stop-and-go jams.

[yungyuc]

This is an interesting problem. But we should start with a better problem statement.