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.

[1]:

import numpy as np
import matplotlib.pyplot as plt

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()

[1]:

<matplotlib.legend.Legend at 0x7fbe9c9840a0>
../../_images/_collections_tutorials_tutorial_gas_expansion_sph_1_1.png

Step 1: Import APIs and the ViscousEulerSPH Model#

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

[2]:

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

# import model
from struphy.models import ViscousEulerSPH

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

[3]:

ViscousEulerSPH.pde()
PDEs solved by model:

Continuity:

∂ₜ ρ + ∇ · (ρ 𝐮) = 0

Momentum:

ρ (∂ₜ 𝐮 + 𝐮 · ∇ 𝐮) = -∇ ( ρ² ∂𝒰(ρ, S)∂ρ ) - ∇ · π

Entropy:

∂ₜ S + 𝐮 · ∇ S = 0

where S denotes the entropy per unit mass and π is the viscous stress tensor.

The viscous stress tensor for a Newtonian fluid is given by:

σ = -μ ( ∇ 𝐮 + (∇ 𝐮)ᵀ - 23 (∇ · 𝐮) 𝐈 )

where μ is the dynamic (shear) viscosity and 𝐈 is the identity tensor.

The internal energy per unit mass can be defined in two ways:

isothermal:    𝒰(ρ, S) = κ(S) \log ρ

polytropic:    𝒰(ρ, S) = κ(S) ργ - 1γ - 1

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.

[4]:

# 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.

[5]:

# 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.

[6]:

# 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.

[7]:

# 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.

[8]:

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.

[9]:

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.

[10]:

# 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.

[11]:

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

[12]:

out = sim.run()

Time stepping: 100%|██████████| 40/40 [00:21<00:00,  1.85step/s]

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.

[13]:

# 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]


No post-processed data in /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1, processing with default options (call out.pproc(...) to choose them)

Post-processing path /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1
Skipping euler_fluid.var backgrounds: Cannot restore initial-condition function: lambdas, nested functions, and closures are not supported.

No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 4047 marker orbits for euler_fluid

100%|██████████| 41/41 [00:00<00:00, 73.27it/s]

Evaluation of distribution functions for euler_fluid

0 starting post-processing of distribution functions for /home/runner/work/struphy/struphy/doc/_collections/tutorials/sim_1/post_processing/kinetic_data/euler_fluid ...
100%|██████████| 1/1 [00:00<00:00, 1281.49it/s]
  0%|          | 0/1 [00:00<?, ?it/s]rank = 0 ----------------------------
self._pproc_rank =0 with xp.sum(data) =np.float64(76.78005716943179) and xp.sum(data_df) =np.float64(76.78005716943179)
self._pproc_rank =0 with xp.sum(data) =np.float64(76.78005716943179) and xp.sum(data_df) =np.float64(76.78005716943179)
self._pproc_rank =0 done.
100%|██████████| 1/1 [00:00<00:00, 21.92it/s]

Evaluation of sph density for euler_fluid


100%|██████████| 1/1 [00:00<00:00, 20.93it/s]

[14]:

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()

[14]:

<matplotlib.colorbar.Colorbar at 0x7fbe33b3ee90>
../../_images/_collections_tutorials_tutorial_gas_expansion_sph_25_1.png 
[15]:

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

i=0
i=4
i=8
i=12
i=16
i=20
i=24
i=28
../../_images/_collections_tutorials_tutorial_gas_expansion_sph_26_1.png