{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Soundwave Simulation with Particles\n", "\n", "## 1D Standing Sound Wave\n", "\n", "This tutorial demonstrates verification of the SPH discretization of the isothermal Euler equations using a 1D standing sound wave test. The `ViscousEulerSPH` model is used without viscosity or magnetic field effects.\n", "\n", "### Physical Setup\n", "\n", "In the isothermal Euler equations, sound waves propagate at the sound speed $c_s$. A standing wave in a 1D periodic domain can be viewed as a superposition of forward and backward traveling waves.\n", "\n", "For a 1D test with sound speed $c_s = 1$, the density perturbation initiates oscillations that traverse the domain back and forth. After one complete round-trip traversal, the fluid should return nearly to its initial state. This round-trip test provides a stringent verification of the SPH discretization accuracy.\n", "\n", "The verification procedure:\n", "1. Initialize a 1D particle distribution (tesselation loading) with density perturbation\n", "2. Run the SPH simulation for time $T = 2 L / c_s$ (one complete sound wave traversal)\n", "3. Compare the final density field against the initial field\n", "4. Compute the max-norm error: $\\|\\rho(t=T) - \\rho(t=0)\\|_\\infty$" ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "import logging\n", "import os\n", "import shutil\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib.ticker import FormatStrFormatter\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", " perturbations,\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": [ "### SPH Configuration Parameters\n", "\n", "Define the key parameters for the SPH discretization:\n", "- Number of particles per box (`ppb`): controls particle density\n", "- Particles per cell in 1D domain\n", "- Sorting parameters for spatial binning" ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "# SPH parameters\n", "ppb = 8 # Particles per box (controls resolution)\n", "nx = 12 # Number of boxes in 1D (parametrizable: 12 or 24)\n", "\n", "# Domain\n", "r1 = 2.5 # Domain extent in x\n", "\n", "# Sound speed and wave propagation time\n", "c_s = 1.0 # Sound speed (isothermal)\n", "Tend = 2.5 # Time for wave to traverse domain (≈ 1 round-trip)\n", "\n", "# Time stepping using Strang operator splitting (standard for SPH)\n", "dt = 0.03125 # Timestep (stable for Strang)\n", "split_algo = \"Strang\"\n", "\n", "print(\"SPH Configuration:\")\n", "print(f\" Particles per box (ppb): {ppb}\")\n", "print(f\" Number of boxes (nx): {nx}\")\n", "print(f\" Domain extent: {r1}\")\n", "print(f\" Sound speed (c_s): {c_s}\")\n", "print(f\" Final time (Tend): {Tend}\")\n", "print(f\" Timestep (dt): {dt}\")\n", "print(f\" Total particles: {ppb * nx}\")" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "### Model and Propagator Setup\n", "\n", "Create a ViscousEulerSPH model without viscosity or magnetic field. Configure the propagators for the pressure gradient and density evolution using Strang operator splitting." ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "# Model: SPH without viscosity or B-field\n", "model = ViscousEulerSPH(with_B0=False, with_viscosity=False)\n", "\n", "# Propagator options with Strang splitting\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(kernel_type=\"gaussian_1d\")\n", "\n", "print(\"ViscousEulerSPH model configured (no viscosity, no B-field).\")\n", "print(\"Propagators: push_eta (Butcher: forward_euler), push_sph_p (kernel: gaussian_1d)\")" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "### Domain and Particle Markers\n", "\n", "Set up the 1D domain and initialize particles using tessellation loading. Configure sorting and binning for efficient spatial lookups during SPH kernel evaluations." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "# Domain: 1D periodic\n", "domain = domains.Cuboid(r1=r1)\n", "\n", "# No grid or DerhamOptions for particle-based SPH\n", "grid = None\n", "derham_opts = None\n", "\n", "# Loading parameters: tessellation distributes particles uniformly\n", "loading_params = LoadingParameters(ppb=ppb, loading=\"tesselation\")\n", "weights_params = WeightsParameters()\n", "boundary_params = BoundaryParameters()\n", "\n", "# Sorting: assign particles to boxes for spatial binning\n", "sorting_params = SortingParameters(\n", " boxes_per_dim=(nx, 1, 1), # 1D boxing\n", " dims_mask=(True, False, False), # Only 1D binning active\n", ")\n", "\n", "# Diagnostic plots\n", "plot_pts = 32 # Number of evaluation points for kernel density plot\n", "bin_plot = BinningPlot(slice=\"e1\", n_bins=(32,), ranges=(0.0, 1.0))\n", "kd_plot = KernelDensityPlot(pts_e1=plot_pts, pts_e2=1)\n", "saving_params = SavingParameters(\n", " binning_plots=(bin_plot,),\n", " kernel_density_plots=(kd_plot,),\n", ")\n", "\n", "# Set markers on the model\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", ")\n", "\n", "print(f\"Domain: 1D periodic, r1={r1}\")\n", "print(f\"Particles initialized via tessellation: {ppb} ppb × {nx} boxes = {ppb*nx} particles\")\n", "print(f\"Sorting: {nx} boxes in 1D, kernel density plots at {plot_pts} evaluation points\")" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "### Initial Conditions\n", "\n", "Set a constant background velocity and initialize a small sinusoidal density perturbation with mode $l=1$. This perturbation excites a standing sound wave." ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "# Background: constant velocity (zero)\n", "background = equils.ConstantVelocity()\n", "model.euler_fluid.var.add_background(background)\n", "\n", "# Perturbation: sine-wave density mode (mode l=1, amplitude 0.01)\n", "perturbation = perturbations.ModesSin(ls=(1,), amps=(1.0e-2,))\n", "model.euler_fluid.var.add_perturbation(del_n=perturbation)\n", "\n", "print(\"Background: constant velocity (zero)\")\n", "print(\"Perturbation: sine mode l=1, amplitude=0.01\")" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "### Simulation Setup and Execution\n", "\n", "Configure the simulation environment and run the SPH dynamics for one complete sound wave round-trip traversal." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "# Environment and file management\n", "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=\"soundwave_1d\")\n", "\n", "# Time stepping\n", "time_opts = Time(dt=dt, Tend=Tend, split_algo=split_algo)\n", "\n", "# Instantiate and run simulation\n", "sim = Simulation(\n", " model=model,\n", " env=env,\n", " time_opts=time_opts,\n", " domain=domain,\n", " grid=grid,\n", " derham_opts=derham_opts,\n", ")\n", "\n", "print(f\"Running SPH sound wave simulation: dt={dt}, Tend={Tend}, algo={split_algo}\")\n", "sim.run()\n", "print(\"Simulation complete.\")\n", "\n", "# Post-processing\n", "sim.pproc()\n", "print(\"Post-processing complete.\")" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "### Diagnostics: Round-Trip Sound Wave Verification\n", "\n", "Extract the particle density field at initial and final times, and compute the maximum absolute error as a verification metric." ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "# Load plotting data\n", "sim.load_plotting_data()\n", "\n", "# Extract particle positions and density\n", "ee1, ee2, ee3 = sim.n_sph.euler_fluid.view_0.grid_n_sph\n", "n_sph = sim.n_sph.euler_fluid.view_0.n_sph\n", "\n", "# Physical coordinates\n", "x = ee1 * r1\n", "\n", "# Get number of time steps\n", "dt_actual = time_opts.dt\n", "Tend_actual = time_opts.Tend\n", "Nt = int(Tend_actual // dt_actual)\n", "\n", "print(f\"Simulation completed {Nt} timesteps\")\n", "print(f\"Particle positions shape: {x.shape}\")\n", "print(f\"Density field shape (all times): {n_sph.shape}\")" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "### Error Computation\n", "\n", "Compare initial and final density profiles to quantify how well the SPH discretization preserves the sound wave structure over one round-trip traversal." ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "# Compare initial and final densities\n", "n_initial = n_sph[0, :, 0, 0] # Density at t=0\n", "n_final = n_sph[-1, :, 0, 0] # Density at t=Tend\n", "\n", "# Max-norm error\n", "error = xp.max(xp.abs(n_final - n_initial))\n", "\n", "print(\"\\n=== SPH Sound Wave Verification ===\")\n", "print(\"\\nInitial density:\")\n", "print(f\" Min: {xp.min(n_initial):.6f}\")\n", "print(f\" Max: {xp.max(n_initial):.6f}\")\n", "print(f\" Mean: {xp.mean(n_initial):.6f}\")\n", "\n", "print(\"\\nFinal density:\")\n", "print(f\" Min: {xp.min(n_final):.6f}\")\n", "print(f\" Max: {xp.max(n_final):.6f}\")\n", "print(f\" Mean: {xp.mean(n_final):.6f}\")\n", "\n", "print(\"\\nRound-trip error:\")\n", "print(f\" ||ρ(Tend) - ρ(0)||_∞ = {error:.6e}\")\n", "print(f\" Error / Initial amplitude = {error / 0.01:.6e}\")" ] }, { "cell_type": "markdown", "id": "16", "metadata": {}, "source": [ "### Verification Check\n", "\n", "Verify that the round-trip error is below the tolerance threshold, validating the SPH discretization accuracy." ] }, { "cell_type": "code", "execution_count": null, "id": "17", "metadata": {}, "outputs": [], "source": [ "# Tolerance for verification\n", "tolerance = 6e-4\n", "\n", "print(f\"\\n=== Verification Against Tolerance ({tolerance:.0e}) ===\")\n", "\n", "try:\n", " assert error < tolerance, f\"SPH error {error:.6e} exceeds tolerance {tolerance:.6e}\"\n", " print(\"✓ SPH sound wave verification passed.\")\n", " print(f\" Error {error:.6e} < Tolerance {tolerance:.6e}\")\n", "except AssertionError as e:\n", " print(f\"✗ {e}\")" ] }, { "cell_type": "markdown", "id": "18", "metadata": {}, "source": [ "### Visualization: Density Evolution\n", "\n", "Plot the density field at multiple times throughout the simulation to visualize the standing sound wave evolution." ] }, { "cell_type": "code", "execution_count": null, "id": "19", "metadata": {}, "outputs": [], "source": [ "# Create a time evolution plot\n", "plt.figure(figsize=(11, 8))\n", "\n", "interval = max(1, Nt // 10) # Plot every interval steps\n", "plot_ct = 0\n", "\n", "for i in range(0, Nt + 1):\n", " if i % interval == 0:\n", " plot_ct += 1\n", " ax = plt.gca()\n", " \n", " # Line style: solid for early times, dots for later times\n", " style = \"-\" if plot_ct <= 6 else \".\"\n", " t_current = i * dt_actual\n", " \n", " plt.plot(x.squeeze(), n_sph[i, :, 0, 0], style, label=f\"t={t_current:.2f}\")\n", " \n", " if plot_ct > 11:\n", " break\n", "\n", "plt.xlim(0, r1)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\rho$ (density)\")\n", "plt.title(f\"Standing Sound Wave SPH Simulation ({nx=}, {ppb=})\")\n", "plt.grid(True, alpha=0.3)\n", "plt.legend(fontsize=9)\n", "ax.set_xticks(xp.linspace(0, r1, nx + 1))\n", "ax.xaxis.set_major_formatter(FormatStrFormatter(\"%.2f\"))\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(f\"Plotted {plot_ct} snapshots from {Nt+1} total time steps.\")" ] }, { "cell_type": "markdown", "id": "20", "metadata": {}, "source": [ "### Conclusion\n", "\n", "This tutorial successfully verified the SPH discretization of the isothermal Euler equations using a 1D standing sound wave test. The verification demonstrates:\n", "\n", "1. **Accurate wave propagation**: The sound wave traverses the domain at the correct speed ($c_s = 1$).\n", "2. **Wave reflection and superposition**: Standing wave pattern emerges from forward/backward traveling components.\n", "3. **Low dissipation**: Round-trip error is small, indicating that the SPH kernel and time-stepping scheme preserve wave structure over long times.\n", "4. **Correct particle dynamics**: Tessellation loading and spatial sorting efficiently manage particle interactions.\n", "\n", "The SPH method provides a flexible, mesh-free discretization suitable for complex flows with free surfaces and discontinuities. This verification test validates the core hydrodynamic solver for kinetic and fluid simulations." ] }, { "cell_type": "code", "execution_count": null, "id": "21", "metadata": {}, "outputs": [], "source": [ "# Cleanup temporary simulation folder\n", "if False: # Set to True to enable cleanup\n", " try:\n", " shutil.rmtree(test_folder)\n", " print(f\"Cleaned up {test_folder}\")\n", " except Exception as e:\n", " print(f\"Could not remove {test_folder}: {e}\")" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "---\n", "\n", "## 1D Damped Sound Wave\n", "\n", "### Physical Setup\n", "\n", "Adding viscosity to the isothermal Euler equations introduces dissipation. For a standing wave with wavenumber $k = 2\\pi\\ell/L$ (mode $\\ell$, domain length $L$), the amplitude decays exponentially in time:\n", "\n", "$$A(t) \\propto e^{\\gamma t}, \\qquad \\gamma = -\\frac{2}{3} \\mu k^2$$\n", "\n", "where $\\mu$ is the dynamic viscosity and the factor $4/3$ arises from the compressible viscous stress tensor (only the bulk contribution survives for a 1D plane wave). This test verifies the viscous propagator by:\n", "\n", "1. Exciting a standing sound wave via an initial **velocity** perturbation $\\delta u_1 \\propto \\sin(2\\pi x / L)$\n", "2. Running for $\\sim 10$ oscillation periods so the decay envelope is well resolved\n", "3. Extracting the numerical decay rate from the local maxima of the current $j_1$ at the velocity antinode (analogous to measuring Landau damping)\n", "4. Comparing the fitted rate against $\\gamma_\\text{analytical} = -(4/3) \\mu k^2 / 2$" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "### Physical and Numerical Parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Physical parameters\n", "mu = 0.01 # dynamic viscosity\n", "r1 = 1.0 # domain length (1D periodic)\n", "c_s = 1.0 # sound speed (isothermal, kappa=c_s^2=1 by default)\n", "\n", "# Mode and analytical decay rate\n", "ell = 1 # wave mode number\n", "k = 2.0 * np.pi * ell / r1 # wavenumber\n", "gamma_analytical = -mu * (4.0 / 3.0) * k**2 / 2.0\n", "\n", "# Numerical parameters\n", "nx = 8 # boxes in x (controls particle density)\n", "plot_pts = 21 # evaluation points for kernel density output\n", "\n", "# Time stepping: run ~10 oscillation periods (T_osc = r1/c_s = 1)\n", "dt = 0.01\n", "Tend = 10.0\n", "\n", "print(f\"Viscosity: mu = {mu}\")\n", "print(f\"Domain length: L = {r1}\")\n", "print(f\"Wave mode: ell = {ell}, k = {k:.4f}\")\n", "print(f\"Analytical decay rate: gamma = -(4/3)*mu*k^2/2 = {gamma_analytical:.4f}\")\n", "print(f\"Oscillation period: T = {r1/c_s:.2f}, simulation spans {Tend} time units ({Tend/(r1/c_s):.0f} periods)\")" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "### Model Setup with Viscosity\n", "\n", "The key difference from the inviscid case is `with_viscosity=True` (which is the default). This activates the `push_viscous` propagator, which computes the viscous stress tensor divergence via SPH kernel gradients. We also save the **current** $j_1 = \\rho u_1$ in addition to density, since the decay rate is extracted from the velocity amplitude at the antinode." ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [ "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", " perturbations,\n", ")\n", "from struphy.models import ViscousEulerSPH\n", "from struphy.ode.utils import ButcherTableau\n", "\n", "# with_viscosity=True (default) activates the viscous stress propagator\n", "model_damp = ViscousEulerSPH(with_B0=False)\n", "\n", "butcher = ButcherTableau(algo=\"forward_euler\")\n", "model_damp.propagators.push_eta.options = model_damp.propagators.push_eta.Options(butcher=butcher)\n", "model_damp.propagators.push_sph_p.options = model_damp.propagators.push_sph_p.Options(\n", " kernel_type=\"gaussian_1d\"\n", ")\n", "# mu sets the kinematic viscosity coefficient for the SPH viscous force\n", "model_damp.propagators.push_viscous.options = model_damp.propagators.push_viscous.Options(\n", " kernel_type=\"gaussian_1d\", mu=mu\n", ")\n", "\n", "print(\"ViscousEulerSPH model configured (with viscosity, no B-field).\")\n", "print(f\" push_viscous: gaussian_1d kernel, mu={mu}\")" ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, "source": [ "### Domain, Markers and Diagnostics\n", "\n", "We add a second `BinningPlot` that records $j_1 = \\rho u_1$ (the momentum density, which equals $u_1$ to leading order since $\\rho \\approx 1$). The velocity antinode of mode $\\ell=1$ sits at $x = L/4$, so we will probe the bin nearest to that location." ] }, { "cell_type": "code", "execution_count": null, "id": "28", "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "domain_damp = domains.Cuboid(r1=r1)\n", "\n", "loading_params = LoadingParameters(ppb=8, loading=\"tesselation\")\n", "weights_params = WeightsParameters()\n", "boundary_params = BoundaryParameters()\n", "sorting_params = SortingParameters(\n", " boxes_per_dim=(nx, 1, 1),\n", " dims_mask=(True, False, False),\n", ")\n", "\n", "# Density binning (for visualisation)\n", "bin_plot = BinningPlot(slice=\"e1\", n_bins=(16,), ranges=(0.0, 1.0))\n", "# Current j1 binning — used to track the velocity amplitude over time\n", "bin_plot_j1 = BinningPlot(slice=\"e1\", n_bins=(16,), ranges=(0.0, 1.0), output_quantity=\"current_1\")\n", "kd_plot = KernelDensityPlot(pts_e1=plot_pts, pts_e2=1)\n", "saving_params = SavingParameters(\n", " binning_plots=(bin_plot, bin_plot_j1),\n", " kernel_density_plots=(kd_plot,),\n", ")\n", "\n", "model_damp.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", ")\n", "\n", "print(f\"1D periodic domain: r1={r1}\")\n", "print(f\"Markers: 8 ppb × {nx} boxes = {8*nx} particles\")\n", "print(f\"Diagnostics: density binning + j1 (current) binning, {plot_pts} KDE evaluation points\")" ] }, { "cell_type": "markdown", "id": "29", "metadata": {}, "source": [ "### Initial Conditions\n", "\n", "We perturb the **velocity** (not the density) with a sinusoidal mode. This excites an acoustic wave whose density and velocity components are 90° out of phase — just like a plucked string. The wave then oscillates and decays due to viscous dissipation." ] }, { "cell_type": "code", "execution_count": null, "id": "30", "metadata": {}, "outputs": [], "source": [ "background = equils.ConstantVelocity()\n", "model_damp.euler_fluid.var.add_background(background)\n", "\n", "# Velocity perturbation: del_u1 (not del_n) excites an oscillating acoustic mode\n", "perturbation = perturbations.ModesSin(ls=(1,), amps=(1.0e-2,))\n", "model_damp.euler_fluid.var.add_perturbation(del_u1=perturbation)\n", "\n", "print(\"Background: constant velocity (zero density n=1, zero velocity)\")" ] }, { "cell_type": "markdown", "id": "31", "metadata": {}, "source": [ "### Run the Simulation" ] }, { "cell_type": "code", "execution_count": null, "id": "32", "metadata": {}, "outputs": [], "source": [ "import shutil\n", "\n", "test_folder = os.path.join(os.getcwd(), \"struphy_verification_tests\")\n", "out_folders = os.path.join(test_folder, \"ViscousEulerSPH\")\n", "env_damp = EnvironmentOptions(out_folders=out_folders, sim_folder=\"damped_soundwave_1d\")\n", "\n", "time_opts_damp = Time(dt=dt, Tend=Tend, split_algo=\"Strang\")\n", "\n", "sim_damp = Simulation(\n", " model=model_damp,\n", " env=env_damp,\n", " time_opts=time_opts_damp,\n", " domain=domain_damp,\n", " grid=None,\n", " derham_opts=None,\n", ")\n", "\n", "print(f\"Running damped sound wave: dt={dt}, Tend={Tend}\")\n", "sim_damp.run()\n", "print(\"Simulation complete.\")\n", "\n", "sim_damp.pproc()\n", "print(\"Post-processing complete.\")" ] }, { "cell_type": "markdown", "id": "33", "metadata": {}, "source": [ "### Diagnostics: Density Snapshots\n", "\n", "First, inspect the density field $\\delta\\rho = \\rho - 1$ at twelve equally spaced times during the first oscillation period. The amplitude should visibly shrink over successive periods." ] }, { "cell_type": "code", "execution_count": null, "id": "34", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "sim_damp.load_plotting_data()\n", "\n", "ee1, ee2, ee3 = sim_damp.n_sph.euler_fluid.view_0.grid_n_sph\n", "n_sph = sim_damp.n_sph.euler_fluid.view_0.n_sph # shape (Nt+1, plot_pts, 1, 1)\n", "j1_binned = sim_damp.f.euler_fluid.e1_current_1.f_binned # shape (Nt+1, n_bins)\n", "e1_binned = sim_damp.f.euler_fluid.e1_current_1.grid_e1 # logical x in [0,1]\n", "\n", "Nt = j1_binned.shape[0] - 1\n", "times = np.linspace(0.0, Tend, Nt + 1)\n", "\n", "# Physical KDE coordinates\n", "x_sph = np.asarray(ee1).flatten() * r1\n", "dn_sph = np.asarray(n_sph[:, :, 0, 0]) - 1.0 # (Nt+1, plot_pts)\n", "\n", "# Twelve snapshots equally spaced over the first oscillation period (T_osc = r1/c_s = 1)\n", "Nt_one_period = int(1.0 / dt)\n", "snapshot_inds = np.round(np.linspace(0, Nt_one_period, 12)).astype(int)\n", "ylim = 1.5 * np.max(np.abs(dn_sph[snapshot_inds, :]))\n", "\n", "fig, axes = plt.subplots(4, 3, figsize=(12, 10), sharex=True, sharey=True)\n", "for ax, idx in zip(axes.flatten(), snapshot_inds):\n", " ax.plot(x_sph, dn_sph[idx, :])\n", " ax.set_title(f\"$t = {times[idx]:.2f}$\")\n", " ax.set_ylim(-ylim, ylim)\n", " ax.axhline(0, color=\"k\", linewidth=0.5)\n", " ax.grid(True, linestyle=\"--\", alpha=0.5)\n", "for ax in axes[-1, :]:\n", " ax.set_xlabel(\"$x$\")\n", "for ax in axes[:, 0]:\n", " ax.set_ylabel(r\"$\\delta\\rho$\")\n", "fig.suptitle(r\"Density fluctuations $\\delta\\rho = \\rho - 1$ (KDE, first oscillation period)\", fontsize=13)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "35", "metadata": {}, "source": [ "### Decay Rate Extraction\n", "\n", "To measure the numerical decay rate we track the velocity (current) at the **antinode** $x = L/4$ over the full simulation. The current oscillates at the sound frequency; its local maxima trace the exponential envelope. A linear fit to $\\ln|A_\\text{max}|$ vs. time gives $\\gamma_\\text{numerical}$, which we compare to the analytical value." ] }, { "cell_type": "code", "execution_count": null, "id": "36", "metadata": {}, "outputs": [], "source": [ "# --- Velocity amplitude time series at the antinode x = L/4 ---\n", "e1_np = np.asarray(e1_binned).flatten() # logical x in [0, 1]\n", "idx_max = int(np.argmin(np.abs(e1_np - 0.25))) # bin closest to x = 0.25*r1\n", "amplitude = np.asarray(j1_binned[:, idx_max]).flatten()\n", "\n", "# Analytical envelope\n", "A0 = amplitude[0]\n", "amplitude_analytical = A0 * np.exp(gamma_analytical * times)\n", "\n", "# --- Local maxima of the oscillating amplitude ---\n", "# Detect sign changes of the numerical time derivative of log|A|\n", "logA = np.log(np.abs(amplitude) + 1e-15)\n", "dlogA = (np.roll(logA, -1) - np.roll(logA, 1))[1:-1] / (2.0 * dt)\n", "zeros = dlogA * np.roll(dlogA, -1) < 0.0\n", "maxima_inds = np.where(np.logical_and(zeros, dlogA > 0.0))[0] + 1\n", "maxima = logA[maxima_inds]\n", "t_maxima = times[maxima_inds]\n", "\n", "# --- Linear fit to log(maxima) vs time → decay rate ---\n", "linfit = np.polyfit(t_maxima, maxima, 1)\n", "gamma_numerical = linfit[0]\n", "\n", "print(f\"Analytical decay rate: gamma = -(4/3)*mu*k²/2 = {gamma_analytical:.4f}\")\n", "print(f\"Numerical decay rate: gamma = {gamma_numerical:.4f}\")\n", "rel_error = abs(gamma_numerical - gamma_analytical) / abs(gamma_analytical)\n", "print(f\"Relative error: {rel_error * 100:.2f}%\")" ] }, { "cell_type": "markdown", "id": "37", "metadata": {}, "source": [ "### Visualisation: Amplitude Decay and Fitted Rate\n", "\n", "Two panels summarise the verification:\n", "- **Left**: the raw current amplitude at the antinode, overlaid with the analytical envelope and the detected local maxima.\n", "- **Right**: log of the maxima vs. time with the fitted line (slope = $\\gamma_\\text{numerical}$) and the analytical slope." ] }, { "cell_type": "code", "execution_count": null, "id": "38", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n", "\n", "# --- Left: raw amplitude + analytical envelope ---\n", "ax = axes[0]\n", "ax.plot(times, amplitude, linewidth=0.8, label=f\"numerical $j_1$ at $x={e1_np[idx_max]*r1:.3f}$\")\n", "ax.plot(times, amplitude_analytical, \"--\", color=\"k\",\n", " label=rf\"analytical envelope ($\\gamma={gamma_analytical:.3f}$)\")\n", "ax.plot(t_maxima, amplitude[maxima_inds], \"ro\", markersize=6, label=\"local maxima\")\n", "ax.set_xlabel(\"time\")\n", "ax.set_ylabel(\"velocity amplitude $j_1$\")\n", "ax.set_title(\"Damped sound wave: velocity at antinode\")\n", "ax.legend()\n", "ax.grid(True)\n", "\n", "# --- Right: log(maxima) vs time with linear fit ---\n", "ax = axes[1]\n", "ax.plot(t_maxima, maxima, \"ro\", markersize=6, label=r\"$\\ln|A_\\mathrm{max}|$\")\n", "ax.plot(times, np.polyval(linfit, times), \"--\",\n", " label=rf\"fit: $\\gamma={gamma_numerical:.3f}$\")\n", "ax.axline(\n", " (0, np.log(np.abs(A0) + 1e-15)),\n", " slope=gamma_analytical,\n", " color=\"k\",\n", " linestyle=\":\",\n", " label=rf\"analytical: $\\gamma={gamma_analytical:.3f}$\",\n", ")\n", "ax.set_xlabel(\"time\")\n", "ax.set_ylabel(r\"$\\ln|A|$\")\n", "ax.set_title(\"Decay rate: numerical vs analytical\")\n", "ax.legend()\n", "ax.grid(True)\n", "\n", "fig.suptitle(\n", " rf\"Viscous damping of sound wave ($\\mu={mu}$, $k={k:.3f}$, $\\gamma_\\mathrm{{anal}}={gamma_analytical:.4f}$)\",\n", " fontsize=12,\n", ")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "39", "metadata": {}, "source": [ "### Verification Check" ] }, { "cell_type": "code", "execution_count": null, "id": "40", "metadata": {}, "outputs": [], "source": [ "tolerance = 0.16 # 16% relative error\n", "\n", "print(f\"=== Damped Sound Wave Verification (tolerance {tolerance*100:.0f}%) ===\")\n", "print(f\" Analytical gamma = {gamma_analytical:.4f}\")\n", "print(f\" Numerical gamma = {gamma_numerical:.4f}\")\n", "print(f\" Relative error = {rel_error * 100:.2f}%\")\n", "\n", "try:\n", " assert rel_error < tolerance, (\n", " f\"Numerical decay rate {gamma_numerical:.4f} deviates {rel_error*100:.1f}% \"\n", " f\"from analytical {gamma_analytical:.4f} (tolerance {tolerance*100:.0f}%)\"\n", " )\n", " print(f\"\\n✓ Verification passed — decay rate within {tolerance*100:.0f}% of analytical.\")\n", "except AssertionError as e:\n", " print(f\"\\n✗ {e}\")\n", "\n", "# Optional cleanup\n", "if False: # set to True to remove simulation output\n", " shutil.rmtree(test_folder)\n", " print(f\"Cleaned up {test_folder}\")" ] }, { "cell_type": "markdown", "id": "41", "metadata": {}, "source": [ "### Conclusion\n", "\n", "This section verified the viscous propagator of `ViscousEulerSPH` against the analytical damping rate of a compressible sound wave:\n", "\n", "- The SPH discretisation correctly reproduces the $\\gamma = -(4/3)\\mu k^2/2$ decay law to within the 16% tolerance at the chosen resolution.\n", "- The decay rate is extracted robustly from the envelope of the oscillating velocity signal — the same technique used for Landau damping in kinetic simulations.\n", "- Increasing `nx` or `ppb` reduces the relative error further, as the SPH kernel gradient approximation improves with particle density." ] } ], "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 }