Gas Expansion with ViscousEulerSPH

This tutorial solves Euler-type fluid dynamics with SPH using the ViscousEulerSPH model. The continuum equations are

\[\begin{split}\begin{aligned} \partial_t \rho + \nabla \cdot (\rho \mathbf u) &= 0\,, \\[2mm] \rho(\partial_t \mathbf u + \mathbf u \cdot \nabla \mathbf u) &= - \nabla \left(\rho^2 \frac{\partial \mathcal U(\rho, S)}{\partial \rho} \right)\,, \\[2mm] \partial_t S + \mathbf u \cdot \nabla S &= 0\,. \end{aligned}\end{split}\]

Here, \(S\) is entropy per unit mass and the internal energy model is

\[\mathcal U(\rho, S) = \kappa(S)\log\rho\,.\]

After SPH discretization, particle trajectories satisfy

\[\begin{split}\begin{aligned} \dot{\mathbf x}_p &= \mathbf v_p\,,\qquad && \mathbf x_p(0) = \mathbf x_{p0}\,, \\[2mm] \dot{\mathbf v}_p &= -\kappa_p(0) \sum_{i=1}^N w_i \left(\frac{1}{\rho^{N,h}(\mathbf x_p)} + \frac{1}{\rho^{N,h}(\mathbf x_i)} \right) \nabla W_h(\mathbf x_p-\mathbf x_i)\,, \end{aligned}\end{split}\]

with smoothed density

\[\rho^{N,h}(\mathbf x)=\sum_{j=1}^N w_j W_h(\mathbf x-\mathbf x_j)\,.\]

In Struphy, the pressure-driven velocity update is handled by PushVinSPHpressure.

Before launching the full gas-expansion case, let us quickly inspect a few available smoothing kernels.

import matplotlib.pyplot as plt
import numpy as np

from struphy.pic.sph_smoothing_kernels import gaussian_uni, linear_uni, trigonometric_uni

x = np.linspace(-1, 1, 200)
out1 = np.zeros_like(x)
out2 = np.zeros_like(x)
out3 = np.zeros_like(x)

for i, xi in enumerate(x):
    out1[i] = trigonometric_uni(xi, 1.0)
    out2[i] = gaussian_uni(xi, 1.0)
    out3[i] = linear_uni(xi, 1.0)
plt.plot(x, out1, label="trigonometric")
plt.plot(x, out2, label="gaussian")
plt.plot(x, out3, label="linear")
plt.title("Some smoothing kernels")
plt.legend()

Step 1: Import APIs and the ViscousEulerSPH Model

Load all Struphy components and plotting utilities needed for this standalone gas-expansion tutorial.

from struphy import (
    BinningPlot,
    BoundaryParameters,
    DerhamOptions,
    EnvironmentOptions,
    KernelDensityPlot,
    LoadingParameters,
    SavingParameters,
    Simulation,
    SortingParameters,
    Time,
    WeightsParameters,
    domains,
    equils,
    grids,
)

# import model
from struphy.models import ViscousEulerSPH

Let us verfiy the model equations by calling the pde method:

ViscousEulerSPH.pde()

Step 2: Define Environment and Domain

Set time step, final time, and simulation domain. These settings are tuned for a 2D expansion embedded in a cuboid geometry.

# environment options
env = EnvironmentOptions()

# time stepping
time_opts = Time(dt=0.04, Tend=1.6, split_algo="Strang")

# geometry
l1 = -3.0
r1 = 3.0
l2 = -3.0
r2 = 3.0
l3 = 0.0
r3 = 1.0
domain = domains.Cuboid(l1=l1, r1=r1, l2=l2, r2=r2, l3=l3, r3=r3)

Step 3: Build the Initial Density Blob

Define a Gaussian density profile in the \((x,y)\) plane and wrap it as a fluid equilibrium object. This serves as the initial condition for the expansion.

# gaussian initial blob
import numpy as np

T_h = 0.2
gamma = 5 / 3
n_fun = lambda x, y, z: np.exp(-(x**2 + y**2) / T_h) / 35

blob = equils.GenericCartesianFluidEquilibrium(n_xyz=n_fun)

Step 4: Set Grid and de Rham Discretization

Choose the spline grid and de Rham options used by the field projections in this SPH run.

# fluid equilibrium (can be used as part of initial conditions)
equil = None

# grid
grid = grids.TensorProductGrid(num_elements=(64, 64, 1))

# derham options
bcs = (("free", "free"), ("free", "free"), None)
derham_opts = DerhamOptions(degree=(3, 3, 1), bcs=bcs)

Step 5: Instantiate ViscousEulerSPH and Simulation

Create a lightweight model with magnetic-background and viscosity terms disabled (with_B0=False, with_viscosity=False), then build the simulation object.

# light-weight model instance
model = ViscousEulerSPH(with_B0=False, with_viscosity=False)

sim = Simulation(
    model,
    env=env,
    time_opts=time_opts,
    domain=domain,
    equil=equil,
    grid=grid,
    derham_opts=derham_opts,
)

Step 6: Configure Marker Sampling

Set dense particle-per-box loading and enable reject_weights to discard very small particle weights. This reduces cost while preserving the dominant density structure.

loading_params = LoadingParameters(ppb=400)
weights_params = WeightsParameters(reject_weights=True, threshold=3e-8)
boundary_params = BoundaryParameters()
nx = 16
ny = 16
sorting_params = SortingParameters(boxes_per_dim=(nx, ny, 1))

Step 7: Add Diagnostics for Visualization

Save both a binning plot and an SPH kernel-density plot so you can compare coarse gridded density against kernel-based reconstruction.

bin_plot = BinningPlot(
    slice="e1_e2",
    n_bins=(64, 64),
    ranges=((0.0, 1.0), (0.0, 1.0)),
    divide_by_jac=False,
)
pts_e1 = 100
pts_e2 = 90
kd_plot = KernelDensityPlot(pts_e1=pts_e1, pts_e2=pts_e2, pts_e3=1)

saving_params = SavingParameters(
    n_markers=1.0,
    binning_plots=(bin_plot,),
    kernel_density_plots=(kd_plot,),
)

model.euler_fluid.set_markers(
    loading_params=loading_params,
    weights_params=weights_params,
    boundary_params=boundary_params,
    sorting_params=sorting_params,
    saving_params=saving_params,
)

Step 8: Select the SPH Kernel for Pressure Updates

For the pressure propagator, choose gaussian_2d as smoothing kernel. You can swap this setting later to study kernel sensitivity.

# propagator options
from struphy import ButcherTableau

butcher = ButcherTableau(algo="forward_euler")
model.propagators.push_eta.options = model.propagators.push_eta.Options(butcher=butcher)

model.propagators.push_sph_p.options = model.propagators.push_sph_p.Options(kernel_type="gaussian_2d")

Step 9: Initialize, Run, and Load Results

Attach the Gaussian background and execute the standard workflow: out = sim.run().

Performance note: early time steps are typically slower because particles are highly concentrated; runtime usually improves as the cloud expands. Running from a console script (especially with MPI/GPU support) is faster than notebook execution.

# background, perturbations and initial conditions
model.euler_fluid.var.add_background(blob)

out = sim.run()

Step 10: Inspect Stored Outputs

After post-processing, use the printed paths to locate simulation artifacts. The following cells compare analytical initialization, marker sampling, SPH reconstruction, and binned density diagnostics.

# analytical functions
n_xyz = blob.n_xyz
n3 = blob.n3

# grids
x = np.linspace(l1, r1, pts_e1)
y = np.linspace(l2, r2, pts_e2)
xx, yy = np.meshgrid(x, y, indexing="ij")
density = out.evaluate("euler_fluid/n", dataset="view_0/n")
ee1, ee2, ee3 = np.meshgrid(density.eta1, density.eta2, density.eta3, indexing="ij")
eta1 = ee1[:, 0, 0]
eta2 = ee2[0, :, 0]
f_e1e2 = out.evaluate("euler_fluid/f", dataset="e1_e2_density/f")
bc_x = np.asarray(f_e1e2["eta1"])
bc_y = np.asarray(f_e1e2["eta2"])

# markers: one (t, marker) variable per saved quantity (x, y, z, v1, v2, v3, weight)
orbits = out.evaluate("euler_fluid/orbits")
positions = np.stack([orbits.x[0], orbits.y[0], orbits.z[0]], axis=-1)
weights = orbits.weight[0].values

# binning and sph eval
n_sph = np.asarray(density)[0]
f_bin = np.asarray(f_e1e2)[0]

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 15))

# plots
plt.subplot(3, 2, 1)
plt.pcolor(xx, yy, n_fun(xx, yy, 0))
plt.axis("square")
plt.title("n_xyz initial")
plt.colorbar()

plt.subplot(3, 2, 2)
plt.pcolor(eta1, eta2, n3(eta1, eta2, 0, squeeze_out=True).T)
plt.axis("square")
plt.title("$\hat{n}^{\t{vol}}$ initial (volume form)")
plt.colorbar()

make_scatter = True
if make_scatter:
    plt.subplot(3, 2, 3)
    ax = plt.gca()
    ax.set_xticks(np.linspace(l1, r1, nx + 1))
    ax.set_yticks(np.linspace(l2, r2, ny + 1))
    plt.tick_params(labelbottom=False)
    coloring = weights
    plt.scatter(positions[:, 0], positions[:, 1], c=coloring, s=0.25)
    plt.grid(c="k")
    plt.axis("square")
    plt.title("$\hat{n}^{\t{vol}}$ initial scatter (random)")
    plt.xlim(l1, r1)
    plt.ylim(l2, r2)
    plt.colorbar()

plt.subplot(3, 2, 4)
ax = plt.gca()
ax.set_xticks(np.linspace(0, 1, nx + 1))
ax.set_yticks(np.linspace(0, 1.0, ny + 1))
plt.tick_params(labelbottom=False)
plt.pcolor(ee1[:, :, 0], ee2[:, :, 0], n_sph[:, :, 0])
plt.grid()
plt.axis("square")
plt.title("n_sph initial (random)")
plt.colorbar()

plt.subplot(3, 2, 5)
ax = plt.gca()
plt.pcolor(bc_x, bc_y, f_bin)
plt.axis("square")
plt.title("n_binned initial (random)")
plt.colorbar()

dt = time_opts.dt
Nt = out.time.size - 1

positions = np.stack([orbits.x, orbits.y, orbits.z], axis=-1)

interval = Nt / 10
plot_ct = 0

plt.figure(figsize=(12, 24))
for i in range(Nt):
    if i % interval == 0:
        print(f"{i=}")
        plot_ct += 1
        plt.subplot(4, 2, plot_ct)
        ax = plt.gca()
        coloring = weights
        plt.scatter(positions[i, :, 0], positions[i, :, 1], c=coloring, s=0.25)
        plt.axis("square")
        plt.title("n0_scatter")
        plt.xlim(l1, r1)
        plt.ylim(l2, r2)
        plt.colorbar()
        plt.title(f"Gas at t={i * dt}")
    if plot_ct == 8:
        break