{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Shock Propagation in Pressureless SPH\n",
    "\n",
    "In this tutorial, you set up and run a **1D Riemann shock benchmark** with the pressureless SPH model in Struphy, comparing numerical shock propagation to the analytical Rankine-Hugoniot weak solution.\n",
    "\n",
    "The notebook focuses on the `PressureLessSPH` model for coupled **density and momentum** evolution:\n",
    "\n",
    "$$\n",
    "\\partial_t \\rho + \\partial_x(\\rho u) = 0,\n",
    "\\qquad\n",
    "\\partial_t(\\rho u) + \\partial_x(\\rho u^2) = 0.\n",
    "$$\n",
    "\n",
    "A short contextual note: in near-constant-density limits, the momentum dynamics are often discussed in relation to inviscid Burgers behavior. Here, the benchmark of interest is the full pressureless Euler Riemann shock.\n",
    "\n",
    "The workflow demonstrates core Struphy concepts: environment setup, particle initialization, time integration, and diagnostics in an interactive format so parameters can be tuned and rerun quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Mathematical Background: Pressureless Euler Shock System\n",
    "\n",
    "The pressureless continuity and momentum equations form a hyperbolic conservation-law system:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  \\partial_t \\rho + \\partial_x(\\rho u) &= 0 \\quad \\text{(Continuity)} \\\\\n",
    "  \\partial_t(\\rho u) + \\partial_x(\\rho u^2) &= 0 \\quad \\text{(Momentum without pressure)}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "For discontinuous initial data, the system develops weak solutions with shocks.\n",
    "\n",
    "### Riemann Problem (Shock Benchmark)\n",
    "\n",
    "We solve the Riemann problem with initial conditions:\n",
    "\n",
    "$$\n",
    "(\\rho, u)(x, 0) = \\begin{cases}\n",
    "  (\\rho_L, u_L) & \\text{if } x < 0.5 \\\\\n",
    "  (\\rho_R, u_R) & \\text{if } x > 0.5\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The **Rankine-Hugoniot shock speed** is:\n",
    "\n",
    "$$\n",
    "s = \\frac{u_L + u_R}{2}.\n",
    "$$\n",
    "\n",
    "This provides a standard and well-established 1D benchmark for validating shock propagation in particle methods."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Step 0: Import Struphy Components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import cunumpy as xp\n",
    "\n",
    "# Struphy imports\n",
    "from struphy import (\n",
    "    BaseUnits,\n",
    "    EnvironmentOptions,\n",
    "    Time,\n",
    "    domains,\n",
    "    equils,\n",
    "    grids,\n",
    "    DerhamOptions,\n",
    "    BoundaryParameters,\n",
    "    LoadingParameters,\n",
    "    WeightsParameters,\n",
    "    SortingParameters,\n",
    "    SavingParameters,\n",
    "    BinningPlot,\n",
    "    Simulation,\n",
    ")\n",
    "from struphy.models import PressureLessSPH\n",
    "from struphy.initial.base import Perturbation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Step 1: Create Environment, Time Integrator, and 1D Domain\n",
    "\n",
    "Set up a 1D periodic domain (stretched into 3D as required by Struphy) with Strang splitting for time integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ====== CONFIGURABLE PARAMETERS ======\n",
    "# Time stepping\n",
    "dt = 1.0e-2      # time step\n",
    "Tend = 1.5       # end time\n",
    "\n",
    "# Domain (1D periodic in logical coords eta1 ∈ [0, 1])\n",
    "l1, r1 = 0.0, 1.0  # eta1 range\n",
    "l2, r2 = 0.0, 1.0  # eta2 range (minimal extent)\n",
    "l3, r3 = 0.0, 1.0  # eta3 range (minimal extent)\n",
    "# =====================================\n",
    "\n",
    "# Environment options\n",
    "env = EnvironmentOptions(sim_folder=\"sim_shock_large\")\n",
    "\n",
    "# Time stepping with Strang splitting\n",
    "time_opts = Time(dt=dt, Tend=Tend, split_algo=\"Strang\")\n",
    "\n",
    "# Geometry: 1D periodic, extended to 3D cuboid\n",
    "domain = domains.Cuboid(l1=l1, r1=r1, l2=l2, r2=r2, l3=l3, r3=r3)\n",
    "\n",
    "print(f\"Domain: eta1 ∈ [{l1}, {r1}], eta2 ∈ [{l2}, {r2}], eta3 ∈ [{l3}, {r3}]\")\n",
    "print(f\"Time stepping: dt={dt}, Tend={Tend}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Step 2: Define Riemann Step for Shock Initial Condition\n",
    "\n",
    "We implement a smooth approximation to the Riemann jump using a tanh transition. This avoids sharp discontinuities in the initial particle distribution while maintaining the shock structure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Comparing Two Riemann Shock Test Cases\n",
    "\n",
    "This notebook includes two distinct test cases to explore pressureless Euler shock dynamics:\n",
    "\n",
    "**Test Case 1: Contact Discontinuity (Rankine-Hugoniot Exact)**  \n",
    "- Initial conditions: $\\rho_L = 1.0, u_L = 0.5$ | $\\rho_R = 0.5, u_R = 0.5$\n",
    "- Shock speed: $s = 0$ (stationary, since $u_L = u_R$)\n",
    "- Validates analytical solution: piecewise-constant profile with **no propagation**\n",
    "\n",
    "**Test Case 2: Delta-Shock (Density Pile-Up)**  \n",
    "- Initial conditions: $\\rho_L = 1.0, u_L = 2.0$ | $\\rho_R = 0.125, u_R = 0.0$\n",
    "- Shock speed: naive estimate $s \\approx 1.0$ (fast-moving compression)\n",
    "- Demonstrates **physical density concentration** at the shock front (SPH forms delta-shock)\n",
    "\n",
    "**To switch test cases**: In the code cell below, uncomment/comment the `test_case` variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "class RiemannStep(Perturbation):\n",
    "    \"\"\"Smooth approximation of a 1D Riemann jump in logical coordinate eta1.\n",
    "    \n",
    "    Args:\n",
    "        left: value on the left side  (x < eta0)\n",
    "        right: value on the right side (x > eta0)\n",
    "        eta0: location of the jump (default 0.5 = middle of domain)\n",
    "        width: width of tanh transition (smaller = sharper transition)\n",
    "        given_in_basis: basis for given variable (\"0\" for density, \"v\" for velocity)\n",
    "        comp: component index for vector quantities\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, left: float, right: float, eta0: float = 0.5, \n",
    "                 width: float = 0.01, given_in_basis: str=\"0\", comp: int = 0):\n",
    "        self.left = left\n",
    "        self.right = right\n",
    "        self.eta0 = eta0\n",
    "        self.width = width\n",
    "        self.given_in_basis = given_in_basis\n",
    "        self.comp = comp\n",
    "\n",
    "    def __call__(self, e1, e2, e3):\n",
    "        \"\"\"Evaluate the smooth step at logical coordinates (e1, e2, e3).\"\"\"\n",
    "        avg = 0.5 * (self.left + self.right)\n",
    "        half_jump = 0.5 * (self.left - self.right)\n",
    "        return avg - half_jump * xp.tanh((e1 - self.eta0) / self.width)\n",
    "\n",
    "\n",
    "# ========== TEST CASE SELECTION ==========\n",
    "# Choose which test case to run:\n",
    "# test_case = \"contact_discontinuity\"  # Satisfies Rankine-Hugoniot exactly\n",
    "test_case = \"delta_shock\"             # Shows density pile-up at shock front\n",
    "test_case = \"contact_discontinuity\"\n",
    "# =========================================\n",
    "\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    # TEST CASE 1: Contact Discontinuity (satisfies Rankine-Hugoniot)\n",
    "    # This is a discontinuity where velocity is equal but density jumps.\n",
    "    # Verifies RH: s = 0.5, mass flux: 0.25, momentum: 0.125 on both sides ✓\n",
    "    rho_L, u_L = 1.0, 0.5\n",
    "    rho_R, u_R = 0.5, 0.5\n",
    "    \n",
    "    shock_speed_theory = (rho_R * u_R - rho_L * u_L) / (rho_R - rho_L)\n",
    "    \n",
    "    print(\"\\n\" + \"=\"*70)\n",
    "    print(\"TEST CASE 1: Contact Discontinuity (Rankine-Hugoniot satisfying)\")\n",
    "    print(\"=\"*70)\n",
    "    print(f\"  Left state:  ρ_L = {rho_L}, u_L = {u_L}, momentum m_L = {rho_L * u_L}\")\n",
    "    print(f\"  Right state: ρ_R = {rho_R}, u_R = {u_R}, momentum m_R = {rho_R * u_R}\")\n",
    "    print(f\"  Shock speed (from RH): s = {shock_speed_theory:.4f}\")\n",
    "    print(f\"  Note: u_L = u_R = {u_L} (contact discontinuity, no compression)\")\n",
    "    print(\"=\"*70)\n",
    "    \n",
    "else:  # delta_shock\n",
    "    # TEST CASE 2: Delta-Shock (high-to-low density with velocity change)\n",
    "    # Initial conditions do NOT satisfy RH initially, but SPH develops\n",
    "    # a delta-shock solution with density pile-up at the shock front.\n",
    "    rho_L, u_L = 1.0, 2.0\n",
    "    rho_R, u_R = 0.125, 0.0\n",
    "    \n",
    "    # Naive shock speed (does NOT satisfy RH for these parameters)\n",
    "    naive_speed = (u_L + u_R) / 2\n",
    "    \n",
    "    print(\"\\n\" + \"=\"*70)\n",
    "    print(\"TEST CASE 2: Delta-Shock (density pile-up at shock front)\")\n",
    "    print(\"=\"*70)\n",
    "    print(f\"  Left state:  ρ_L = {rho_L}, u_L = {u_L}, momentum m_L = {rho_L * u_L}\")\n",
    "    print(f\"  Right state: ρ_R = {rho_R}, u_R = {u_R}, momentum m_R = {rho_R * u_R}\")\n",
    "    print(f\"  Naive shock speed: s ≈ {naive_speed:.4f}\")\n",
    "    print(\"  Note: These parameters lead to SPH forming a delta-shock solution\")\n",
    "    print(\"        where mass concentrates at the shock front (density pile-up).\")\n",
    "    print(\"=\"*70)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "## Step 3: Configure Grid and de Rham Discretization\n",
    "\n",
    "Set up the tensor-product grid and de Rham options for the 1D problem. Since this is truly 1D physics in a 3D logical domain, we use minimal resolution in non-1D directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grid: fine resolution in eta1, minimal in others\n",
    "grid = grids.TensorProductGrid(num_elements=(64, 1, 1))\n",
    "\n",
    "# de Rham options: periodic boundaries in all directions\n",
    "derham_opts = DerhamOptions(degree=(3, 1, 1))\n",
    "\n",
    "print(\"Grid elements: (64, 1, 1)\")\n",
    "print(\"Boundary conditions: periodic in all directions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## Step 4: Instantiate the PressureLessSPH Model\n",
    "\n",
    "Create a model instance for the Riemann shock benchmark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model instance for the pressureless Euler shock system\n",
    "model = PressureLessSPH()\n",
    "\n",
    "# Simulation object\n",
    "sim = Simulation(\n",
    "    model,\n",
    "    env=env,\n",
    "    time_opts=time_opts,\n",
    "    domain=domain,\n",
    "    equil=None,  # No background equilibrium needed\n",
    "    grid=grid,\n",
    "    derham_opts=derham_opts,\n",
    ")\n",
    "\n",
    "print(\"PressureLessSPH model instantiated (no external field).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## Step 5: Configure Particle Markers with SPH Diagnostics\n",
    "\n",
    "Set up particle loading, weights, boundaries, and binning diagnostics for the Riemann shock benchmark:\n",
    "\n",
    "- **Np = 10000** particles by default (can be tuned)\n",
    "- **Periodic boundaries** in all directions (1D domain)\n",
    "- **64 bins** along $\\eta_1$\n",
    "- Binned outputs for both **density** and **current_1** (used to reconstruct velocity snapshots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ====== CONFIGURABLE PARTICLE PARAMETERS ======\n",
    "Np_particles = 20000  # Number of particles\n",
    "n_bins = 64        # Number of bins for diagnostics\n",
    "# ==============================================\n",
    "\n",
    "loading_params = LoadingParameters(Np=Np_particles)\n",
    "weights_params = WeightsParameters()\n",
    "\n",
    "# Periodic boundaries in all directions (1D periodic domain)\n",
    "boundary_params = BoundaryParameters(\n",
    "    bc=(\"periodic\", \"periodic\", \"periodic\"),\n",
    "    bc_sph=(\"periodic\", \"periodic\", \"periodic\")\n",
    ")\n",
    "\n",
    "# Sorting: n_bins boxes in eta1 direction, minimal in others\n",
    "sorting_params = SortingParameters(\n",
    "    boxes_per_dim=(n_bins, 1, 1),\n",
    "    dims_mask=(True, False, False)  # Only sort in eta1\n",
    ")\n",
    "\n",
    "# Binning diagnostics for density and current in eta1\n",
    "bin_plot_density = BinningPlot(\n",
    "    slice=\"e1\",\n",
    "    n_bins=(n_bins,),\n",
    "    ranges=(0.0, 1.0),\n",
    "    output_quantity=\"density\"\n",
    ")\n",
    "\n",
    "bin_plot_current_1 = BinningPlot(\n",
    "    slice=\"e1\",\n",
    "    n_bins=(n_bins,),\n",
    "    ranges=(0.0, 1.0),\n",
    "    output_quantity=\"current_1\"\n",
    ")\n",
    "\n",
    "saving_params = SavingParameters(binning_plots=(bin_plot_density, bin_plot_current_1))\n",
    "\n",
    "model.cold_fluid.set_markers(\n",
    "    loading_params=loading_params,\n",
    "    weights_params=weights_params,\n",
    "    boundary_params=boundary_params,\n",
    "    sorting_params=sorting_params,\n",
    "    saving_params=saving_params,\n",
    ")\n",
    "\n",
    "print(f\"Particle setup: Np={Np_particles}, {n_bins} bins for diagnostics\")\n",
    "print(\"Binned outputs: e1_density and e1_current_1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Step 6: Set Propagator Options (No External Field)\n",
    "\n",
    "Configure the time integrators for particle position and velocity updates. **No external potential is applied** ($\\phi = 0$), so the benchmark is the unforced pressureless Euler shock system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "from struphy import ButcherTableau\n",
    "\n",
    "# Forward Euler time integration for both position and velocity\n",
    "butcher = ButcherTableau(algo=\"forward_euler\")\n",
    "model.propagators.push_eta.options = model.propagators.push_eta.Options(butcher=butcher)\n",
    "\n",
    "# No external field (phi=None means zero forcing)\n",
    "\n",
    "print(\"Propagators configured: forward_euler, no external field\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "# Riemann Shock Benchmark\n",
    "\n",
    "In this section, we set up and run the large-amplitude Riemann shock, then compare numerical and analytical shock propagation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Step 7: Apply Riemann Shock Initial Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constant background (zero velocity, unit density)\n",
    "background = equils.ConstantVelocity(ux=0.0, uy=0.0, uz=0.0, n=1.0, p0=0.0)\n",
    "model.cold_fluid.var.add_background(background)\n",
    "\n",
    "# Perturbations: density and velocity jumps (both smooth tanh transitions)\n",
    "del_n = RiemannStep(\n",
    "    left=rho_L - 1.0,      # density perturbation on left\n",
    "    right=rho_R - 1.0,     # density perturbation on right\n",
    "    eta0=0.5,\n",
    "    width=0.01\n",
    ")\n",
    "\n",
    "del_u1 = RiemannStep(\n",
    "    left=u_L,          # velocity on left\n",
    "    right=u_R,         # velocity on right\n",
    "    eta0=0.5,\n",
    "    width=0.01,\n",
    "    given_in_basis=\"v\"  # given in velocity basis\n",
    ")\n",
    "\n",
    "model.cold_fluid.var.add_perturbation(del_n=del_n, del_u1=del_u1)\n",
    "\n",
    "print(\"Riemann shock initial conditions applied:\")\n",
    "print(\"  Background: ρ=1.0, u=(0, 0, 0)\")\n",
    "print(\"  Perturbation: smooth tanh transition with width=0.01\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## Step 8: Run Riemann Shock Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n\" + \"=\"*60)\n",
    "print(\"Running LARGE-AMPLITUDE RIEMANN SHOCK simulation...\")\n",
    "print(\"=\"*60)\n",
    "out = sim.run()\n",
    "print(\"Simulation completed.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## Step 9: Access Results\n",
    "\n",
    "Field data is post-processed lazily, the first time it is accessed through the `out` object returned by `run`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "## Step 10: Visualize Density and Velocity Evolution\n",
    "\n",
    "Plot density and reconstructed velocity snapshots at multiple times to track shock propagation in both conserved fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract binned outputs\n",
    "rho_binned = np.asarray(out.evaluate(\"cold_fluid/f\", dataset=\"e1_density/f\"))\n",
    "current1_binned = np.asarray(out.evaluate(\"cold_fluid/f\", dataset=\"e1_current_1/f\"))\n",
    "t_grid = out.time\n",
    "eta1_bins = np.linspace(0, 1, n_bins + 1)[:-1]  # bin centers\n",
    "\n",
    "# Reconstruct velocity from binned current and density: u1 = j1 / rho\n",
    "rho_floor = 1.0e-12\n",
    "u1_binned = np.divide(current1_binned, np.maximum(rho_binned, rho_floor))\n",
    "\n",
    "# Plot 6 snapshots with density (left) and velocity (right), one per row\n",
    "fig, axes = plt.subplots(6, 2, figsize=(12, 16), sharex=True)\n",
    "\n",
    "Nt = t_grid.size\n",
    "snapshot_indices = np.linspace(0, Nt - 1, 6, dtype=int)\n",
    "\n",
    "u1_lim = max(np.max(np.abs(u1_binned)), 1.0e-6)\n",
    "\n",
    "# Analytical shock speed (depends on test case)\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    shock_speed_analytical = (rho_R * u_R - rho_L * u_L) / (rho_R - rho_L)\n",
    "else:\n",
    "    shock_speed_analytical = (u_L + u_R) / 2\n",
    "\n",
    "for row, t_idx in enumerate(snapshot_indices):\n",
    "    rho_profile = rho_binned[t_idx, :]\n",
    "    u1_profile = u1_binned[t_idx, :]\n",
    "    \n",
    "    # Calculate analytical shock position at this time step\n",
    "    shock_pos_t = (0.5 + shock_speed_analytical * t_grid[t_idx]) % 1.0\n",
    "    shock2_pos_t = (shock_pos_t - 0.5) % 1.0\n",
    "\n",
    "    # Density snapshots (left column)\n",
    "    ax_rho = axes[row, 0]\n",
    "    ax_rho.plot(eta1_bins, rho_profile, color=\"tab:blue\", linewidth=2, label=r\"$\\rho$\")\n",
    "    ax_rho.axvline(shock_pos_t, color=\"orange\", linestyle=\":\", linewidth=2.5, alpha=0.8, label=\"Analytical shock\")\n",
    "    ax_rho.axvline(shock2_pos_t, color=\"orange\", linestyle=\":\", linewidth=2.5, alpha=0.8)\n",
    "    ax_rho.set_title(f\"t = {t_grid[t_idx]:.4f}\")\n",
    "    ax_rho.grid(True, alpha=0.3)\n",
    "    ax_rho.set_ylim([0, 1.2])\n",
    "    ax_rho.set_ylabel(r\"$\\rho$\")\n",
    "    if row == len(snapshot_indices) - 1:\n",
    "        ax_rho.set_xlabel(r\"$\\eta_1$ (logical coordinate)\")\n",
    "    ax_rho.legend(loc=\"upper right\")\n",
    "\n",
    "    # Velocity snapshots (right column)\n",
    "    ax_u = axes[row, 1]\n",
    "    ax_u.plot(eta1_bins, u1_profile, color=\"tab:red\", linewidth=2, label=r\"$u_1 = j_1/\\rho$\")\n",
    "    ax_u.axvline(shock_pos_t, color=\"orange\", linestyle=\":\", linewidth=2.5, alpha=0.8, label=\"Analytical shock\")\n",
    "    ax_u.axvline(shock2_pos_t, color=\"orange\", linestyle=\":\", linewidth=2.5, alpha=0.8)\n",
    "    ax_u.grid(True, alpha=0.3)\n",
    "    ax_u.set_ylim([-1.1 * u1_lim, 1.1 * u1_lim])\n",
    "    if row == len(snapshot_indices) - 1:\n",
    "        ax_u.set_xlabel(r\"$\\eta_1$ (logical coordinate)\")\n",
    "    ax_u.set_ylabel(r\"$u_1$\")\n",
    "    ax_u.legend(loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.suptitle(\"Riemann Shock Evolution: Density and Velocity Snapshots\", fontsize=14, y=0.995)\n",
    "plt.show()\n",
    "\n",
    "print(f\"\\nCaptured snapshots from t=0 to t={t_grid[-1]:.4f}\")\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    print(\"Contact discontinuity: expected zero shock speed (s = 0)\")\n",
    "else:\n",
    "    print(f\"Delta-shock: naive shock speed estimate s ≈ {(u_L + u_R)/2:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25",
   "metadata": {},
   "source": [
    "# Analytical Solutions and Comparison\n",
    "\n",
    "Now we compute analytical solutions and compare them to the numerical results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "## Analytical Solutions: Rankine-Hugoniot Conditions\n",
    "\n",
    "For pressureless Euler, the Rankine-Hugoniot jump conditions across a shock moving with speed $s$ are:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  s(\\rho_R - \\rho_L) &= \\rho_R u_R - \\rho_L u_L \\quad \\text{(mass conservation)} \\\\\n",
    "  s(\\rho_R u_R - \\rho_L u_L) &= \\rho_R u_R^2 - \\rho_L u_L^2 \\quad \\text{(momentum conservation)}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "### Test Case 1: Contact Discontinuity\n",
    "\n",
    "When $u_L = u_R$, the discontinuity is a *contact discontinuity* with zero shock speed:\n",
    "\n",
    "$$\n",
    "s = \\frac{\\rho_R u_R - \\rho_L u_L}{\\rho_R - \\rho_L}\n",
    "$$\n",
    "\n",
    "The weak solution is simply:\n",
    "\n",
    "$$\n",
    "(\\rho, u)(x, t) = \\begin{cases}\n",
    "  (\\rho_L, u_L) & \\text{if } x < x_0 \\\\\n",
    "  (\\rho_R, u_R) & \\text{if } x > x_0\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "with no propagation ($s = 0$ when $u_L = u_R$).\n",
    "\n",
    "### Test Case 2: Delta-Shock (Compression)\n",
    "\n",
    "When $u_L > u_R$ and $\\rho_L < \\rho_R$, the true Riemann solution is a **delta-shock**: a moving discontinuity where mass concentrates at the shock front. The weak solution exhibits density pile-up (higher density near the shock than either state). The simple estimate $s \\approx (u_L + u_R)/2$ is approximate; the SPH method will capture the physical density compression at the discontinuity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute shock position at final time\n",
    "t_final = t_grid[-1]\n",
    "x0_shock = 0.5\n",
    "\n",
    "# Shock speed depends on test case\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    # Contact discontinuity: shock speed from Rankine-Hugoniot\n",
    "    shock_speed = (rho_R * u_R - rho_L * u_L) / (rho_R - rho_L)\n",
    "else:  # delta_shock\n",
    "    # Delta-shock: velocity average (naive estimate, true solution is more complex)\n",
    "    shock_speed = (u_L + u_R) / 2\n",
    "\n",
    "shock_pos_final = x0_shock + shock_speed * t_final\n",
    "# Account for periodicity\n",
    "shock_pos_final = shock_pos_final % 1.0\n",
    "\n",
    "print(\"\\nAnalytical Shock Solution:\")\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    print(f\"  Shock speed (RH formula): s = (ρ_R u_R - ρ_L u_L) / (ρ_R - ρ_L) = {shock_speed:.4f}\")\n",
    "    print(\"  Note: Contact discontinuity has u_L = u_R, so velocity is continuous.\")\n",
    "else:\n",
    "    print(f\"  Shock speed (naive estimate): s ≈ (u_L + u_R)/2 = {shock_speed:.4f}\")\n",
    "    print(\"  Note: True pressureless Euler solution is a delta-shock with density pile-up.\")\n",
    "print(f\"  Initial shock position: x₀ = {x0_shock}\")\n",
    "print(f\"  Shock position at t={t_final:.4f}: x_shock ≈ {shock_pos_final:.4f}\")\n",
    "print(f\"  Distance traveled: {shock_speed * t_final:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "## Compare Numerical Solution with Analytical Prediction\n",
    "\n",
    "**For Contact Discontinuity (Test Case 1)**: The analytical solution is exact—velocity is continuous across the jump, and density follows a perfect step.\n",
    "\n",
    "**For Delta-Shock (Test Case 2)**: The analytical solution shown is a piecewise-constant reference only. The true pressureless Euler solution is a *delta-shock* where mass concentrates at the discontinuity (density pile-up), not captured by the simple step function. The SPH method will show this effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get final density profile from simulation\n",
    "rho_numerical = rho_binned[-1, :]\n",
    "\n",
    "# Construct analytical shock solution at final time\n",
    "eta1_fine = np.linspace(0, 1, 1000)\n",
    "shock_pos = x0_shock + shock_speed * t_final\n",
    "rho_analytical = np.where(eta1_fine < shock_pos, rho_L, rho_R)\n",
    "\n",
    "# Determine test case for title and labels\n",
    "if test_case == \"contact_discontinuity\":\n",
    "    title_prefix = \"Contact Discontinuity\"\n",
    "    info_text = \"Expected: stationary density jump at x = 0.5\"\n",
    "else:\n",
    "    title_prefix = \"Delta-Shock\"\n",
    "    info_text = \"Note: SPH may show density pile-up at shock front\"\n",
    "\n",
    "# Plot comparison\n",
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "# Numerical solution (binned)\n",
    "ax.step(eta1_bins, rho_numerical, where='mid', label='Numerical (SPH)', linewidth=2, color='blue')\n",
    "\n",
    "# Analytical weak solution\n",
    "ax.plot(eta1_fine, rho_analytical, '--', label='Analytical (piecewise constant)', linewidth=2, color='red')\n",
    "\n",
    "# Annotations\n",
    "ax.axvline(shock_pos_final, color='orange', linestyle=':', linewidth=2, alpha=0.7, label=f'Shock front at t={t_final:.4f}')\n",
    "ax.set_xlabel(r'$\\eta_1$ (logical coordinate)', fontsize=12)\n",
    "ax.set_ylabel(r'$\\rho$ (density)', fontsize=12)\n",
    "ax.set_title(f'{title_prefix}: Numerical vs Analytical at t={t_final:.4f}\\n({info_text})', fontsize=13)\n",
    "ax.grid(True, alpha=0.3)\n",
    "ax.legend(fontsize=11)\n",
    "ax.set_ylim([0, max(rho_analytical)*1.3] if max(rho_analytical) > 0 else [0, 1.2])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"\\nComparison at t={t_final:.4f}:\")\n",
    "print(f\"  Numerical shock position (steepest gradient): ≈ {eta1_bins[np.argmin(np.diff(rho_numerical))]}\")\n",
    "print(f\"  Analytical shock position: {shock_pos_final:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "## Compare Velocity at Final Time (Numerical vs Analytical)\n",
    "\n",
    "Using the same weak-solution shock states, we compare the reconstructed numerical velocity $u_1 = j_1/\\rho$ to the piecewise analytical profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical velocity at final time from binned current/density\n",
    "u1_numerical = u1_binned[-1, :]\n",
    "\n",
    "# Analytical velocity weak solution at final time\n",
    "u1_analytical = np.where(eta1_fine < shock_pos, u_L, u_R)\n",
    "\n",
    "# Plot comparison\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "\n",
    "ax.step(eta1_bins, u1_numerical, where='mid', label='Numerical $u_1$ (SPH)', linewidth=2, color='tab:red')\n",
    "ax.plot(eta1_fine, u1_analytical, '--', label='Analytical $u_1$ (Rankine-Hugoniot)', linewidth=2, color='black')\n",
    "\n",
    "ax.axvline(shock_pos_final, color='orange', linestyle=':', linewidth=2, alpha=0.7, label=f'Shock front at t={t_final:.4f}')\n",
    "ax.set_xlabel(r'$\\eta_1$ (logical coordinate)', fontsize=12)\n",
    "ax.set_ylabel(r'$u_1$ (velocity)', fontsize=12)\n",
    "ax.set_title(f'Large-Amplitude Shock Velocity: Numerical vs Analytical at t={t_final:.4f}', fontsize=13)\n",
    "ax.grid(True, alpha=0.3)\n",
    "ax.legend(fontsize=11)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"\\nVelocity comparison at t={t_final:.4f}:\")\n",
    "print(f\"  Numerical left/right plateau (approx): {np.mean(u1_numerical[:5]):.4f} / {np.mean(u1_numerical[-5:]):.4f}\")\n",
    "print(f\"  Analytical left/right states: {u_L:.4f} / {u_R:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "# Discussion: Shock Dynamics and Density Evolution\n",
    "\n",
    "In this Riemann-shock benchmark, density changes are intentionally large (from 1.0 to 0.125). This follows directly from the coupled pressureless Euler system:\n",
    "\n",
    "- The continuity equation $\\partial_t \\rho + \\partial_x(\\rho u) = 0$ couples density evolution to velocity gradients.\n",
    "- Across a shock, both $\\rho$ and $u$ jump and satisfy the Rankine-Hugoniot condition.\n",
    "- The SPH method captures the weak solution through particle redistribution and kernel-based remapping.\n",
    "\n",
    "The particle method naturally captures shocks through:\n",
    "\n",
    "1. **Particle compression**: particles accumulate near the discontinuity, increasing local density.\n",
    "2. **SPH kernel smoothing**: kernel averaging reconstructs a resolved shock layer from particle data.\n",
    "3. **Consistent momentum update**: velocity and density evolve together under the same conservation-law dynamics.\n",
    "\n",
    "## Extension\n",
    "\n",
    "To further explore:\n",
    "\n",
    "- Increase `Np_particles` for higher resolution near the shock front.\n",
    "- Tune kernel parameters in `WeightsParameters()` for different smoothing behavior.\n",
    "- Compare multiple time-stepping schemes via `ButcherTableau` options.\n",
    "- Use `PressureLessSPH/params_riemann_shock.py` as a standalone script counterpart to this notebook benchmark."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env (3.12.3)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}