Hybrid MHD-kinetic energy exchange

Energetic ions drifting near the Alfvén speed exchange energy with a shear-Alfvén wave through Struphy’s pressure-coupling scheme, while the total energy stays conserved. The energetic-ion Maxwellian is specified by $n_h=0.1,\qquad u_{z,h}=1,\qquad v_{\mathrm{th},z,h}=0.3,$ with transverse thermal speeds $v_{\mathrm{th},x,h}=v_{\mathrm{th},y,h}=1$ . The bulk has $\mathbf{B}_0=\mathbf{e}_z$ , $n_0=1$ and transverse velocity coefficient-noise amplitude $10^{-3}$ .

Field and fluid energy oscillating as they exchange energy with the energetic ion population — Actual output from the script below. Interactive controls are available in the plot above.Open graph in a new tab ↗

Physical problem

Two-way kinetic-MHD coupling

Model: LinearMHDVlasovPC
Domain: Cuboid (r3=20)
Grid: 1 × 1 × 32
FEEC degree: 1 × 1 × 3
Integrator: Analytic
Evolution: 2,000 steps
Visualization: Plotly

PDEs solved by model:

MHD continuity:

\frac{\partial ρ ~}{\partial t} + \nabla \cdot (ρ_{0} U ~) = 0

MHD momentum:

ρ_{0} \frac{\partial U ~}{\partial t} + \nabla p ~ + \frac{A _{h}}{A _{b}} \nabla \cdot P ~_{h, ⊥} = (\nabla \times B ~) \times B_{0} + J_{0} \times B ~

J_{0} = \nabla \times B_{0}

MHD pressure:

\frac{\partial p ~}{\partial t} + \nabla \cdot (p_{0} U ~) + \frac{2}{3} p_{0} \nabla \cdot U ~ = 0

MHD induction:

\frac{\partial B ~}{\partial t} - \nabla \times (U ~ \times B_{0}) = 0

Energetic-particle Vlasov equation:

\frac{\partial f _{h}}{\partial t} + (v + U ~_{⊥}) \cdot \nabla f_{h} + [\frac{1}{ϵ} v \times (B_{0} + B ~) - \nabla U ~_{⊥} \cdot v] \cdot \frac{\partial f _{h}}{\partial v} = 0

Perpendicular pressure tensor:

P ~_{h, ⊥} = \int v_{⊥} v_{⊥}^{⊤} f_{h} d v

Simulation geometry

Domain used in this run

The actual Cuboid geometry reconstructed from this run's metadata. Switch between the 3D view and physical-coordinate slices.

Domain

Cuboid

Documentation ↗

Loading Cuboid…

Generated and exported by Struphy.Drag to rotate · Scroll to zoom · Shift-drag to pan

Parameters

Parameter	Value	Description
Loading parameters…

Performance

Profiling

Download profile .h5

Struphy's built-in scope-profiler instruments every propagator, pusher and solver call — these are the real per-call timings from the run above.

Durations

Gantt (rank 0)

Region statistics

Region	Count	Avg (s)	Min (s)	Max (s)	Total (s)

Complete source

Simulation script

Download .py

Run it locally

1. Install and compile

Use Python 3.10 or newer with Struphy's system prerequisites installed. Run these commands in your Python environment:

pip install "struphy[pproc]>=3.4.0"
pip install "plasma-plots[plotly]>=0.1.1"
struphy compile

PNG exports need Chrome or Chromium. If Chrome is not installed, run kaleido_get_chrome.

2. Run the script

Download the .py file above, open a terminal in its directory, then run:

python hybrid-alfven-ion-coupling.py

Figures are saved in the current directory. Add --show to display them before saving.

"""Hybrid MHD-kinetic energy exchange with Struphy's LinearMHDVlasovPC model.

A population of energetic ions, drifting close to the Alfvén speed along a
uniform background field, is evolved together with a linear-MHD shear-Alfvén
wave via the pressure-coupling scheme. The two subsystems continuously trade
energy back and forth through the coupling term while the discretization
conserves the total -- the defining property of a structure-preserving hybrid
scheme.

Install dependencies and compile the kernels (Python 3.10 or newer):

    pip install "struphy[pproc]>=3.4.0"
    pip install "plasma-plots[plotly]>=0.1.1"
    struphy compile

PNG exports require Chrome or Chromium. If Chrome is not installed, run:

    kaleido_get_chrome

Save this file and run it from the directory where you want the output:

    python hybrid-alfven-ion-coupling.py

Figures are saved in the current directory; add --show to display them before saving.
"""

import argparse

import numpy as np
import plotly.graph_objects as go

from struphy import (
    BoundaryParameters,
    DerhamOptions,
    EnvironmentOptions,
    FieldsBackground,
    LoadingParameters,
    SavingParameters,
    Simulation,
    SortingParameters,
    Time,
    WeightsParameters,
    domains,
    equils,
    grids,
    maxwellians,
    perturbations,
)
from struphy.models import LinearMHDVlasovPC
from plasma_plots import save_figure


def create_simulation() -> Simulation:
    model = LinearMHDVlasovPC()
    model.em_fields.b_field.save_data = True

    # A periodic domain aligned with the background field, so energetic ions can
    # free-stream along B0 and resonate with the Alfvén wave.
    domain = domains.Cuboid(r3=20.0)
    equil = equils.HomogenSlab(B0z=1.0, n0=1.0, beta=0.1)
    grid = grids.TensorProductGrid(num_elements=(1, 1, 32))
    derham_opts = DerhamOptions(degree=(1, 1, 3))
    time_opts = Time(dt=0.02, Tend=40.0)

    # A modest population of energetic ions, sorted into boxes for efficient
    # accumulation of the pressure-coupling term onto the field grid.
    model.energetic_ions.set_markers(
        loading_params=LoadingParameters(Np=8000, seed=1),
        weights_params=WeightsParameters(),
        boundary_params=BoundaryParameters(),
        sorting_params=SortingParameters(boxes_per_dim=(1, 1, 8), do_sort=True),
        saving_params=SavingParameters(),
    )

    model.propagators.push_eta_pc.options = model.propagators.push_eta_pc.Options()
    model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
    model.propagators.pc6d.options = model.propagators.pc6d.Options()
    model.propagators.shearalfven.options = model.propagators.shearalfven.Options()
    model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options()

    # A small-amplitude shear-Alfvén seed, and energetic ions drifting near the
    # resonant parallel velocity u3 = v_A = 1 with a finite thermal spread.
    model.mhd.velocity.add_background(FieldsBackground())
    model.mhd.velocity.add_perturbation(perturbations.Noise(amp=1e-3, comp=0, seed=42))
    model.energetic_ions.var.add_background(
        maxwellians.Maxwellian3D(n=(0.1, None), u3=(1.0, None), vth3=(0.3, None)),
    )

    env = EnvironmentOptions(
        out_folders="struphy_gallery_runs",
        sim_folder="hybrid_alfven_ion_coupling",
    )
    sim = Simulation(
        model=model,
        name="Hybrid MHD-kinetic energy exchange",
        description=(
            "Energetic ions drifting near the Alfvén speed exchange energy with a "
            "shear-Alfvén wave through Struphy’s pressure-coupling scheme, while "
            "the total energy stays conserved."
            r" The energetic-ion Maxwellian is specified by $$n_h=0.1,\qquad u_{z,h}=1,\qquad v_{\mathrm{th},z,h}=0.3,$$"
            r" with transverse thermal speeds :math:`v_{\mathrm{th},x,h}=v_{\mathrm{th},y,h}=1`."
            r" The bulk has :math:`\mathbf{B}_0=\mathbf{e}_z`, :math:`n_0=1` and transverse velocity coefficient-noise amplitude :math:`10^{-3}`."
        ),
        env=env,
        time_opts=time_opts,
        domain=domain,
        equil=equil,
        grid=grid,
        derham_opts=derham_opts,
    )
    return sim


def pproc(sim: Simulation, show: bool = False):
    output = sim.output
    output.pproc()  # on every rank; under MPI the output is not processed implicitly

    # LinearMHDVlasovPC tracks each subsystem's energy as a scalar every step:
    # en_B/en_U (field/fluid), en_f (kinetic energetic ions), en_tot (total).
    scalars = output.scalars
    time = np.asarray(output.time)
    en_B = np.asarray(scalars["en_B"])
    en_U = np.asarray(scalars["en_U"])
    en_f_change = np.asarray(scalars["en_f"].plasma.analysis.drift())
    total_drift = np.asarray(scalars["en_tot"].plasma.analysis.drift() / scalars["en_tot"].isel(t=0))

    relative_drift = float(scalars["en_tot"].plasma.analysis.relative_error().max())
    print(f"Max relative drift in total energy (should be ~0): {relative_drift:.2e}")

    figure = go.Figure(
        data=[
            go.Scatter(x=time, y=en_B, mode="lines", name="Field energy (en_B)", line={"color": "#168aad", "width": 2.5}),
            go.Scatter(x=time, y=en_U, mode="lines", name="Fluid kinetic energy (en_U)", line={"color": "#f77f00", "width": 2.5}),
            go.Scatter(x=time, y=en_f_change, mode="lines", name="Energetic-ion energy change (en_f − en_f₀)", line={"color": "#d62828", "width": 2.5}),
            go.Scatter(
                x=time,
                y=total_drift,
                mode="lines",
                name="Total energy drift (relative)",
                line={"color": "#6a4c93", "width": 2, "dash": "dot"},
                yaxis="y2",
            ),
        ],
    )
    figure.update_layout(
        title="Hybrid MHD-kinetic energy exchange",
        xaxis_title="t [a.u.]",
        yaxis_title="Energy [a.u.]",
        yaxis2={
            "title": "Relative total-energy drift",
            "overlaying": "y",
            "side": "right",
            "showgrid": False,
        },
        template="plotly_white",
        autosize=True,
        legend={"x": 0.02, "y": 0.98, "bgcolor": "rgba(255,255,255,0.82)"},
        margin={"l": 70, "r": 70, "t": 80, "b": 60},
    )

    save_figure(figure, "hybrid-alfven-ion-coupling", show=show)


if __name__ == "__main__":
    argparser = argparse.ArgumentParser(description="Run the hybrid alfven ion coupling example.")
    argparser.add_argument(
        "--pproc-only",
        action="store_true",
        help="Run post-processing on an existing simulation instead of running a new one.",
    )
    argparser.add_argument("--show", action="store_true", help="Show the figures before saving them.")
    args = argparser.parse_args()

    simulation = create_simulation()
    if not args.pproc_only:
        simulation.run()
    pproc(simulation, show=args.show)

See the post-processing guide for more ways to inspect the result.