diff --git a/15_Sod_Shock_Tube/CASE/DATA/run.cfg b/15_Sod_Shock_Tube/CASE/DATA/run.cfg
new file mode 100644
index 0000000..5d182f4
--- /dev/null
+++ b/15_Sod_Shock_Tube/CASE/DATA/run.cfg
@@ -0,0 +1,5 @@
+[job_defaults]
+
+n_procs: 1
+n_threads: 1
+
diff --git a/15_Sod_Shock_Tube/CASE/DATA/setup.xml b/15_Sod_Shock_Tube/CASE/DATA/setup.xml
new file mode 100644
index 0000000..3de3879
--- /dev/null
+++ b/15_Sod_Shock_Tube/CASE/DATA/setup.xml
@@ -0,0 +1,242 @@
+
+
+
+
+
+
+
+
+
+
+ 436
+
+
+ 0
+ 0.00045662100456621
+
+
+
+ X0
+ X1
+ Y0
+ Y1
+ Z0
+ Z1
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ density = 1;
+ density = 0.125;
+ 1.17862
+
+
+ 0.0
+
+
+
+
+ 0.0
+
+
+
+
+ 1005
+
+
+
+
+ 0
+
+
+
+
+ 0
+
+
+
+ 0.028966
+ 1
+ 293.15
+
+
+ 0
+ 0
+ 0
+
+
+
+ 0
+ 0
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ all[]
+ x < 0
+ x >= 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 1
+ 1e-08
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 1.0
+ off
+
+
+
+ 0
+
+
+
+
+
+
+
+
+
+
+ 0
+
+
+
+ 1e+12
+ -1e+12
+
+ 1
+ 1e-08
+ SGDH
+
+
+
+
+
+
+
+
+
+
+ 1
+ 0
+
+
+
+ velocity[0] = 0.;
+velocity[1] = 0.;
+velocity[2] = 0.;
+ velocity[0] = 0.;
+velocity[1] = 0.;
+velocity[2] = 0.;
+
+
+
+
+
+
+
+
+
+
+
+
+ pressure = 1;
+ pressure = 0.1;
+ 2
+ 1e-12
+
+
+ 0
+
+
+ 1
+ 1e-08
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/15_Sod_Shock_Tube/POST/exactSodTube.py b/15_Sod_Shock_Tube/POST/exactSodTube.py
new file mode 100644
index 0000000..f0b7867
--- /dev/null
+++ b/15_Sod_Shock_Tube/POST/exactSodTube.py
@@ -0,0 +1,461 @@
+#!/usr/bin/env python3
+"""
+Exact solution of the 1D Riemann problem for the Euler equations
+of a perfect gas.
+
+This script computes the analytical solution of the Riemann problem
+(shock and/or rarefaction) for given left/right initial conditions.
+The solution is sampled on a uniform 1D mesh and exported in VTU/PVD
+format for visualization (ParaView).
+
+Features:
+- Perfect gas (constant gamma)
+- Exact resolution of pressure and velocity in the star region
+- Temporal sampling at constant time step
+- ParaView-compatible export
+
+The exact Riemann solver formulation for the 1D Euler equations of
+a perfect gas follows:
+
+[1] E. F. Toro,
+ "Riemann Solvers and Numerical Methods for Fluid Dynamics",
+ 3rd Edition, Springer, 2009.
+ Chapters 4 and 5 (exact solver, shocks and rarefactions).
+
+Usage:
+ pvpython exactSodTube.py
+
+Author: fhermet
+"""
+import os
+import math
+import numpy as np
+import xml.etree.ElementTree as ET
+
+# ------------------------------------------------------------
+# Physical and numerical parameters
+# ------------------------------------------------------------
+gamma = 1.4
+
+xmin, xmax = -1.0, 1.0
+nx = 32000
+x0 = 0.0 # discontinuity position
+
+# time (must match setup.xml: 436 iterations, dt = 4.5662e-4 s)
+dt = 0.00045662100456621
+n_iter = 436
+t_final = dt * n_iter
+output_stride = n_iter # write a VTU every 'output_stride' steps (set to 1 for all)
+output_dir = "."
+
+n_steps = int(round(t_final / dt))
+t_final_effective = n_steps * dt
+print(f"n_steps = {n_steps}, t_final_effective = {t_final_effective}")
+
+# ------------------------------------------------------------
+# Exact Riemann functions (perfect gas)
+# ------------------------------------------------------------
+def sound_speed(rho, p, gamma=gamma):
+ """
+ Speed of sound for a perfect gas.
+ """
+ return math.sqrt(gamma * p / rho)
+
+def f_shock(p, rhoK, pK, gamma=gamma):
+ """
+ Function f(p) for a shock wave in the exact Riemann solver.
+
+ Parameters
+ ----------
+ p : float
+ Tested pressure.
+ rhoK : float
+ Initial state density (left or right).
+ pK : float
+ Initial state pressure.
+ gamma : float
+ Ratio of specific heats.
+
+ Returns
+ -------
+ float
+ Shock contribution to the compatibility function.
+ """
+ AK = 2.0 / ((gamma + 1.0) * rhoK)
+ BK = (gamma - 1.0) / (gamma + 1.0) * pK
+ return (p - pK) * math.sqrt(AK / (p + BK))
+
+def f_rare(p, rhoK, pK, gamma=gamma):
+ """
+ Function f(p) for a rarefaction wave in the exact Riemann solver.
+ """
+ aK = sound_speed(rhoK, pK, gamma)
+ return (2.0 * aK / (gamma - 1.0)) * ((p / pK) ** ((gamma - 1.0) / (2.0 * gamma)) - 1.0)
+
+def f_total(p, rhol, ul, pl, rhor, ur, pr, gamma=gamma):
+ """
+ Total compatibility function for the Riemann problem.
+
+ Corresponds to:
+ f_L(p) + f_R(p) + (u_R - u_L)
+
+ Used to determine the star pressure p*.
+ """
+ if p > pl:
+ fL = f_shock(p, rhol, pl, gamma)
+ else:
+ fL = f_rare(p, rhol, pl, gamma)
+ if p > pr:
+ fR = f_shock(p, rhor, pr, gamma)
+ else:
+ fR = f_rare(p, rhor, pr, gamma)
+ return fL + fR + (ur - ul)
+
+def df_dp_side(p, rhoK, pK, gamma=gamma):
+ """
+ Derivative of f(p) with respect to p for one side (left or right).
+ """
+ if p > pK:
+ AK = 2.0 / ((gamma + 1.0) * rhoK)
+ BK = (gamma - 1.0) / (gamma + 1.0) * pK
+ sqrt_term = math.sqrt(AK / (p + BK))
+ return sqrt_term * (1.0 - 0.5 * (p - pK) / (p + BK))
+ else:
+ aK = sound_speed(rhoK, pK, gamma)
+ return (1.0 / (rhoK * aK)) * (p / pK) ** (-(gamma + 1.0) / (2.0 * gamma))
+
+def df_total(p, rhol, ul, pl, rhor, ur, pr, gamma=gamma):
+ """
+ Derivative of the total compatibility function with respect to pressure.
+ """
+ return df_dp_side(p, rhol, pl, gamma) + df_dp_side(p, rhor, pr, gamma)
+
+def pressure_guess_pvrs(rhol, ul, pl, rhor, ur, pr, gamma=gamma):
+ """
+ Initial estimate of the star pressure p* using the PVRS approximation
+ (Pressure-Velocity Riemann Solver), cf. Toro (2009), Section 4.6.
+
+ Returns
+ -------
+ float
+ Initial pressure for Newton iteration.
+ """
+ aL = sound_speed(rhol, pl, gamma)
+ aR = sound_speed(rhor, pr, gamma)
+ pPV = 0.5 * (pl + pr) - 0.5 * (ur - ul) * (rhol * aL + rhor * aR) / (rhol + rhor)
+ return max(1e-6, pPV)
+
+def solve_p_star(rhol, ul, pl, rhor, ur, pr, gamma=gamma,
+ tol=1e-10, max_iter=50):
+ """
+ Solve for the star pressure p* using Newton's method, cf. Toro (2009), Section 4.4.
+
+ Parameters
+ ----------
+ tol : float
+ Convergence tolerance.
+ max_iter : int
+ Maximum number of iterations.
+
+ Returns
+ -------
+ float
+ Pressure in the star region.
+ """
+ p = pressure_guess_pvrs(rhol, ul, pl, rhor, ur, pr, gamma)
+ for _ in range(max_iter):
+ f = f_total(p, rhol, ul, pl, rhor, ur, pr, gamma)
+ df = df_total(p, rhol, ul, pl, rhor, ur, pr, gamma)
+ if abs(df) < 1e-20:
+ break
+ p_new = p - f / df
+ if p_new < 1e-6:
+ p_new = 1e-6
+ if abs(p_new - p) < tol:
+ return p_new
+ p = p_new
+ return p
+
+def compute_u_star(pstar, rhol, ul, pl, rhor, ur, pr, gamma=gamma):
+ """
+ Compute the velocity in the star region u* from p*.
+ """
+ if pstar > pl:
+ fL = f_shock(pstar, rhol, pl, gamma)
+ else:
+ fL = f_rare(pstar, rhol, pl, gamma)
+ if pstar > pr:
+ fR = f_shock(pstar, rhor, pr, gamma)
+ else:
+ fR = f_rare(pstar, rhor, pr, gamma)
+ return 0.5 * (ul + ur) + 0.5 * (fR - fL)
+
+def star_region_densities(pstar, rhol, pl, rhor, pr, gamma=gamma):
+ """
+ Compute left and right densities in the star region.
+
+ Returns
+ -------
+ rhoL_star, rhoR_star : tuple(float, float)
+ Left and right star densities.
+ """
+ if pstar > pl:
+ p_ratio = pstar / pl
+ rhoL_star = rhol * (p_ratio + (gamma - 1.0) / (gamma + 1.0)) / \
+ ((gamma - 1.0) / (gamma + 1.0) * p_ratio + 1.0)
+ else:
+ rhoL_star = rhol * (pstar / pl) ** (1.0 / gamma)
+ if pstar > pr:
+ p_ratio = pstar / pr
+ rhoR_star = rhor * (p_ratio + (gamma - 1.0) / (gamma + 1.0)) / \
+ ((gamma - 1.0) / (gamma + 1.0) * p_ratio + 1.0)
+ else:
+ rhoR_star = rhor * (pstar / pr) ** (1.0 / gamma)
+ return rhoL_star, rhoR_star
+
+def sample_riemann(x, t, rhol, ul, pl, rhor, ur, pr,
+ gamma=gamma, x0=0.0):
+ """
+ Sample the exact Riemann problem solution at time t, cf. Toro (2009), Section 4.5.
+
+ Parameters
+ ----------
+ x : array_like
+ Spatial positions.
+ t : float
+ Time.
+ x0 : float
+ Initial discontinuity position.
+
+ Returns
+ -------
+ rho, u, p : ndarray
+ Density, velocity and pressure fields.
+ """
+ x = np.asarray(x)
+
+ if t == 0.0:
+ rho = np.where(x < x0, rhol, rhor)
+ u = np.where(x < x0, ul, ur)
+ p = np.where(x < x0, pl, pr)
+ return rho, u, p
+
+ # solve once for this problem (independent of t)
+ pstar = solve_p_star(rhol, ul, pl, rhor, ur, pr, gamma)
+ ustar = compute_u_star(pstar, rhol, ul, pl, rhor, ur, pr, gamma)
+ rhoL_star, rhoR_star = star_region_densities(pstar, rhol, pl, rhor, pr, gamma)
+
+ aL = sound_speed(rhol, pl, gamma)
+ aR = sound_speed(rhor, pr, gamma)
+
+ # wave speeds
+ if pstar > pl:
+ sL = ul - aL * math.sqrt((gamma + 1.0)/(2.0*gamma) * (pstar/pl - 1.0) + 1.0)
+ headL = tailL = None # not used
+ aL_star = None
+ else:
+ headL = ul - aL
+ aL_star = aL * (pstar / pl) ** ((gamma - 1.0)/(2.0*gamma))
+ tailL = ustar - aL_star
+ sL = None
+
+ if pstar > pr:
+ sR = ur + aR * math.sqrt((gamma + 1.0)/(2.0*gamma) * (pstar/pr - 1.0) + 1.0)
+ headR = tailR = None
+ aR_star = None
+ else:
+ headR = ur + aR
+ aR_star = aR * (pstar / pr) ** ((gamma - 1.0)/(2.0*gamma))
+ tailR = ustar + aR_star
+ sR = None
+
+ xi = (x - x0) / t
+
+ rho = np.zeros_like(xi)
+ u = np.zeros_like(xi)
+ p = np.zeros_like(xi)
+
+ for i, s in enumerate(xi):
+ if s < ustar:
+ # left side
+ if pstar > pl:
+ # left shock
+ if s < sL:
+ rho[i] = rhol
+ u[i] = ul
+ p[i] = pl
+ else:
+ rho[i] = rhoL_star
+ u[i] = ustar
+ p[i] = pstar
+ else:
+ # left rarefaction
+ if s < headL:
+ rho[i] = rhol
+ u[i] = ul
+ p[i] = pl
+ elif s > tailL:
+ rho[i] = rhoL_star
+ u[i] = ustar
+ p[i] = pstar
+ else:
+ u_i = (2.0/(gamma+1.0)) * (aL + 0.5*(gamma-1.0)*ul + s)
+ a_i = aL - 0.5*(gamma-1.0)*(u_i - ul)
+ rho[i] = rhol * (a_i / aL) ** (2.0/(gamma-1.0))
+ u[i] = u_i
+ p[i] = pl * (a_i / aL) ** (2.0*gamma/(gamma-1.0))
+ else:
+ # right side
+ if pstar > pr:
+ # right shock
+ if s > sR:
+ rho[i] = rhor
+ u[i] = ur
+ p[i] = pr
+ else:
+ rho[i] = rhoR_star
+ u[i] = ustar
+ p[i] = pstar
+ else:
+ # right rarefaction
+ if s > headR:
+ rho[i] = rhor
+ u[i] = ur
+ p[i] = pr
+ elif s < tailR:
+ rho[i] = rhoR_star
+ u[i] = ustar
+ p[i] = pstar
+ else:
+ u_i = (2.0/(gamma+1.0)) * (-aR + 0.5*(gamma-1.0)*ur + s)
+ a_i = aR + 0.5*(gamma-1.0)*(u_i - ur)
+ rho[i] = rhor * (a_i / aR) ** (2.0/(gamma-1.0))
+ u[i] = u_i
+ p[i] = pr * (a_i / aR) ** (2.0*gamma/(gamma-1.0))
+
+ return rho, u, p
+
+# ------------------------------------------------------------
+# VTU / PVD writing
+# ------------------------------------------------------------
+def write_vtu(filename, x, rho, u, p, E):
+ """
+ Write a 1D VTU file containing the physical fields.
+
+ Written fields:
+ - rho: density
+ - velocity: velocity
+ - pressure: pressure
+ - energy: total energy
+ - temperature: temperature (perfect gas)
+
+ Parameters
+ ----------
+ filename : str
+ Path to the VTU file.
+ """
+ npts = len(x)
+ ncells = npts - 1
+
+ vtk = ET.Element("VTKFile", type="UnstructuredGrid",
+ version="0.1", byte_order="LittleEndian")
+ grid = ET.SubElement(vtk, "UnstructuredGrid")
+ piece = ET.SubElement(grid, "Piece",
+ NumberOfPoints=str(npts),
+ NumberOfCells=str(ncells))
+
+ point_data = ET.SubElement(piece, "PointData", Scalars="rho")
+
+ def add_array(name, data):
+ arr = ET.SubElement(point_data, "DataArray",
+ type="Float64", Name=name, format="ascii")
+ arr.text = " ".join(map(str, data))
+
+ add_array("rho", rho)
+ add_array("velocity", u)
+ add_array("pressure", p)
+ add_array("energy", E)
+ add_array("temperature", p/(rho*287.1))
+
+
+ points = ET.SubElement(piece, "Points")
+ coords = ET.SubElement(points, "DataArray",
+ type="Float64", NumberOfComponents="3",
+ format="ascii")
+ coords.text = " ".join(f"{xi} 0 0" for xi in x)
+
+ cells = ET.SubElement(piece, "Cells")
+
+ conn = ET.SubElement(cells, "DataArray",
+ type="Int32", Name="connectivity", format="ascii")
+ conn.text = " ".join(f"{i} {i+1}" for i in range(ncells))
+
+ offs = ET.SubElement(cells, "DataArray",
+ type="Int32", Name="offsets", format="ascii")
+ offs.text = " ".join(str(2 * (i+1)) for i in range(ncells))
+
+ types = ET.SubElement(cells, "DataArray",
+ type="UInt8", Name="types", format="ascii")
+ types.text = " ".join("3" for _ in range(ncells)) # 3 = VTK_LINE
+
+ os.makedirs(os.path.dirname(filename), exist_ok=True)
+ tree = ET.ElementTree(vtk)
+ tree.write(filename)
+ print(f"[VTU] Written: {filename}")
+
+def write_pvd(filename, files, times):
+ """
+ Write a PVD file referencing a time series of VTU files.
+
+ Parameters
+ ----------
+ files : list of str
+ VTU file names.
+ times : list of float
+ Associated times.
+ """
+ vtk = ET.Element("VTKFile", type="Collection",
+ version="0.1", byte_order="LittleEndian")
+ coll = ET.SubElement(vtk, "Collection")
+
+ for f, t in zip(files, times):
+ ET.SubElement(coll, "DataSet",
+ timestep=str(t), group="", part="0", file=f)
+
+ os.makedirs(os.path.dirname(filename), exist_ok=True)
+ tree = ET.ElementTree(vtk)
+ tree.write(filename)
+ print(f"[PVD] Written: {filename}")
+
+# ------------------------------------------------------------
+# Main program
+# ------------------------------------------------------------
+os.makedirs(output_dir, exist_ok=True)
+
+dx = (xmax - xmin) / nx
+x = np.linspace(xmin + 0.5*dx, xmax - 0.5*dx, nx)
+
+# Initial conditions
+rhol, ul, pl = 1.0, 0.0, 1
+rhor, ur, pr = 0.125, 0.0, 0.1
+
+vtu_files = []
+time_values = []
+
+# time loop at constant dt (sampling)
+for step in range(1, n_steps + 1):
+ t = step * dt
+
+ rho, u, p = sample_riemann(x, t, rhol, ul, pl, rhor, ur, pr, gamma, x0=x0)
+ E = p/(gamma - 1.0) + 0.5*rho*u**2
+
+ if (step % output_stride == 0) or (step == n_steps):
+ base = f"exact_{step:06d}.vtu"
+ fullpath = os.path.join(output_dir, base)
+ write_vtu(fullpath, x, rho, u, p, E)
+ vtu_files.append(base)
+ time_values.append(t)
+
+# PVD file
+pvd_path = os.path.join(output_dir, "exact.pvd")
+write_pvd(pvd_path, vtu_files, time_values)
diff --git a/15_Sod_Shock_Tube/POST/postprocessSodTube.py b/15_Sod_Shock_Tube/POST/postprocessSodTube.py
new file mode 100644
index 0000000..0de7d77
--- /dev/null
+++ b/15_Sod_Shock_Tube/POST/postprocessSodTube.py
@@ -0,0 +1,436 @@
+#!/usr/bin/env pvpython
+# -*- coding: utf-8 -*-
+"""
+ParaView post-processing (pvpython): 1D extraction and multi-resolution comparison.
+
+This script opens several results (code_saturne .case files and/or reference
+solutions in VTK .pvd/.vtu format), selects the available time step closest to
+a user-specified time t_user, and samples the physical fields along a segment
+P1->P2 using the ParaView PlotOverLine filter.
+
+The extracted 1D profiles (density rho, pressure p, velocity u_x) are then
+plotted with Matplotlib to compare different resolutions and/or an exact
+solution. Zoom insets and arrow annotations highlight the shock tube structures
+(rarefaction fan, contact surface, shock wave).
+
+Inputs:
+- code_saturne cases: CASE/RESU/*/postprocessing/RESULTS_FLUID_DOMAIN.case
+- Exact reference: .pvd/.vtu file provided in CASES via the "path" key
+
+Outputs:
+- High-resolution PNG figure
+
+Notes:
+- Must be run with pvpython (ParaView environment) for paraview.simple.
+- Explicit cleanup of ParaView objects (Delete) to limit memory usage.
+
+Usage:
+ pvpython postprocessSodTube.py
+
+Author: fhermet
+"""
+
+
+import os, glob
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+from paraview.simple import OpenDataFile, PlotOverLine, GetTimeKeeper, Delete
+from paraview.simple import servermanager
+from vtk.util.numpy_support import vtk_to_numpy
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes, mark_inset
+
+# ---------------------------
+# Matplotlib
+# ---------------------------
+mpl.rcParams.update({
+ "text.usetex": True,
+ "font.family": "serif",
+})
+
+# ---------------------------
+# Parameters
+# ---------------------------
+t_user = 0.2
+
+# sampling segment
+P1 = (-1.0, 0.0, 0.0)
+P2 = ( 1.0, 0.0, 0.0)
+
+ZOOMS = {
+ "rarefaction": (-0.1, 0.1),
+ "contact": (0.10, 0.25),
+ "contact_v": (0.10, 0.30),
+}
+
+ARROWS = {
+ "rho": dict(xy=(0.2, 0.32), xytext=(-0.3, 0.15), text=r"Contact surface", fontsize=12),
+ "p": dict(xy=(0.38, 0.2), xytext=(0.5, 0.3), text=r"Shock wave", fontsize=12),
+ "u": dict(xy=(-0.2, 0.20), xytext=(-0.75, 0.35),text=r"Rarefaction fan", fontsize=12),
+}
+
+CASES = [
+ dict(case_dir="../CASE", resolution=2000, label=r"$N=2000$", color="navy", ls="-", lw=1.4),
+ dict(path="exact.pvd", resolution=32000, label=r"Exact", color="black", ls="--", lw=1.4),
+]
+
+FIELD_NAMES = {
+ "P": ["Pressure", "pressure"],
+ "rho": ["Density", "rho"],
+ "V": ["Velocity", "velocity"],
+}
+
+# ---------------------------
+# Utilities
+# ---------------------------
+def nearest_time(timesteps, t):
+ """
+ Return the available time step closest to the requested time.
+
+ Parameters
+ ----------
+ timesteps : iterable of float
+ Available time steps from the ParaView reader.
+ t : float
+ User-requested time.
+
+ Returns
+ -------
+ float
+ Actual time used (nearest neighbor).
+ """
+ ts = np.asarray(list(timesteps), dtype=float)
+ return float(ts[np.argmin(np.abs(ts - float(t)))])
+
+def reader_timesteps(reader):
+ """
+ Retrieve the list of time steps from a ParaView reader.
+
+ Works with both:
+ - code_saturne files (.case)
+ - VTK files (.pvd, .vtu)
+
+ Parameters
+ ----------
+ reader : paraview.simple.Proxy
+ Already opened ParaView reader.
+
+ Returns
+ -------
+ list of float
+ Available time steps.
+ """
+ if getattr(reader, "TimestepValues", None):
+ return list(reader.TimestepValues)
+ return list(GetTimeKeeper().TimestepValues)
+
+def resolve_saturne_case(case_dir, case_filename="RESULTS_FLUID_DOMAIN.case"):
+ """
+ Automatically resolve the path to a code_saturne .case file.
+
+ Assumption:
+ case_dir/RESU/*/postprocessing/RESULTS_FLUID_DOMAIN.case
+
+ Parameters
+ ----------
+ case_dir : str
+ Root directory of the code_saturne computation.
+ case_filename : str
+ Name of the .case file (default: standard code_saturne name).
+
+ Returns
+ -------
+ str
+ Full path to the unique .case file found.
+
+ Raises
+ ------
+ FileNotFoundError
+ If no .case file is found.
+ RuntimeError
+ If multiple .case files are found.
+ """
+ pattern = os.path.join(case_dir, "RESU", "*", "postprocessing", case_filename)
+ matches = sorted(glob.glob(pattern))
+ if not matches:
+ raise FileNotFoundError(f"No .case file found via:\n {pattern}")
+ if len(matches) > 1:
+ raise RuntimeError("Multiple .case files found:\n " + "\n ".join(matches))
+ return matches[0]
+
+def get_np_any(dataset, names):
+ """
+ Extract a VTK field as a NumPy array.
+
+ Searches successively in:
+ - PointData
+ - CellData
+
+ Parameters
+ ----------
+ dataset : vtkDataSet
+ Data retrieved via servermanager.Fetch().
+ names : list of str
+ Possible field names (e.g. ["Pressure", "pressure"]).
+
+ Returns
+ -------
+ numpy.ndarray
+ NumPy array of the field found.
+
+ Raises
+ ------
+ RuntimeError
+ If no matching field is found.
+ """
+ pd, cd = dataset.GetPointData(), dataset.GetCellData()
+ for name in names:
+ arr = pd.GetArray(name) or cd.GetArray(name)
+ if arr is not None:
+ return vtk_to_numpy(arr)
+ raise RuntimeError(f"Array not found among: {names}")
+
+def sample_line(reader, line, t):
+ """
+ Sample physical fields along a line at a given time.
+
+ Extracted fields:
+ - abscissa (arc_length or x coordinate)
+ - velocity Ux
+ - pressure P
+ - density rho
+
+ Compatible with:
+ - code_saturne (.case)
+ - exact solutions (.pvd / .vtu)
+
+ Parameters
+ ----------
+ reader : paraview.simple.Proxy
+ ParaView reader for the case.
+ line : paraview.simple.Proxy
+ Configured PlotOverLine filter.
+ t : float
+ Time at which to sample.
+
+ Returns
+ -------
+ x : numpy.ndarray
+ Abscissa along the line.
+ Ux : numpy.ndarray
+ X-component of velocity.
+ P : numpy.ndarray
+ Pressure.
+ rho : numpy.ndarray
+ Density.
+ """
+ reader.UpdatePipeline(time=t)
+ line.UpdatePipeline(time=t)
+ data = servermanager.Fetch(line)
+
+ arc = data.GetPointData().GetArray("Point_X")
+ x = vtk_to_numpy(arc) if arc is not None else vtk_to_numpy(data.GetPoints().GetData())[:, 0]
+
+ P = get_np_any(data, FIELD_NAMES["P"])
+ rho = get_np_any(data, FIELD_NAMES["rho"])
+ V = get_np_any(data, FIELD_NAMES["V"])
+ Ux = V[:, 0] if (hasattr(V, "ndim") and V.ndim == 2) else V
+
+ return x, Ux, P, rho
+
+def add_zoom(ax, x_ref, y_ref, xlim, ylim=None, loc="upper right",
+ width="38%", height="38%", borderpad=1.0):
+ """
+ Add a zoom inset on an existing matplotlib axis.
+
+ All curves already plotted on the main axis are automatically
+ copied into the inset with the same style.
+
+ Parameters
+ ----------
+ ax : matplotlib.axes.Axes
+ Main axis.
+ x_ref : numpy.ndarray
+ Reference abscissa (used for auto-ylim).
+ y_ref : numpy.ndarray
+ Reference field for computing zoom bounds.
+ xlim : tuple(float, float)
+ X-bounds of the zoom.
+ ylim : tuple(float, float) or None
+ Y-bounds of the zoom (auto if None).
+ loc : str
+ Inset position.
+ width, height : str
+ Relative size of the inset.
+ borderpad : float
+ Padding between inset and main axis.
+
+ Returns
+ -------
+ matplotlib.axes.Axes
+ Created inset axis.
+ """
+ axins = inset_axes(ax, width=width, height=height, loc=loc, borderpad=borderpad)
+
+ # copy all curves already plotted (same styles)
+ for line in ax.get_lines():
+ axins.plot(line.get_xdata(), line.get_ydata(),
+ color=line.get_color(), ls=line.get_linestyle(), lw=line.get_linewidth())
+
+ axins.set_xlim(*xlim)
+
+ if ylim is None:
+ m = (x_ref >= xlim[0]) & (x_ref <= xlim[1])
+ yy = y_ref[m]
+ if yy.size:
+ axins.set_ylim(yy.min() - 0.01, yy.max() + 0.01)
+ else:
+ axins.set_ylim(*ylim)
+
+ axins.grid(True, alpha=0.25)
+ axins.tick_params(labelsize=8)
+ mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.3", lw=0.8)
+ return axins
+
+def add_arrow(ax, xy, xytext, text, fontsize=12):
+ """
+ Add an arrow annotation on a plot.
+
+ Parameters
+ ----------
+ ax : matplotlib.axes.Axes
+ Target axis.
+ xy : tuple(float, float)
+ Point the arrow points to.
+ xytext : tuple(float, float)
+ Text position.
+ text : str
+ Displayed text.
+ fontsize : int
+ Font size.
+ """
+ ax.annotate(
+ text, xy=xy, xytext=xytext, textcoords="data", xycoords="data",
+ ha="left", va="center", fontsize=fontsize,
+ arrowprops=dict(arrowstyle="->", lw=1.2),
+ bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.8),
+ )
+
+# ---------------------------
+# Extraction (case loop)
+# ---------------------------
+results = []
+
+for k, c in enumerate(CASES):
+ # ----------------------------------------------------------------
+ # 1) Resolve the file path to open
+ # - either an exact file (.pvd/.vtu) provided directly
+ # - or a code_saturne computation (.case) resolved automatically
+ # ----------------------------------------------------------------
+ filepath = c.get("path") or resolve_saturne_case(c["case_dir"])
+ print(f"[{k}] Opening: {filepath}")
+
+ # ----------------------------------------------------------------
+ # 2) Create ParaView reader and query metadata
+ # ----------------------------------------------------------------
+ reader = OpenDataFile(filepath)
+ reader.UpdatePipelineInformation()
+
+ # ----------------------------------------------------------------
+ # 3) Retrieve available time steps
+ # ----------------------------------------------------------------
+ ts = reader_timesteps(reader)
+ if not ts:
+ raise RuntimeError(f"No time steps detected in {filepath}")
+
+ # ----------------------------------------------------------------
+ # 4) Select the time step closest to the requested one
+ # (important if t_user does not match an output exactly)
+ # ----------------------------------------------------------------
+ t = nearest_time(ts, t_user)
+ print(f"Requested time={t_user} -> used={t}")
+
+ # ----------------------------------------------------------------
+ # 5) Create PlotOverLine filter
+ # - line defined by P1 -> P2
+ # - resolution controlled by the case (mesh / reference)
+ # ----------------------------------------------------------------
+ line = PlotOverLine(Input=reader,
+ Point1=list(P1),
+ Point2=list(P2),
+ Resolution=int(c["resolution"]))
+
+ # ----------------------------------------------------------------
+ # 6) Sample physical fields along the line
+ # - abscissa x (arc length or coordinate)
+ # - velocity Ux
+ # - pressure P
+ # - density rho
+ # ----------------------------------------------------------------
+ x, Ux, P, rho = sample_line(reader, line, t)
+
+ # ----------------------------------------------------------------
+ # 7) Store results for later plotting
+ # (plotting is intentionally separated from extraction)
+ # ----------------------------------------------------------------
+ results.append(dict(case=c, t=t, x=x, Ux=Ux, P=P, rho=rho))
+
+ # ----------------------------------------------------------------
+ # 8) Explicit cleanup of ParaView objects
+ # (important to avoid memory leaks with pvpython)
+ # ----------------------------------------------------------------
+ Delete(line)
+ Delete(reader)
+
+# Reference for auto-ylim of zooms
+x_ref = results[0]["x"]
+Ux_ref = results[0]["Ux"]
+P_ref = results[0]["P"]
+rho_ref = results[0]["rho"]
+
+# ---------------------------
+# Plot
+# ---------------------------
+fig, (axR, axP, axU) = plt.subplots(3, 1, sharex=True, figsize=(10, 14))
+
+for r in results:
+ c = r["case"]
+ style = dict(color=c.get("color"), ls=c.get("ls", "-"), lw=c.get("lw", 1.4), label=c.get("label", "case"))
+ axU.plot(r["x"], r["Ux"], **style)
+ axP.plot(r["x"], r["P"], **style)
+ axR.plot(r["x"], r["rho"],**style)
+
+# Zooms
+add_zoom(axR, x_ref, rho_ref, ZOOMS["contact"], loc="upper right", borderpad=1.5)
+add_zoom(axP, x_ref, P_ref, ZOOMS["rarefaction"], loc="lower left", width="22%", height="22%", borderpad=2.5)
+add_zoom(axU, x_ref, Ux_ref, ZOOMS["contact_v"], loc="upper right", width="22%", height="22%", borderpad=2.5)
+
+# Arrows
+add_arrow(axR, **ARROWS["rho"])
+add_arrow(axP, **ARROWS["p"])
+add_arrow(axU, **ARROWS["u"])
+
+# Labels
+axU.set_ylabel(r"$u_x$", fontsize=16)
+axP.set_ylabel(r"$p$", fontsize=16)
+axR.set_ylabel(r"$\rho$",fontsize=16)
+
+for ax in (axR, axP, axU):
+ ax.set_xlabel(r"$x$", fontsize=16)
+ ax.tick_params(labelsize=11)
+ ax.grid(True, alpha=0.3)
+
+# Single legend without duplicates (on axP)
+handles, labels = axP.get_legend_handles_labels()
+uniq = {}
+for h, lab in zip(handles, labels):
+ uniq.setdefault(lab, h)
+axP.legend(uniq.values(), uniq.keys(), fancybox=True, frameon=True, framealpha=0.9, fontsize=12)
+
+fig.suptitle(rf"Comparison at $t = {results[0]['t']:.2f}\,\mathrm{{s}}$", fontsize=20, y=0.96)
+fig.tight_layout(rect=[0, 0, 1, 0.96])
+
+out = f"results_t{t_user:.3f}.png"
+fig.savefig(out, dpi=300, bbox_inches="tight")
+print(f"Figure saved: {out}")
diff --git a/15_Sod_Shock_Tube/README.md b/15_Sod_Shock_Tube/README.md
new file mode 100644
index 0000000..5458936
--- /dev/null
+++ b/15_Sod_Shock_Tube/README.md
@@ -0,0 +1,48 @@
+Sod Shock Tube
+==============
+
+Classical 1D Riemann problem (Sod, 1978) solved with the compressible
+module of code_saturne (constant gamma, Euler equations).
+
+A diaphragm at x = 0 separates a high-pressure region (left) from a
+low-pressure region (right). Removing the diaphragm produces three
+characteristic wave structures: a left-moving rarefaction fan, a contact
+discontinuity, and a right-moving shock wave.
+
+Setup
+-----
+
+- Compressible model: constant gamma (gamma = 1.4)
+- Turbulence: off (inviscid)
+- Mesh: 2000 x 1 x 1 Cartesian (generated at runtime), x in [-1, 1]
+- Boundary conditions: symmetry on all faces
+- Time stepping: fixed dt, CFL ~ 0.54, 436 iterations (t_final ~ 0.2 s)
+
+Initial conditions:
+
+| Region | Density | Pressure | Velocity |
+|--------|---------|----------|----------|
+| Left | 1.0 | 1.0 | 0.0 |
+| Right | 0.125 | 0.1 | 0.0 |
+
+The user-defined file `CASE/SRC/cs_user_physical_properties.cpp` forces
+molecular viscosity, volume viscosity, and thermal conductivity to zero
+for a true inviscid (Euler) simulation.
+
+Post-processing
+---------------
+
+Scripts in `POST/` require pvpython (ParaView environment):
+
+- `exactSodTube.py`: generates the exact Riemann solution in VTU/PVD format.
+- `postprocessSodTube.py`: extracts 1D profiles from code_saturne results
+ and compares them against the exact solution.
+
+Reference
+---------
+
+G. A. Sod, "A survey of several finite difference methods for systems of
+nonlinear hyperbolic conservation laws", J. Comput. Phys., 27(1):1-31, 1978.
+
+E. F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics",
+3rd Edition, Springer, 2009.
diff --git a/15_Sod_Shock_Tube/smgr.xml b/15_Sod_Shock_Tube/smgr.xml
new file mode 100644
index 0000000..6453908
--- /dev/null
+++ b/15_Sod_Shock_Tube/smgr.xml
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/16_Sedov_Blast_Wave/CASE/DATA/run.cfg b/16_Sedov_Blast_Wave/CASE/DATA/run.cfg
new file mode 100644
index 0000000..89ecabe
--- /dev/null
+++ b/16_Sedov_Blast_Wave/CASE/DATA/run.cfg
@@ -0,0 +1,5 @@
+[job_defaults]
+
+n_procs: 1
+n_threads: 1
+
diff --git a/16_Sedov_Blast_Wave/CASE/DATA/setup.xml b/16_Sedov_Blast_Wave/CASE/DATA/setup.xml
new file mode 100644
index 0000000..ed39457
--- /dev/null
+++ b/16_Sedov_Blast_Wave/CASE/DATA/setup.xml
@@ -0,0 +1,238 @@
+
+
+
+
+
+
+
+
+ 1
+ 0.5
+
+
+
+
+
+
+ 1
+ 10
+ 0.1
+ 5e-05
+ 0.1
+
+
+
+ X0
+ X1
+ Y0
+ Y1
+ Z0
+ Z1
+
+
+
+
+
+ 0
+
+
+
+
+
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ density = 1;
+ 1.17862
+
+
+ 0.0
+
+
+
+
+ 0
+
+
+
+
+ 1005
+
+
+
+
+ 0
+
+
+
+
+ 0
+
+
+
+ 0.028966
+ 101325
+ 0.00348310693138
+
+
+ 0
+ 0
+ 0
+
+
+
+ 0
+ 0
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ all[]
+
+
+
+
+
+
+
+
+
+
+
+
+ 1
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 1.0
+ off
+
+
+
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+ 1e+12
+ -1e+12
+ 1
+ SGDH
+
+
+
+
+
+
+
+
+ 1
+ 0
+
+
+
+ velocity[0] = 0.;
+velocity[1] = 0.;
+velocity[2] = 0.;
+
+
+
+
+
+
+
+
+
+
+
+
+ pressure = 1;
+ 2
+
+
+
+ 1
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/16_Sedov_Blast_Wave/CASE/SRC/cs_user_initialization.cpp b/16_Sedov_Blast_Wave/CASE/SRC/cs_user_initialization.cpp
new file mode 100644
index 0000000..dff4d6f
--- /dev/null
+++ b/16_Sedov_Blast_Wave/CASE/SRC/cs_user_initialization.cpp
@@ -0,0 +1,130 @@
+/*============================================================================
+ * User initialization for the Sedov blast wave problem.
+ *
+ * A high-pressure region (r < r0) surrounded by a low-pressure ambient
+ * creates a cylindrical shock wave propagating radially outward.
+ *============================================================================*/
+
+/* code_saturne version 9.1 */
+
+/*
+ This file is part of code_saturne, a general-purpose CFD tool.
+
+ Copyright (C) 1998-2025 EDF S.A.
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
+ Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+/*----------------------------------------------------------------------------*/
+
+#include "cs_defs.h"
+
+/*----------------------------------------------------------------------------
+ * Standard C library headers
+ *----------------------------------------------------------------------------*/
+
+#include
+#include
+
+#if defined(HAVE_MPI)
+#include
+#endif
+
+/*----------------------------------------------------------------------------
+ * PLE library headers
+ *----------------------------------------------------------------------------*/
+
+#include
+
+/*----------------------------------------------------------------------------
+ * Local headers
+ *----------------------------------------------------------------------------*/
+
+#include "cs_headers.h"
+
+/*----------------------------------------------------------------------------*/
+
+BEGIN_C_DECLS
+
+/*----------------------------------------------------------------------------*/
+/*!
+ * \file cs_user_initialization.cpp
+ *
+ * \brief Initialization prior to solving time steps.
+ */
+/*----------------------------------------------------------------------------*/
+
+/*============================================================================
+ * User function definitions
+ *============================================================================*/
+
+/*----------------------------------------------------------------------------*/
+/*!
+ * \brief Define initial conditions for variables.
+ *
+ * This function is not yet called by code_saturne, use
+ * `cs_user_initialization` for the time being
+ * This function is called before the beginning of the computation
+ * allowing an overload of the GUI defined initialization (called just
+ * after cs_gui_initial_conditions).
+ *
+ * \param[in, out] domain pointer to a cs_domain_t structure
+ */
+/*----------------------------------------------------------------------------*/
+
+void
+cs_user_initialization(cs_domain_t *domain)
+{
+ const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
+ const cs_real_3_t *cell_cen =
+ (const cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;
+
+ /* Access field pointers */
+ cs_real_t *rho = CS_F_(rho)->val;
+ cs_real_t *p = CS_F_(p)->val;
+ cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
+
+ /* Sedov blast wave parameters:
+ * rho0 - uniform ambient density
+ * pinfty - ambient pressure (very low)
+ * phot - high pressure inside the blast region
+ * r0 - initial blast radius */
+ const cs_real_t rho0 = 1.0;
+ const cs_real_t pinfty = 0.01;
+ const cs_real_t phot = 100.0;
+ const cs_real_t r0 = 0.005;
+
+ /* Loop over all cells: set uniform density and zero velocity,
+ * then apply a pressure discontinuity at radius r0 */
+ for (cs_lnum_t c = 0; c < n_cells; c++) {
+ double x = cell_cen[c][0];
+ double y = cell_cen[c][1];
+ double r = sqrt(x*x + y*y);
+
+ rho[c] = rho0;
+ vel[c][0] = 0.0;
+ vel[c][1] = 0.0;
+ vel[c][2] = 0.0;
+
+ if (r < r0)
+ p[c] = phot; /* hot region: energy source */
+ else
+ p[c] = pinfty; /* ambient: near-vacuum pressure */
+ }
+}
+
+/*----------------------------------------------------------------------------*/
+
+END_C_DECLS
diff --git a/16_Sedov_Blast_Wave/README.md b/16_Sedov_Blast_Wave/README.md
new file mode 100644
index 0000000..5a31c7c
--- /dev/null
+++ b/16_Sedov_Blast_Wave/README.md
@@ -0,0 +1,38 @@
+Sedov Blast Wave
+================
+
+2D Sedov blast wave problem solved with the compressible module of
+code_saturne (constant gamma, Euler equations).
+
+A point-like energy release at the origin creates a cylindrical shock
+wave that propagates radially. The shock radius follows the self-similar
+Sedov-Taylor solution: R(t) proportional to t^(2/5) in 2D (t^(1/2) for
+the cylindrical symmetry used here with a quarter-domain).
+
+Setup
+-----
+
+- Compressible model: constant gamma (gamma = 1.4)
+- Turbulence: off (inviscid)
+- Mesh: 200 x 200 x 1 Cartesian (generated at runtime), domain [0, 0.2]^2
+- Boundary conditions:
+ - x0, y0: symmetry (quarter-domain exploiting axial symmetry)
+ - x1, y1: supersonic outlet
+ - z0, z1: symmetry (2D)
+- Time stepping: adaptive (CFL = 1), dt_ref = 5e-5 s, t_final = 0.5 s
+
+Initial conditions (set in `cs_user_initialization.cpp`):
+
+| Region | Density | Pressure | Velocity |
+|-----------|---------|----------|----------|
+| r < 0.005 | 1.0 | 100.0 | 0.0 |
+| r >= 0.005| 1.0 | 0.01 | 0.0 |
+
+Reference
+---------
+
+L. I. Sedov, "Similarity and Dimensional Methods in Mechanics",
+Academic Press, 1959.
+
+G. I. Taylor, "The formation of a blast wave by a very intense explosion",
+Proc. R. Soc. Lond. A, 201(1065):159-186, 1950.
diff --git a/16_Sedov_Blast_Wave/smgr.xml b/16_Sedov_Blast_Wave/smgr.xml
new file mode 100644
index 0000000..9deb963
--- /dev/null
+++ b/16_Sedov_Blast_Wave/smgr.xml
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+