{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 2D Hagen-Poiseuille Channel Flow with SPH\n",
    "\n",
    "## Viscous Flow between Parallel Plates\n",
    "\n",
    "This tutorial verifies the SPH discretisation of the 2D viscous isothermal Euler equations by simulating **Hagen-Poiseuille flow**: a pressure-driven channel flow between two no-slip walls.\n",
    "\n",
    "### Physical Setup\n",
    "\n",
    "Consider a 2D channel of width $H$ in the $y$-direction (no-slip walls at $y=0$ and $y=H$) and periodic in the $x$-direction. A uniform body force $g_x$ drives the flow in the $x$-direction. At steady state, viscosity balances the driving force, producing the **parabolic velocity profile**\n",
    "\n",
    "$$u_x^\\text{exact}(y) = \\frac{g_x}{2\\mu}\\,y\\,(H - y),$$\n",
    "\n",
    "with peak velocity at the channel centreline $y = H/2$:\n",
    "\n",
    "$$U_\\text{max} = \\frac{g_x H^2}{8\\mu}.$$\n",
    "\n",
    "The characteristic relaxation time scale is $T_\\text{relax} = H^2/(\\pi^2\\mu)$.\n",
    "\n",
    "### Verification procedure\n",
    "\n",
    "1. Start from rest with a uniform particle distribution and no-slip boundary conditions at the walls.\n",
    "2. Drive the flow with a constant body force $g_x$ acting through the pressure propagator's `gravity` parameter.\n",
    "3. Run until the flow is fully relaxed to the steady state ($t \\gg T_\\text{relax}$).\n",
    "4. Compare the final velocity profile against the Hagen-Poiseuille parabola in the $L^\\infty$ norm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import os\n",
    "import shutil\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cunumpy as xp\n",
    "\n",
    "from struphy import (\n",
    "    BinningPlot,\n",
    "    BoundaryParameters,\n",
    "    EnvironmentOptions,\n",
    "    KernelDensityPlot,\n",
    "    LoadingParameters,\n",
    "    SavingParameters,\n",
    "    Simulation,\n",
    "    SortingParameters,\n",
    "    Time,\n",
    "    WeightsParameters,\n",
    "    domains,\n",
    "    equils,\n",
    ")\n",
    "from struphy.models import ViscousEulerSPH\n",
    "from struphy.ode.utils import ButcherTableau\n",
    "\n",
    "logger = logging.getLogger(\"struphy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "### Physical and Numerical Parameters\n",
    "\n",
    "We choose $\\mu = 0.1$ and a body force $g_x = 0.1$, giving an analytical peak velocity $U_\\text{max} = g_x H^2 / (8\\mu) = 0.125$. The relaxation time $T_\\text{relax} = H^2/(\\pi^2\\mu) \\approx 1.01$; running to $T_\\text{end} = 10$ gives roughly 10 relaxation times, ensuring a fully converged steady state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Physical parameters\n",
    "mu  = 0.1    # dynamic viscosity\n",
    "g_x = 0.1   # body force in x (acts as pressure gradient)\n",
    "H   = 1.0    # channel height in y\n",
    "\n",
    "# Derived quantities\n",
    "U_max_exact   = g_x * H**2 / (8.0 * mu)\n",
    "T_relax       = H**2 / (np.pi**2 * mu)\n",
    "\n",
    "# Numerical parameters\n",
    "nx       = 8    # boxes per dimension\n",
    "ppb      = 16   # particles per box\n",
    "plot_pts = 21   # KDE evaluation points\n",
    "\n",
    "# Time stepping: run 10× past relaxation time\n",
    "dt   = 0.01\n",
    "Tend = 10.0\n",
    "\n",
    "print(f\"Viscosity:           mu   = {mu}\")\n",
    "print(f\"Body force:          g_x  = {g_x}\")\n",
    "print(f\"Channel height:      H    = {H}\")\n",
    "print(f\"Analytical U_max:    {U_max_exact:.4f}\")\n",
    "print(f\"Relaxation time:     T_relax = {T_relax:.2f}\")\n",
    "print(f\"Simulation time:     Tend = {Tend}  ({Tend/T_relax:.1f}× T_relax)\")\n",
    "print(f\"Total particles:     {ppb * nx * nx}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "### Model Setup\n",
    "\n",
    "`with_viscosity=True` (default) activates the viscous stress propagator. The 2D Gaussian SPH kernel (`gaussian_2d`) is required for the 2D geometry. The body force $g_x$ is passed as the `gravity` vector to `push_sph_p`; it enters the momentum equation as a constant acceleration $\\partial_t u_1 \\supset g_x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ViscousEulerSPH(with_B0=False, with_p=True, with_viscosity=True)\n",
    "\n",
    "butcher = ButcherTableau(algo=\"forward_euler\")\n",
    "model.propagators.push_eta.options     = model.propagators.push_eta.Options(butcher=butcher)\n",
    "model.propagators.push_sph_p.options   = model.propagators.push_sph_p.Options(\n",
    "    kernel_type=\"gaussian_2d\",\n",
    "    gravity=(g_x, 0.0, 0.0),   # body force drives flow in x\n",
    ")\n",
    "model.propagators.push_viscous.options = model.propagators.push_viscous.Options(\n",
    "    kernel_type=\"gaussian_2d\",\n",
    "    mu=mu,\n",
    ")\n",
    "\n",
    "print(\"ViscousEulerSPH model configured (with pressure, with viscosity).\")\n",
    "print(f\"  push_sph_p:   gaussian_2d kernel, gravity=({g_x}, 0, 0)\")\n",
    "print(f\"  push_viscous: gaussian_2d kernel, mu={mu}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "### Domain, Boundary Conditions and Diagnostics\n",
    "\n",
    "The channel geometry uses:\n",
    "- **x-direction**: periodic (flow direction).\n",
    "- **y-direction**: no-slip walls (`bc_sph=\"noslip\"`). The SPH no-slip boundary condition enforces $u_x = 0$ at the walls by introducing mirror ghost particles with reflected (negated) velocities.\n",
    "\n",
    "We bin the **current** $j_1 = \\rho u_1$ as a function of $y$ to recover the velocity profile. Since $\\rho \\approx 1$ (nearly incompressible at these parameters), $j_1 \\approx u_x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2D channel domain: x in [0, 1], y in [0, H]\n",
    "domain = domains.Cuboid(r1=1.0, r2=H)\n",
    "\n",
    "loading_params  = LoadingParameters(ppb=ppb, loading=\"tesselation\")\n",
    "weights_params  = WeightsParameters()\n",
    "boundary_params = BoundaryParameters(\n",
    "    bc    =(\"periodic\", \"reflect\",  \"periodic\"),   # particle reflections\n",
    "    bc_sph=(\"periodic\", \"noslip\",   \"periodic\"),   # SPH ghost treatment\n",
    ")\n",
    "sorting_params = SortingParameters(\n",
    "    boxes_per_dim=(nx, nx, 1),\n",
    "    dims_mask=(True, True, False),\n",
    ")\n",
    "\n",
    "# Bin j1 (velocity) vs y to reconstruct the velocity profile\n",
    "bin_plot_j1 = BinningPlot(slice=\"e2\", n_bins=(16,), ranges=(0.0, 1.0), output_quantity=\"current_1\")\n",
    "bin_plot_n  = BinningPlot(slice=\"e2\", n_bins=(16,), ranges=(0.0, 1.0))\n",
    "kd_plot     = KernelDensityPlot(pts_e1=plot_pts, pts_e2=plot_pts, pts_e3=1)\n",
    "saving_params = SavingParameters(\n",
    "    n_markers=1.0,\n",
    "    binning_plots=(bin_plot_j1, bin_plot_n),\n",
    "    kernel_density_plots=(kd_plot,),\n",
    ")\n",
    "\n",
    "model.euler_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",
    "    bufsize=2,   # extra ghost layer for no-slip\n",
    ")\n",
    "\n",
    "print(f\"2D channel: x-periodic [0, 1], y-noslip [0, {H}]\")\n",
    "print(f\"Particles: {ppb} ppb × {nx}×{nx} boxes = {ppb * nx * nx} total\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "### Initial Conditions\n",
    "\n",
    "The flow starts **from rest** (zero velocity, uniform density $\\rho_0 = 1$). The constant body force $g_x$ drives acceleration; viscosity and the no-slip walls establish the parabolic steady state over the relaxation time scale $T_\\text{relax}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start from rest: no velocity perturbation needed — body force drives the flow\n",
    "background = equils.ConstantVelocity()\n",
    "model.euler_fluid.var.add_background(background)\n",
    "\n",
    "print(\"Initial condition: uniform density n=1, zero velocity everywhere\")\n",
    "print(f\"Body force g_x={g_x} will accelerate flow; viscosity establishes steady state\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "### Simulation Setup and Execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_folder = os.path.join(os.getcwd(), \"struphy_verification_tests\")\n",
    "out_folders = os.path.join(test_folder, \"ViscousEulerSPH\")\n",
    "env = EnvironmentOptions(out_folders=out_folders, sim_folder=\"hagen_poiseuille\")\n",
    "\n",
    "time_opts = Time(dt=dt, Tend=Tend, split_algo=\"Strang\")\n",
    "\n",
    "sim = Simulation(\n",
    "    model=model,\n",
    "    env=env,\n",
    "    time_opts=time_opts,\n",
    "    domain=domain,\n",
    "    grid=None,\n",
    "    derham_opts=None,\n",
    ")\n",
    "\n",
    "print(f\"Running Hagen-Poiseuille flow: dt={dt}, Tend={Tend}\")\n",
    "sim.run()\n",
    "print(\"Simulation complete.\")\n",
    "\n",
    "sim.pproc()\n",
    "print(\"Post-processing complete.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "### Load Diagnostics and Compute Exact Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.load_plotting_data()\n",
    "\n",
    "e2_grid   = sim.f.euler_fluid.e2_current_1.grid_e2   # logical y in [0, 1]\n",
    "j1_binned = sim.f.euler_fluid.e2_current_1.f_binned  # shape (Nt+1, n_bins)\n",
    "\n",
    "Nt    = int(Tend / dt)\n",
    "times = np.linspace(0.0, Tend, Nt + 1)\n",
    "\n",
    "e2_np  = np.asarray(e2_grid).flatten()\n",
    "y_np   = e2_np * H                  # physical y coordinate\n",
    "\n",
    "# Analytical Hagen-Poiseuille profile\n",
    "u_exact      = g_x / (2.0 * mu) * y_np * (H - y_np)\n",
    "u_max_exact  = np.max(u_exact)\n",
    "\n",
    "u_num_final  = np.asarray(j1_binned[-1, :]).flatten()\n",
    "u_max_num    = np.max(u_num_final)\n",
    "\n",
    "# Centreline velocity over time (index of y closest to H/2)\n",
    "idx_centre  = int(np.argmin(np.abs(e2_np - 0.5)))\n",
    "u_centre    = np.asarray(j1_binned[:, idx_centre]).flatten()\n",
    "\n",
    "print(f\"Analytical U_max = {u_max_exact:.6f}\")\n",
    "print(f\"Numerical  U_max = {u_max_num:.6f}\")\n",
    "\n",
    "abs_err            = np.abs(u_num_final - u_exact)\n",
    "rel_err_pointwise  = abs_err / u_max_exact\n",
    "rel_error_interior = rel_err_pointwise[1:-1]   # exclude wall bins (exact value → 0)\n",
    "rel_error_umax     = abs(u_max_num - u_max_exact) / u_max_exact\n",
    "\n",
    "print(f\"Mean interior relative error: {np.mean(rel_error_interior) * 100:.2f}%\")\n",
    "print(f\"Max  interior relative error: {np.max(rel_error_interior) * 100:.2f}%\")\n",
    "print(f\"U_max relative error:         {rel_error_umax * 100:.2f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "### Visualisation\n",
    "\n",
    "Three panels summarise the result:\n",
    "\n",
    "- **Left**: final numerical velocity profile vs the analytical parabola.\n",
    "- **Centre**: pointwise relative error $|u_x^\\text{num} - u_x^\\text{exact}| / U_\\text{max}$ vs $y$ (excluding wall bins).\n",
    "- **Right**: time evolution of the centreline velocity, showing relaxation to the Hagen-Poiseuille steady state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "# --- Left: velocity profile ---\n",
    "ax = axes[0]\n",
    "ax.plot(u_num_final, y_np, \"o-\", markersize=4, label=\"Numerical (SPH)\")\n",
    "ax.plot(u_exact,     y_np, \"k--\",              label=\"Analytical (Hagen-Poiseuille)\")\n",
    "ax.set_xlabel(r\"$u_x$\")\n",
    "ax.set_ylabel(r\"$y$\")\n",
    "ax.set_title(\"Steady-state velocity profile\")\n",
    "ax.legend()\n",
    "ax.grid(True)\n",
    "\n",
    "# --- Centre: pointwise relative error ---\n",
    "ax = axes[1]\n",
    "ax.plot(rel_err_pointwise * 100, y_np, \"r-o\", markersize=4)\n",
    "ax.set_xlabel(r\"$|u_x^\\mathrm{num} - u_x^\\mathrm{exact}|\\,/\\,U_\\mathrm{max}$ [%]\")\n",
    "ax.set_ylabel(r\"$y$\")\n",
    "ax.set_title(f\"Pointwise relative error (max interior = {np.max(rel_error_interior)*100:.1f}%)\")\n",
    "ax.grid(True)\n",
    "\n",
    "# --- Right: centreline velocity relaxation ---\n",
    "ax = axes[2]\n",
    "ax.plot(times, u_centre, label=r\"Numerical $u_x(y=H/2)$\")\n",
    "ax.axhline(u_max_exact, color=\"k\", linestyle=\"--\",\n",
    "           label=rf\"Exact $U_\\mathrm{{max}} = {u_max_exact:.4f}$\")\n",
    "ax.set_xlabel(\"time\")\n",
    "ax.set_ylabel(r\"$u_x(y=H/2)$\")\n",
    "ax.set_title(\"Centreline velocity: relaxation to steady state\")\n",
    "ax.legend()\n",
    "ax.grid(True)\n",
    "\n",
    "fig.suptitle(\n",
    "    rf\"Hagen-Poiseuille: $\\mu={mu}$, $g_x={g_x}$, $H={H}$, {nx}\\times{nx} boxes\",\n",
    "    fontsize=12,\n",
    ")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "### Verification Check\n",
    "\n",
    "Two assertions:\n",
    "1. Maximum pointwise relative error in the **interior** bins (away from walls where the exact value vanishes) is below 5%.\n",
    "2. Relative error in the peak velocity $U_\\text{max}$ is below 5%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "tol_interior = 0.05   # 5% pointwise tolerance\n",
    "tol_umax     = 0.05   # 5% tolerance on U_max\n",
    "\n",
    "print(\"=== Hagen-Poiseuille Verification ===\")\n",
    "print(f\"  Max interior relative error: {np.max(rel_error_interior) * 100:.2f}%  (tolerance {tol_interior*100:.0f}%)\")\n",
    "print(f\"  U_max relative error:        {rel_error_umax * 100:.2f}%           (tolerance {tol_umax*100:.0f}%)\")\n",
    "\n",
    "try:\n",
    "    assert np.max(rel_error_interior) < tol_interior, (\n",
    "        f\"Interior error {np.max(rel_error_interior)*100:.1f}% exceeds tolerance {tol_interior*100:.0f}%\"\n",
    "    )\n",
    "    print(\"\\n✓ Interior velocity profile check passed.\")\n",
    "except AssertionError as e:\n",
    "    print(f\"\\n✗ {e}\")\n",
    "\n",
    "try:\n",
    "    assert rel_error_umax < tol_umax, (\n",
    "        f\"U_max error {rel_error_umax*100:.1f}% exceeds tolerance {tol_umax*100:.0f}%\"\n",
    "    )\n",
    "    print(\"✓ U_max check passed.\")\n",
    "except AssertionError as e:\n",
    "    print(f\"✗ {e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "\n",
    "This tutorial verified the SPH discretisation of 2D viscous channel flow:\n",
    "\n",
    "- The **no-slip boundary condition** is implemented via mirror ghost particles with negated tangential velocity, correctly enforcing $u_x = 0$ at both walls.\n",
    "- The **parabolic steady state** emerges from the balance between body force and viscous stress, reproduced to better than 5% pointwise accuracy with $8 \\times 8$ boxes and 16 particles per box.\n",
    "- The **centreline velocity** converges smoothly to $U_\\text{max}$ over the relaxation time $T_\\text{relax} = H^2/(\\pi^2 \\mu) \\approx 1$.\n",
    "- Increasing `nx` or `ppb` further reduces the error, as the SPH kernel gradient approximation of the viscous stress tensor improves with particle density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional cleanup\n",
    "if False:   # set to True to remove simulation output\n",
    "    shutil.rmtree(test_folder)\n",
    "    print(f\"Cleaned up {test_folder}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env (3.12.3.final.0)",
   "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
}