API Reference

Abstract type for scalar fields (temperature, composition, etc.)

AdvancedThreadManager

Advanced thread management with CPU affinity, NUMA awareness, and work-stealing.

AsyncCommManager{T}

Advanced asynchronous communication manager for overlapping computation and communication.

CPUParallelizer{T}

Advanced CPU parallelization system with SIMD, NUMA, and task-based parallelism.

CombinerConfig

Configuration structure for field combination process.

DynamicLoadBalancer

Dynamic load balancing system that adapts to computational heterogeneity.

EAB2ALUCacheEntry{T}

Type-stable cache entry for EAB2 (Exponential Adams-Bashforth 2) method. Stores banded matrices and their LU factorizations for each spherical harmonic degree l.

ERK2Cache{T}

Cached data structure for Exponential 2nd Order Runge-Kutta method. Stores precomputed matrix exponentials and φ functions for each spherical harmonic mode.

EnergyTracker

Structure for tracking energy conservation during simulation.

FieldCombiner{T}

Structure for combining distributed GeoDynamo.jl output files using the consistent field structures and parameter system.

HybridParallelizer{T}

Coordinates MPI and threads for maximum CPU parallelization.

MasterParallelizer{T}

Comprehensive parallelization system combining all techniques.

MemoryOptimizer{T}

Advanced memory management with cache optimization and NUMA awareness.

NaNDetectionConfig

Configuration for NaN/Inf detection during simulation.

ParallelIOOptimizer{T}

Advanced parallel I/O optimization with asynchronous writes and data staging.

PerformanceMonitor

Comprehensive performance monitoring for parallel efficiency analysis.

Phi2ConditioningMonitor

Global structure for monitoring φ₂ function conditioning during ERK2 integration. Tracks worst conditioning, series expansion usage, and numerical issues.

RadialDomain

Specification of the radial discretization for the spherical shell.

Radial Grid

The simulation domain is a spherical shell between inner radius ri (e.g., inner core boundary) and outer radius ro (e.g., core-mantle boundary). The radial direction is discretized using Chebyshev collocation for spectral accuracy.

Fields

• N: Number of radial points
• local_range: Indices of radial points local to this MPI process
• r: Radial coordinate values [1, N]
• dr_matrices: Radial derivative matrices (1st, 2nd order, etc.)
• radial_laplacian: Matrix for radial part of Laplacian
• integration_weights: Quadrature weights for radial integration

SHTnsKitConfig <: AbstractSHTnsConfig

Main configuration structure for spherical harmonic transforms using SHTnsKit. This struct encapsulates all parameters needed for transforms and parallelization.

Fields

Core SHTnsKit Configuration

• sht_config: The underlying SHTnsKit.SHTConfig object

Grid Parameters

• nlat::Int: Number of latitude (theta) points (Gauss-Legendre grid)
• nlon::Int: Number of longitude (phi) points (uniform grid)
• lmax::Int: Maximum spherical harmonic degree
• mmax::Int: Maximum spherical harmonic order (≤ lmax)
• nlm::Int: Total number of (l,m) spectral mode pairs

Parallelization Infrastructure

• pencils::NamedTuple: Collection of PencilArrays Pencil objects for different data orientations (:theta, :phi, :r, :spec, :mixed)
• fft_plans::Dict{Symbol,Any}: Precomputed FFTW plans for FFT operations
• transpose_plans::Dict{Symbol,Any}: Plans for data redistribution between pencils

Auxiliary Data

• memory_estimate::String: Human-readable memory usage estimate
• l_values::Vector{Int}: Spherical harmonic degree for each spectral index
• m_values::Vector{Int}: Spherical harmonic order for each spectral index
• theta_grid::Vector{Float64}: Latitude values (Gauss-Legendre nodes)
• phi_grid::Vector{Float64}: Longitude values (uniform spacing)
• gauss_weights::Vector{Float64}: Gauss-Legendre quadrature weights

Internal

• _buffer_cache::Dict{Symbol,Any}: Reusable work arrays to reduce allocations

Usage

config = create_shtnskit_config(lmax=32, mmax=32, nlat=64, nlon=128)

SHTnsPhysicalField{T}

A scalar field represented on the physical (θ, φ, r) grid.

Storage Layout

Field values stored in a 3D array:

• Dimension 1: Latitude index (1 to nlat), Gauss-Legendre points
• Dimension 2: Longitude index (1 to nlon), uniform spacing
• Dimension 3: Radial level index

Pencil Orientation

The pencil field specifies how data is distributed across MPI processes. Different pencil orientations are optimal for different operations:

• phi pencil: Optimal for FFTs (all longitudes local)
• theta pencil: Optimal for Legendre transforms (all latitudes local)
• r pencil: Optimal for radial operations (all radii local)

SHTnsSpectralField{T}

A scalar field represented in spectral space as spherical harmonic coefficients.

Storage Layout

Coefficients are stored as separate real and imaginary parts in 3D arrays:

• Dimension 1: Combined (l,m) spectral index (1 to nlm)
• Dimension 2: Dummy dimension (size 1) for PencilArrays compatibility
• Dimension 3: Radial level index

Boundary Conditions

Each (l,m) mode can have independent boundary conditions at the inner and outer radial boundaries. Common types:

• DIRICHLET: Specified value at boundary
• NEUMANN: Specified derivative at boundary

Fields

• config: SHTnsKit configuration (grid and transform parameters)
• nlm: Total number of spectral modes
• data_real, data_imag: Real/imaginary parts of coefficients
• pencil: PencilArrays pencil defining data distribution
• bc_type_inner/outer: Boundary condition type for each mode
• boundary_values: Boundary values [2, nlm] for inner/outer

SHTnsTorPolField{T}

Vector field in toroidal-poloidal spectral representation.

The Toroidal-Poloidal Decomposition

Any solenoidal (divergence-free) vector field can be written as: v = ∇ × (T r̂) + ∇ × ∇ × (P r̂)

where T and P are scalar functions called the toroidal and poloidal potentials.

Why This Representation?

1. Automatically satisfies ∇·v = 0 (mass conservation for incompressible flow)
2. Decouples certain physical processes in the equations
3. Reduces storage: 2 scalars instead of 3 vector components
4. Natural for spherical geometry and spectral methods

Fields

• toroidal: Toroidal potential T(l,m,r) in spectral space
• poloidal: Poloidal potential P(l,m,r) in spectral space

SHTnsVectorField{T}

A vector field represented by its three spherical components in physical space.

Components:

• r_component: Radial component v_r(θ, φ, r)
• θ_component: Latitudinal component v_θ(θ, φ, r)
• φ_component: Longitudinal component v_φ(θ, φ, r)

Each component is a SHTnsPhysicalField, potentially with different pencil orientations for optimal computation of different operations.

SIMDOptimizer{T}

Advanced SIMD vectorization for mathematical operations with AVX/NEON support.

SimulationState{T}

Unified simulation state with comprehensive parallelization and diagnostics. Type-stable version with concrete cache types for optimal performance.

SolenoidalMonitor

Structure for monitoring solenoidal constraint ∇·v = 0 for velocity and magnetic fields.

SpectralToPhysicalConverter{T}

Structure for converting spectral field data to physical space using the GeoDynamo.jl infrastructure (PencilArrays + SHTns).

TaskGraph

Represents computation as a directed acyclic graph for optimal scheduling.

TaskNode

Represents a computation task in a dependency graph.

ThreadingAccelerator{T} (Backward Compatibility)

Basic CPU threading acceleration for existing code.

_get_influence_cache(domain::RadialDomain, dr_matrix)

Get or create cached influence functions and matrix for given domain. Note: drmatrix should be a BandedMatrix (defined in linearalgebra.jl)

_get_tau_cache(domain::RadialDomain)

Get or create cached tau polynomials and derivatives for given domain.

_shtns_make_transpose(pair)

Create a PencilArrays transpose plan between two pencil configurations. Used internally to set up efficient data redistribution operations.

Arguments

• pair: A Pair of source and destination Pencil objects (src => dest)

Returns

• A Transposition object that can be used with mul! for data redistribution

add_internal_sources_local!(field::AbstractScalarField{T}, domain::RadialDomain) where T

Add volumetric sources (completely local operation). This works for any scalar field with radial source profile.

analyze_load_balance(pencil::Pencil)

Analyze and report load balance for a given pencil decomposition.

apply_banded_full!(out, B, v)

Apply banded matrix to full vector.

apply_composition_boundary_conditions!(field; time_index=nothing)

Refresh composition boundary values from the BoundaryConditions subsystem and enforce Dirichlet data in spectral space.

apply_enhanced_implicit_step!(state, dt)

Apply enhanced implicit time integration with advanced solvers.

apply_exponential_filter!(config::SHTnsKitConfig, alm::Matrix{ComplexF64};
                           order::Int=16, cutoff::Float64=0.65,
                           alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Apply an exponential spectral filter for dealiasing. filter(l) = exp(-α * (l/lmax)^order) where α is chosen so filter(cutoff*lmax) = 0.5

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• order: Filter order (higher = sharper cutoff, default 16)
• cutoff: Fraction of lmax where filter = 0.5 (default 0.65)
• alm_out: Optional output array

apply_flux_bc_direct!(spec_real, spec_imag, local_lm, lm_idx,
                     field::AbstractScalarField, domain, r_range)

Direct modification of boundary values to approximately satisfy flux BC. This is the simplest but least accurate method.

apply_flux_bc_influence_matrix!(spec_real, spec_imag, local_lm, lm_idx,
                               apply_inner, apply_outer, field::AbstractScalarField,
                               domain, r_range)

Apply flux boundary conditions using the influence matrix method.

apply_flux_bc_spectral_complete!(temp_field, domain)

Complete implementation of flux boundary conditions in spectral space. This modifies the spectral coefficients to satisfy ∂T/∂r = prescribed_flux.

apply_flux_bc_tau!(spec_real, spec_imag, local_lm, lm_idx, 
                   apply_inner, apply_outer, field::AbstractScalarField, 
                   domain, r_range)

Apply flux boundary conditions using the tau method. This is the generalized version that works with any scalar field.

apply_geometric_factors_spectral!(field::AbstractScalarField{T}, domain::RadialDomain) where T

Apply geometric factors (1/r, 1/(r sin θ)) in spectral space. For gradients in spherical coordinates. This is generic.

apply_horizontal_laplacian!(config::SHTnsKitConfig, alm::Matrix{ComplexF64};
                             alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Apply the horizontal (angular) Laplacian to spectral coefficients. ∇²h Yl^m = -l(l+1) Y_l^m

This scales each coefficient by -l(l+1), which is the eigenvalue of the horizontal Laplacian on the unit sphere.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• alm_out: Optional output array (if nothing, modifies alm in-place)

Returns

The Laplacian-transformed coefficients

apply_inverse_horizontal_laplacian!(config::SHTnsKitConfig, alm::Matrix{ComplexF64};
                                     alm_out::Union{Matrix{ComplexF64},Nothing}=nothing,
                                     regularize_l0::Bool=true)

Apply the inverse horizontal Laplacian to spectral coefficients. This scales each coefficient by -1/(l(l+1)), which is useful for solving Poisson equations.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• alm_out: Optional output array
• regularize_l0: If true, set l=0 mode to zero (since 1/0 is undefined)

Returns

The inverse-Laplacian-transformed coefficients

apply_magnetic_boundary_conditions!(fields; time_index=nothing)

Refresh magnetic boundary data from the BoundaryConditions subsystem and enforce the corresponding Dirichlet values in spectral space.

apply_scalar_flux_bc_spectral!(field::AbstractScalarField{T}, domain::RadialDomain;
                               method::Symbol=:tau) where T

Apply flux boundary conditions to a scalar field in spectral space. Methods available: :tau (most robust), :influence_matrix, :direct (simplest).

apply_spectral_filter!(config::SHTnsKitConfig, alm::Matrix{ComplexF64},
                       filter_func::Function;
                       alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Apply a custom spectral filter to coefficients. The filter function takes (l, m) and returns a scaling factor.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• filter_func: Function (l, m) -> scale_factor
• alm_out: Optional output array

Example

# Exponential filter for dealiasing
exp_filter(l, m) = exp(-(l/lmax)^16)
apply_spectral_filter!(config, alm, exp_filter)

apply_tau_correction!(profile::Vector{T}, tau_coeffs, domain::RadialDomain) where T

Add tau correction to the radial profile.

apply_temperature_boundary_conditions!(field; time_index=nothing)

Refresh cached boundary values through the BoundaryConditions subsystem and enforce Dirichlet data directly in spectral space when appropriate.

apply_velocity_boundary_conditions!(fields; time_index=nothing)

Refresh velocity boundary data from the BoundaryConditions subsystem and immediately enforce Dirichlet constraints in spectral space.

apply_velocity_component_flux_bc!(field::SHTnsSpectralField{T},
                                   domain::RadialDomain,
                                   dr_matrix::BandedMatrix{T},
                                   method::Symbol) where T

Apply stress-free boundary conditions to a single velocity component (toroidal or poloidal).

All methods now correctly enforce ∂T/∂r = T/r for stress-free boundaries.

Methods

• :tau - High-order accurate tau method (recommended, default)
• :direct - First-order finite difference approximation
• :physical_stress - First-order finite difference (equivalent to :direct)

Performance

Uses pre-allocated workspace buffers if available (set via setvelocityworkspace!). Otherwise allocates temporary arrays (slower).

apply_velocity_flux_bc_direct!(spec_real, spec_imag, local_lm, lm_idx,
                               apply_inner, apply_outer, domain, r_range)

Apply stress-free boundary conditions using direct substitution.

Stress-Free Boundary Condition Physics

For the toroidal velocity potential T in spherical coordinates, the stress-free (free-slip) boundary condition requires zero tangential stress at the boundary.

The tangential stress tensor component is: σrθ ∝ r ∂/∂r(vθ/r) + (1/r)∂v_r/∂θ

For toroidal flow (vr = 0), this reduces to: σrθ ∝ r ∂/∂r(vθ/r) = ∂vθ/∂r - v_θ/r

Setting σrθ = 0 gives: ∂vθ/∂r = v_θ/r

Since v_θ ∝ T for toroidal flow, this becomes: ∂T/∂r = T/r

This is NOT the same as simple Neumann (∂T/∂r = 0)!

Implementation

Using first-order finite difference at boundary: (T[2] - T[1])/Δr = T[1]/r[1]

Solving for T[1]: T[1] = T[2] / (1 + Δr/r[1])

Similarly for outer boundary: (T[nr] - T[nr-1])/Δr = T[nr]/r[nr] T[nr] = T[nr-1] / (1 - Δr/r[nr])

apply_velocity_flux_bc_direct_ws!(spec_real, spec_imag, local_lm, lm_idx,
                                  apply_inner, apply_outer,
                                  domain, r_range, ws, tid)

Workspace-based version of direct stress-free BC (zero allocation). Uses pre-allocated buffers from workspace ws for thread tid.

Stress-Free Boundary Condition Physics

For the toroidal velocity potential T in spherical coordinates, the stress-free (free-slip) boundary condition requires zero tangential stress at the boundary.

The tangential stress tensor component is: σrθ ∝ r ∂/∂r(vθ/r) + (1/r)∂v_r/∂θ

For toroidal flow (vr = 0), this reduces to: σrθ ∝ r ∂/∂r(vθ/r) = ∂vθ/∂r - v_θ/r

Setting σrθ = 0 gives: ∂vθ/∂r = v_θ/r

Since v_θ ∝ T for toroidal flow, this becomes: ∂T/∂r = T/r

This is NOT the same as simple Neumann (∂T/∂r = 0)!

Implementation

Using first-order finite difference at boundary: (T[2] - T[1])/Δr = T[1]/r[1]

Solving for T[1]: T[1] = T[2] / (1 + Δr/r[1])

Similarly for outer boundary: (T[nr] - T[nr-1])/Δr = T[nr]/r[nr] T[nr] = T[nr-1] / (1 - Δr/r[nr])

apply_velocity_flux_bc_physical_stress!(spec_real, spec_imag, local_lm, lm_idx,
                                        apply_inner, apply_outer, domain, r_range)

Apply flux boundary conditions for proper stress-free boundaries. Enforces ∂T/∂r = T/r at boundaries, which corresponds to zero tangential stress.

Physical Justification

For stress-free boundaries: τ = ∂vtan/∂r - vtan/r = 0 In spectral form with v_tan = T(r) × f(θ,φ): ∂(T×f)/∂r - (T×f)/r = 0 (∂T/∂r - T/r) × f = 0 => ∂T/∂r = T/r

This is the CORRECT condition for stress-free boundaries in spherical coordinates.

apply_velocity_flux_bc_physical_stress_ws!(spec_real, spec_imag, local_lm, lm_idx,
                                           apply_inner, apply_outer,
                                           domain, r_range, ws, tid)

Workspace-based version of physical stress BC (zero allocation). Enforces ∂T/∂r = T/r for proper stress-free boundaries.

apply_velocity_flux_bc_spectral!(fields::SHTnsVelocityFields{T}, domain::RadialDomain;
                                 method::Symbol=:tau) where T

Apply stress-free boundary conditions to velocity field components in spectral space.

This function enforces the correct stress-free condition for the toroidal velocity potential. For stress-free boundaries:

• Poloidal (radial) component: P = 0 (Dirichlet, vr = 0, handled by enforcevelocityboundaryvalues!)
• Toroidal (tangential) component: ∂T/∂r = T/r (stress-free condition, handled here)

Arguments

• fields: Velocity field structure with toroidal and poloidal components
• domain: Radial domain information
• method: Method for applying BCs (default :tau recommended for accuracy)
  • :tau - High-order accurate tau method (recommended)
  • :direct - First-order finite difference
  • :physical_stress - First-order finite difference (equivalent to :direct)

Physical Interpretation

For stress-free boundaries, the zero tangential stress condition requires: σrθ = r ∂/∂r(vθ/r) = ∂vθ/∂r - vθ/r = 0

For the toroidal potential T where v_θ ∝ T, this becomes: ∂T/∂r = T/r at boundaries

Note: This is NOT the same as simple Neumann (∂T/∂r = 0)!

apply_velocity_flux_bc_tau!(spec_real, spec_imag, local_lm, lm_idx,
                            apply_inner, apply_outer, dr_matrix, domain, r_range)

Apply stress-free boundary conditions using the tau method for velocity components.

Stress-Free Boundary Condition Physics

For the toroidal velocity potential T in spherical coordinates, the stress-free (free-slip) boundary condition requires zero tangential stress at the boundary.

The correct condition is: ∂T/∂r = T/r

This means the target flux at each boundary is NOT zero, but T/r:

• Inner boundary: target_flux = T[1]/r[1]
• Outer boundary: target_flux = T[nr]/r[nr]

The tau method adds a correction polynomial to enforce this condition exactly.

apply_velocity_flux_bc_tau_ws!(spec_real, spec_imag, local_lm, lm_idx,
                               apply_inner, apply_outer, dr_matrix,
                               domain, r_range, ws, tid)

Workspace-based version of tau method stress-free BC (zero allocation). Uses pre-allocated buffers for profiles, derivatives, and corrections.

Stress-Free Boundary Condition Physics

For the toroidal velocity potential T in spherical coordinates, the stress-free (free-slip) boundary condition requires zero tangential stress at the boundary.

The correct condition is: ∂T/∂r = T/r

This means the target flux at each boundary is NOT zero, but T/r:

• Inner boundary: target_flux = T[1]/r[1]
• Outer boundary: target_flux = T[nr]/r[nr]

The tau method adds a correction polynomial to enforce this condition exactly.

batch_convert_directory(input_dir::String, output_dir::String = "";
                       pattern::String = "combined_global_time_",
                       precision::Type{T} = Float64,
                       compute_diagnostics::Bool = true) where T

Convert all spectral files matching a pattern in a directory.

batch_magnetic_transforms!(mag_fields::SHTnsMagneticFields{T}) where T

Perform batched transforms for better cache efficiency using shtnskit_transforms.jl

batch_shtnskit_transforms!(specs::Vector{SHTnsSpectralField{T}},
                          physs::Vector{SHTnsPhysicalField{T}}) where T

Batch process multiple transforms using SHTnsKit with PencilArrays.

batch_spectral_to_physical!(specs, physs)

Compatibility wrapper that calls batch_shtnskit_transforms! for batched spectral→physical transforms using SHTnsKit with PencilArrays/MPI.

batch_velocity_transforms!(fields::SHTnsVelocityFields{T}) where T

Perform batched transforms for better cache efficiency using shtnskit_transforms.jl

batch_write_spectral_fields!(nc_file, fields::Dict, config::OutputConfig, field_info::FieldInfo)

Write multiple spectral fields in batched operations for better I/O performance

build_banded_A(T, domain, diffusivity, l) -> BandedMatrix{T}

Construct banded A = ν*(d²/dr² + (2/r)d/dr − l(l+1)/r²) in banded storage.

build_influence_matrix(G_inner, G_outer, dr_matrix, domain)

Construct a 2×2 matrix mapping influence amplitudes to boundary flux errors. Rows correspond to (inner, outer) boundaries; columns to (inner, outer) influence functions.

build_rhs_cnab2!(rhs, un, nl, nl_prev, dt, matrices)

Build RHS for CNAB2 IMEX: rhs = un/dt + (1-θ)·L·un + (3/2)·nl − (1/2)·nl_prev, where θ = matrices.theta and L is the diffusivity-scaled linear operator.

build_theta_derivative_matrix(::Type{T}, config::SHTnsKitConfig) where T

Build sparse matrix for θ-derivatives in spectral space. This matrix couples different l modes with the same m.

calculate_linear_neighbors(rank::Int, dim::Int, proc_dims::Tuple, boundaries::Symbol)

Calculate neighbor ranks for non-Cartesian topologies.

check_field_for_nan(field_data::AbstractArray, field_name::String, config::NaNDetectionConfig, step::Int)

Check a field for NaN or Inf values. Returns (hasnan, hasinf, nancount, infcount).

check_simulation_state_for_nan(state::SimulationState, step::Int;
                                 config::NaNDetectionConfig=DEFAULT_NAN_CONFIG)

Comprehensive NaN/Inf check across all simulation fields. Returns true if any NaN/Inf detected.

check_solenoidal_constraint!(state::SimulationState)

Check and record solenoidal constraint for velocity and magnetic fields.

check_spectral_field_for_nan(field::SHTnsSpectralField, field_name::String,
                              config::NaNDetectionConfig, step::Int)

Check both real and imaginary parts of a spectral field.

combine_distributed_time(output_dir::String, time::Float64; 
                        config::CombinerConfig = create_combiner_config(),
                        precision::Type{T} = Float64) where T

Main function to combine distributed files for a specific time.

combine_magnetic_spectral!(combiner::FieldCombiner{T}, lm_mapping::Dict) where T

Combine magnetic spectral coefficients from all distributed files.

combine_physical_fields!(combiner::FieldCombiner{T}) where T

Combine physical space temperature field from distributed files.

combine_spectral_fields!(combiner::FieldCombiner{T}) where T

Combine spectral field data from distributed files.

combine_time_series(output_dir::String, time_range::Tuple{Float64,Float64}; 
                   config::CombinerConfig = create_combiner_config(),
                   precision::Type{T} = Float64) where T

Combine multiple time points into a time series.

combine_to_global!(combiner::FieldCombiner{T}) where T

Combine distributed field data into global field structures.

combine_velocity_spectral!(combiner::FieldCombiner{T}, lm_mapping::Dict) where T

Combine velocity spectral coefficients from all distributed files.

compute_all_gradients_spectral!(field::AbstractScalarField{T}, domain::RadialDomain) where T

Compute all gradient components (θ, φ, r) in spectral space. This is the main driver function that works for any scalar field.

compute_all_nonlinear_terms!(fields, temp_field, comp_field, mag_field, domain)

Compute all nonlinear forcing terms for the momentum equation.

Governing Equation (Dimensional)

In a rotating reference frame with angular velocity Ω:

∂u/∂t = -u×ω - ∇p/ρ + 2Ω(ẑ×u) + ν∇²u + (αg/ρ)T·r̂ + (Δρg/ρ)C·r̂ + (1/μ₀ρ)(∇×B)×B

where:

• u: velocity, ω = ∇×u: vorticity
• Ω: rotation rate, ẑ: rotation axis
• ν: kinematic viscosity
• α: thermal expansion coefficient, g: gravity
• T: temperature perturbation, C: composition perturbation
• B: magnetic field, μ₀: permeability

Non-Dimensionalization (Magnetic Diffusion Time Scaling)

Length scale: d (shell thickness or ball radius) Time scale: τ = d²/η (magnetic diffusion time) Velocity scale: U = η/d Magnetic field scale: B₀ Temperature scale: ΔT

Dimensionless parameters:

• E = ν/(Ω d²): Ekman number
• Pm = ν/η: Magnetic Prandtl number
• Pr = ν/κ: Prandtl number
• Sc = ν/D: Schmidt number
• Ra = (αgΔT d³)/(νκ): Rayleigh number
• Ra_C = (Δρg d³)/(ρνD): Compositional Rayleigh number

Non-Dimensional Momentum Equation

(E/Pm) ∂ũ/∂τ + (E/Pm)(∇×ũ)×ũ + ẑ×ũ = -∇p̃ + (Pm/Pr)Ra·T̃·r̂ + (Pm/Sc)Ra_C·C̃·r̂ + (∇×B̃)×B̃ + E∇²ũ

where tilde denotes dimensionless quantities.

Implementation Notes

After dividing Eq. (1) by E/Pm, the explicit RHS entering the time integrator is: RHS = ũ×ω̃ - (Pm/E)(ẑ×ũ) + (Pm/E)(Pm/Pr)Ra·T̃·r̂ + (Pm/E)(Pm/Sc)Ra_C·C̃·r̂ + (Pm/E)(∇×B̃)×B̃

• The advection term ũ×ω̃ already matches −(∇×ũ)×ũ without extra scaling.
• Coriolis, buoyancy, and Lorentz forces carry the (Pm/E) prefactor (=rossby_factor).
• Viscous diffusion is treated implicitly with coefficient Pm (passed via diffusivity).

The time derivative retains the (E/Pm) mass term and is handled by the integrator.

compute_boundary_fluxes(profile::Vector{T}, dr_matrix, domain::RadialDomain) where T

Compute flux (dT/dr) at both boundaries using the derivative matrix. Note: drmatrix should be a BandedMatrix{T} (defined in linearalgebra.jl)

compute_chebyshev_polynomial(n::Int, domain::RadialDomain)

Compute Chebyshev polynomial T_n on the radial grid.

compute_divergence_spectral(tor_spec::SHTnsSpectralField, pol_spec::SHTnsSpectralField,
                            domain::RadialDomain) where T

Compute divergence of a vector field in spectral space. For toroidal-poloidal decomposition: v = ∇×(T r̂) + ∇×∇×(P r̂) The divergence is: ∇·v = ∇·(∇×∇×(P r̂)) = 0 theoretically But numerically we check: div_max = max |∇·v|

compute_enstrophy(config::SHTnsKitConfig, Tlm::Matrix{ComplexF64}; real_field::Bool=true)

Compute enstrophy (mean square vorticity) from toroidal coefficients. Enstrophy is related to the rotational part of the kinetic energy.

compute_field_energy(field_data::Array{T,3}) where T

Compute L2 energy of a 3D field: E = ½ ∫ |f|² dV ≈ ½ Σ |f_ijk|²

compute_global_diagnostics(combiner::FieldCombiner{T}) where T

Compute global diagnostics for the combined fields.

compute_global_diagnostics(converter::SpectralToPhysicalConverter{T}) where T

Compute global field diagnostics using MPI reductions.

compute_horizontal_gradient_magnitude(config::SHTnsKitConfig, alm::Matrix{ComplexF64})

Compute the magnitude of the horizontal gradient |∇h f|² in spectral space. |∇h f|² = l(l+1) |f_lm|² (summed over all modes)

This is useful for computing gradient energy or penalty terms.

Arguments

• config: SHTnsKit configuration
• alm: Spectral coefficients

Returns

Scalar value of the integrated horizontal gradient magnitude squared

compute_influence_functions_flux(oc_domain::RadialDomain)

Compute influence functions for flux BCs. These are solutions to the homogeneous equation with specific BCs.

compute_inner_tau_function(domain::RadialDomain)

Compute tau function for inner boundary only (zero derivative at outer).

compute_magnetic_helicity(mag_fields::SHTnsMagneticFields{T}) where T

Compute magnetic helicity using enhanced spectral integration

compute_nonlinear!(parallelizer, temp_field, vel_fields, domain)

Advanced CPU computation with SIMD, NUMA, and task-based parallelism.

compute_outer_tau_function(domain::RadialDomain)

Compute tau function for outer boundary only (zero derivative at inner).

compute_phi1_function(A, expA)

Compute φ1(A) = (exp(A) - I) / A efficiently with comprehensive error handling. Uses series expansion for small ||A|| to avoid numerical issues.

compute_phi2_function(A, expA; l=0)

Compute φ2(A) = (exp(A) - I - A) / A² efficiently with comprehensive error handling. Uses series expansion for small ||A|| to avoid numerical issues. Tracks conditioning statistics when monitoring is enabled.

compute_phi_gradient_spectral!(field::AbstractScalarField{T}) where T

Compute ∂field/∂φ using spherical harmonic properties (local operation). This is generic and works for any scalar field.

compute_radial_gradient_spectral!(field::AbstractScalarField{T}, domain::RadialDomain) where T

Compute ∂field/∂r using banded matrix derivative operator in spectral space (local operation). This uses the pre-computed derivative matrices from the field for optimal accuracy and efficiency.

compute_scalar_advection_local!(field::AbstractScalarField{T}, vel_fields) where T

Compute -u·∇field in physical space (completely local operation). This works for any scalar field.

compute_scalar_energy(field::AbstractScalarField{T}, oc_domain::RadialDomain) where T

Compute energy ∫ field² dV

compute_scalar_energy_spectrum(config::SHTnsKitConfig, alm::Matrix{ComplexF64}; real_field::Bool=true)

Compute the energy spectrum per spherical harmonic degree l using SHTnsKit v1.1.15.

Returns

Vector of length lmax+1 with energy at each degree l.

compute_scalar_rms(field::AbstractScalarField{T}, oc_domain::RadialDomain) where T

Compute RMS value of scalar field.

compute_spectral_energy_diagnostics!(diagnostics, component, real_part, imag_part, field_info)

Compute spectral energy distribution diagnostics using SHTns configuration

compute_tangential_stress_components(v_theta, v_phi, dv_theta_dr, dv_phi_dr, r_val)

Compute tangential stress components from velocity and derivatives.

For incompressible flow with vr = 0: τrθ = μ(∂vθ/∂r - vθ/r) τrφ = μ(∂vφ/∂r - v_φ/r)

Arguments

• v_theta, v_phi: Tangential velocity components [nlat, nlon]
• dv_theta_dr, dv_phi_dr: Radial derivatives [nlat, nlon]
• r_val: Radial position

Returns

• tau_theta, tau_phi: Stress components [nlat, nlon] with μ = 1
• max_stress: Maximum stress magnitude

compute_tau_coefficients_both(flux_error_inner::T, flux_error_outer::T, domain::RadialDomain) where T

Compute tau polynomial coefficients for both boundaries. Uses highest two Chebyshev modes as tau functions.

compute_tau_coefficients_inner(flux_error::T, domain::RadialDomain) where T

Compute tau coefficient for inner boundary only. Uses a single tau function that doesn't affect outer boundary.

compute_tau_coefficients_outer(flux_error::T, domain::RadialDomain) where T

Compute tau coefficient for outer boundary only. Uses a single tau function that doesn't affect inner boundary.

compute_tau_correction_both_boundaries!(correction::Vector{T}, ...)

In-place version that writes to pre-allocated correction buffer.

compute_tau_correction_inner_boundary!(correction::Vector{T}, ...)

In-place version that writes to pre-allocated correction buffer.

compute_tau_correction_outer_boundary!(correction::Vector{T}, ...)

In-place version that writes to pre-allocated correction buffer.

compute_theta_gradient_spectral!(field::AbstractScalarField{T}) where T

Compute ∂field/∂θ using spherical harmonic recurrence relations (local operation). This is generic and works for any scalar field.

compute_theta_recurrence_coefficients(::Type{T}, config::SHTnsKitConfig) where T

Pre-compute recurrence coefficients for θ-derivatives. Store as [nlm, 2] matrix: [:, 1] for l-1 coupling, [:, 2] for l+1 coupling

compute_total_energy!(state::SimulationState)

Compute and record total energy of the simulation state.

compute_total_scalar_energy(config::SHTnsKitConfig, alm::Matrix{ComplexF64}; real_field::Bool=true)

Compute total energy of a scalar field from its spectral coefficients.

compute_total_vector_energy(config::SHTnsKitConfig, Slm::Matrix{ComplexF64}, Tlm::Matrix{ComplexF64}; real_field::Bool=true)

Compute total kinetic energy of a vector field from spheroidal/toroidal coefficients.

compute_vector_energy(v_r::Array{T,3}, v_theta::Array{T,3}, v_phi::Array{T,3}) where T

Compute kinetic/magnetic energy: E = ½ ∫ |v|² dV = ½ ∫ (vr² + vθ² + v_φ²) dV

compute_vector_energy_spectrum(config::SHTnsKitConfig, Slm::Matrix{ComplexF64}, Tlm::Matrix{ComplexF64}; real_field::Bool=true)

Compute the kinetic energy spectrum per spherical harmonic degree l for a vector field decomposed into spheroidal (Slm) and toroidal (Tlm) components.

Returns

Vector of length lmax+1 with kinetic energy at each degree l.

convert_spectral_file(input_filename::String, output_filename::String = "";
                     precision::Type{T} = Float64, 
                     compute_diagnostics::Bool = true) where T

Main function to convert a single spectral data file to physical space.

convert_to_physical!(converter::SpectralToPhysicalConverter{T}) where T

Convert all loaded spectral data to physical space using SHTns transforms.

create_boundary_slice(size_arr::Tuple, dim::Int, side::Symbol, width::Int)

Create array slice for boundary data extraction.

create_combiner_config(; kwargs...)

Create combiner configuration with keyword arguments.

create_computation_pencils(topology, dims, config)

Create specialized pencils for different stages of computation.

create_enhanced_field_variables!(nc_file, field_info::FieldInfo, config::OutputConfig, 
                                available_fields::Vector{String})

Enhanced field variable setup that leverages SHTns configuration

create_enhanced_output_config(state)

Create output configuration with enhanced I/O settings.

create_erk2_cache(config, domain, diffusivity, dt; use_krylov=false, m=20, tol=1e-8)

Create ERK2 cache with precomputed matrix functions for all spherical harmonic modes.

create_erk2_config(; lmax, mmax, nlat, nlon, optimize_for_erk2=true)

Create an SHTnsKit configuration for ERK2 timestepping.

create_etd_cache(config, domain, diffusivity, dt) -> ETDCache

Build per-l exponential cache for the linear operator Al = diffusivity*(d²/dr² + (2/r)d/dr − l(l+1)/r²). Computes exp(dt Al) and phi1(dt A_l) via dense methods. Single-rank recommended.

create_field_combiner(output_dir::String, time::Float64; 
                     precision::Type{T} = Float64) where T

Create a field combiner by reading configuration from distributed files.

create_halo_slice(size_arr::Tuple, dim::Int, side::Symbol, width::Int)

Create array slice for halo region insertion.

create_memory_efficient_output_buffer(data_size::Int, precision::DataType)

Create appropriately sized buffer for output operations

create_parameter_template(filename::String)

Create a template parameter file for users to customize.

create_pencil_array(::Type{T}, pencil::Pencil; init=:zero) where T

Create a PencilArray with specified initialization.

create_pencil_decomposition_shtnskit(nlat, nlon, nr, sht_config, comm, optimize)

Create PencilArrays decomposition optimized for spherical harmonic transforms.

The Pencil Decomposition Strategy

In a pencil decomposition, 3D data is distributed across MPI processes such that one dimension is always fully local (not split across processes). This enables efficient operations along that dimension without MPI communication.

Physical Space Pencils (nlat × nlon × nr)

theta pencil:  [θ local] × [φ distributed] × [r distributed]
               → Best for operations along latitude (Legendre transforms)

phi pencil:    [θ distributed] × [φ local] × [r distributed]
               → Best for FFTs along longitude

r pencil:      [θ distributed] × [φ distributed] × [r local]
               → Best for radial operations (derivatives, boundary conditions)

Spectral Space Pencil (nlm × 1 × nr)

The spectral pencil stores spherical harmonic coefficients indexed by a combined (l,m) index. The middle dimension is 1 (dummy) for compatibility.

Arguments

• nlat, nlon, nr: Grid dimensions
• sht_config: SHTnsKit configuration (provides nlm)
• comm: MPI communicator
• optimize: Whether to optimize process topology for load balance

Returns

NamedTuple with pencil configurations: (:theta, :θ, :phi, :φ, :r, :spec, :mixed)

create_pencil_fft_plans(pencils, dims) -> Dict{Symbol,Any}

Create precomputed FFTW plans for efficient FFT operations.

Why Precomputed Plans?

FFTW achieves best performance when FFT plans are created once and reused. Plan creation involves:

1. Analyzing the input size and memory layout
2. Choosing optimal algorithm (Cooley-Tukey, Bluestein, etc.)
3. Possibly benchmarking different strategies

Plan Types Created

• :phi_forward / :phi_backward: FFT/IFFT along longitude (dimension 2)
• :theta_forward / :theta_backward: FFT/IFFT for theta pencil orientation

Technical Notes

• Plans operate on the parent() array of PencilArrays (the underlying Julia array)
• We use FFTW directly rather than PencilFFTPlan because we need single-dimension transforms, not full multi-dimensional distributed FFTs
• The plans are tied to specific array sizes and memory layouts

Arguments

• pencils: NamedTuple of Pencil configurations
• dims: Tuple (nlat, nlon, nr) of grid dimensions

Returns

Dict mapping plan names to FFTW plan objects. Contains :fallback => true on error.

create_pencil_topology(shtns_config; optimize=true)

Create enhanced pencil decomposition for SHTns grids. Supports both 1D and 2D decompositions with automatic optimization. Accepts an object with fields nlat, nlon, and nlm (e.g., SHTnsKitConfig).

create_shtns_aware_output_config(shtns_config::SHTnsKitConfig, pencils::NamedTuple; 
                                base_config::OutputConfig = default_config())

Create output configuration that integrates with SHTns configuration and pencil decomposition

create_shtns_physical_field(T, config, oc_domain, pencil) -> SHTnsPhysicalField{T}

Create a new physical space field initialized to zero.

Arguments

• T: Element type (typically Float64)
• config: SHTnsKit configuration providing grid dimensions
• oc_domain: RadialDomain (for consistency with spectral field API)
• pencil: PencilArrays Pencil defining the data distribution

Returns

A new SHTnsPhysicalField with all grid values initialized to zero.

create_shtns_spectral_field(T, config, oc_domain, pencil_spec) -> SHTnsSpectralField{T}

Create a new spectral field initialized to zero with default Dirichlet BCs.

Arguments

• T: Element type (typically Float64)
• config: SHTnsKit configuration providing nlm and other parameters
• oc_domain: RadialDomain specifying the radial discretization
• pencil_spec: PencilArrays Pencil defining the data distribution

Returns

A new SHTnsSpectralField with:

• All spectral coefficients initialized to zero
• Dirichlet boundary conditions on all modes
• Zero boundary values

create_shtns_vector_field(T, config, oc_domain, pencils) -> SHTnsVectorField{T}

Create a new vector field with three physical space components.

Arguments

• T: Element type
• config: SHTnsKit configuration
• oc_domain: RadialDomain specification
• pencils: Either a NamedTuple with :theta/:θ, :phi/:φ, :r keys, or a tuple (pencilθ, pencilφ, pencil_r)

Returns

A new SHTnsVectorField with all components initialized to zero. Each component uses a potentially different pencil orientation for optimal computation of different operations.

create_shtnskit_config(; lmax, mmax, nlat, nlon, nr, optimize_decomp) -> SHTnsKitConfig

Create and initialize a complete SHTnsKit configuration for spherical harmonic transforms with MPI parallelization.

Keyword Arguments

• lmax::Int: Maximum spherical harmonic degree (required)
• mmax::Int=lmax: Maximum spherical harmonic order (≤ lmax)
• nlat::Int: Number of latitude points. Must be ≥ lmax+1 for numerical accuracy. Defaults to max(lmax+2, parameter system value)
• nlon::Int: Number of longitude points. Must be ≥ 2*mmax+1 for alias-free transforms. Powers of 2 are preferred for FFT efficiency.
• nr::Int=i_N: Number of radial points (from parameter system)
• optimize_decomp::Bool=true: Whether to optimize MPI process topology

Returns

• SHTnsKitConfig: Fully initialized configuration ready for transforms

Algorithm

1. Create base SHTnsKit configuration with Gauss-Legendre grid
2. Set up pencil decomposition for MPI parallelization
3. Create FFT plans for longitude transforms
4. Create transpose plans for data redistribution
5. Initialize grid coordinates and quadrature weights

Example

# Create config for lmax=63 simulation
config = create_shtnskit_config(lmax=63, nlat=96, nlon=192)

create_shtnskit_transpose_plans(pencils) -> Dict{Symbol,Any}

Create transpose plans for redistributing data between pencil orientations.

What is a Transpose?

In pencil decomposition, a "transpose" is an MPI all-to-all communication that redistributes data so a different dimension becomes local. For example:

• theta → phi transpose: makes longitude local (for FFTs)
• phi → r transpose: makes radius local (for radial derivatives)

Why Plans?

Like FFT plans, transpose plans encode the communication pattern and can optimize buffer allocation and MPI operations.

Transpose Constraints

PencilArrays requires that source and destination pencils differ in at most one distributed dimension. If they differ in two dimensions, we need a multi-step transpose (e.g., phi → theta → r).

Plan Keys Created

• :theta_to_phi, :phi_to_theta: Switch between theta and phi local
• :theta_to_r, :r_to_theta: Switch between theta and r local
• :phi_to_r, :r_to_phi: May be direct or marked as :multi_step

Usage

# Transpose data from theta-local to phi-local orientation
mul!(dest_array, transpose_plans[:theta_to_phi], src_array)

create_spectral_converter(filename::String; precision::Type{T} = Float64) where T

Create a spectral to physical converter by reading configuration from a NetCDF file.

create_transpose_plans(pencils)

Create enhanced transpose plans between pencil orientations. Includes caching and communication optimization.

eab2_update!(u, nl, nl_prev, etd, config)

Apply EAB2 update per (l,m): u^{n+1} = E u^n + dtphi1(3/2 nl^n − 1/2 nl^{n−1}).

eab2_update_krylov!(u, nl, nl_prev, domain, diffusivity, config, dt; m=20, tol=1e-8)

EAB2 update using Krylov exp/φ1 actions and banded LU for φ1.

eab2_update_krylov_cached!(u, nl, nl_prev, alu_map, domain, ν, config, dt; m=20, tol=1e-8)

Same as eab2updatekrylov!, but reuses cached banded A and LU per l.

enforce_composition_boundary_values!(field)

Anchor spectral boundary modes to stored Dirichlet values when applicable.

enforce_magnetic_boundary_values!(fields)

Anchor magnetic toroidal/poloidal spectral data to cached Dirichlet boundary values on the inner and outer radial surfaces.

enforce_temperature_boundary_values!(field)

Anchor spectral boundary coefficients to stored Dirichlet values when needed.

enforce_velocity_boundary_values!(fields)

Anchor toroidal and poloidal spectral coefficients to the currently cached Dirichlet boundary values on the inner and outer radial surfaces.

erk2_stage_residual_stats(buffers) -> NamedTuple

Compute diagnostic statistics for the difference between stage nonlinear terms and the base-step nonlinear terms.

estimate_field_count() -> Int

Estimate the number of field arrays typically allocated simultaneously. Used for memory usage estimation. Returns 6 as a reasonable default covering velocity (3 components), temperature, composition, and magnetic field.

estimate_memory_usage(pencils, field_count::Int, precision::Type)

Estimate memory usage for given pencil configuration.

estimate_memory_usage_shtnskit(nlat, nlon, lmax, field_count, T)

Estimate memory usage for SHTnsKit-based transforms with PencilArrays.

eval_with_context(expr, param_dict::Dict{Symbol, Any})

Evaluate an expression by substituting symbols with values from param_dict.

evaluate_chebyshev_derivative(n::Int, r::T, domain::RadialDomain) where T

Evaluate derivative of Chebyshev polynomial at a point.

evaluate_tau_derivative_inner(tau::Vector, domain::RadialDomain)

Evaluate tau function derivative at inner boundary.

evaluate_tau_derivative_outer(tau::Vector, domain::RadialDomain)

Evaluate tau function derivative at outer boundary.

exchange_dimension_halos!(data::Array, pencil::Pencil, dim::Int, 
                          halo_width::Int, boundaries::Symbol, comm::MPI.Comm)

Perform halo exchange along a specific dimension using MPI point-to-point communication.

exp_action_krylov(Aop!, v, dt; m=20, tol=1e-8) -> y ≈ exp(dt A) v

Simple Arnoldi-based approximation of the exponential action.

extract_coefficients_for_shtnskit!(coeffs_buffer, spec_real, spec_imag, r_local, config)

Extract spectral coefficients into SHTnsKit's expected matrix format.

Data Format Conversion

Our internal format: Linear array indexed by combined (l,m) index SHTnsKit format: Matrix indexed by [l+1, m+1]

Arguments

• coeffs_buffer: Pre-allocated (lmax+1) × (mmax+1) complex matrix (output)
• spec_real, spec_imag: Real and imaginary parts of spectral coefficients
• r_local: Radial level index (1-based)
• config: SHTnsKit configuration

Threading

Uses @threads for parallel filling of the coefficient matrix.

extract_coefficients_for_shtnskit(spec_real, spec_imag, r_local, config) -> Matrix{ComplexF64}

High-level coefficient extraction with automatic buffer management and MPI gathering.

This is the main entry point for preparing spectral coefficients for SHTnsKit. It handles:

1. Buffer allocation/reuse from cache
2. Local coefficient extraction
3. MPI Allreduce to combine coefficients from all processes

Why MPI Gathering?

Spectral coefficients may be distributed across MPI processes. SHTnsKit needs the complete coefficient matrix, so we sum partial contributions from all processes.

Returns

A complete (lmax+1) × (mmax+1) coefficient matrix ready for SHTnsKit.synthesis()

extract_divergence_coefficients(config::SHTnsKitConfig, Slm::Matrix{ComplexF64})

Extract divergence spectral coefficients from spheroidal potential. The divergence field in spectral space is related to Slm by the horizontal Laplacian.

Returns

Matrix of divergence coefficients.

extract_local_radial_profile!(profile, data, local_lm, nr, r_range)

In-place version to avoid allocations; writes the local radial line into profile for the given local_lm using the provided r_range.

extract_physical_slice_generic!(slice_buffer, phys_data, r_local, config)

Generic extraction for any pencil orientation using pre-allocated buffer.

extract_physical_slice_phi_local!(slice_buffer, phys_data, r_local, config)

Extract physical slice when in phi-local pencil using pre-allocated buffer.

extract_spectral_subset(field::SHTnsSpectralField{T}; l_max::Int = 10) where T

Extract a subset of spectral modes for analysis.

extract_vector_component_generic!(component_buffer, v_data, r_local, config)

Generic extraction for vector components using pre-allocated buffer.

extract_vorticity_coefficients(config::SHTnsKitConfig, Tlm::Matrix{ComplexF64})

Extract vorticity spectral coefficients from toroidal potential. The vorticity field in spectral space is related to Tlm by the horizontal Laplacian.

Returns

Matrix of vorticity coefficients.

find_distributed_files(output_dir::String, time::Float64, filename_prefix::String = "geodynamo")

Find all distributed files for a specific simulation time.

find_package_root()

Find the root directory of the GeoDynamo.jl package.

get_comm()

Get MPI communicator, initializing MPI if needed. Provides thread-safe lazy initialization.

get_default_nlat() -> Int

Get the default number of latitude points for the grid. Tries to read from the parameter system first, falls back to 64. Used during precompilation when parameters may not be available.

get_default_nlon() -> Int

Get the default number of longitude points for the grid. Tries to read from the parameter system first, falls back to 128. Power of 2 is preferred for efficient FFT operations.

get_dimension_neighbors(pencil::Pencil, dim::Int, boundaries::Symbol) -> (Int, Int)

Get left and right neighbor process ranks for a given dimension. Returns (leftneighbor, rightneighbor) where MPI.MPIPROCNULL indicates no neighbor.

get_eab2_alu_cache!(caches, key, ν, T, domain) -> Dict{Int,Tuple{BandedMatrix{T},BandedLU{T}}}

Retrieve or initialize a cache mapping l -> (Abanded, LU(Abanded)) for EAB2. Reinitializes if ν or nr changed. Type-stable version using EAB2ALUCacheEntry.

get_erk2_cache!(caches, key, diffusivity, config, domain, dt; use_krylov=false)

Retrieve or create ERK2 cache with automatic invalidation when parameters change. Type-stable version using concrete ERK2Cache{T} type.

get_flux_value(lm_idx::Int, boundary::Int, field::AbstractScalarField)

Get prescribed flux value for given mode and boundary from scalar field. This is a generic version that works with any scalar field.

get_nprocs()

Get total number of MPI processes.

get_parameters()

Get the current global parameters. If not set, loads default parameters. Thread-safe lazy initialization with type-stable return.

get_pencil_orientation(pencil) -> Symbol

Determine which dimension(s) are fully local in a pencil decomposition.

Returns

• :theta_phi: Both angular dimensions local (serial or single-node)
• :theta: Latitude (θ) dimension is local
• :phi: Longitude (φ) dimension is local
• :r: Radial dimension is local

Usage

Used to choose optimal transform strategies based on data layout.

get_process_dimensions(topology) -> Tuple

Extract process grid dimensions from topology.

get_rank()

Get MPI rank of current process.

get_shtnskit_performance_stats()

Get performance statistics for SHTnsKit transforms with PencilArrays. Returns information about the v1.1.15+ features being used.

get_shtnskit_version_info()

Get version and capability information about the SHTnsKit installation.

get_velocity_workspace(::Type{T})::Union{VelocityWorkspace{T}, Nothing} where T

Get the global velocity workspace if available and matches type T. Returns nothing if not set or type mismatch.

hybrid_compute_nonlinear!(parallelizer, temp_field, vel_fields, domain)

Compute nonlinear terms using hybrid MPI+threading parallelism.

index_to_lm_shtnskit(idx, lmax, mmax) -> (l, m)

Convert linear spectral index to spherical harmonic degree (l) and order (m).

Index Ordering

The linear index follows the convention where m varies fastest within each l:

• idx=1: (l=0, m=0)
• idx=2: (l=1, m=0)
• idx=3: (l=1, m=1)
• idx=4: (l=2, m=0)
• idx=5: (l=2, m=1)
• idx=6: (l=2, m=2)
• ...

Performance Note

This function uses a linear search. For performance-critical code with many lookups, consider precomputing the lvalues and mvalues arrays (stored in SHTnsKitConfig).

infer_pencils_from_transpose_name(name::Symbol, plan) -> (src_pencil, dest_pencil)

Infer source and destination pencils from transpose operation name and plan. For TransposeOperator, we need to extract pencil info from context since the operator itself doesn't expose src/dest pencils directly.

initialize_enhanced_simulation(::Type{T}=Float64; kwargs...)

Initialize simulation with all parallelization features enabled.

initialize_parameters(config_file::String = "")

Initialize the global parameter system.

initialize_master_simulation(::Type{T}=Float64; kwargs...)

Initialize simulation with comprehensive CPU parallelization.

install_erk2_cache_bundle!(target, bundle)

Copy ERK2 caches from bundle into the target cache dictionary.

list_available_times(output_dir::String, filename_prefix::String = "geodynamo")

List all available simulation times in the output directory.

load_distributed_fields!(combiner::FieldCombiner{T}) where T

Load field data from all distributed files into the combiner structure.

load_erk2_cache_bundle!(target, path) -> metadata

Load caches from path and install them into target, returning metadata.

load_erk2_cache_bundle(path) -> (caches, metadata)

Load ERK2 caches and associated metadata from disk.

load_parameters(config_file::String = "")

Load parameters from a configuration file. If no file is specified, loads from the default config/default_params.jl file.

Arguments

• config_file::String: Path to parameter file (optional)

Returns

• GeoDynamoParameters: Loaded parameters

load_parameters_from_file(config_file::String)

Load parameters from a Julia file containing parameter definitions.

load_physical_temperature!(field::SHTnsPhysicalField{T}, nc_file, var_name::String) where T

Load physical space temperature data from NetCDF file.

load_single_file!(combiner::FieldCombiner{T}, filename::String) where T

Load data from a single file (no combination needed).

load_single_spectral_field!(field::SHTnsSpectralField{T}, nc_file, 
                           real_var::String, imag_var::String) where T

Load a single spectral field from NetCDF file.

load_spectral_coefficients!(field::SHTnsSpectralField{T}, nc_file, 
                            real_var_name::String, imag_var_name::String,
                            l_values::Vector{Int}, m_values::Vector{Int}, 
                            r_spectral::Vector{Float64}) where T

Load spectral coefficients from NetCDF file into a spectral field structure.

load_spectral_data!(converter::SpectralToPhysicalConverter{T}, filename::String) where T

Load spectral field data from NetCDF file into GeoDynamo field structures.

main_batch_convert(input_dir::String; kwargs...)

Main entry point for batch conversion of files.

main_combine_time(output_dir::String, time::Float64; kwargs...)

Main entry point for combining distributed files for a single time.

main_combine_time_series(output_dir::String, time_range::Tuple{Float64,Float64}; kwargs...)

Main entry point for combining time series.

main_convert_file(filename::String; kwargs...)

Main entry point for converting a single file.

map_spectral_coefficients!(global_tor_real, global_tor_imag, global_pol_real, global_pol_imag,
                          local_tor_real, local_tor_imag, local_pol_real, local_pol_imag,
                          local_l, local_m, lm_mapping)

Map local spectral coefficients to global arrays using (l,m) index mapping.

maybe_log_erk2_stage_residual!(label, buffers, step)

Emit a diagnostic log entry when ERK2 diagnostics are enabled.

modify_for_flux_inner!(spec_real, spec_imag, local_lm, prescribed_flux,
                      dr_matrix, domain, r_range)

Modify coefficients near inner boundary to approximate flux condition. Uses low-order extrapolation.

modify_for_flux_outer!(spec_real, spec_imag, local_lm, prescribed_flux,
                      dr_matrix, domain, r_range)

Modify coefficients near outer boundary to approximate flux condition. Uses low-order extrapolation.

normalize_influence_function!(G::Vector{T}, domain::RadialDomain, which_boundary::Int) where T

Normalize influence function to have unit flux at the specified boundary.

optimize_communication_order(plans::Dict)

Determine optimal order for transpose operations to minimize communication.

optimize_erk2_transforms!(config::SHTnsKitConfig)

Optimize SHTnsKit transforms for ERK2 timestepping with PencilFFTs. This function pre-warms transform plans and optimizes memory layout.

optimize_fft_performance!(config::SHTnsKitConfig)

Optimize FFT performance by warming up FFTW plans and checking efficiency.

optimize_field_data_layout!(field_data::Dict, field_info::FieldInfo)

Optimize data layout for output based on pencil decomposition

optimize_magnetic_memory_layout!(mag_fields::SHTnsMagneticFields{T}) where T

Optimize memory layout for better cache performance using pencil topology

optimize_process_topology(nprocs::Int, dims::Tuple{Int,Int,Int})

Find optimal 2D process grid for given number of processes and problem dimensions. Minimizes communication volume.

optimize_process_topology_shtnskit(nprocs, nlat, nlon) -> Tuple{Int,Int}

Find optimal 2D MPI process grid for theta-phi parallelization.

The Optimization Problem

Given nprocs MPI processes, we need to factor it as nprocs = p_theta × p_phi such that:

1. Each process gets at least 2 grid points in each direction
2. Load imbalance (due to non-divisibility) is minimized
3. The decomposition is valid (exact factorization)

Algorithm

Iterates through all valid factorizations and scores them by load imbalance. The load imbalance for a dimension is: |gridsize mod processes| / gridsize

Example

For nprocs=12, nlat=64, nlon=128:

• Possible factorizations: (1,12), (2,6), (3,4), (4,3), (6,2), (12,1)
• (4,3) gives nlat/4=16 points/proc in θ, nlon/3≈43 points/proc in φ
• This might be optimal if it minimizes total imbalance

Arguments

• nprocs: Total number of MPI processes
• nlat: Number of latitude points
• nlon: Number of longitude points

Returns

• (p_theta, p_phi): Optimal process grid dimensions

optimize_velocity_memory_layout!(fields::SHTnsVelocityFields{T}) where T

Optimize memory layout for better cache performance using pencil topology

perform_analysis_direct!(phys, spec, config)

Direct analysis without transpose (fallback).

perform_analysis_from_phi_pencil!(phys_phi, spec, config)

Perform analysis from phi-pencil data.

perform_analysis_phi_local!(phys, spec, config)

Perform analysis when physical field is in phi-pencil orientation.

This is the most efficient analysis path because:

1. The phi (longitude) dimension is fully local on each process
2. SHTnsKit's FFT operates entirely in local memory
3. No MPI communication needed during the transform itself

perform_analysis_with_transpose!(phys, spec, config, to_phi_plan)

Perform analysis with transpose to phi-pencil.

perform_synthesis_direct!(spec, phys, config)

Direct synthesis method that handles any pencil orientation.

This is the default/fallback method. It works regardless of the physical field's pencil orientation by using a generic storage function that handles the index mapping appropriately.

Note

For phi-pencil physical fields, perform_synthesis_phi_local! is more efficient as it can use optimized storage.

perform_synthesis_phi_local!(spec, phys, config)

Perform synthesis when physical field is in phi-pencil orientation.

This is the most efficient synthesis path because:

1. The phi (longitude) dimension is fully local on each process
2. SHTnsKit's FFT operates entirely in local memory
3. No MPI communication needed during the transform itself

Algorithm for each radial level:

1. Extract spectral coefficients into SHTnsKit's (lmax+1, mmax+1) matrix format
2. Call SHTnsKit.synthesis() which does Legendre transform + FFT
3. Store the resulting (nlat, nlon) physical field slice

perform_synthesis_to_phi_pencil!(spec, phys_phi, config)

Core synthesis routine that writes directly to a phi-pencil array.

This is the workhorse function called by other synthesis routines. It assumes the destination array is already in phi-pencil orientation.

perform_synthesis_with_transpose!(spec, phys, config, back_plan)

Perform synthesis when physical field is NOT in phi-pencil orientation.

Strategy

When the target physical field is in a non-phi pencil (e.g., theta or r pencil), we can't directly use SHTnsKit because it requires all longitude points to be local.

Solution:

1. Create a temporary phi-pencil array
2. Perform synthesis to the temporary array (longitude local)
3. Transpose the result to the target pencil orientation

This involves one extra MPI all-to-all communication (the transpose) but ensures the FFT can operate on contiguous local data.

phi1_action_krylov(BA, LU_A, v, dt; m=20, tol=1e-8) -> y ≈ φ1(dt A) v

Compute φ1(dt A) v = A^{-1}[(exp(dt A) − I) v]/dt using Krylov exp(action) and banded solve.

print_geodynamo_banner(config, nprocs::Int, nthreads::Int)

Print the GeoDynamo startup banner with configuration information.

print_output_verification_report(output_dir::String, hist_numbers::Vector{Int}, nprocs::Int;
                                geometry::String="shell")

Print a comprehensive report verifying output files from all ranks across multiple history snapshots.

print_pencil_axes(pencils)

Print the axes_local tuple for each pencil, showing the local index ranges for all three axes. This helps verify which axes are distributed (those with nontrivial subranges across ranks) and which axis is contiguous locally.

print_pencil_info(pencils)

Print detailed information about pencil decomposition.

print_shtnskit_config_summary(nlat, nlon, lmax, mmax, nlm, nprocs, memory_estimate)

Print configuration summary for SHTnsKit setup.

print_transpose_statistics()

Print accumulated transpose timing statistics.

read_distributed_metadata(filename::String)

Read metadata from a distributed NetCDF file to determine global configuration.

read_file_metadata(filename::String)

Read metadata from NetCDF file to determine grid configuration.

report_energy_conservation(step::Int; interval::Int=100)

Report energy conservation statistics periodically.

report_phi2_conditioning(step::Int; interval::Int=100)

Report φ₂ conditioning statistics periodically during simulation.

report_solenoidal_constraint(step::Int; interval::Int=100)

Report solenoidal constraint violations periodically.

reset_energy_tracker!()

Reset energy tracking data.

reset_phi2_monitor!()

Reset φ₂ conditioning monitor statistics.

reset_solenoidal_monitor!()

Reset solenoidal constraint monitoring data.

rotate_field_90x!(config::SHTnsKitConfig, alm::Matrix{ComplexF64};
                  alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Rotate a scalar field by 90 degrees around the x-axis. Equivalent to: Z(-π/2) * Y(π/2) * Z(π/2)

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• alm_out: Optional output array

Returns

The rotated coefficients

rotate_field_90y!(config::SHTnsKitConfig, alm::Matrix{ComplexF64};
                  alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Rotate a scalar field by 90 degrees around the y-axis. This is a special case with optimized Wigner d-matrix values.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• alm_out: Optional output array

Returns

The rotated coefficients

rotate_field_euler!(config::SHTnsKitConfig, alm::Matrix{ComplexF64},
                    alpha::Real, beta::Real, gamma::Real;
                    alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Rotate a scalar field by Euler angles (ZYZ convention) in spectral space. The rotation is: R = Rz(gamma) * Ry(beta) * Rz(alpha)

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• alpha, beta, gamma: Euler angles in radians (ZYZ convention)
• alm_out: Optional output array

Returns

The rotated coefficients

rotate_field_y!(config::SHTnsKitConfig, alm::Matrix{ComplexF64}, beta::Real;
                alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Rotate a scalar field around the y-axis by angle beta in spectral space. Uses Wigner d-matrices (small Wigner rotation matrices).

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• beta: Rotation angle in radians
• alm_out: Optional output array

Returns

The rotated coefficients

rotate_field_z!(config::SHTnsKitConfig, alm::Matrix{ComplexF64}, alpha::Real;
                alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Rotate a scalar field around the z-axis by angle alpha in spectral space. This is a pure phase rotation: alm[l,m] -> alm[l,m] * exp(-imalpha)

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients (modified in-place if alm_out is nothing)
• alpha: Rotation angle in radians
• alm_out: Optional output array (if nothing, modifies alm in-place)

Returns

The rotated coefficients (alm_out if provided, otherwise alm)

run_enhanced_simulation!(state::SimulationState{T})

Run geodynamo simulation with all parallelization optimizations.

run_master_simulation!(state::SimulationState{T})

Run geodynamo simulation with maximum CPU parallelization optimizations.

save_combined_fields(combiner::FieldCombiner{T}, output_filename::String; 
                    config::CombinerConfig = create_combiner_config()) where T

Save combined fields to NetCDF file using GeoDynamo I/O system.

save_combined_time_series(combiners::Vector{FieldCombiner{T}}, output_dir::String, 
                         filename_prefix::String; 
                         config::CombinerConfig = create_combiner_config()) where T

Save combined time series to a single NetCDF file.

save_erk2_cache_bundle(path, caches; metadata=Dict())

Persist a dictionary of ERK2 caches to disk along with optional metadata.

save_parameters(params::GeoDynamoParameters, filename::String)

Save parameters to a Julia file.

save_physical_fields(converter::SpectralToPhysicalConverter{T}, output_filename::String) where T

Save converted physical space fields to NetCDF file using parallel I/O.

set_parameters!(params::GeoDynamoParameters; validate::Bool=true, strict::Bool=false)

Set the global parameters (thread-safe) with optional validation.

Arguments

• params::GeoDynamoParameters: Parameters to set
• validate::Bool=true: Whether to validate parameters before setting
• strict::Bool=false: If true, throw error on invalid parameters; if false, proceed with warnings

set_shtnskit_threads(num_threads::Int)

Configure the number of threads used by SHTnsKit transforms. Uses SHTnsKit.shtnsusethreads when available.

set_velocity_workspace!(ws::Union{VelocityWorkspace{T}, Nothing}) where T

Register a global VelocityWorkspace to be used by velocity kernels when available. Pass nothing to disable and fall back to internal buffers. Type-stable version that only accepts VelocityWorkspace or nothing.

setup_shtns_metadata!(nc_file, shtns_config::SHTnsKitConfig, pencils::NamedTuple)

Add SHTns-specific metadata to NetCDF file

shtnskit_analysis_inplace!(config::SHTnsKitConfig, f::Matrix{Float64},
                            alm_out::Matrix{ComplexF64})

In-place physical-to-spectral analysis using SHTnsKit v1.1.15. Writes result directly to alm_out, avoiding allocation of temporary arrays.

Arguments

• config: SHTnsKit configuration
• f: Input physical field (nlat × nlon)
• alm_out: Output spectral coefficients (lmax+1 × mmax+1), modified in-place

shtnskit_physical_to_spectral!(phys, spec)

Transform physical space values to spectral coefficients (analysis).

The Analysis Operation

Given physical field values f(θ,φ), compute spherical harmonic coefficients: fl^m = ∫∫ f(θ,φ) × Yl^m*(θ,φ) sin(θ) dθ dφ

The integral is computed numerically using:

• Gauss-Legendre quadrature for the θ integral
• FFT for the φ integral (exploiting periodicity)

Implementation

Uses SHTnsKit's analysis function which:

1. Performs the FFT along longitude (extracting Fourier modes)
2. Performs the Legendre transform (computing l coefficients for each m)

Arguments

• phys::SHTnsPhysicalField: Source physical field values
• spec::SHTnsSpectralField: Destination spectral field (modified in-place)

Side Effects

Modifies spec.data_real and spec.data_imag with the computed coefficients

shtnskit_qst_to_spatial!(config::SHTnsKitConfig, Qlm, Slm, Tlm, vr, vtheta, vphi)

Convert QST spectral coefficients to full 3D spatial vector field using SHTnsKit v1.1.15.

This is more efficient than separate scalar + vector synthesis as it handles all three components in a single call.

Arguments

• Qlm: Q (radial) spectral coefficients
• Slm: S (spheroidal/poloidal) spectral coefficients
• Tlm: T (toroidal) spectral coefficients
• vr, vtheta, vphi: Output spatial components (modified in-place)

shtnskit_spatial_to_qst!(config::SHTnsKitConfig, vr, vtheta, vphi, Qlm, Slm, Tlm)

Convert full 3D spatial vector field to QST spectral coefficients using SHTnsKit v1.1.15.

Arguments

• vr, vtheta, vphi: Input spatial components
• Qlm: Output Q (radial) spectral coefficients (modified in-place)
• Slm: Output S (spheroidal/poloidal) spectral coefficients (modified in-place)
• Tlm: Output T (toroidal) spectral coefficients (modified in-place)

shtnskit_spectral_to_physical!(spec, phys)

Transform spectral coefficients to physical space values (synthesis).

The Synthesis Operation

Given spherical harmonic coefficients fl^m, compute the physical field: f(θ,φ) = Σ{l,m} fl^m × Yl^m(θ,φ)

where Y_l^m are the spherical harmonics.

Implementation

Uses SHTnsKit's synthesis function which:

1. Performs the Legendre transform (summing over l for each m)
2. Performs the FFT along longitude (summing over m)

Arguments

• spec::SHTnsSpectralField: Source spectral field with coefficients
• phys::SHTnsPhysicalField: Destination physical field (modified in-place)

Side Effects

Modifies phys.data with the synthesized field values

shtnskit_synthesis_inplace!(config::SHTnsKitConfig, alm::Matrix{ComplexF64},
                             f_out::Matrix{Float64})

In-place spectral-to-physical synthesis using SHTnsKit v1.1.15. Writes result directly to f_out, avoiding allocation of temporary arrays.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients (lmax+1 × mmax+1)
• f_out: Output physical field (nlat × nlon), modified in-place

shtnskit_vector_analysis!(vec_phys::SHTnsVectorField{T},
                         tor_spec::SHTnsSpectralField{T},
                         pol_spec::SHTnsSpectralField{T}) where T

Vector analysis using SHTnsKit with PencilArrays.

Toroidal-Poloidal Analysis

Decomposes a 3D velocity field into toroidal and poloidal scalars.

For a solenoidal vector field v (∇·v = 0): v = ∇×(T r̂) + ∇×∇×(P r̂)

Mathematical Note on Analysis

SHTnsKit.spattoSHsphtor takes (vθ, vφ) and returns (P, T).

This is mathematically valid for solenoidal fields because:

1. The solenoidal constraint ∇·v = 0 couples vr to (vθ, v_φ)
2. The decomposition into T (rotational) and P (potential) parts is unique
3. The radial component v_r = l(l+1)/r * P is implicitly determined

However, this assumes EXACT solenoidality. In numerical simulations with finite precision, v_r may not exactly satisfy ∇·v = 0.

Alternative: Use Full 3-Component Analysis

For better numerical accuracy, one could use: Qcoeffs = analysis(vr * r / l(l+1)) # Recover P from vr S, T = spattoSHsphtor(vθ, v_φ) # Decompose tangential

Then check: Qcoeffs ≈ Scoeffs (should match for solenoidal field)

Current implementation uses 2-component analysis which is standard practice for solenoidal MHD simulations.

shtnskit_vector_synthesis!(tor_spec::SHTnsSpectralField{T},
                          pol_spec::SHTnsSpectralField{T},
                          vec_phys::SHTnsVectorField{T};
                          domain::Union{RadialDomain,Nothing}=nothing) where T

Vector synthesis using SHTnsKit spheroidal-toroidal decomposition with PencilArrays.

Toroidal-Poloidal Decomposition

For a solenoidal vector field v (∇·v = 0): v = ∇×(T r̂) + ∇×∇×(P r̂)

where T = toroidal scalar, P = poloidal scalar.

In spherical components: vr = l(l+1)/r * P * Ylm (from poloidal only) vθ, vφ from both T and P (computed by SHTnsKit.SHsphtortospat)

CRITICAL: SHTnsKit.SHsphtortospat returns ONLY tangential components. The radial component v_r MUST be computed separately from the poloidal scalar.

Parameters

• domain: Optional RadialDomain needed for computing vr with correct radial scaling. If not provided, vr will be set to zero (suitable for tests).

spectral_gradient!(config::SHTnsKitConfig, Slm::Matrix{ComplexF64},
                   grad_theta::Matrix{Float64}, grad_phi::Matrix{Float64})

Compute the horizontal gradient of a scalar field in spectral space. Uses SHTnsKit.SHtograd_spat for efficient computation.

Arguments

• Slm: Spectral coefficients of the scalar field
• grad_theta: Output θ-component of gradient (modified in-place)
• grad_phi: Output φ-component of gradient (modified in-place)

store_coefficients_from_shtnskit!(spec_real, spec_imag, coeffs_matrix, r_local, config)

Convert SHTnsKit coefficient matrix format back to our linear spectral storage.

This is the inverse of extract_coefficients_for_shtnskit!:

• Input: SHTnsKit's (lmax+1) × (mmax+1) coefficient matrix
• Output: Our linear-indexed spectral arrays (real and imaginary parts separate)

Physical Constraint

For real-valued physical fields, the m=0 coefficients must be purely real. This function enforces this by zeroing the imaginary part for m=0 modes.

store_physical_slice_generic!(phys_data, phys_slice, r_local, config)

Generic storage for any pencil orientation.

store_physical_slice_phi_local!(phys_data, phys_slice, r_local, config)

Copy a 2D physical field slice into a 3D array at a specific radial level.

Optimized for Phi-Local Layout

When the physical field is in phi-pencil orientation, the (θ, φ) indices correspond directly to the array's first two dimensions, making this a straightforward copy operation.

Arguments

• phys_data: 3D destination array (θ × φ × r)
• phys_slice: 2D source array from SHTnsKit (θ × φ)
• r_local: Radial index to write to
• config: Configuration for grid dimensions

store_scalar_component_generic!(v_component, field, r_local, config)

Store a scalar field into a component array for any pencil orientation. Used for storing the radial component v_r from synthesized field.

store_vector_components_generic!(v_theta, v_phi, vt_field, vp_field, r_local, config)

Store vector components for any pencil orientation.

store_zero_component_generic!(v_component, r_local, config)

Set a component to zero at a given radial level. Used at r=0 (ball geometry) where v_r must be zero for regularity.

synchronize_halos!(arr::PencilArray; halo_width::Int=1, boundaries::Symbol=:all)

Synchronize ghost/halo regions for parallel finite difference computations.

Parameters

• arr: PencilArray to synchronize
• halo_width: Width of halo region (number of ghost cells), default 1
• boundaries: Boundary handling (:all, :periodic, :nonperiodic), default :all

Notes

This function implements MPI point-to-point communication to exchange boundary data between neighboring processes in each decomposed dimension. It's essential for maintaining accuracy in finite difference stencil operations near subdomain boundaries.

For spectral methods, explicit halo exchange may not be needed as the global nature of spectral transforms handles boundary coupling through transpose operations.

synchronize_pencil_data!(field)

Synchronize PencilArray data across MPI processes to ensure consistency.

synchronize_pencil_transforms!(field::SHTnsSpectralField{T}) where T

Ensure all pending PencilFFTs operations are completed and data is consistent across processes.

test_halo_exchange(pencil::Pencil, ::Type{T}=Float64; halo_width::Int=1, verbose::Bool=true) where T

Test halo exchange functionality by creating a test array with rank-based values. Returns true if halo exchange is working correctly.

Example

pencils = create_pencil_topology(shtns_config)
test_halo_exchange(pencils.θ, Float64; halo_width=1, verbose=true)

transform_field_and_gradients_to_physical!(field::AbstractScalarField{T}) where T

Transform scalar field and all gradient components to physical space in a single batched operation to minimize communication.

transpose_with_timer!(dest, src, plan, label="")

Perform transpose with optional timing and statistics.

truncate_spectral_modes!(config::SHTnsKitConfig, alm::Matrix{ComplexF64},
                         lmax_new::Int, mmax_new::Int=lmax_new;
                         alm_out::Union{Matrix{ComplexF64},Nothing}=nothing)

Truncate spectral coefficients to a lower resolution. Sets all modes with l > lmaxnew or m > mmaxnew to zero.

Arguments

• config: SHTnsKit configuration
• alm: Input spectral coefficients
• lmax_new: New maximum degree
• mmax_new: New maximum order (default = lmax_new)
• alm_out: Optional output array

update_derived_parameters!(params::GeoDynamoParameters)

Update derived parameters based on primary parameters.

update_global_parameters!()

Update all global parameter variables with values from the current parameter struct.

validate_magnetic_configuration(mag_fields::SHTnsMagneticFields{T}, config::SHTnsKitConfig) where T

Validate magnetic field configuration consistency with SHTns setup

validate_mpi_consistency!(field::SHTnsSpectralField{T}) where T

Check that spectral field data is consistent across MPI processes after time stepping.

validate_output_compatibility(field_info::FieldInfo, shtns_config::SHTnsKitConfig)

Validate that output field information is compatible with SHTns configuration

validate_parameters(params::GeoDynamoParameters; strict::Bool=false)

Validate all simulation parameters for physical correctness and numerical stability.

Arguments

• params::GeoDynamoParameters: Parameters to validate
• strict::Bool=false: If true, throw errors on invalid params; if false, issue warnings

Returns

• (is_valid::Bool, errors::Vector{String}, warnings::Vector{String})

validate_pencil_decomposition(config::SHTnsKitConfig)

Validate that pencil decomposition is optimal for the problem size and MPI configuration.

validate_stress_free_boundary(v_r, v_theta, v_phi, r_val, theta, phi; tolerance=0.05)

Validate that velocity field satisfies stress-free boundary condition.

For stress-free boundaries with vr = 0, the condition is: ∂vθ/∂r - vθ/r = 0 (zero tangential stress in θ direction) ∂vφ/∂r - v_φ/r = 0 (zero tangential stress in φ direction)

Arguments

• v_r, v_theta, v_phi: Velocity components at boundary [nlat, nlon]
• r_val: Radial position of boundary
• theta, phi: Coordinate arrays
• tolerance: Maximum allowed |stress| / |v|

Returns

• is_valid: Boolean indicating if condition is satisfied
• max_violation: Maximum relative stress violation

validate_velocity_configuration(fields::SHTnsVelocityFields{T}, config::SHTnsKitConfig) where T

Validate velocity field configuration consistency with SHTns setup

verify_all_ranks_wrote(output_dir::String, hist_number::Int, nprocs::Int;
                      file_type::String="hist", geometry::String="shell")

Verify that all MPI ranks wrote their output files for a given history number.

Returns a tuple (success::Bool, missingranks::Vector{Int}, fileinfo::Dict)

write_enhanced_coordinates!(nc_file, field_info::FieldInfo, config::OutputConfig)

Write coordinates with enhanced precision and metadata from SHTns config

write_grid_file!(config::OutputConfig, field_info::FieldInfo,
                shtns_config::Union{SHTnsKitConfig,Nothing},
                metadata::Dict{String,Any})

Write a separate grid file containing all coordinate and grid information. This is written only once by rank 0 at the start of the simulation.

write_single_spectral_component!(nc_file, component, field_data, config)

Write a single spectral component with enhanced memory handling

write_spectral_data_parallel!(nc_file, real_name, imag_name, real_data, imag_data, field_info)

Write spectral data using parallel I/O strategies based on pencil decomposition

zero_scalar_work_arrays!(field::AbstractScalarField{T}) where T

Efficiently zero all work arrays for a scalar field.

AbstractBoundaryCondition{T}

Abstract base type for all boundary conditions.

BoundaryCache{T}

Cache structure for storing processed boundary data.

BoundaryConditionSet{T}

Complete set of boundary conditions for inner and outer boundaries.

BoundaryData{T}

Common data structure for boundary condition data from any source.

BoundaryLocation

Enumeration for boundary locations.

BoundaryType

Enumeration for boundary condition types.

FieldType

Enumeration for different physical field types.

_get_cached_bc_shtns_config(lmax, mmax, nlat, nlon)

Get or create a cached SHTnsKit configuration for boundary transforms. Reuses configurations to avoid repeated setup/teardown overhead.

add_noise_to_boundary(boundary_data::BoundaryData, noise_amplitude::Real, 
                     noise_type::Symbol=:gaussian)

Add noise to existing boundary data.

apply_boundary_conditions!(field, field_type::FieldType, solver_state)

Apply boundary conditions during solver operations.

This function integrates boundary conditions with the timestepping and solving process.

apply_boundary_conditions_to_rhs!(rhs, state, field_type::FieldType)

Apply boundary conditions to right-hand side during timestepping.

This function modifies the RHS vector to enforce boundary conditions during implicit and explicit timestepping methods.

apply_composition_bc_to_rhs!(rhs, comp_field)

Apply composition boundary conditions to RHS vector.

apply_composition_boundaries!(comp_field, boundary_set::BoundaryConditionSet; time::Float64=0.0)

Apply composition boundary conditions from a BoundaryConditionSet to a composition field. This is a convenience wrapper that loads the boundary set and applies it.

apply_composition_boundary_conditions!(comp_field, time_index::Int=1)

Apply composition boundary conditions to the field.

apply_gaussian_smoothing!(data::AbstractMatrix, theta::Vector, phi::Vector, radius::Real)

Apply Gaussian smoothing kernel to 2D data array.

apply_magnetic_bc_to_rhs!(rhs, magnetic_field)

Apply magnetic field boundary conditions to RHS vector.

apply_magnetic_boundary_conditions!(magnetic_field, time_index::Int=1)

Apply magnetic field boundary conditions to the field.

apply_temperature_bc_to_rhs!(rhs, temp_field)

Apply temperature boundary conditions to RHS vector.

apply_temperature_boundaries!(temp_field, boundary_set::BoundaryConditionSet; time::Float64=0.0)

Apply temperature boundary conditions from a BoundaryConditionSet to a temperature field. This is a convenience wrapper that loads the boundary set and applies it.

apply_temperature_boundary_conditions!(temp_field, time_index::Int=1)

Apply temperature boundary conditions to the field.

apply_velocity_bc_to_rhs!(rhs, velocity_field)

Apply velocity boundary conditions to RHS vector.

apply_velocity_boundary_conditions!(velocity_field, time_index::Int=1)

Apply velocity boundary conditions to the field.

bilinear_interpolate(data::Matrix{T}, theta_idx::Tuple{Int, Int}, phi_idx::Tuple{Int, Int},
                    theta_weights::Tuple{T, T}, phi_weights::Tuple{T, T}) where T

Perform bilinear interpolation on a 2D data array. Handles periodic boundary conditions in phi direction.

cache_boundary_data!(cache::BoundaryCache{T}, key::String, data::Array{T}) where T

Cache processed boundary data.

calculate_potential_field_boundary(coeffs::Dict, theta::Vector, phi::Vector, lmax::Int)

Calculate magnetic field from spherical harmonic coefficients of the potential.

check_interpolation_bounds(boundary_data::BoundaryData, target_theta::Vector{T}, 
                          target_phi::Vector{T}) where T

Check if target grid is within the bounds of the source grid and warn if extrapolation is needed.

clear_bc_shtns_config_cache!()

Clear the cached SHTnsKit configurations for boundary transforms. Call this when grid parameters change or to free memory.

clear_boundary_cache!(cache::BoundaryCache)

Clear all cached boundary data.

combine_programmatic_patterns(patterns::Vector{Tuple{Symbol, Real}}, config; 
                             parameters::Vector{Dict}=Dict[], description::String="")

Combine multiple programmatic patterns with different amplitudes.

compute_boundary_condition_residual(field, field_type::FieldType)

Compute residual to check how well boundary conditions are satisfied.

Returns a measure of how much the current field violates the boundary conditions. Useful for monitoring solution quality and debugging.

compute_radial_curl_component(B_r, B_theta, B_phi, theta, phi)

Compute the radial component of ∇×B to verify insulating boundary conditions.

For insulating boundaries, (∇×B)r should be zero. In spherical coordinates: (∇×B)r = (1/(r sin θ))[∂(Bφ sin θ)/∂θ - ∂Bθ/∂φ]

Arguments

• B_r, B_theta, B_phi: Magnetic field components on boundary [nlat, nlon]
• theta, phi: Coordinate arrays

Returns

• curl_r: Radial component of curl [nlat, nlon]
• max_curl: Maximum absolute value of (∇×B)_r

copy_boundary_conditions!(dest_field, src_field, field_type::FieldType)

Copy boundary conditions from one field to another.

create_boundary_data(values::Array{T}, field_type::String; 
                    theta=nothing, phi=nothing, time=nothing,
                    units="", description="", file_path="programmatic") where T

Create a BoundaryData structure from raw data.

create_composition_interpolation_cache(boundary_set::BoundaryConditionSet, config)

Create interpolation cache for composition boundaries.

create_hybrid_composition_boundaries(file_spec::String, prog_spec::Tuple, config; swap_boundaries=false)

Create hybrid composition boundaries (one from file, one programmatic).

create_hybrid_magnetic_boundaries(file_spec::String, prog_spec::Tuple, config; swap_boundaries=false)

Create hybrid magnetic boundaries (one from file, one programmatic).

create_hybrid_temperature_boundaries(file_spec::String, prog_spec::Tuple, config; swap_boundaries=false)

Create hybrid temperature boundaries (one from file, one programmatic).

create_hybrid_velocity_boundaries(file_spec::String, prog_spec::Tuple, config; swap_boundaries=false)

Create hybrid velocity boundaries (one from file, one programmatic).

create_interpolation_cache(boundary_data::BoundaryData, target_theta::Vector{T}, 
                          target_phi::Vector{T}) where T

Create interpolation cache for efficient repeated interpolations.

create_layered_composition_boundary(config, layer_specs::Vector{Tuple{Real, Real, Real}})

Create layered composition boundary conditions.

Arguments

• layer_specs: Vector of (colatitudestart, colatitudeend, composition) tuples

Example

# Create three-layer composition structure
layer_specs = [
    (0.0, π/3, 0.8),     # High composition in top layer
    (π/3, 2π/3, 0.4),    # Medium composition in middle layer  
    (2π/3, π, 0.1)       # Low composition in bottom layer
]

create_layered_temperature_boundary(config, layer_specs::Vector{Tuple{Real, Real, Real}})

Create layered temperature boundary conditions.

Arguments

• layer_specs: Vector of (colatitudestart, colatitudeend, temperature) tuples

Example

# Create three-layer temperature structure
layer_specs = [
    (0.0, π/3, 1000.0),    # High temperature in top layer
    (π/3, 2π/3, 500.0),   # Medium temperature in middle layer  
    (2π/3, π, 100.0)      # Low temperature in bottom layer
]

create_magnetic_interpolation_cache(boundary_set::BoundaryConditionSet, config)

Create interpolation cache for magnetic boundaries.

create_programmatic_boundary(pattern::Symbol, config, amplitude::Real=1.0;
                            parameters::Dict=Dict(), description::String="")

Create programmatically generated boundary conditions.

Available patterns:

• :uniform - Uniform value
• :y11 - Y₁₁ spherical harmonic
• :plume - Gaussian plume pattern
• :hemisphere - Hemispherical pattern
• :dipole - Dipolar pattern (Y₁₀)
• :quadrupole - Quadrupolar pattern
• :custom - User-defined function

create_programmatic_composition_boundaries(inner_spec::Tuple, outer_spec::Tuple, config)

Create fully programmatic composition boundaries.

create_programmatic_magnetic_boundaries(inner_spec::Tuple, outer_spec::Tuple, config)

Create fully programmatic magnetic boundaries.

create_programmatic_magnetic_boundary(pattern::Symbol, config, amplitude::Real=1.0; 
                                    parameters::Dict=Dict())

Create programmatically generated magnetic boundary conditions.

Available patterns:

• :insulating - Insulating boundary (Br = 0, ∂Btan/∂r = 0)
• :perfect_conductor - Perfect conductor (B_tan = 0)
• :dipole - Dipolar magnetic field pattern
• :quadrupole - Quadrupolar magnetic field pattern
• :potential_field - Potential field from spherical harmonic coefficients
• :uniform_field - Uniform magnetic field
• :custom - User-defined magnetic field function

create_programmatic_temperature_boundaries(inner_spec::Tuple, outer_spec::Tuple, config)

Create fully programmatic temperature boundaries.

create_programmatic_velocity_boundaries(inner_spec::Tuple, outer_spec::Tuple, config)

Create fully programmatic velocity boundaries.

create_programmatic_velocity_boundary(pattern::Symbol, config, amplitude::Real=1.0; 
                                    parameters::Dict=Dict())

Create programmatically generated velocity boundary conditions.

Available patterns:

• :no_slip - Zero velocity at boundary (amplitude ignored)
• :stress_free - Zero stress at boundary (amplitude ignored)
• :uniform_rotation - Uniform rotation with angular velocity amplitude
• :differential_rotation - Differential rotation pattern
• :zonal_flow - Zonal (east-west) flow pattern
• :meridional_flow - Meridional (north-south) flow pattern
• :custom - User-defined velocity function

create_temperature_interpolation_cache(boundary_set::BoundaryConditionSet, config)

Create interpolation cache for temperature boundaries.

create_time_dependent_programmatic_boundary(pattern::Symbol, config, time_span::Tuple{Real, Real}, 
                                           ntime::Int; amplitude::Real=1.0, 
                                           parameters::Dict=Dict(), description::String="",
                                           field_type::String="temperature")

Create time-dependent programmatically generated boundary conditions.

create_velocity_interpolation_cache(boundary_set::BoundaryConditionSet, config)

Create interpolation cache for velocity boundaries.

determine_field_type_from_name(field_name::String) -> FieldType

Determine field type from field name string.

enforce_boundary_conditions_in_solution!(solution, state, field_type::FieldType)

Enforce boundary conditions in the solution vector after timestepping.

This function ensures that the solution satisfies the boundary conditions after each timestep, which may be necessary for certain discretization schemes.

enforce_composition_bc_in_solution!(solution, comp_field)

Enforce composition boundary conditions in solution vector.

enforce_composition_boundary_constraints!(comp_field, bc_spec::Dict)

Enforce specific composition boundary constraints based on user specification.

Arguments

• comp_field: Composition field structure
• bc_spec: Dictionary specifying boundary condition types

BC Specification Format

bc_spec = Dict(
    :inner => :dirichlet,  # or :neumann, :flux
    :outer => :dirichlet,  # or :neumann, :flux
    :inner_value => 1.0,   # for Dirichlet (0.0-1.0)
    :outer_value => 0.0,   # for Dirichlet (0.0-1.0)
    :inner_flux => 0.01,   # for Neumann (∂C/∂r)
    :outer_flux => 0.0     # for Neumann (∂C/∂r)
)

Boundary Condition Types

• :dirichlet - Fixed composition: C = C₀ at boundary
• :neumann or :flux - Fixed compositional flux: ∂C/∂r = q at boundary

Physical Interpretation

• Dirichlet (Fixed C): Typical for prescribed composition boundaries

  • Inner (ICB): Fixed light element concentration at inner core
  • Outer (CMB): Fixed composition at core-mantle boundary

• Neumann (Fixed flux): Typical for compositional flux boundaries

  • Inner: Compositional release/freezing (∂C/∂r = q)
  • Outer: Zero flux (no exchange, ∂C/∂r = 0)

Examples

# Fixed composition at both boundaries
enforce_composition_boundary_constraints!(comp_field, Dict(
    :inner => :dirichlet, :inner_value => 1.0,
    :outer => :dirichlet, :outer_value => 0.0
))

# Compositional release from below, sealed top
enforce_composition_boundary_constraints!(comp_field, Dict(
    :inner => :flux, :inner_flux => 0.01,
    :outer => :flux, :outer_flux => 0.0
))

enforce_magnetic_bc_in_solution!(solution, magnetic_field)

Enforce magnetic field boundary conditions in solution vector.

enforce_magnetic_boundary_constraints!(magnetic_field, bc_type::Symbol)

Enforce magnetic boundary condition constraints based on magnetohydrodynamic physics.

Boundary condition types:

• :insulating - Insulating boundary (σ = 0): No normal current (J_n = 0)

  • Physical interpretation: (∇×B)r = μ₀Jr = 0 at boundary
  • Implementation: Potential field matching (B continuous, tangential current-free)
  • Spectral: Typically Dirichlet BC for both components matching external potential

• :perfect_conductor - Perfect conductor (σ → ∞): Zero tangential magnetic field

  • Physical interpretation: Btangential = 0, Br free (from ∇·B = 0)
  • Implementation: Tangential components zero, radial component determined by matching
  • Spectral: Toroidal (T) = 0 Dirichlet, Poloidal (Q) matches

• :potential_field - General potential field matching

  • Implementation: Match external potential field at boundary
  • Spectral: Dirichlet BC for both toroidal and poloidal components

Notes:

For most geodynamo applications:

• Inner boundary (ICB): Often insulating or potential field
• Outer boundary (CMB): Often insulating (matching Earth's mantle)

enforce_temperature_bc_in_solution!(solution, temp_field)

Enforce temperature boundary conditions in solution vector.

enforce_temperature_boundary_constraints!(temp_field, bc_spec::Dict)

Enforce specific temperature boundary constraints based on user specification.

Arguments

• temp_field: Temperature field structure
• bc_spec: Dictionary specifying boundary condition types

BC Specification Format

bc_spec = Dict(
    :inner => :dirichlet,  # or :neumann, :flux
    :outer => :dirichlet,  # or :neumann, :flux
    :inner_value => 1000.0,  # for Dirichlet
    :outer_value => 300.0,   # for Dirichlet
    :inner_flux => 0.1,      # for Neumann (∂T/∂r)
    :outer_flux => 0.0       # for Neumann (∂T/∂r)
)

Boundary Condition Types

• :dirichlet - Fixed temperature: T = T₀ at boundary
• :neumann or :flux - Fixed heat flux: ∂T/∂r = q at boundary
• :mixed - Mixed/Robin: α T + β ∂T/∂r = γ (future)

Physical Interpretation

• Dirichlet (Fixed T): Typical for prescribed temperature boundaries

  • Inner (CMB): Fixed temperature from core evolution
  • Outer (surface): Fixed temperature from atmospheric coupling

• Neumann (Fixed flux): Typical for heat flux boundaries

  • Inner: Uniform heating from below (∂T/∂r = q)
  • Outer: Zero flux (insulated, ∂T/∂r = 0)

Examples

# Fixed temperature at both boundaries
enforce_temperature_boundary_constraints!(temp_field, Dict(
    :inner => :dirichlet, :inner_value => 1000.0,
    :outer => :dirichlet, :outer_value => 300.0
))

# Uniform heating from below, insulated top
enforce_temperature_boundary_constraints!(temp_field, Dict(
    :inner => :flux, :inner_flux => 0.1,
    :outer => :flux, :outer_flux => 0.0
))

enforce_velocity_bc_in_solution!(solution, velocity_field)

Enforce velocity boundary conditions in solution vector.

enforce_velocity_boundary_constraints!(velocity_field, bc_type::Symbol=:no_slip)

Enforce specific velocity boundary constraints based on boundary condition type.

NOTE: This function is for programmatically setting boundary constraints. For boundary conditions loaded from files or other sources, use loadvelocityboundary_conditions!() instead.

Arguments

• velocity_field: Velocity field structure with toroidal and poloidal components
• bc_type: Type of boundary condition (:noslip, :stressfree, :impermeable)

Boundary Condition Mapping (for solenoidal/incompressible flows):

• No-slip: vr = vθ = v_φ = 0 → Q = T = 0 at boundaries (Dirichlet)
• Stress-free: v_r = 0 (Dirichlet), tangential stress = 0 → Q = 0 (Dirichlet), ∂T/∂r = 0 (Neumann)
• Impermeable: v_r = 0 → Q = 0 at boundaries (Dirichlet), T unconstrained

Field naming convention:

• velocity_field.poloidal actually stores Q (radial component)
• velocity_field.toroidal actually stores T (tangential toroidal component)
• S component (spheroidal) is zero for solenoidal flows

estimate_interpolation_error(boundary_data::BoundaryData, target_theta::Vector{T}, 
                            target_phi::Vector{T}) where T

Estimate interpolation error based on grid resolution differences.

find_composition_boundary_time_index(boundary_set::BoundaryConditionSet, current_time::Float64)

Find the appropriate time index for the current simulation time.

find_grid_indices(coords::Vector{T}, target::T; is_periodic::Bool=false) where T

Find the two surrounding indices in a coordinate array for interpolation. Handles periodic coordinates (e.g., longitude) if is_periodic=true.

find_temperature_boundary_time_index(boundary_set::BoundaryConditionSet, current_time::Float64)

Find the appropriate time index for the current simulation time.

get_boundary_condition_summary(field, field_type::FieldType)

Get a summary of the current boundary condition state.

get_boundary_data(field, field_type::FieldType)

Get boundary data from field using unified interface with fallback support.

get_boundary_statistics(boundary_data::BoundaryData)

Get statistical information about boundary data.

get_cached_data(cache::BoundaryCache{T}, key::String) where T

Retrieve cached boundary data.

get_comm()

Get MPI communicator (with fallback).

get_composition_boundary_data(comp_field)

Get boundary data from field or fallback cache.

get_composition_time_index(comp_field)

Get current time index from field or fallback cache.

get_current_boundaries(field, field_type::FieldType)

Get current boundary values for any field type.

get_current_composition_boundaries(comp_field)

Get current composition boundary conditions.

get_current_magnetic_boundaries(magnetic_field)

Get current magnetic field boundary conditions.

get_current_simulation_time(solver_state)

Extract current simulation time from solver state.

get_current_temperature_boundaries(temp_field)

Get current temperature boundary conditions.

get_current_velocity_boundaries(velocity_field)

Get current velocity boundary conditions.

get_default_boundary_type(field_type::FieldType, location::BoundaryLocation) -> BoundaryType

Get default boundary condition type for a field at a specific location.

get_default_units(field_type::FieldType) -> String

Get default units for a field type.

get_field_from_state(state, field_type::FieldType)

Extract the appropriate field from simulation state.

get_interpolation_statistics(boundary_data::BoundaryData, interpolated_data::Array{T}) where T

Compute statistics comparing original and interpolated data.

get_interpolation_weights(coords::Vector{T}, target::T, indices::Tuple{Int, Int}) where T

Calculate interpolation weights for linear interpolation.

get_magnetic_boundary_data(magnetic_field)

Return (boundary_set, cache) for the magnetic field, falling back to a module-level cache when the field struct lacks boundary storage.

get_magnetic_time_index(magnetic_field)

Fetch the currently active boundary time index, honoring legacy cache storage.

get_netcdf_file_info(filename::String)

Get information about a NetCDF boundary condition file.

get_nprocs()

Get number of MPI processes (with fallback).

get_rank()

Get MPI rank (with fallback).

get_temperature_boundary_data(temp_field)

Get boundary data from field or fallback cache.

get_temperature_time_index(temp_field)

Get current time index from field or fallback cache.

get_time_index(field, field_type::FieldType)

Get current time index from field using unified interface with fallback support.

get_velocity_boundary_data(velocity_field)

Return (boundary_set, cache) for the provided velocity field, falling back to module-level storage when the struct does not carry boundary metadata.

get_velocity_time_index(velocity_field)

Fetch the current boundary time index for a velocity field, honoring fallback storage.

infer_composition_bc_type(boundary::BoundaryData)

Infer boundary condition type (Dirichlet or Neumann) from boundary metadata.

Checks the description field for keywords:

• "flux", "neumann" → NEUMANN (∂C/∂r specified)
• Otherwise → DIRICHLET (C specified)

infer_temperature_bc_type(boundary::BoundaryData)

Infer boundary condition type (Dirichlet or Neumann) from boundary metadata.

Checks the description field for keywords:

• "flux", "neumann", "heat flux" → NEUMANN (∂T/∂r specified)
• Otherwise → DIRICHLET (T specified)

infer_velocity_bc_type(boundary::BoundaryData)

Infer boundary condition type (Dirichlet or Neumann) from boundary metadata.

initialize_boundary_conditions!(field, field_type::FieldType, config)

Initialize boundary condition support for a field structure.

Adds the necessary boundary condition fields to existing field structures without breaking compatibility with existing code.

interpolate_boundary_to_grid(boundary_data::BoundaryData, target_theta::Vector{T}, 
                            target_phi::Vector{T}, time_index::Int=1) where T

Interpolate boundary data to a target grid using bilinear interpolation.

interpolate_with_cache(boundary_data::BoundaryData, cache::Dict, time_index::Int=1)

Perform interpolation using pre-computed cache for efficiency.

load_composition_boundaries_from_files(inner_file::String, outer_file::String, config)

Load composition boundary conditions from NetCDF files.

load_composition_boundary_conditions!(comp_field, boundary_specs::Dict)

Load composition boundary conditions from various sources.

Arguments

• comp_field: SHTnsCompositionField structure
• boundary_specs: Dictionary specifying boundary sources

Examples

# NetCDF files for both boundaries
boundary_specs = Dict(
    :inner => "cmb_composition.nc",
    :outer => "surface_composition.nc"
)

# Mixed NetCDF and programmatic
boundary_specs = Dict(
    :inner => "cmb_composition.nc", 
    :outer => (:uniform, 0.1)  # 10% concentration
)

load_magnetic_boundaries_from_files(inner_file::String, outer_file::String, config)

Load magnetic field boundary conditions from NetCDF files.

load_magnetic_boundary_conditions!(magnetic_field, boundary_specs::Dict)

Load magnetic field boundary conditions from various sources.

Arguments

• magnetic_field: SHTnsMagneticField structure
• boundary_specs: Dictionary specifying boundary sources

Examples

# Insulating inner, potential field outer
boundary_specs = Dict(
    :inner => (:insulating, 0.0),
    :outer => (:potential_field, "geomagnetic_coefficients.nc")
)

# Perfect conductor boundaries
boundary_specs = Dict(
    :inner => (:perfect_conductor, 0.0),
    :outer => (:perfect_conductor, 0.0)
)

# NetCDF files for both boundaries
boundary_specs = Dict(
    :inner => "cmb_magnetic.nc",
    :outer => "surface_magnetic.nc"
)

load_temperature_boundaries_from_files(inner_file::String, outer_file::String, config)

Load temperature boundary conditions from NetCDF files.

load_temperature_boundary_conditions!(temp_field, boundary_specs::Dict)

Load temperature boundary conditions from various sources.

Arguments

• temp_field: SHTnsTemperatureField structure
• boundary_specs: Dictionary specifying boundary sources

Examples

# NetCDF files for both boundaries
boundary_specs = Dict(
    :inner => "cmb_temperature.nc",
    :outer => "surface_temperature.nc"
)

# Mixed NetCDF and programmatic
boundary_specs = Dict(
    :inner => "cmb_temperature.nc", 
    :outer => (:uniform, 300.0)
)

load_velocity_boundaries_from_files(inner_file::String, outer_file::String, config)

Load velocity boundary conditions from NetCDF files.

load_velocity_boundary_conditions!(velocity_field, boundary_specs::Dict)

Load velocity boundary conditions from various sources.

Arguments

• velocity_field: SHTnsVelocityField structure
• boundary_specs: Dictionary specifying boundary sources

Examples

# No-slip boundaries at both surfaces
boundary_specs = Dict(
    :inner => (:no_slip, 0.0),
    :outer => (:no_slip, 0.0)
)

# Stress-free boundaries
boundary_specs = Dict(
    :inner => (:stress_free, 0.0),
    :outer => (:stress_free, 0.0)
)

# NetCDF file for inner, no-slip for outer
boundary_specs = Dict(
    :inner => "cmb_velocity.nc",
    :outer => (:no_slip, 0.0)
)

log_boundary_condition_status(state, rank::Int=0)

Log the status of all boundary conditions in the simulation state.

magnetic_to_qst_coefficients(B_r, B_theta, B_phi, config)

Convert physical magnetic field components to QST spectral coefficients.

For magnetic fields:

• Q: Radial component coefficients (B_r)
• S: Spheroidal tangential component coefficients
• T: Toroidal tangential component coefficients

print_boundary_data_info(boundary_data::BoundaryData, prefix::String="")

Print detailed information about boundary data.

print_boundary_info(boundary_set::BoundaryConditionSet)

Print comprehensive information about a boundary condition set.

print_boundary_summary(field, field_type::FieldType)

Print a summary of loaded boundary conditions for any field type.

read_netcdf_boundary_data(filename::String; precision::Type{T}=Float64) where T

Read boundary condition data from a NetCDF file.

Returns a BoundaryData structure containing all boundary information.

reset_boundary_conditions!(field, field_type::FieldType)

Reset/clear boundary conditions for a field.

set_programmatic_composition_boundaries!(comp_field, inner_spec::Tuple, outer_spec::Tuple)

Set programmatic composition boundary conditions.

set_programmatic_magnetic_boundaries!(magnetic_field, inner_spec::Tuple, outer_spec::Tuple)

Set programmatic magnetic boundary conditions.

set_programmatic_temperature_boundaries!(temp_field, inner_spec::Tuple, outer_spec::Tuple)

Set programmatic temperature boundary conditions.

set_programmatic_velocity_boundaries!(velocity_field, inner_spec::Tuple, outer_spec::Tuple)

Set programmatic velocity boundary conditions.

shtns_physical_to_spectral(physical_data::Matrix{T}, config; return_complex::Bool=false) where T

Transform physical boundary data to spectral coefficients using SHTnsKit. This is a common utility function used by both thermal and composition modules.

Arguments

• physical_data: 2D array of physical values on (nlat, nlon) grid
• config: SHTnsKit configuration with lmax, mmax, nlm fields
• return_complex: If true, returns complex coefficients; otherwise real part only

Returns

Vector of spectral coefficients of length nlm.

Performance (v1.1.15)

Uses cached SHTnsKit configurations to avoid repeated setup overhead.

shtns_spectral_to_physical(coeffs::Vector, config, nlat::Int, nlon::Int)

Transform spectral coefficients to physical boundary data using SHTnsKit. Inverse of shtnsphysicalto_spectral.

Arguments

• coeffs: Vector of spectral coefficients of length nlm
• config: SHTnsKit configuration with lmax, mmax fields
• nlat, nlon: Output grid dimensions

Returns

Matrix of physical values on (nlat, nlon) grid.

smooth_boundary_data(boundary_data::BoundaryData, smoothing_radius::Real)

Apply spatial smoothing to boundary data.

spherical_harmonic(l::Int, m::Int, theta::Real, phi::Real)

Compute spherical harmonic Y_l^m(θ, φ).

spherical_harmonic_phi_derivative(l::Int, m::Int, theta::Real, phi::Real)

Compute ∂Y_l^m/∂φ.

spherical_harmonic_theta_derivative(l::Int, m::Int, theta::Real, phi::Real)

Compute ∂Y_l^m/∂θ.

update_boundary_conditions_for_timestep!(state, current_time::Float64)

Update all boundary conditions for the current timestep.

This function should be called at the beginning of each timestep to ensure all fields have up-to-date boundary conditions.

update_composition_time_index!(comp_field, time_index::Int)

Update time index in field or fallback cache.

update_magnetic_time_index!(magnetic_field, time_index)

Persist the provided time index to both the magnetic field structure and the module-level cache when available.

update_temperature_time_index!(temp_field, time_index::Int)

Update time index in field or fallback cache.

update_time_dependent_boundaries!(field, field_type::FieldType, current_time::Float64)

Update time-dependent boundary conditions for any field type.

update_time_dependent_composition_boundaries!(comp_field, current_time::Float64)

Update time-dependent composition boundary conditions.

update_time_dependent_magnetic_boundaries!(magnetic_field, current_time::Float64)

Update time-dependent magnetic boundary conditions.

update_time_dependent_temperature_boundaries!(temp_field, current_time::Float64)

Update time-dependent temperature boundary conditions.

update_time_dependent_velocity_boundaries!(velocity_field, current_time::Float64)

Update time-dependent velocity boundary conditions.

update_velocity_time_index!(velocity_field, time_index)

Update cached time indices for velocity boundary conditions in both the field and the module-level fallback cache (when present).

validate_boundary_compatibility(inner::BoundaryData, outer::BoundaryData, field_name::String)

Validate that inner and outer boundary data are compatible.

validate_boundary_files(field_type::FieldType, boundary_specs::Dict, config)

Validate boundary condition files for any field type.

validate_composition_boundary_files(boundary_specs::Dict, config)

Validate composition boundary condition files.

validate_composition_range(boundary_data::BoundaryData)

Validate that composition values are in the physical range [0, 1].

validate_field_boundary_compatibility(field, field_type::FieldType, boundary_set::BoundaryConditionSet)

Validate that a field structure is compatible with boundary conditions.

validate_flux_boundary_values(values::Array, field_type::String;
                              typical_value_range::Tuple=(0.0, 1e10),
                              check_units::Bool=true)

Validate that boundary values are physically reasonable for flux boundary conditions.

For heat flux boundaries:

• Values should be in W/m² or dimensionless heat flux
• Typical range: 0.01 - 1000 W/m² for Earth
• Check that values aren't accidentally in Kelvin (which would be 100-1000× larger)

Arguments

• values: Boundary data array
• field_type: Type of field ("temperature", "composition")
• typical_value_range: Expected range for flux values
• check_units: Whether to perform unit sanity checks

Returns

• is_valid: Boolean
• warnings: Array of warning messages

validate_insulating_boundary(B_r, B_theta, B_phi, theta, phi; tolerance=1e-2)

Validate that a magnetic field satisfies the insulating boundary condition (∇×B)_r = 0.

Arguments

• B_r, B_theta, B_phi: Magnetic field components [nlat, nlon]
• theta, phi: Coordinate arrays
• tolerance: Maximum allowed |(∇×B)_r| / |B|

Returns

• is_valid: Boolean indicating if condition is satisfied
• max_violation: Maximum |(∇×B)_r| / |B|

validate_interpolation_grids(src_theta::Vector, src_phi::Vector, 
                            tgt_theta::Vector, tgt_phi::Vector)

Validate that source and target grids are compatible for interpolation.

validate_magnetic_boundary_files(boundary_specs::Dict, config)

Validate magnetic field boundary condition files.

validate_netcdf_boundary_file(filename::String, required_vars::Vector{String}=[])

Validate that a NetCDF boundary condition file has the required structure.

validate_temperature_boundary_files(boundary_specs::Dict, config)

Validate temperature boundary condition files.

validate_temperature_flux_boundary(boundary::BoundaryData, bc_type::Int)

Validate that a temperature boundary condition has correct data type for its BC type.

For Neumann (flux) BCs: ensures values represent heat flux, not temperature For Dirichlet BCs: ensures values represent temperature

validate_temperature_range(boundary_data::BoundaryData)

Validate that temperature values are in a reasonable physical range.

validate_velocity_boundary_files(boundary_specs::Dict, config)

Validate velocity boundary condition files.

velocity_to_qst_coefficients(v_r, v_theta, v_phi, config)

Convert physical velocity components (vr, vθ, v_φ) to QST spectral coefficients by using proper SHTnsKit decomposition.

The QST decomposition used by SHTnsKit:

• Q: Radial component coefficients (transforms like scalar field)
• S: Spheroidal horizontal component coefficients (curl-free part)
• T: Toroidal horizontal component coefficients (divergence-free part)

Arguments

• v_r: Radial velocity component [nlat, nlon]
• v_theta: Colatitude velocity component [nlat, nlon]
• v_phi: Azimuthal velocity component [nlat, nlon]
• config: SHTnsKit configuration

Returns

• (Q_coeffs, S_coeffs, T_coeffs): QST spectral coefficients

write_netcdf_boundary_data(filename::String, boundary_data::BoundaryData)

Write boundary condition data to a NetCDF file.

InitialConditions

Module for loading and generating initial conditions for geodynamo simulations. Supports loading from NetCDF files, generating random fields, and setting prescribed analytical patterns.

generate_random_composition!(comp_field, amplitude, modes_range)

Generate random composition initial conditions.

generate_random_initial_conditions!(field, field_type::Symbol;
                                   amplitude=1.0, modes_range=1:10,
                                   seed=nothing)

Generate random initial conditions for any field type.

Arguments

• field: Field structure to initialize
• field_type: Type of field (:temperature, :magnetic, :velocity, :composition)
• amplitude: Overall amplitude of random perturbations
• modes_range: Range of spherical harmonic modes to excite
• seed: Random seed for reproducibility (optional)

Examples

# Random temperature field
generate_random_initial_conditions!(temp_field, :temperature, amplitude=0.1)

# Random magnetic field with specific modes
generate_random_initial_conditions!(mag_field, :magnetic,
                                   amplitude=0.01, modes_range=1:20, seed=42)

generate_random_magnetic!(mag_field, amplitude, modes_range)

Generate random magnetic initial conditions.

generate_random_temperature!(temp_field, amplitude, modes_range)

Generate random temperature initial conditions.

generate_random_velocity!(vel_field, amplitude, modes_range)

Generate random velocity initial conditions.

load_composition_initial_conditions!(comp_field, file_path::String)

Load composition initial conditions from NetCDF file.

load_initial_conditions!(field, field_type::Symbol, file_path::String)

Load initial conditions from NetCDF file for any field type.

Arguments

• field: Field structure (temperature, magnetic, velocity, or composition type)
• field_type: Field type (:temperature, :magnetic, :velocity, :composition)
• file_path: Path to NetCDF file containing initial conditions

File Format

NetCDF files should contain:

• For scalar fields: spectral coefficients array
• For vector fields: toroidal and poloidal spectral coefficients
• Coordinate arrays: lm indices, radial grid

load_magnetic_initial_conditions!(mag_field, file_path::String)

Load magnetic initial conditions from NetCDF file.

load_temperature_initial_conditions!(temp_field, file_path::String)

Load temperature initial conditions from NetCDF file.

load_velocity_initial_conditions!(vel_field, file_path::String)

Load velocity initial conditions from NetCDF file.

randomize_magnetic_field!(field; amplitude, lmax, domain=nothing)

Populate magnetic toroidal/poloidal fields with random perturbations.

randomize_scalar_field!(field; amplitude, lmax, domain=nothing)

Populate a scalar spectral field (temperature/composition) with random perturbations up to degree lmax. If a radial domain is provided and includes r=0, ball regularity is enforced.

randomize_vector_field!(field; amplitude, lmax, domain=nothing)

Populate velocity-like toroidal/poloidal fields with random perturbations up to degree lmax.

save_initial_conditions(field, field_type::Symbol, file_path::String)

Save current field state as initial conditions to NetCDF file.

This function is useful for saving generated or computed initial conditions for later use in simulations.

set_analytical_composition!(comp_field, pattern, amplitude; parameters...)

Set analytical composition patterns.

set_analytical_initial_conditions!(field, field_type::Symbol, pattern::Symbol;
                                  amplitude=1.0, parameters...)

Set analytical initial conditions based on predefined patterns.

Patterns

• :conductive - Conductive temperature profile
• :dipole - Dipolar magnetic field
• :convective - Small convective velocity pattern
• :stratified - Stratified composition profile

Examples

# Conductive temperature profile
set_analytical_initial_conditions!(temp_field, :temperature, :conductive)

# Earth-like dipolar magnetic field
set_analytical_initial_conditions!(mag_field, :magnetic, :dipole, amplitude=1.0)

set_analytical_magnetic!(mag_field, pattern, amplitude; parameters...)

Set analytical magnetic field patterns.

set_analytical_temperature!(temp_field, pattern, amplitude; parameters...)

Set analytical temperature patterns.

set_analytical_velocity!(vel_field, pattern, amplitude; parameters...)

Set analytical velocity patterns.

apply_shell_composition_boundaries!(comp_field, boundary_set; time=0)

apply_shell_temperature_boundaries!(temp_field, boundary_set; time=0)

Wrapper around core boundary application for shell geometry.

create_shell_composition_field(T, cfg::ShellConfig; nr=GeoDynamo.i_N)

create_shell_hybrid_composition_boundaries(inner_spec, outer_spec, cfg::ShellConfig; precision=Float64)

Create composition boundaries for shell geometry. Dispatches to appropriate function based on spec types:

• Both Tuple: fully programmatic boundaries
• String + Tuple: hybrid (one from file, one programmatic)
• Both String: both from files

create_shell_hybrid_temperature_boundaries(inner_spec, outer_spec, cfg::ShellConfig; precision=Float64)

Create temperature boundaries for shell geometry. Dispatches to appropriate function based on spec types:

• Both Tuple: fully programmatic boundaries
• String + Tuple: hybrid (one from file, one programmatic)
• Both String: both from files

create_shell_magnetic_fields(T, cfg::ShellConfig; nr_oc=GeoDynamo.i_N, nr_ic=GeoDynamo.i_Nic)

create_shell_pencils(cfg::ShellConfig; optimize=true)

Create a shell-oriented pencil decomposition using the core topology helper.

create_shell_physical_field(T, cfg::ShellConfig, domain::GeoDynamo.RadialDomain, pencil)

create_shell_radial_domain(nr=i_N) -> RadialDomain

Create a radial domain suitable for a spherical shell using the global parameter d_rratio for inner radius ratio. This is a thin wrapper over GeoDynamo.create_radial_domain.

create_shell_spectral_field(T, cfg::ShellConfig, domain::GeoDynamo.RadialDomain, pencil)

create_shell_temperature_field(T, cfg::ShellConfig; nr=GeoDynamo.i_N)

create_shell_vector_field(T, cfg::ShellConfig, domain::GeoDynamo.RadialDomain, pencils)

create_shell_velocity_fields(T, cfg::ShellConfig; nr=GeoDynamo.i_N)

apply_ball_composition_regularity!(comp_field)

apply_ball_temperature_regularity!(temp_field)

Convenience to enforce scalar regularity on the temperature spectral field. Call after assembling or updating temp_field.spectral.

ball_physical_to_spectral!(phys::GeoDynamo.SHTnsPhysicalField,
                           spec::GeoDynamo.SHTnsSpectralField)

Wrapper for transforms in a solid sphere that enforces scalar regularity at r=0 after analysis. Use this for scalar fields (temperature, composition, etc.).

ball_vector_analysis!(vec::GeoDynamo.SHTnsVectorField,
                      tor::GeoDynamo.SHTnsSpectralField,
                      pol::GeoDynamo.SHTnsSpectralField)

Wrapper for vector analysis in a solid sphere; enforces vector regularity at r=0 after transforming to spectral toroidal/poloidal.

create_ball_composition_field(T, cfg::BallConfig; nr=GeoDynamo.i_N)

create_ball_hybrid_composition_boundaries(inner_spec, outer_spec, cfg::BallConfig; precision=Float64)

create_ball_hybrid_temperature_boundaries(inner_spec, outer_spec, cfg::BallConfig; precision=Float64)

create_ball_magnetic_fields(T, cfg::BallConfig; nr=GeoDynamo.i_N)

Create magnetic fields for a solid sphere. Since a "core" split is not used in a ball, we pass the same domain for both oc and ic to reuse the core implementation.

create_ball_pencils(cfg::BallConfig; optimize=true)

create_ball_physical_field(T, cfg::BallConfig, domain::GeoDynamo.RadialDomain, pencil)

create_ball_radial_domain(nr=i_N) -> RadialDomain

Create a radial domain for a solid sphere (inner radius = 0). Uses a cosine-stretched grid similar to the shell for compatibility, but sets the inner radius to zero and adjusts the coordinate columns to match expectations of downstream operators.

create_ball_spectral_field(T, cfg::BallConfig, domain::GeoDynamo.RadialDomain, pencil)

create_ball_temperature_field(T, cfg::BallConfig; nr=GeoDynamo.i_N)

create_ball_vector_field(T, cfg::BallConfig, domain::GeoDynamo.RadialDomain, pencils)

create_ball_velocity_fields(T, cfg::BallConfig; nr=GeoDynamo.i_N)

enforce_ball_scalar_regularity!(spec::GeoDynamo.SHTnsSpectralField)

Enforce scalar regularity at r=0 for solid sphere: for l>0, the scalar amplitude must vanish at r=0. Sets inner radial plane to zero for all nonzero l modes (both real and imaginary parts).

enforce_ball_vector_regularity!(tor_spec::GeoDynamo.SHTnsSpectralField,
                                pol_spec::GeoDynamo.SHTnsSpectralField)

Enforce vector-field regularity at r=0 for solid sphere. For smooth fields, both toroidal and poloidal potentials behave like r^{l+1}, so they vanish at r=0 for all l ≥ 1. Zeros the inner radial plane for l≥1.