Solovev Configuration

Note

For general information about finite element discretization, basis functions, mesh parameters, polynomial degrees, boundary conditions, and matrix/operator dimensions, see Overview.

This script runs a magnetic relaxation simulation for a Solovev configuration. The script is located at scripts/config_scripts/solovev.py.

Problem Statement

The Solovev configuration solves the magnetohydrostatic (MHS) equilibrium equation:

\[\mathbf{J} \times \mathbf{B} = \nabla p \quad \text{in } \Omega\]

with boundary conditions:

\[\mathbf{B} \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega\]

and the constraint:

\[\nabla \cdot \mathbf{B} = 0 \quad \text{in } \Omega\]

where: - \(\mathbf{B}: \Omega \to \mathbb{R}^3\) is the magnetic field (2-form, area form) - \(\mathbf{J} = \nabla \times \mathbf{B}\) is the current density (1-form) - \(p: \Omega \to \mathbb{R}\) is the plasma pressure (0-form, scalar) - \(\mathbf{n}\) is the outward unit normal vector on the boundary - \(\Omega\) is the toroidal plasma domain

For the force-free case (\(\nabla p = 0\)), the equation becomes:

\[\mathbf{J} \times \mathbf{B} = 0 \quad \text{in } \Omega\]

which implies \(\mathbf{J} \parallel \mathbf{B}\) (force-free condition).

Magnetic Relaxation Method

The relaxation method solves the time-dependent MHD equations:

\[\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times \mathbf{E}\]

where the electric field \(\mathbf{E}\) is given by:

\[\mathbf{E} = \mathbf{u} \times \mathbf{B} - \eta \mathbf{J}\]

and: - \(\mathbf{u}\) is the velocity field (divergence-free) - \(\eta\) is the resistivity parameter

The velocity is chosen to drive the system toward equilibrium:

\[\mathbf{u} = \mathbf{f} = \mathbf{J} \times \mathbf{B} - \nabla p\]

where \(\mathbf{f}\) is the force imbalance. The velocity is projected onto the divergence-free space using the Leray projection \(P_{\mathrm{Leray}}\).

Time-Stepping Scheme

The time-stepping uses a midpoint rule (Crank-Nicolson):

\[\mathbf{B}^{n+1} = \mathbf{B}^n + \Delta t \nabla \times \mathbf{E}(\mathbf{B}^{\mathrm{mid}})\]

where \(\mathbf{B}^{\mathrm{mid}} = (\mathbf{B}^n + \mathbf{B}^{n+1})/2\) is evaluated at the midpoint.

The system is solved using a Picard iteration scheme with adaptive time-stepping.

Usage:

python scripts/config_scripts/solovev.py run_name=test_run boundary_type=rotating_ellipse n_r=16 n_theta=16 n_zeta=8

The script accepts various parameters:

  • run_name: Name for the output files

  • boundary_type: Type of boundary (e.g., rotating_ellipse, cerfon)

  • n_r, n_theta, n_zeta: Number of splines in each direction

  • p_r, p_theta, p_zeta: Polynomial degrees

  • maxit: Maximum number of iterations

  • And many more (see script for full list)

The script generates HDF5 files with simulation data and PDF plots showing:

  • Force trace over iterations

  • Pressure contours

  • Magnetic field structure

Discretization Parameters

This script uses a 3D DeRham sequence with mesh parameters \(n_r, n_\theta, n_\zeta\), polynomial degrees \(p_r, p_\theta, p_\zeta\), and boundary conditions (typically clamped in radial direction, periodic in poloidal and toroidal directions). For details on computing the number of DOFs \(N_0, N_1, N_2, N_3\) from these parameters, see Overview.

Finite Element Spaces

The magnetic field \(\mathbf{B}\) is represented as a 2-form (area form) in the finite element space:

\[V_2 = \text{span}\{\Lambda_2^i\}_{i=1}^{N_2}\]

The current density \(\mathbf{J} = (\nabla \times)_h \mathbf{B}\) is a 1-form, where \((\nabla \times)_h: V_2 \to V_1\) is the curl operator (weak form).

Matrix and Operator Dimensions

The 2-form mass matrix \(M_2 \in \mathbb{R}^{N_2 \times N_2}\) and 1-form mass matrix \(M_1 \in \mathbb{R}^{N_1 \times N_1}\) are used. The Leray projection \(P_{\mathrm{Leray}} \in \mathbb{R}^{N_2 \times N_2}\) ensures divergence-free velocity fields. Cross-product projections compute \(\mathbf{J} \times \mathbf{B}\) in the finite element space.

Initial Magnetic Field

The initial magnetic field is constructed from a harmonic eigenmode by solving:

\[(M_2 \Delta_2) \mathbf{v} = \lambda \mathbf{v}\]

The eigenvector \(\mathbf{v}\) corresponding to the smallest eigenvalue is extracted and normalized to unit L2 norm.

Code Walkthrough

Block 1: Imports and Configuration (lines 1-17)

Imports JAX, NumPy, and MRX modules for DeRham sequences, I/O, mappings, relaxation, and plotting. Enables 64-bit precision for numerical accuracy.

Block 2: Main Function (lines 19-41)

The main function parses command-line arguments using parse_args() and updates the default configuration. It then calls run() with the updated configuration. Command-line parameters override defaults using the format parameter_name=value.

Block 3: Setup and Initialization (lines 44-91)

The run() function creates output directory structure, initializes trace dictionary, and selects boundary mapping based on boundary_type:

  • tokamak: Uses cerfon_map (circular cross-section)

  • helix: Uses helical_map (helical boundary)

  • rotating_ellipse: Uses rotating_ellipse_map (rotating elliptical cross-section)

Sets up DeRham sequence with polar coordinates and Dirichlet boundary conditions, assembles all necessary matrices (mass, derivative, cross-product projections, Leray projection), and initializes magnetic field from harmonic eigenmode:

\[\mathbf{B}_{\mathrm{harm}} = \text{eigh}(M_2 \Delta_2)[1][:, 0]\]

The field is normalized to unit L2 norm: \(\hat{\mathbf{B}} = \mathbf{B}_{\mathrm{harm}} / \|\mathbf{B}_{\mathrm{harm}}\|_2\)

Block 4: Magnetic Relaxation Loop (lines 93-94)

Calls run_relaxation_loop() which performs the time-stepping iteration to relax the magnetic field toward equilibrium. This solves the MHD force balance equation:

\[\mathbf{J} \times \mathbf{B} = \nabla p\]

or for the force-free case:

\[\mathbf{J} \times \mathbf{B} = 0\]

where \(\mathbf{J} = \nabla \times \mathbf{B}\) is the current density. The solver uses an implicit time-stepping scheme with Picard iteration.

Block 5: Post-processing and Output (lines 96-122)

After the relaxation loop completes:

  • Computes final pressure field using diagnostics.pressure()

  • Saves all trace data (iterations, forces, energies, etc.) to HDF5 file

  • Generates plots using generate_solovev_plots() which creates: * Pressure contour plots * Force and energy convergence traces * Magnetic field visualizations (if save_B=True)

The Solovev configuration provides a simple analytical equilibrium that serves as a test case for the relaxation solver. The script demonstrates how to set up and run magnetic relaxation simulations for various boundary geometries.