From 4a938993d222913734cb428dae48a547df56817e Mon Sep 17 00:00:00 2001 From: fhermet Date: Tue, 14 Apr 2026 12:57:26 +0200 Subject: [PATCH 1/2] 15_Sod_Shock_Tube: add compressible Sod shock tube tutorial --- 15_Sod_Shock_Tube/CASE/DATA/run.cfg | 5 + 15_Sod_Shock_Tube/CASE/DATA/setup.xml | 242 ++++++++++ 15_Sod_Shock_Tube/POST/exactSodTube.py | 461 +++++++++++++++++++ 15_Sod_Shock_Tube/POST/postprocessSodTube.py | 436 ++++++++++++++++++ 15_Sod_Shock_Tube/README.md | 48 ++ 15_Sod_Shock_Tube/smgr.xml | 7 + 6 files changed, 1199 insertions(+) create mode 100644 15_Sod_Shock_Tube/CASE/DATA/run.cfg create mode 100644 15_Sod_Shock_Tube/CASE/DATA/setup.xml create mode 100644 15_Sod_Shock_Tube/POST/exactSodTube.py create mode 100644 15_Sod_Shock_Tube/POST/postprocessSodTube.py create mode 100644 15_Sod_Shock_Tube/README.md create mode 100644 15_Sod_Shock_Tube/smgr.xml 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 @@ + + + + + + + + + all[] + + + + + all[] + + + + 1 + + + + + + + + + + + + + + + 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 From 698a6e7b5ffb923864de64c8ddcf969ca49ac5c2 Mon Sep 17 00:00:00 2001 From: fhermet Date: Tue, 14 Apr 2026 14:58:16 +0200 Subject: [PATCH 2/2] 16_Sedov_Blast_Wave: add compressible Sedov blast wave tutorial --- 16_Sedov_Blast_Wave/CASE/DATA/run.cfg | 5 + 16_Sedov_Blast_Wave/CASE/DATA/setup.xml | 238 ++++++++++++++++++ .../CASE/SRC/cs_user_initialization.cpp | 130 ++++++++++ 16_Sedov_Blast_Wave/README.md | 38 +++ 16_Sedov_Blast_Wave/smgr.xml | 7 + 5 files changed, 418 insertions(+) create mode 100644 16_Sedov_Blast_Wave/CASE/DATA/run.cfg create mode 100644 16_Sedov_Blast_Wave/CASE/DATA/setup.xml create mode 100644 16_Sedov_Blast_Wave/CASE/SRC/cs_user_initialization.cpp create mode 100644 16_Sedov_Blast_Wave/README.md create mode 100644 16_Sedov_Blast_Wave/smgr.xml 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 @@ + + + + + + + all[] + + + + + all[] + + + + 1 + + + + + + 10 + + + + + + + + + + 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 @@ + + + + + + +