Verification of MHD Wave Dispersion Relations#

1D Slab Waves in a Homogeneous Magnetized Plasma#

This tutorial demonstrates verification of the dispersion relations for the shear Alfvén, slow magnetosonic, and fast magnetosonic waves in a homogeneous magnetized plasma slab using the LinearMHD model.

Physical Setup#

We consider a one-dimensional equilibrium with constant density \(n_0\), pressure \(p_0\), and a uniform background magnetic field

\[\mathbf{B}_0 = (B_{0x}, B_{0y}, B_{0z}).\]

Small-amplitude perturbations are launched with wave vector \(\mathbf{k} = k\,\hat{\mathbf{z}}\), so the three linear MHD branches reduce to the standard slab-wave dispersion relations. In Struphy units, the Alfvén and sound speeds are

\[v_A^2 = \frac{|\mathbf{B}_0|^2}{n_0}, \qquad c_S^2 = \gamma\,\frac{p_0}{n_0}.\]

The three wave branches are:

\[\text{shear Alfv\'en}: \qquad \omega = v_A\,k\,\frac{B_{0z}}{|\mathbf{B}_0|},\]
\[\text{slow magnetosonic}: \qquad \omega = k\sqrt{\frac{1}{2}(c_S^2 + v_A^2)\left(1 - \sqrt{1 - \delta}\right)},\]
\[\text{fast magnetosonic}: \qquad \omega = k\sqrt{\frac{1}{2}(c_S^2 + v_A^2)\left(1 + \sqrt{1 - \delta}\right)},\]

with

\[\delta = \frac{4 B_{0z}^2 c_S^2 v_A^2}{(c_S^2 + v_A^2)^2 |\mathbf{B}_0|^2}.\]

We verify these wave speeds by:

  1. Initializing small velocity perturbations (noise)

  2. Running a transient simulation

  3. Extracting dominant wave frequencies via FFT power spectrum

  4. Comparing fitted wave speeds against analytical predictions

[1]:

import logging
import os
import shutil

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import cunumpy as xp

from struphy import (
    DerhamOptions,
    EnvironmentOptions,
    Simulation,
    Time,
    domains,
    equils,
    grids,
    perturbations,
)
from struphy.models import LinearMHD

logger = logging.getLogger("struphy")

Model and Equilibrium Parameters#

Define the magnetic field configuration, plasma density, and pressure (via plasma beta) that sets up the equilibrium. These parameters determine the wave speeds we will verify.

[2]:

# Magnetic field components
B0x = 0.0
B0y = 1.0
B0z = 1.0

# Plasma parameters
beta = 3.0  # ratio of thermal to magnetic pressure
n0 = 0.7    # reference density

# Compute thermal pressure from beta and magnetic pressure
Bsquare = B0x**2 + B0y**2 + B0z**2
p0 = beta * Bsquare / 2.0

print(f"Magnetic field: B0 = ({B0x}, {B0y}, {B0z})")
print(f"Alfvén speed: v_A = {np.sqrt(Bsquare / n0):.4f}")
print(f"Sound speed: c_s = {np.sqrt(5/3 * p0 / n0):.4f}")
print(f"Beta (plasma): {beta}")

Magnetic field: B0 = (0.0, 1.0, 1.0)
Alfvén speed: v_A = 1.6903
Sound speed: c_s = 2.6726
Beta (plasma): 3.0

Domain and Numerical Discretization#

Set up a 1D periodic domain in the \(z\)-direction with a tensor-product grid. Since the wave propagates along \(z\), we keep \(x\) and \(y\) as single-element.

[3]:

# 1D domain (extended in z)
domain = domains.Cuboid(r3=60.0)  # z ∈ [0, 60)

# Grid: 1D in z with 64 elements (fine enough to resolve waves)
grid = grids.TensorProductGrid(num_elements=(1, 1, 64))

# Derham options for 1D:
# - degree=(1,1,3): cubic edges in z for smooth wave representation
derham_opts = DerhamOptions(degree=(1, 1, 3))

print(f"Domain: z ∈ [0, {domain.params['r3']})")
print(f"Grid elements: {grid.num_elements}")
print(f"Derham degree: {derham_opts.degree}")

Domain: z ∈ [0, 60.0)
Grid elements: (1, 1, 64)
Derham degree: (1, 1, 3)

Model Instantiation and Propagator Options#

Create a LinearMHD model and configure the propagators for the shear Alfvén and magnetosonic branches. We use the implicit time-stepping algorithm for stability over the long simulation window.

[4]:

# Model instance
model = LinearMHD()

# Choose implicit or explicit time-stepping
algo = "implicit"  # Alternative: "explicit"

# Propagator options
model.propagators.shear_alf.options = model.propagators.shear_alf.Options(algo=algo)

Initial Conditions#

Add small random perturbations to velocity components. These perturbations excite a broadband spectrum of waves that we will analyze via FFT. The background equilibrium (magnetic field, density, pressure) is set via the equil object in the Simulation.

[5]:

# Set velocity perturbations (small noise to excite waves)
model.mhd.velocity.add_perturbation(perturbations.Noise(amp=0.1, comp=0, seed=123))
model.mhd.velocity.add_perturbation(perturbations.Noise(amp=0.1, comp=1, seed=123))
model.mhd.velocity.add_perturbation(perturbations.Noise(amp=0.1, comp=2, seed=123))

# Build equilibrium object
equil = equils.HomogenSlab(B0x=B0x, B0y=B0y, B0z=B0z, beta=beta, n0=n0)

Simulation Setup and Execution#

Configure the simulation environment, time-stepping parameters, and run the transient dynamics. We run for a long enough time to collect sufficient wave cycles for accurate FFT analysis.

[6]:

# Environment and file management
test_folder = os.path.join(os.getcwd(), "struphy_verification_tests")
out_folders = os.path.join(test_folder, "LinearMHD")
env = EnvironmentOptions(out_folders=out_folders, sim_folder="slab_waves_1d")

# Time-stepping: dt=0.15 is stable for implicit stepping
time_opts = Time(dt=0.15, Tend=180.0)

# Instantiate and run simulation
sim = Simulation(
    model=model,
    env=env,
    time_opts=time_opts,
    domain=domain,
    grid=grid,
    derham_opts=derham_opts,
    equil=equil,
)

print(f"Running simulation: dt={time_opts.dt}, Tend={time_opts.Tend}")
out = sim.run()
print("Simulation complete.")

Running simulation: dt=0.15, Tend=180.0
Time stepping: 100%|██████████| 1200/1200 [01:12<00:00, 16.55step/s]
Simulation complete.

Diagnostics: FFT Analysis and Wave Speed Verification#

Extract velocity and pressure data over time, compute power spectra, and fit dominant frequencies to extract wave speeds. We compare the fitted speeds against the theoretical predictions for shear Alfvén, slow and fast magnetosonic waves.

Shear Alfvén wave: expected speed \(\approx {0.707} \, v_A\)

Slow and fast magnetosonic waves: derived from dispersion relation involving both \(v_A\) and \(c_s\).

Computing the Dispersion Relation#

A wave \(\propto e^{i(kz - \omega t)}\) shows up as a peak at \((\omega, k)\) in the Fourier transform of the field over time and space. We therefore

  1. take the field along \(z\) at the first \(x\) and \(y\) grid point, giving data \(f(t, z)\) on uniform grids,

  2. compute the space-time power spectrum \(|\hat f(\omega, k)|\) with a 2D FFT and keep the non-negative frequencies and wave numbers,

  3. at each \(k\) (in the band between 1/8 and 1/2 of the resolved wave numbers), locate the local maxima in \(\omega\) that rise above a fraction noise_level of that column’s maximum; each peak belongs to one wave branch,

  4. fit a line \(\omega = v\,k + b\) to each branch. The slope \(v\) is the phase velocity of that wave.

[7]:

from scipy.fft import fft2, fftfreq
from scipy.signal import argrelextrema


def power_spectrum(field, component=0):
    """Space-time power spectrum |F(omega, k)| of a field along z, at the first x and y grid point."""
    if "component" in field.dims:
        field = field.isel(component=component)
    data = field.isel(eta1=0, eta2=0).transpose("t", "eta3")
    time, z = data.t.values, data.Z.values
    nt, nz = data.shape
    power = (2.0 / nt) * (2.0 / nz) * np.abs(fft2(data.values))[: nt // 2, : nz // 2]
    omega = 2 * np.pi * fftfreq(nt, time[1] - time[0])[: nt // 2]
    k = 2 * np.pi * fftfreq(nz, z[1] - z[0])[: nz // 2]
    return omega, k, power


def fit_branches(omega, k, power, n_branches, noise_level, order=10):
    """Fit omega = v * k + b to each of the n_branches spectral peaks; returns [(v, b), ...] sorted by omega."""
    k_fit, omega_fit = [], [[] for _ in range(n_branches)]
    for i in range(k.size // 8, k.size // 2):
        column = power[:, i]
        maxima = argrelextrema(column, np.greater, order=order)[0]
        peaks = sorted(j for j in maxima if column[j] > noise_level * column.max())
        if not peaks:
            continue
        assert len(peaks) == n_branches, (
            f"Found {len(peaks)} branches at k={k[i]:.3f}, expected {n_branches}. "
            "Try another noise_level or order."
        )
        k_fit.append(k[i])
        for branch, j in zip(omega_fit, peaks):
            branch.append(omega[j])
    return [np.polyfit(k_fit, branch, deg=1) for branch in omega_fit]


def plot_spectrum(omega, k, power, fits, theory, title):
    """Normalized power spectrum with the fitted branches (dotted) and the theoretical ones, omega = v * k (dashed)."""
    fig, ax = plt.subplots(figsize=(7, 6))
    normalized = np.maximum(power**2 / (power**2).max(), 1e-15)
    levels = np.logspace(-15, 0, 31)
    mappable = ax.contourf(k, omega, normalized, levels=levels, norm=LogNorm(), cmap="plasma")
    fig.colorbar(mappable, ax=ax, ticks=[1e-12, 1e-9, 1e-6, 1e-3, 1e0], format="%.0e")
    for n, (slope, intercept) in enumerate(fits):
        ax.plot(k, slope * k + intercept, "w:", lw=2, label=f"fit {n + 1}: v = {slope:.4f}")
    for label, speed in theory.items():
        ax.plot(k, speed * k, "--", label=f"{label}: v = {speed:.4f}")
    ax.set(xlabel="$k$", ylabel=r"$\omega$", title=title, xlim=(0, k[-1]), ylim=(0, omega[-1]))
    ax.legend(loc="upper left")
    plt.show()

[8]:

# Field data is post-processed lazily on first access through `out`
# Extract velocity and pressure time-series
u_of_t = out.evaluate("mhd/velocity")
p_of_t = out.evaluate("mhd/pressure")

gamma = 5 / 3  # Adiabatic index

# 1. Shear Alfvén wave analysis from the x-component of the velocity
print("\n=== Shear Alfvén Wave Analysis ===")
omega, k, power = power_spectrum(u_of_t, component=0)
fits_alfven = fit_branches(omega, k, power, n_branches=1, noise_level=0.5)

# Theoretical Alfvén speed
vA = xp.sqrt(Bsquare / n0)
v_alfven_theory = vA * B0z / xp.sqrt(Bsquare)
v_alfven_fit = float(fits_alfven[0][0])

plot_spectrum(omega, k, power, fits_alfven, {"shear Alfvén": v_alfven_theory}, title="$u_x$ power spectrum")

print(f"Théoretical Alfvén speed: {v_alfven_theory:.6f}")
print(f"Fitted Alfvén speed:      {v_alfven_fit:.6f}")
print(f"Relative error:           {abs(v_alfven_fit - v_alfven_theory) / v_alfven_theory * 100:.2f}%")

rel_error_alfven = xp.abs(v_alfven_fit - v_alfven_theory) / v_alfven_theory
assert rel_error_alfven < 0.07, f"Alfvén wave speed relative error {rel_error_alfven:.4f} exceeds tolerance"
print("✓ Alfvén wave verification passed.\n")


=== Shear Alfvén Wave Analysis ===
../../_images/_collections_tutorials_tutorial_linear_mhd_slab_waves_1d_15_1.png 
Théoretical Alfvén speed: 1.195229
Fitted Alfvén speed:      1.145898
Relative error:           4.13%
✓ Alfvén wave verification passed.


[9]:

# 2. Magnetosonic waves analysis from pressure
print("=== Slow and Fast Magnetosonic Wave Analysis ===")
omega, k, power = power_spectrum(p_of_t)
fits_sonic = fit_branches(omega, k, power, n_branches=2, noise_level=0.4)

# Theoretical magnetosonic speeds
cS = xp.sqrt(gamma * p0 / n0)
delta = (4 * B0z**2 * cS**2 * vA**2) / ((cS**2 + vA**2) ** 2 * Bsquare)
v_slow_theory = xp.sqrt(0.5 * (cS**2 + vA**2) * (1.0 - xp.sqrt(1.0 - delta)))
v_fast_theory = xp.sqrt(0.5 * (cS**2 + vA**2) * (1.0 + xp.sqrt(1.0 - delta)))

v_slow_fit = float(fits_sonic[0][0])
v_fast_fit = float(fits_sonic[1][0])

plot_spectrum(
    omega, k, power, fits_sonic,
    {"slow magnetosonic": v_slow_theory, "fast magnetosonic": v_fast_theory},
    title="$p$ power spectrum",
)

print("\nSlow Magnetosonic Wave:")
print(f"  Théoretical speed: {v_slow_theory:.6f}")
print(f"  Fitted speed:      {v_slow_fit:.6f}")
print(f"  Relative error:    {abs(v_slow_fit - v_slow_theory) / v_slow_theory * 100:.2f}%")

print("\nFast Magnetosonic Wave:")
print(f"  Théoretical speed: {v_fast_theory:.6f}")
print(f"  Fitted speed:      {v_fast_fit:.6f}")
print(f"  Relative error:    {abs(v_fast_fit - v_fast_theory) / v_fast_theory * 100:.2f}%")

rel_error_slow = xp.abs(v_slow_fit - v_slow_theory) / v_slow_theory
rel_error_fast = xp.abs(v_fast_fit - v_fast_theory) / v_fast_theory

assert rel_error_slow < 0.05, f"Slow wave speed relative error {rel_error_slow:.4f} exceeds tolerance"
assert rel_error_fast < 0.19, f"Fast wave speed relative error {rel_error_fast:.4f} exceeds tolerance"
print("\n✓ Magnetosonic wave verification passed.")

=== Slow and Fast Magnetosonic Wave Analysis ===
../../_images/_collections_tutorials_tutorial_linear_mhd_slab_waves_1d_16_1.png 

Slow Magnetosonic Wave:
  Théoretical speed: 1.073990
  Fitted speed:      1.052736
  Relative error:    1.98%

Fast Magnetosonic Wave:
  Théoretical speed: 2.974314
  Fitted speed:      3.221092
  Relative error:    8.30%

✓ Magnetosonic wave verification passed.

Conclusion#

This tutorial successfully verified the dispersion relations for three branches of MHD waves in a homogeneous magnetized plasma:

  1. Shear Alfvén wave propagates perpendicular to the ambient magnetic field at the magnetic tension-driven speed.

  2. Slow magnetosonic wave is a hybrid mode combining acoustic and magnetic effects, with lower phase velocity.

  3. Fast magnetosonic wave is a magnetically-dominated mode with higher phase velocity, approaching \(\sqrt{c_s^2 + v_A^2}\) at high frequencies.

The FFT-based power spectrum analysis provides a robust method to extract wave speeds from simulation data and verify the model’s fidelity to the underlying MHD theory. This verification validates both the physical model implementation and the numerical discretization.

[10]:

# Cleanup temporary simulation folder
if False:  # Set to True to enable cleanup
    try:
        shutil.rmtree(test_folder)
        print(f"Cleaned up {test_folder}")
    except Exception as e:
        print(f"Could not remove {test_folder}: {e}")