Magnetohydrodynamics (MHD)

Solve fluid-electromagnetic coupling problems in plasma physics and astrophysics

Overview

Magnetohydrodynamics (MHD) studies electrically conducting fluids where fluid flow and magnetic fields interact. iTensor provides numerical solvers designed specifically for MHD problems in astrophysics, fusion research, and plasma physics.

The MHD module combines fluid dynamics equations with Maxwell's electromagnetic equations to model systems like stellar magnetospheres, accretion disks, and laboratory plasmas.

iTensor's MHD solver seamlessly connects with its tensor infrastructure, allowing calculations in arbitrary coordinate systems and curved spacetime metrics.

Core MHD Equations

The MHD module implements the following coupled partial differential equations:

Mass Conservation

∂ρ/∂t + ∇·(ρv) = 0

Where ρ is mass density and v is fluid velocity. This equation tracks the flow of mass in the system.

Momentum Equation

ρ(∂v/∂t + v·∇v) = -∇p + J×B + ρg

This equation balances fluid acceleration against pressure gradients (∇p), Lorentz forces (J×B), and gravitational forces (ρg).

Induction Equation

∂B/∂t = ∇×(v×B - η∇×B)

Describes magnetic field evolution. In ideal MHD (η = 0), field lines are "frozen" into the fluid.

Magnetic Divergence Constraint

∇·B = 0

States that magnetic field lines have no sources or sinks (no magnetic monopoles).

MHD Solver Types

iTensor implements several types of MHD solvers, each optimized for different regimes:

Ideal MHD

Assumes perfect conductivity (η = 0). Field lines are frozen into the fluid, making this suitable for many large-scale astrophysical phenomena.

MHDSolver(type='ideal')

Resistive MHD

Includes finite resistivity (η > 0), allowing magnetic reconnection and Ohmic dissipation. Essential for studying magnetic energy conversion to heat.

MHDSolver(type='resistive', eta=1e-3)

Hall MHD

Includes the Hall term for systems where ion and electron dynamics decouple. Important for small-scale plasma phenomena where ion inertia matters.

MHDSolver(type='hall', ion_skin_depth=0.1)

Relativistic MHD

For systems with relativistic fluid velocities, like jets from black holes or neutron stars. Handles Lorentz transformations automatically.

MHDSolver(type='relativistic')

Numerical Methods

iTensor's MHD solvers use specialized numerical techniques to maintain stability and accuracy:

Finite Volume Method - Ensures conservation of mass, momentum, and energy
Constrained Transport - Maintains ∇·B = 0 to machine precision
Riemann Solvers - Computes accurate fluxes at cell interfaces
- HLLD solver for MHD-specific shock capturing
- Roe solver for accurate treatment of all wave modes
Adaptive Mesh Refinement - Concentrates resolution where needed
Strong Stability Preserving Time Integration - RK3-SSP for robustness

Method Selection

Choose different numerical methods based on your problem requirements:

from itensor.mhd import MHDSolver # Creating a solver with specific numerical methods solver = MHDSolver( spatial_method='finite_volume', riemann_solver='hlld', reconstruction='weno5', time_integrator='ssp_rk3', div_cleaning='constrained_transport' )

Examples

Orszag-Tang Vortex

A standard test problem for MHD codes that develops complex turbulent flows from simple initial conditions:

from itensor.mhd import MHDSolver, UniformGrid2D import numpy as np # Create a 2D grid grid = UniformGrid2D(nx=512, ny=512, xmin=0, xmax=1, ymin=0, ymax=1) # Initialize solver solver = MHDSolver(grid=grid, gamma=5/3) # gamma is adiabatic index # Set initial conditions for Orszag-Tang vortex x, y = grid.get_coordinates() X, Y = np.meshgrid(x, y) # Velocity field vx = -np.sin(2*np.pi*Y) vy = np.sin(2*np.pi*X) # Magnetic field Bx = -np.sin(2*np.pi*Y) By = np.sin(4*np.pi*X) # Density and pressure rho = 25/(36*np.pi) p = 5/(12*np.pi) # Set the initial state solver.set_state( density=rho, velocity_x=vx, velocity_y=vy, magnetic_x=Bx, magnetic_y=By, pressure=p ) # Run simulation solver.evolve( final_time=1.0, cfl=0.4, # Courant number for time step control output_interval=0.1 ) # Visualize final state solver.plot_field('density', title='Orszag-Tang Density at t=1.0') solver.plot_field('magnetic_pressure', title='Magnetic Pressure at t=1.0') solver.plot_field_lines('magnetic', background='vorticity')

Initial Density

[Density plot: t=0.0]

Intermediate Stage

[Density plot: t=0.5]

Final State

[Density plot: t=1.0]

Magnetic Reconnection

Simulation of magnetic field line reconnection in a current sheet:

from itensor.mhd import MHDSolver, UniformGrid2D import numpy as np # Create a grid focused on the reconnection region grid = UniformGrid2D(nx=400, ny=200, xmin=-10, xmax=10, ymin=-5, ymax=5) # Initialize resistive MHD solver solver = MHDSolver(grid=grid, type='resistive', eta=1e-3) # Set up a Harris current sheet with perturbation x, y = grid.get_coordinates() X, Y = np.meshgrid(x, y) # Magnetic field with hyperbolic tangent profile B0 = 1.0 L = 0.5 # Current sheet width Bx = B0 * np.tanh(Y/L) By = 0.1 * B0 * np.sin(2*np.pi*X/10) # Small perturbation # Pressure balance p0 = 0.5 p = p0 + (B0**2 - Bx**2)/2 # Increased pressure in current sheet # Set the initial state solver.set_state( density=1.0, # Uniform density velocity_x=0.0, # Initially at rest velocity_y=0.0, magnetic_x=Bx, # Current sheet profile magnetic_y=By, # Perturbation pressure=p # Pressure balance ) # Run the simulation solver.evolve( final_time=50.0, cfl=0.4, output_interval=5.0 ) # Visualize results solver.animate_field('current_density', interval=200) solver.plot_time_series('reconnection_rate')

Boundary Conditions

iTensor supports multiple boundary condition types for MHD simulations:

Physical Boundaries

outflow - Zero-gradient for all variables
reflecting - Mirrors velocity, preserves normal B
conducting_wall - Enforces B·n = 0 (perfect conductor)
fixed - Holds variables at specified values

Periodic Boundaries

Wraps the computational domain, connecting opposite edges. Useful for studying instabilities and turbulence without boundary artifacts.

solver.set_boundary('periodic', ['x', 'y'])

# Setting different boundary conditions for each edge solver.set_boundary('fixed', 'left', density=1.0, velocity_x=1.0, magnetic_y=2.0) solver.set_boundary('outflow', 'right') solver.set_boundary('reflecting', ['top', 'bottom'])

Analysis Tools

The MHD module includes specialized diagnostics for plasma analysis:

Field Quantities

current_density - J = ∇×B
vorticity - ω = ∇×v
magnetic_pressure - p_m = |B|²/2
beta - β = p/p_m (plasma beta)
alfven_speed - v_A = |B|/√ρ

Integral Diagnostics

total_energy() - Sum of kinetic, thermal, magnetic energy
magnetic_helicity() - H_m = ∫A·B dV
cross_helicity() - H_c = ∫v·B dV
reconnection_rate() - Rate of flux reconnection

# Tracking conservation laws energy_history = [] magnetic_helicity_history = [] for step in range(100): solver.evolve(dt=0.1) energy_history.append(solver.total_energy()) magnetic_helicity_history.append(solver.magnetic_helicity()) # Plot energy conservation if step % 10 == 0: solver.plot_conservation()

Use these tools to verify that your simulation is behaving physically, conserving appropriate quantities, and correctly resolving important MHD phenomena.

Applications

Solar Physics

Model solar flares, coronal mass ejections, and magnetic field evolution in the solar atmosphere.

Solar Flare Example →

Accretion Disks

Investigate magnetorotational instability and angular momentum transport around compact objects.

Accretion Disk Example →

Fusion Research

Simulate plasma instabilities in tokamak and stellarator configurations for fusion energy research.

Tokamak Instability Example →

Astrophysical Jets

Study the formation and propagation of relativistic jets from active galactic nuclei and stellar objects.

Jet Propagation Example →

Previous: Differential Equations Next: Visualization Tools