API Reference

_BUFFER_CACHE_LOCK

Global ReentrantLock for thread-safe access to SHTnsKitConfig buffer caches. All access to config.buffercache should be protected by this lock.

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.

ERK2BoundarySide{T}

Encodes how to enforce a boundary condition at one boundary (inner or outer) for the ERK2 exponential integrator.

All BC types are expressed as a linear constraint: sumj stencil[j] * u[j] + lcorrection * u[b] = value

which is solved for u[b]: u[b] = (value - sum{j≠b} stencil[j] * u[j]) / (stencil[b] + lcorrection)

BC types and their encoding:

• :dirichlet → stencil = delta{b}, value = BCval (u[b] = BC_val)
• :neumann → stencil = d1[b,:], value = q (du/dr = q)
• :stressfreetor → stencil = d1[b,:], lcorrection = -1/rb (dT/dr - T/r = 0)
• :noslip_pol → stencil = d1[b,:], value = 0 (dP/dr = 0)
• :stressfreepol → stencil = d2[b,:], value = 0 (d²P/dr² = 0)
• :insulatinginner → stencil = d1[b,:], lcorrection = -l/r_b (d/dr - l/r)P = 0
• :insulatingouter → stencil = d1[b,:], lcorrection = +(l+1)/r_b (d/dr + (l+1)/r)P = 0

ERK2BoundarySpec{T}

Complete boundary condition specification for both boundaries.

inner_mode_values and outer_mode_values are optional per-mode value overrides (indexed by lmidx). When provided, they override `bcside.value` for specific modes. Used e.g. for rotating inner core: T(l=1,m=0) = rotomega * rinner.

ERK2Cache{T}

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

ERK2InfluenceMatrix{T}

Precomputed influence matrix data for enforcing P = 0 at boundaries while using derivative boundary conditions in the exponential operator.

This matches the Fortran DD_2DCODE influence matrix method for poloidal velocity.

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.

GradientWorkspace{T}

Shared workspace for transient gradient spectral fields. Reused across scalar field computations that run sequentially (temperature, composition). These arrays are zeroed and recomputed every timestep, so sharing is safe.

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)

SHTnsPhysField{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)

SHTnsSpecField{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 SHTnsPhysField, 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.

_copy_spectral_to_field!(field::SHTnsSpecField{T}, real_data::AbstractArray, imag_data::AbstractArray)

Copy restart spectral data into a SHTnsSpecField. Handles the dimension mismatch between the restart file format [nlm, nr] and the field storage format [nlm, 1, nr].

_get_influence_cache(domain::RadialDomain, ∂r)

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

_get_tau_cache(domain::RadialDomain)

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

_load_file_bcs!(state::SimulationState)

Load file-based boundary conditions from paths specified in parameters. Modifies the field's mutable boundaryinterpolationcache in place.

_load_restart_file(filepath, tracker, config; pencils=nothing)

Load simulation state from a specific restart NetCDF file path using parallel I/O. All ranks open the file collectively and read their local slices.

_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

_spectral_curl_torpol!(dst_tor_r, dst_tor_i, dst_pol_r, dst_pol_i,
                       src_tor_r, src_tor_i, src_pol_r, src_pol_i,
                       ℓ_factors, d1_matrix, d²_matrix, domain, config, T)

Shared spectral curl for toroidal-poloidal fields: (∇×V)tor = [l(l+1)/r² - d²/dr² - 2/r d/dr] Vpol (∇×V)pol = -l(l+1)/r² Vtor

Used by both current density (j = ∇×B) and induction curl (∇×(u×B)).

add_internal_sources_local!(𝔽::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_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,
                     𝔽::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, 𝔽::AbstractScalarField,
                               domain, r_range, owns_mode::Bool)

Apply flux boundary conditions using the influence matrix method.

MPI Safety

The owns_mode parameter indicates whether this process owns the lm mode. All processes must call this function for each mode to ensure Allreduce is called the same number of times by all processes (prevents deadlock).

apply_flux_bc_tau!(spec_real, spec_imag, local_lm, lm_idx,
                   apply_inner, apply_outer, 𝔽::AbstractScalarField,
                   domain, r_range, owns_mode::Bool)

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

MPI Safety

The owns_mode parameter indicates whether this process owns the lm mode. All processes must call this function for each mode to ensure Allreduce is called the same number of times by all processes (prevents deadlock).

apply_geometric_factors_spectral!(ws::GradientWorkspace{T}, 𝔽::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_influence_matrix_correction!(result, influence, bc_inner_val, bc_outer_val)

Apply influence matrix correction to enforce P = 0 at boundaries. This subtracts the Green's function response that would give non-zero boundary values.

After this correction, result[1] ≈ bcinnerval and result[nr] ≈ bcouterval (typically both zero for no-penetration).

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_scalar_flux_bc_spectral!(𝔽::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).

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call MPI collectives the same number of times, preventing deadlock with uneven lm distribution.

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_velocity_poloidal_influence_correction!(field, influence_matrices, config)

Apply influence matrix correction to enforce P = 0 at boundaries for velocity poloidal. This is called after erk2finalizefield! to enforce the no-penetration condition while preserving the derivative BC constraints used in the exponential operator.

Matches DD_2DCODE's influence matrix method for velocity poloidal.

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!(ℬ::SHTnsMagneticFields{T}) where T

Perform batched transforms for better cache efficiency using shtnskit_transforms.jl

batch_shtnskit_transforms!(specs::Vector{SHTnsSpecField{T}},
                          physs::Vector{SHTnsPhysField{T}}) where T

Batch process multiple transforms using SHTnsKit with PencilArrays.

MPI Safety

This function processes transforms sequentially to avoid calling MPI collectives (Barrier) from multiple threads simultaneously. MPI collective operations must be called from the same thread on all processes to avoid deadlock.

Note: The individual transforms themselves are still efficient as SHTnsKit performs optimized Legendre transforms and FFTs internally.

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!(𝒰::SHTnsVelocityFields{T}) where T

Perform batched transforms for better cache efficiency using shtnskit_transforms.jl

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_erk2_bc_spec_magnetic_pol(T, domain) -> ERK2BoundarySpec{T}

Build ERK2 boundary spec for magnetic poloidal field (insulating BCs, l-dependent). Inner: (d/dr - l/r)P = 0. Outer: (d/dr + (l+1)/r)P = 0.

build_erk2_bc_spec_magnetic_tor(T, nr) -> ERK2BoundarySpec{T}

Build ERK2 boundary spec for magnetic toroidal field (insulating: BT = 0 at both).

build_erk2_bc_spec_scalar(T, domain, i_bc) -> ERK2BoundarySpec{T}

Build ERK2 boundary spec for a scalar field (temperature or composition). i_bc encoding: 1=DD, 2=DN, 3=ND, 4=NN.

build_erk2_bc_spec_velocity_pol(T, domain, i_vel_bc) -> ERK2BoundarySpec{T}

Build ERK2 boundary spec for the velocity poloidal component.

Uses proper derivative boundary conditions matching DD_2DCODE:

• No-slip: ∂P/∂r = 0 (first derivative)
• Stress-free: ∂²P/∂r² = 0 (second derivative)

The no-penetration condition (P = 0 at boundaries) is enforced separately via the influence matrix method, which is applied after the exponential operator in erk2finalizevelocity_poloidal!.

ivelbc: 1=no-slip both, 2=no-slip inner/stress-free outer, 3=stress-free inner/no-slip outer, 4=stress-free both.

build_erk2_bc_spec_velocity_tor(T, domain, i_vel_bc; config=nothing, rot_omega=0.0) -> ERK2BoundarySpec{T}

Build ERK2 boundary spec for the velocity toroidal component. ivelbc: 1=no-slip both, 2=no-slip inner/stress-free outer, 3=stress-free inner/no-slip outer, 4=stress-free both.

If rot_omega ≠ 0 and inner BC is no-slip, the inner boundary value for (l=1, m=0) is set to rot_omega * r_inner (rotating inner core).

build_influence_matrix(G_inner, G_outer, ∂r, 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_lm_lookup_tables(lmax, mmax) -> (l_values, m_values)

Build precomputed lookup tables for converting linear indices to (l, m).

Returns

• l_values: Vector where l_values[idx] gives the degree l for linear index idx
• m_values: Vector where m_values[idx] gives the order m for linear index idx

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.

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

build_∂θ(::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_parallel_netcdf_support(comm)

Verify that parallel NetCDF (MPI-IO via HDF5) is available at runtime. Called once at initialization. Errors immediately if not supported.

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::SHTnsSpecField, field_name::String,
                              config::NaNDetectionConfig, step::Int)

Check both real and imaginary parts of a spectral field.

check_transform_synchronization(config; strict::Bool=false) -> Bool

Verify that the SHTnsKit transform configuration is safe for parallel execution.

Checks:

1. All processes have same number of local radial levels
2. Spectral mode distribution is valid
3. FFT plans are properly initialized

Arguments

• config: SHTnsKit configuration object with pencil information
• strict: If true, throw an error on invalid configuration (recommended for production)

Returns

true if configuration is safe for parallel transforms.

Example

# Validate before running transforms
check_transform_synchronization(config; strict=true)

clear_buffer_cache!(config)

Thread-safe clearing of all cached buffers. Useful when changing configurations or to free memory.

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!(𝔽::AbstractScalarField{T}, domain::RadialDomain, ws::GradientWorkspace{T}) 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!(𝒰, 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 (magnetic diffusion time scaling)

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

where:

• τ = L²/η (magnetic diffusion time scaling)
• E = ν/(2ΩL²) is the Ekman number
• r factor in buoyancy represents linear gravity profile g(r) ∝ r

Implementation Notes

The explicit RHS entering the time integrator is: RHS = -(E/Pm)·(∇×ũ)×ũ - (ẑ×ũ) + (Pm/Pr)·Ra·T̃·r·r̂ + (Pm/Sc)·Ra_C·C̃·r·r̂ + (∇×B̃)×B̃

Coefficients (Fortran DD_2DCODE convention, mass coeff = E on du/dt):

• Advection: E (= d_Ro = Rossby number, equals E for magnetic diffusion time)
• Coriolis: 1 (no scaling)
• Thermal buoyancy: (Pm/Pr) * Ra * r (with radial factor)
• Compositional buoyancy: (Pm/Sc) * Ra_C * r (with radial factor)
• Lorentz: 1/Pm

Viscous diffusion is treated implicitly with coefficient E (Ekman number). Mass coefficient E is applied in the time-stepping matrices.

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

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

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

Compute Chebyshev polynomial T_n on the radial grid.

compute_composition_nonlinear!(𝔽::SHTnsCompositionField{T},
                               vel_fields, 𝒟ᵒᶜ::RadialDomain;
                               geometry::Symbol = get_parameters().geometry) where T

Compute the nonlinear advection term -u·∇C for the composition equation.

This function implements the EXPLICIT part of the composition equation: ∂C/∂t + u·∇C = (Pm/Sc) ∇²C

The advection term -u·∇C is computed using the following optimized workflow:

1. Compute ∇C in spectral space (no MPI communication required)
2. Batched transform of C and ∇C to physical space
3. Compute advection -u·∇C locally in physical space
4. Transform result back to spectral space
5. Apply boundary conditions

Arguments

• 𝔽: Composition field structure to update
• vel_fields: Velocity field structure (provides u for advection)
• 𝒟ᵒᶜ: Radial domain information
• geometry: Geometry type (:shell or :ball) for transform selection

Side Effects

• Updates 𝔽.nonlinear with the computed advection term
• Updates 𝔽.gradient with ∇C in physical space
• Applies boundary conditions to spectral representation

Performance Notes

• Gradients computed in spectral space avoid MPI communication
• Batched transforms minimize data transfer overhead
• Physical space operations are embarrassingly parallel

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

Check solenoidal constraint for toroidal-poloidal decomposition. For v = ∇×(T r̂) + ∇×∇×(P r̂), the divergence ∇·v = 0 is satisfied analytically by construction (divergence of a curl is identically zero). Therefore this function always returns (0.0, 0.0) — the solenoidal constraint is guaranteed by the representation, not by numerical accident.

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(𝒟ᵒᶜ::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(ℬ::SHTnsMagneticFields{T}) where T

Compute magnetic helicity using enhanced spectral integration

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!(𝔽::AbstractScalarField{T}, ws::GradientWorkspace{T}) where T

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

compute_radial_gradient_spectral!(𝔽::AbstractScalarField{T}, domain::RadialDomain, ws::GradientWorkspace{T}) where T

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

Note: Radial data is always local (not MPI distributed). Only spectral (lm) modes are distributed across MPI processes.

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

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

compute_scalar_energy(𝔽::AbstractScalarField{T}, 𝒟ᵒᶜ::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(𝔽::AbstractScalarField{T}, 𝒟ᵒᶜ::RadialDomain) where T

Compute RMS value of scalar field.

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_theta_gradient_spectral!(𝔽::AbstractScalarField{T}, ws::GradientWorkspace{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_composition_matrices(config, domain, diffusivity, dt;
                             i_cmp_bc, theta, T)

Create implicit time-stepping matrices for the composition equation with boundary conditions embedded in the matrix rows (matching Fortran cmpbcC).

Arguments

• config: SHTnsKitConfig with l_values
• domain: RadialDomain with derivative matrices and grid
• diffusivity: Compositional diffusivity (Pm/Sc in magnetic diffusion time)
• dt: Timestep
• i_cmp_bc: BC type (1=DD, 2=DN, 3=ND, 4=NN)
• theta: Implicit parameter (default d_implicit)
• T: Numeric type (default Float64)

Returns

• SHTnsImplicitMatrices{T} with BCs embedded in system matrices

create_computation_pencils(topology, dims, config)

Create specialized pencils for different stages of computation.

create_dirichlet_bc(T, nr, value) -> ERK2BoundarySide{T}

Create a Dirichlet BC side (u = value at boundary).

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, bc_spec=nothing)

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

Boundary rows of A are zeroed only for l=0 (where the spherical Laplacian has a null space, making A singular). For l≥1, the full operator is retained (non-singular due to the -l(l+1)/r² term), enabling accurate LU-based phi1/phi2 computation and O(h³) enforcement corrections per step.

create_erk2_cache_composition(T, config, domain, diffusivity, dt, i_cmp_bc; kwargs...)

Create ERK2 cache for composition field with embedded BCs. Wrapper around createerk2cache_scalar with composition-specific defaults.

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

Create ERK2 cache for magnetic poloidal field with insulating BCs embedded in the matrix A.

This matches DD_2DCODE's approach where the matrix has:

• Inner boundary row: (∂/∂r - l/r) operator → insulating interior
• Outer boundary row: (∂/∂r + (l+1)/r) operator → insulating exterior

Embedding BCs in A ensures exp(dt*A) automatically preserves solutions satisfying the insulating constraints, rather than requiring post-evolution correction.

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

Create ERK2 cache for magnetic toroidal field with insulating BCs embedded in the matrix A.

This matches DD_2DCODE's approach where the matrix has:

• Inner boundary row: identity → BT = 0
• Outer boundary row: identity → BT = 0

For Dirichlet BCs (BT = 0), we zero the boundary rows except for the diagonal (identity), which makes exp(dt*A)|_boundary = identity, preserving BT = 0.

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

Create ERK2 cache for scalar fields (temperature or composition) with boundary rows zeroed in the matrix A.

Boundary condition types (matching DD2DCODE tmpbcT / cmpbc_C):

• i_bc = 1: Dirichlet-Dirichlet (fixed value both boundaries)
• i_bc = 2: Dirichlet-Neumann (fixed value inner, fixed flux outer)
• i_bc = 3: Neumann-Dirichlet (fixed flux inner, fixed value outer)
• i_bc = 4: Neumann-Neumann (fixed flux both, l=0 inner uses Dirichlet)

For exponential integrators, boundary rows are always zeroed:

• exp(dt*A)[boundary,:] = [1, 0, ..., 0] → preserves boundary value
• Actual BC enforcement (Dirichlet/Neumann) is done via post-processing in erk2preparefield! and erk2finalizefield! using enforceerk2bc!

create_erk2_cache_temperature(T, config, domain, diffusivity, dt, i_tmp_bc; kwargs...)

Create ERK2 cache for temperature field with embedded BCs. Wrapper around createerk2cache_scalar with temperature-specific defaults.

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_insulating_inner_bc(T, d1_row, r_inv) -> ERK2BoundarySide{T}

Create an insulating magnetic poloidal inner BC: (d/dr - l/r)P = 0. The l-factor is applied dynamically during enforcement.

create_insulating_outer_bc(T, d1_row, r_inv) -> ERK2BoundarySide{T}

Create an insulating magnetic poloidal outer BC: (d/dr + (l+1)/r)P = 0. The (l+1) factor is split: fixedcorrection = +1/r (the "+1" part), and lsign = +1 with uselcorrection = true (the "l" part). Total correction = lsignlrinv + fixed_correction = (l+1)/r.

create_magnetic_poloidal_matrices(config, domain, diffusivity, dt;
                                  theta, T)

Create implicit time-stepping matrices for the poloidal magnetic component with insulating boundary conditions embedded in the matrix rows (matching Fortran magbcPol).

Boundary conditions (l-dependent):

• Inner: (∂/∂r - l/r) BP = 0 (field matches r^l interior solution)
• Outer: (∂/∂r + (l+1)/r) BP = 0 (field matches r^{-(l+1)} exterior decay)

create_magnetic_toroidal_matrices(config, domain, diffusivity, dt;
                                  theta, T)

Create implicit time-stepping matrices for the toroidal magnetic component with insulating boundary conditions embedded in the matrix rows (matching Fortran magbcTor).

Both boundary rows use identity (BT = 0) for insulating exterior/interior.

create_neumann_bc(T, d1_row, value) -> ERK2BoundarySide{T}

Create a Neumann BC side (du/dr = value at boundary).

create_noslip_pol_bc(T, d1_row) -> ERK2BoundarySide{T}

Create a no-slip poloidal BC side: dP/dr = 0.

create_parallel_netcdf(filename, config, field_info, metadata, comm)

Open a new NetCDF file in parallel mode. All ranks call this collectively. Defines global dimensions and variables based on field_info.

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; nr=i_N, 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_radial_domain(nr=nothing; rratio=nothing, kl=nothing) -> RadialDomain

Create a radial domain for a spherical shell geometry.

Arguments

• nr: Number of radial points (defaults to get_parameters().i_N)
• rratio: Inner/outer radius ratio (defaults to get_parameters().d_rratio)
• kl: Bandwidth for finite differences (defaults to get_parameters().i_KL)

Returns

A RadialDomain with Chebyshev-like radial spacing optimized for spectral accuracy.

create_shtns_composition_field(::Type{T}, config::SHTnsKitConfig,
                               𝒟ᵒᶜ::RadialDomain,
                               pencils=nothing, pencil_spec=nothing) where T

Create a composition field structure with all necessary arrays for spectral computation.

Arguments

• T: Numeric type (e.g., Float64)
• config: SHTnsKitConfig with transform parameters (lmax, mmax, nlat, nlon)
• 𝒟ᵒᶜ: Radial domain information for outer core
• pencils: Optional pencil decomposition (defaults to config.pencils)
• pencil_spec: Optional spectral pencil (defaults to pencils.spec)

Returns

• SHTnsCompositionField{T}: Fully initialized composition field structure

Example

config = create_shtns_config(lmax=32, mmax=32, nlat=64, nlon=128)
domain = create_radial_domain(nr=64, ri=0.35, ro=1.0)
𝔽 = create_shtns_composition_field(Float64, config, domain)

Notes

Default boundary conditions are no-flux (NEUMANN) at both boundaries, appropriate for typical compositional convection problems.

create_shtns_physical_field(T, config, 𝒟ᵒᶜ, pencil) -> SHTnsPhysField{T}

Create a new physical space field initialized to zero.

Arguments

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

Returns

A new SHTnsPhysField with all grid values initialized to zero.

create_shtns_spectral_field(T, config, 𝒟ᵒᶜ, pencil_spec) -> SHTnsSpecField{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
• 𝒟ᵒᶜ: RadialDomain specifying the radial discretization
• pencil_spec: PencilArrays Pencil defining the data distribution

Returns

A new SHTnsSpecField with:

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

create_shtns_vector_field(T, config, 𝒟ᵒᶜ, pencils) -> SHTnsVectorField{T}

Create a new vector field with three physical space components.

Arguments

• T: Element type
• config: SHTnsKit configuration
• 𝒟ᵒᶜ: 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_stress_free_pol_bc(T, d2_row) -> ERK2BoundarySide{T}

Create a stress-free poloidal BC side: d²P/dr² = 0.

create_stress_free_tor_bc(T, d1_row, r_inv) -> ERK2BoundarySide{T}

Create a stress-free toroidal BC side: dT/dr - T/r = 0.

create_temperature_matrices(config, domain, diffusivity, dt;
                             i_tmp_bc, theta, T)

Create implicit time-stepping matrices for the temperature equation with boundary conditions embedded in the matrix rows (matching Fortran tmpbcT).

Arguments

• config: SHTnsKitConfig with l_values
• domain: RadialDomain with derivative matrices and grid
• diffusivity: Thermal diffusivity (Pm/Pr in magnetic diffusion time)
• dt: Timestep
• i_tmp_bc: BC type (1=DD, 2=DN, 3=ND, 4=NN)
• theta: Implicit parameter (default d_implicit)
• T: Numeric type (default Float64)

Returns

• SHTnsImplicitMatrices{T} with BCs embedded in system matrices

create_transpose_plans(pencils)

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

create_velocity_green_matrices(config, domain, diffusivity;
                                theta, T)

Create Green's function matrices for the influence matrix method (poloidal pressure). These use Dirichlet BCs (identity rows) at both boundaries, matching Fortran velbcGre.

create_velocity_poloidal_influence_matrices(T, config, domain, diffusivity, dt, i_vel_bc; theta)

Create influence matrices for poloidal velocity to enforce P = 0 at boundaries while using derivative BCs (∂P/∂r = 0 or ∂²P/∂r² = 0) in the implicit/exponential system.

This matches DD2DCODE's velmatrices routine.

create_velocity_poloidal_matrices(config, domain, diffusivity, dt;
                                  theta, i_vel_bc, T)

Create implicit time-stepping matrices for the poloidal velocity component with boundary conditions embedded in the matrix rows.

The boundary rows of the system matrix are replaced with the BC equations:

• No-slip: first derivative row (∂P/∂r = value)
• Stress-free: second derivative row (∂²P/∂r² = value)

create_velocity_toroidal_matrices(config, domain, diffusivity, dt;
                                  theta, i_vel_bc, T)

Create implicit time-stepping matrices for the toroidal velocity component with boundary conditions embedded in the matrix rows.

The boundary rows of the system matrix are replaced with the BC equations:

• No-slip: identity row (T = value)
• Stress-free: ∂T/∂r - T/r = 0

This ensures the implicit solve enforces BCs exactly rather than applying them as post-processing.

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}).

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

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.

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

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

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

Mass Coefficient

For equations of the form c * du/dt = ν*L*u + NL (e.g. velocity with c=dE), pass `masscoeff=c. The operator becomesA = (ν/c)*Land NL is scaled by1/c`.

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

enforce_erk2_bc!(result, bc_side, boundary_idx, l, nr)

Enforce a boundary condition on the result vector at the given boundary index. Modifies result[boundary_idx] to satisfy the linear constraint encoded in bc_side.

enforce_velocity_boundary_values!(𝒰)

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

erk2_finalize_field!(buffers, u, cache, config, dt)

Finalize ERK2 second stage by applying phi2 correction.

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

erk2_prepare_field!(buffers, u, nl, cache, config, dt)

Prepare ERK2 first stage by computing linear evolution and k1 terms.

MPI Safety

Uses global loop bounds (1:nlm) to ensure all processes call Allreduce the same number of times, preventing deadlock with uneven lm distribution.

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, nr, lmax, field_count, T)

Estimate memory usage for SHTnsKit-based transforms with PencilArrays.

Arguments

• nlat: Number of latitude points
• nlon: Number of longitude points
• nr: Number of radial points
• lmax: Maximum spherical harmonic degree
• field_count: Number of fields to estimate for
• T: Element type (e.g., Float64)

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_all_fields(state::SimulationState)

Extract all field data from the simulation state for output/restart writing. Returns a Dict with spectral field data (real/imag parts) for velocity, magnetic, temperature, and composition fields.

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_coefficients_pair_for_shtnskit(spec1_real, spec1_imag, spec2_real, spec2_imag, r_local, config)

Extract two spectral coefficient matrices efficiently for vector transforms.

This is optimized for the common case in vector synthesis/analysis where we need both toroidal and poloidal coefficients. It avoids one copy operation compared to calling extract_coefficients_for_shtnskit twice.

Returns

Tuple (coeffs1, coeffs2) of coefficient matrices.

extract_dense_row(banded::BandedMatrix{T}, row::Int) -> Vector{T}

Extract a full dense row from a banded matrix.

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.

WARNING: MPI Synchronization

This function contains MPI.Allreduce! which is a collective operation. When called inside a per-radial loop, ALL MPI processes must call this function the same number of times, otherwise deadlock will occur. Ensure even radial distribution or use global loop bounds.

extract_physical_slice_phi_local!(slice_buffer, phys_data, r_local, config)

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

WARNING: MPI Synchronization

This function contains MPI.Allreduce! which is a collective operation. When called inside a per-radial loop, ALL MPI processes must call this function the same number of times, otherwise deadlock will occur. Ensure even radial distribution or use global loop bounds.

extract_spectral_subset(field::SHTnsSpecField{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.

WARNING: MPI Synchronization

This function contains MPI.Allreduce! which is a collective operation. When called inside a per-radial loop, ALL MPI processes must call this function the same number of times, otherwise deadlock will occur. Ensure even radial distribution or use global loop bounds.

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_cached_buffer!(create_func::Function, config, key::Symbol)
get_cached_buffer!(config, key::Symbol) do ... end

Thread-safe accessor for buffer cache. Returns existing buffer if present, otherwise creates a new one using create_func() and caches it.

Note: The function parameter comes FIRST to support Julia's do block syntax. When using do block, Julia desugars it to pass the closure as the first argument.

Arguments

• create_func::Function: Zero-argument function to create buffer if not cached
• config: SHTnsKitConfig object containing the buffer cache
• key::Symbol: Key to look up in the buffer cache

Returns

The cached or newly created buffer.

Example

buffer = get_cached_buffer!(config, :my_buffer) do
    zeros(Float64, nlat, nlon)
end

get_comm()

Get MPI communicator, initializing MPI if needed. Thread-safe via double-checked locking with _MPI_INIT_LOCK.

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_erk2_composition_cache!(caches, diffusivity, T, config, domain, dt, i_cmp_bc; use_krylov=false, m=20, tol=1e-8)

Retrieve or create ERK2 cache for composition field with embedded BCs. Uses BC type specified by icmpbc (1=DD, 2=DN, 3=ND, 4=NN), matching DD2DCODE's cmpbc_C.

get_erk2_influence_matrices!(cache, key, T, config, domain, diffusivity, dt, i_vel_bc; theta)

Get or create ERK2 influence matrices for poloidal velocity. Caches matrices by key for reuse across timesteps. Invalidates and recomputes when dt changes.

get_erk2_magnetic_poloidal_cache!(caches, diffusivity, T, config, domain, dt; use_krylov=false, m=20, tol=1e-8)

Retrieve or create ERK2 cache for magnetic poloidal field with embedded insulating BCs. Uses l-dependent insulating BCs matching DD2DCODE's magbc_Pol:

• Inner boundary: (∂/∂r - l/r)P = 0
• Outer boundary: (∂/∂r + (l+1)/r)P = 0

get_erk2_magnetic_toroidal_cache!(caches, diffusivity, T, config, domain, dt; use_krylov=false, m=20, tol=1e-8)

Retrieve or create ERK2 cache for magnetic toroidal field with embedded insulating BCs. Uses Dirichlet BCs (BT = 0) at both boundaries, matching DD2DCODE's magbc_Tor.

get_erk2_temperature_cache!(caches, diffusivity, T, config, domain, dt, i_tmp_bc; use_krylov=false, m=20, tol=1e-8)

Retrieve or create ERK2 cache for temperature field with embedded BCs. Uses BC type specified by itmpbc (1=DD, 2=DN, 3=ND, 4=NN), matching DD2DCODE's tmpbc_T.

get_flux_value(lm_idx::Int, boundary::Int, 𝔽::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.

index_to_lm_fast(idx, config) -> (l, m)

Fast O(1) conversion from linear spectral index to (l, m) using precomputed tables.

Uses the lvalues and mvalues arrays stored in SHTnsKitConfig during initialization. This is significantly faster than index_to_lm_shtnskit for repeated lookups.

Arguments

• idx: Linear spectral index (1-based)
• config: SHTnsKitConfig containing precomputed lvalues and mvalues

Returns

Tuple (l, m) for the spherical harmonic degree and order.

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 SHTnsKit m-major convention (m varies slowest, l fastest):

• idx=1: (l=0, m=0)
• idx=2: (l=1, m=0)
• idx=3: (l=2, m=0)
• ...
• idx=lmax+1: (l=lmax, m=0)
• idx=lmax+2: (l=1, m=1)
• idx=lmax+3: (l=2, m=1)
• ...

Performance Note

This function uses a linear search. For performance-critical code with many lookups, use index_to_lm_fast with precomputed lookup tables from 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_parameters(config_file::String = "")

Initialize the global parameter system.

initialize_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::SHTnsPhysField{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::SHTnsSpecField{T}, nc_file, 
                           real_var::String, imag_var::String) where T

Load a single spectral field from NetCDF file.

load_spectral_coefficients!(field::SHTnsSpecField{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,
                      ∂r, 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,
                      ∂r, 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_magnetic_memory_layout!(ℬ::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!(𝒰::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.

Pivot singularity threshold: abs(pivot) < eps(T) * PIVOT_SINGULARITY_FACTOR

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

Print the GeoDynamo startup banner with configuration information.

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, nr, lmax, mmax, nlm, nprocs, memory_estimate)

Print configuration summary for SHTnsKit setup.

print_transpose_statistics()

Print accumulated transpose timing statistics.

Reciprocal condition number threshold for matrix function fallback

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.

restore_fields_from_restart!(state::SimulationState{T}, restart_data::Dict{String,Any})

Restore simulation state fields from restart data loaded by read_restart!. Copies spectral coefficients from the restart data Dict into the state's field arrays.

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_simulation!(state::SimulationState{T}; restart_file::String="", restart_dir::String="", restart_time::Float64=0.0)

Run geodynamo simulation with all parallelization optimizations. Supports restarting from a saved NetCDF restart file.

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

Evaluate a parsed expression using only safe arithmetic operations and known parameter values. No arbitrary code execution is possible.

safe_parse_value(value_str::AbstractString, param_dict::Dict{Symbol, Any})

Parse a parameter value string without using eval. Supports:

• Integer and float literals (including scientific notation)
• Boolean literals (true, false)
• Symbol literals (:name)
• String literals ("...")
• Mathematical constants (π)
• Simple arithmetic expressions referencing previously-defined parameters

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.

Series term convergence: opnorm(term) < eps(T) * SERIES_CONVERGENCE_FACTOR

set_composition_rhs_bc!(rhs_real, rhs_imag, local_lm, nr;
                         inner_value=0.0, outer_value=0.0,
                         inner_value_imag=0.0, outer_value_imag=0.0)

Set boundary values in the RHS vector for the composition solve. Matches Fortran cmp_setbc: sets RHS boundary rows to prescribed values.

For standard compositional convection, boundary values are typically zero (homogeneous Neumann: ∂C/∂r=0, or homogeneous Dirichlet: C=0). Non-zero imaginary parts are used when file-based spectral BCs are loaded.

set_magnetic_rhs_bc!(rhs_real, rhs_imag, local_lm, nr;
                      inner_value=0.0, outer_value=0.0)

Set boundary values in the RHS vector for the magnetic field solve. Matches Fortran magsetbcTor / tim_zerobc: sets RHS boundary rows to zero.

For insulating BCs, all boundary values are zero (homogeneous conditions).

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_temperature_rhs_bc!(rhs_real, rhs_imag, local_lm, nr;
                         inner_value=0.0, outer_value=0.0,
                         inner_value_imag=0.0, outer_value_imag=0.0)

Set boundary values in the RHS vector for the temperature solve. Matches Fortran tmp_setbc: sets RHS boundary rows to prescribed values.

For standard geodynamo simulations, boundary values are typically zero (homogeneous Dirichlet: T=0, or homogeneous Neumann: ∂T/∂r=0). Non-zero imaginary parts are used when file-based spectral BCs are loaded.

set_velocity_rhs_bc_poloidal!(rhs_real, rhs_imag, local_lm, nr;
                               inner_value=0.0, outer_value=0.0)

Set boundary values in the RHS vector for the poloidal velocity solve. Matches Fortran: sets RHS boundary rows to zero for homogeneous BCs.

For no-slip: RHS boundary = 0 (∂P/∂r = 0) For stress-free: RHS boundary = 0 (∂²P/∂r² = 0)

set_velocity_rhs_bc_toroidal!(rhs_real, rhs_imag, local_lm, nr;
                               inner_value=0.0, outer_value=0.0)

Set boundary values in the RHS vector for the toroidal velocity solve. Matches Fortran velsetbcTor: sets RHS boundary rows to zero (or prescribed value).

For no-slip: RHS boundary = 0 (or prescribed velocity, e.g., rotating inner core) For stress-free: RHS boundary = 0 (homogeneous condition ∂T/∂r - T/r = 0)

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.

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::SHTnsPhysField: Source physical field values
• spec::SHTnsSpecField: 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::SHTnsSpecField: Source spectral field with coefficients
• phys::SHTnsPhysField: 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::SHTnsSpecField{T},
                         pol_spec::SHTnsSpecField{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::SHTnsSpecField{T},
                          pol_spec::SHTnsSpecField{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).

solve_banded!(x, lu, b)

Solve the banded linear system lu * x = b in-place.

In-place aliasing x === b is explicitly supported and safe. The forward substitution sweep processes rows in ascending order: at step i, it reads x[j] for j < i (already holding the intermediate result y[j]) and b[i] (not yet overwritten when aliased). The back substitution sweep processes rows in descending order with the same safe access pattern.

WARNING: Do not change the loop ordering without verifying aliasing safety.

solve_composition_implicit_step!(solution, rhs, matrices;
                                  bc_inner, bc_outer, bc_inner_imag, bc_outer_imag)

Solve the implicit composition system with boundary conditions embedded in the matrix. Before solving, sets the RHS boundary values appropriately.

This matches the Fortran DD_2DCODE approach where each rank has full radial profiles for its subset of (l,m) modes:

1. Loop over local lm modes
2. Set RHS boundary rows to prescribed values (cmp_setbc)
3. Solve banded system with BC rows in matrix
4. Solution automatically satisfies BCs

Arguments

• solution: Output spectral field
• rhs: Input RHS spectral field (modified in place for BCs)
• matrices: SHTnsImplicitMatrices with composition BCs embedded
• bc_inner: Optional vector of inner boundary values (real part) per mode (default: zeros)
• bc_outer: Optional vector of outer boundary values (real part) per mode (default: zeros)
• bc_inner_imag: Optional vector of inner boundary values (imag part) per mode (default: zeros)
• bc_outer_imag: Optional vector of outer boundary values (imag part) per mode (default: zeros)

solve_magnetic_implicit_step!(solution, rhs, matrices, component;
                               mag_bc_inner=nothing, prev_bc_inner=nothing)

Solve the implicit magnetic system with boundary conditions embedded in the matrix. Before solving, sets the RHS boundary values appropriately.

This matches the Fortran DD_2DCODE approach where each rank has full radial profiles for its subset of (l,m) modes:

1. Loop over local lm modes
2. Set RHS boundary rows to 0 (insulating BCs)
3. Solve banded system with BC rows in matrix
4. Solution automatically satisfies BCs

Arguments

• solution: Output spectral field
• rhs: Input RHS spectral field (modified in place for BCs)
• matrices: SHTnsImplicitMatrices with BCs embedded
• component: :toroidal or :poloidal
• mag_bc_inner: Optional inner boundary values for toroidal (conducting IC case)
• prev_bc_inner: Previous step inner BC values (for incremental form)

solve_temperature_implicit_step!(solution, rhs, matrices;
                                  bc_inner, bc_outer, bc_inner_imag, bc_outer_imag)

Solve the implicit temperature system with boundary conditions embedded in the matrix. Before solving, sets the RHS boundary values appropriately.

This matches the Fortran DD_2DCODE approach where each rank has full radial profiles for its subset of (l,m) modes:

1. Loop over local lm modes
2. Set RHS boundary rows to prescribed values (tmp_setbc)
3. Solve banded system with BC rows in matrix
4. Solution automatically satisfies BCs

Arguments

• solution: Output spectral field
• rhs: Input RHS spectral field (modified in place for BCs)
• matrices: SHTnsImplicitMatrices with temperature BCs embedded
• bc_inner: Optional vector of inner boundary values (real part) per mode (default: zeros)
• bc_outer: Optional vector of outer boundary values (real part) per mode (default: zeros)
• bc_inner_imag: Optional vector of inner boundary values (imag part) per mode (default: zeros)
• bc_outer_imag: Optional vector of outer boundary values (imag part) per mode (default: zeros)

solve_velocity_implicit_step!(solution, rhs, matrices, component;
                               i_vel_bc=1, domain=nothing,
                               rot_omega=0.0, current_field=nothing)

Solve the implicit velocity system with boundary conditions embedded in the matrix. Before solving, sets the RHS boundary values appropriately.

This matches the Fortran DD_2DCODE approach where each rank has full radial profiles for its subset of (l,m) modes:

1. Loop over local lm modes
2. Set RHS boundary rows to BC values (typically 0)
3. Solve banded system with BC rows in matrix
4. Solution automatically satisfies BCs

Arguments

• solution: Output spectral field
• rhs: Input RHS spectral field (modified in place for BCs)
• matrices: SHTnsImplicitMatrices with BCs embedded
• component: :toroidal or :poloidal
• i_vel_bc: Velocity BC type (1-4)
• domain: RadialDomain (needed for rotating IC boundary)
• rot_omega: Inner core rotation rate (for no-slip toroidal l=1,m=0)
• current_field: Current velocity field (for incremental form of rotating IC BC)

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::SHTnsSpecField{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!(𝔽::AbstractScalarField{T}, ws::GradientWorkspace{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, label=:default)

Perform transpose with optional timing and statistics.

Arguments

• dest::PencilArray: Destination array
• src::PencilArray: Source array
• label::Symbol: Label for timing statistics (default: :default)

Notes

Uses PencilArrays.transpose! which creates a Transposition internally. Enable timing with ENABLE_TIMING[] = true.

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.

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

Validate magnetic field configuration consistency with SHTns setup

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

Check that spectral field data is valid (no NaN/Inf) across all MPI processes. Returns the field after validation. Warns if any process has invalid data.

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_radial_distribution(pencils; warn_uneven::Bool=true, strict::Bool=true) -> Bool

Validate that radial dimension has compatible distribution across all pencils.

MPI Synchronization Requirement

The SHTnsKit transforms use MPI.Allreduce inside per-radial-level loops. All processes must have the SAME number of local radial levels, otherwise processes will enter/exit the loop at different times causing MPI DEADLOCK.

Arguments

• pencils: Named tuple of pencil configurations
• warn_uneven: If true, emit warning for uneven distribution
• strict: If true (default), throw an error instead of just warning

Returns

true if distribution is valid (all processes have same local radial count). false if there's a potential synchronization issue.

Default Behavior

Strict mode is enabled by default to prevent MPI deadlock in production runs. Set strict=false only for debugging purposes.

Example

# Recommended: use default strict mode
validate_radial_distribution(pencils)

# For debugging only: disable strict mode
validate_radial_distribution(pencils; strict=false)

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

Validate velocity field configuration consistency with SHTns setup

verify_all_ranks_wrote(output_dir, hist_number; geometry="shell", expected_dims=nothing)

Verify that the single parallel output file exists for a given history number and has correct global dimensions.

write_coordinate_data!(ds, field_info, config)

Write coordinate arrays. Only rank 0 writes coordinates (they are global/shared).

write_field_data!(ds, fields, config, field_info)

Write field data using parallel offset writes. Each rank writes its local pencil slice at the correct global position.

write_grid_file!(config, field_info, shtns_config, metadata)

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

zero_scalar_work_arrays!(𝔽::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.

ProgrammaticBoundarySet{T}

Wrapper around BoundaryConditionSet that also stores BC types for each boundary. Returned by create_programmatic_temperature_boundaries and create_programmatic_composition_boundaries.

SpectralBoundaryCoefficients{T}

Stores spectral boundary condition coefficients for inner and outer boundaries. Row 1 = inner boundary, Row 2 = outer boundary (matching Fortran convention).

Fields

• bc_real::Matrix{T}: Real part of spectral BCs [2, nlm]
• bc_imag::Matrix{T}: Imaginary part of spectral BCs [2, nlm]
• nlm::Int: Number of spectral modes
• source_file::String: Path to the source file
• source_format::Symbol: Format used to load (:physical or :spectral)

Ylm(l::Int, m::Int)

Spherical harmonic pattern specifier for programmatic boundary conditions.

Examples

# Create Y₂¹ pattern with amplitude 0.5
boundary = create_programmatic_boundary(Ylm(2, 1), config, 0.5)

# Create Y₄⁻² pattern
boundary = create_programmatic_boundary(Ylm(4, -2), config, 1.0)

# Time-dependent oscillating Y₃² pattern
boundary = create_time_dependent_programmatic_boundary(Ylm(3, 2), config, (0.0, 10.0), 100;
    amplitude=0.3)

_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. Thread-safe implementation using a lock for the check-then-set pattern.

Load physical-space boundary data from NetCDF, transform to spectral via SHTnsKit.

Load pre-computed spectral boundary coefficients from NetCDF (Fortran-compatible format).

_parse_boundary_spec(spec::Tuple, cfg) -> (values::Matrix, bc_type::BoundaryType)

Parse a boundary specification tuple into physical-space values and BC type.

Supported spec formats:

• (:uniform, value) or (:dirichlet, value): Dirichlet BC with uniform value
• (:neumann, value): Neumann BC with uniform flux value

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

Add noise to existing boundary data.

apply_boundary_conditions!(𝔽, 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, pbs::ProgrammaticBoundarySet; time=0.0)

Apply programmatic composition boundary conditions to a composition field. Same interface as apply_temperature_boundaries!.

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_temperature_bc_to_rhs!(rhs, temp_field)

Apply temperature boundary conditions to RHS vector.

apply_temperature_boundaries!(temp_field, pbs::ProgrammaticBoundarySet; time=0.0)

Apply programmatic temperature boundary conditions to a temperature field.

Transforms physical-space boundary data to spectral coefficients and stores them in the field's boundary_values matrix, along with setting BC types.

Arguments

• temp_field: SHTnsTemperatureField to apply boundaries to
• pbs: ProgrammaticBoundarySet containing boundary data and BC types
• time: Current simulation time (for time-dependent boundaries)

apply_velocity_bc_to_rhs!(rhs, velocity_field)

Apply velocity boundary conditions to RHS vector.

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.

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. Thread-safe implementation.

clear_boundary_cache!(cache::BoundaryCache)

Clear all cached boundary data.

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.

copy_boundary_conditions!(dest_𝔽, src_𝔽, field_type::FieldType)

Copy boundary conditions from one field to another. Both fields must have boundary condition fields already defined.

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_interpolation_cache(boundary_data::BoundaryData, target_theta::Vector{T}, 
                          target_phi::Vector{T}) where T

Create interpolation cache for efficient repeated interpolations.

create_programmatic_boundary(ylm::Ylm, config, amplitude::Real=1.0; kwargs...)

Create a boundary pattern from a spherical harmonic Yₗᵐ using SHTnsKit synthesis. Values are produced on the Gauss-Legendre grid of the given config.

Examples

boundary = create_programmatic_boundary(Ylm(2, 1), config, 0.5)
boundary = create_programmatic_boundary(Ylm(4, -2), config, 1.0)

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

Create a ProgrammaticBoundarySet for composition from programmatic specifications. Same interface as create_programmatic_temperature_boundaries.

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

Create a ProgrammaticBoundarySet for temperature from programmatic specifications.

Arguments

• inner_spec: Tuple specifying inner boundary, e.g. (:uniform, 100.0)
• outer_spec: Tuple specifying outer boundary, e.g. (:uniform, 250.0)
• cfg: SHTnsKitConfig with grid parameters

Returns

A ProgrammaticBoundarySet{Float64} containing boundary data and BC types.

Examples

bset = create_programmatic_temperature_boundaries((:uniform, 1.0), (:uniform, 0.0), cfg)
bset = create_programmatic_temperature_boundaries((:dirichlet, 1.0), (:neumann, 0.0), cfg)

create_time_dependent_programmatic_boundary(ylm::Ylm, config, time_span, ntime; kwargs...)

Create a time-dependent boundary pattern from a spherical harmonic Yₗᵐ with cosine time modulation. Values are produced on the Gauss-Legendre grid.

Example

boundary = create_time_dependent_programmatic_boundary(Ylm(3, 2), config, (0.0, 10.0), 100;
    amplitude=0.5, time_factor=2π)

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_magnetic_bc_in_solution!(solution, magnetic_field)

Enforce magnetic field boundary conditions in solution vector.

enforce_temperature_bc_in_solution!(solution, temp_field)

Enforce temperature boundary conditions in solution vector.

enforce_velocity_bc_in_solution!(solution, velocity_field)

Enforce velocity boundary conditions in solution vector.

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_boundary_time_index(boundary_set::BoundaryConditionSet, current_time::Float64)

Find the appropriate time index for the current simulation time. Used by field-specific boundary condition modules.

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.

get_bc_vectors_from_field(field)

Extract boundary condition vectors from a field's cache for use in implicit solves.

Returns

A named tuple (inner_real, outer_real, inner_imag, outer_imag) where each element is either a Vector{T} of length nlm or nothing if no file BCs are loaded.

get_boundary_condition_summary(𝔽, 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, initializing MPI if needed. Always returns a valid communicator (never nothing).

get_current_boundaries(field, field_type::FieldType)

Get current boundary values for any field type.

get_current_simulation_time(solver_state)

Extract current simulation time from solver state.

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};
                          is_periodic::Bool=false) where T

Calculate interpolation weights for linear interpolation. Handles periodic wrapping when is_periodic=true and indices wrap around (e.g., (n, 1)).

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_time_index(field, field_type::FieldType)

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

initialize_boundary_conditions!(𝔽, field_type::FieldType)

Initialize boundary condition support for a field structure.

IMPORTANT: The field structure must already have the following fields defined:

• boundaryconditionset (can be nothing)
• boundaryinterpolationcache (Dict{String, Any})
• boundarytimeindex (Ref{Int})

For scalar fields (TEMPERATURE, COMPOSITION):

• bctypeinner, bctypeouter (Vector{Int})
• boundary_values (Matrix)

For vector fields (VELOCITY, MAGNETIC):

• toroidal and poloidal components, each with bctypeinner, bctypeouter, boundary_values

Note: The field struct must be a mutable struct since this function modifies its fields.

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_spectral_bc_from_file(filename::String, config; format::Symbol=:physical, T::Type=Float64)

Load boundary conditions from a NetCDF file and return spectral coefficients.

Arguments

• filename: Path to the NetCDF file
• config: SHTnsKitConfig (or any object with lmax, mmax, nlm fields)
• format: :physical for physical-space data, :spectral for pre-computed coefficients
• T: Numeric type for the coefficients (default Float64)

NetCDF File Formats

Physical format (:physical):

• Variables: "inner" [nlat, nlon] and/or "outer" [nlat, nlon]
• Or a single variable with 3rd dimension of size 2: "temperature" nlat, nlon, 2

Spectral format (:spectral, Fortran-compatible):

• Variables: "Cbc_Re" [2, nlm] and "Cbc_Im" [2, nlm]
• Row 1 = inner boundary, Row 2 = outer boundary

Returns

• SpectralBoundaryCoefficients{T} with spectral BC data

log_boundary_condition_status(state, rank::Int=0)

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

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_type::FieldType)

Reset/clear boundary conditions for a field.

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.

store_bc_in_field!(field, bc_coeffs::SpectralBoundaryCoefficients)

Store spectral boundary coefficients into a field's boundary_interpolation_cache. The field must have a boundary_interpolation_cache::Dict{String, Any} attribute.

After calling this, get_bc_vectors_from_field(field) will return the stored vectors.

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_time_dependent_boundaries!(field, field_type::FieldType, current_time::Float64)

Update time-dependent boundary conditions for any field type.

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_field_boundary_compatibility(𝔽, field_type::FieldType, boundary_set::BoundaryConditionSet)

Validate that a field structure is compatible with boundary conditions.

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_netcdf_boundary_file(filename::String, required_vars::Vector{String}=[])

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

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

Write boundary condition data to a NetCDF file.

BoundaryDerivativeCache{T}

Cache of boundary values and radial derivatives for a spectral field. All values are stored for m >= 0 (SHTnsKit triangular indexing) and are expanded to negative m using conjugate symmetry.

GauntTensorCache(lmax, lmax_topo; nth=0, nph=0)
GauntTensorCache{T}

Cache for pre-computed Gaunt tensor integrals used in topography mode coupling.

Gaunt integrals arise from triple products of spherical harmonics:

\[G_{\ell m, \ell' m', LM} = \int Y_\ell^{m*} Y_{\ell'}^{m'} Y_L^{M} d\Omega\]

This cache stores three types of Gaunt tensors:

• G: Basic Gaunt integrals for shift terms
• G_∇: Gradient Gaunt integrals G^{(∇)} for slope terms
• G_cross: Cross Gaunt integrals G^{(×)} for tangential coupling

Arguments

• lmax::Int: Maximum spherical harmonic degree for field variables
• lmax_topo::Int: Maximum degree for topography expansion
• nth::Int=0: Number of θ quadrature points (auto if 0)
• nph::Int=0: Number of φ quadrature points (auto if 0)

Example

gaunt = GauntTensorCache(64, 16)
precompute_gaunt_tensors!(gaunt)

See also: precompute_gaunt_tensors!, get_gaunt_tensor

GauntTensorCache{T}(lmax::Int, lmax_topo::Int; nth::Int=0, nph::Int=0) where T

Create a Gaunt tensor cache with SHTnsKit configuration.

Arguments

• lmax::Int: Maximum spherical harmonic degree for fields
• lmax_topo::Int: Maximum degree for topography
• nth::Int: Number of theta points (default: 2*max(lmax,lmax_topo) + 2)
• nph::Int: Number of phi points (default: 2*nth)

StefanState(; kwargs...)
StefanState{T}

Stefan condition state for tracking ICB (inner core boundary) evolution.

The Stefan condition relates the rate of phase boundary motion to the heat flux discontinuity across the interface:

\[\rho L \frac{\partial h_i}{\partial t} = k_{oc} \frac{\partial T}{\partial n}\bigg|_{oc} - k_{ic} \frac{\partial T}{\partial n}\bigg|_{ic}\]

This structure tracks the ICB topography evolution including latent heat release and optional Clapeyron slope corrections.

Keyword Arguments

• lmax::Int=32: Maximum spherical harmonic degree
• ri::Float64=0.35: Inner core radius (nondimensional)
• k_ic::Float64=1.0: Inner core thermal conductivity
• k_oc::Float64=1.0: Outer core thermal conductivity
• rho::Float64=1.0: Density at ICB
• L::Float64=1.0: Latent heat of fusion
• use_clapeyron::Bool=false: Include Clapeyron slope correction
• clapeyron_slope::Float64=0.0: Clapeyron slope dT_m/dP (K/Pa)

Example

state = StefanState(lmax=64, ri=0.35, k_ic=2.0, k_oc=1.0)

See also: initialize_stefan_state!, update_icb_topography!

StefanState{T}(; lmax=32, ri=0.35, kwargs...) where T

Create an uninitialized Stefan state with default parameters.

Arguments

• lmax::Int=32: Maximum spherical harmonic degree
• ri::Float64: Inner core radius
• k_ic::Float64=1.0: Inner core thermal conductivity (nondimensional)
• k_oc::Float64=1.0: Outer core thermal conductivity (nondimensional)
• rho::Float64=1.0: Density (nondimensional)
• L::Float64=1.0: Latent heat (nondimensional)
• use_clapeyron::Bool=false: Include Clapeyron correction
• clapeyron_slope::Float64=0.0: Clapeyron slope dT_m/dP

TopographyCouplingConfig

Configuration structure for enabling/disabling topography coupling.

Fields

• enabled::Bool: Master switch for all topography coupling (default: false)
• velocity_coupling::Bool: Enable velocity BC topography corrections
• magnetic_coupling::Bool: Enable magnetic BC topography corrections
• thermal_coupling::Bool: Enable thermal BC topography corrections
• stefan_enabled::Bool: Enable Stefan condition for ICB evolution
• include_shift_terms::Bool: Include O(εh) shift terms (more accurate but slower)
• include_slope_terms::Bool: Include O(ε∇h) slope terms (required for coupling)
• epsilon::Float64: Topography amplitude parameter ε (default: 0.01)
• lmax_topo::Int: Maximum degree for topography expansion (default: same as simulation)

Example

config = TopographyCouplingConfig(
    enabled = true,
    velocity_coupling = true,
    magnetic_coupling = true,
    epsilon = 0.02
)

See also: enable_topography!, get_topography_config

TopographyData
TopographyData{T}

Complete topography data container for both ICB and CMB boundaries.

This structure holds the spectral topography fields for inner and outer boundaries, along with pre-computed Gaunt tensor caches for efficient mode coupling calculations.

Fields

• icb::Union{TopographyField, Nothing}: Inner core boundary (ICB) topography
• cmb::Union{TopographyField, Nothing}: Core-mantle boundary (CMB) topography
• gaunt_cache::Union{GauntTensorCache, Nothing}: Pre-computed Gaunt tensors
• epsilon::Float64: Topography amplitude parameter ε
• is_time_dependent::Bool: Whether topography evolves (e.g., Stefan condition)

Example

topo_data = TopographyData{Float64}(icb_field, cmb_field, gaunt, 0.02, false)

See also: TopographyField, create_topography_data, GauntTensorCache

TopographyData{T}() where T

Create empty TopographyData with no topography.

TopographyField(lmax, mmax, radius, location)
TopographyField{T}

Spectral representation of boundary topography h(θ, φ) expanded in spherical harmonics.

The topography is represented as:

\[h(\theta, \phi) = \sum_{L=0}^{L_{max}} \sum_{M=-L}^{L} h_{LM} Y_L^M(\theta, \phi)\]

Arguments

• lmax::Int: Maximum spherical harmonic degree L
• mmax::Int: Maximum azimuthal order M (capped at lmax)
• radius::Float64: Reference boundary radius
• location::BoundaryLocation: INNER_BOUNDARY (ICB) or OUTER_BOUNDARY (CMB)

Fields

• radius: Reference radius of the boundary
• lmax, mmax: Maximum degree and order
• nlm: Total number of spectral coefficients
• coeffs_real, coeffs_imag: Real and imaginary parts of h_{LM}
• location: Boundary location enum
• rms_amplitude, max_amplitude: Topography statistics

Example

topo = TopographyField(32, 32, 0.35, INNER_BOUNDARY)

See also: TopographyData, set_topography_coefficients!

TopographyField{T}(lmax::Int, mmax::Int, radius::T, location::BoundaryLocation) where T

Create an empty TopographyField with zero coefficients.

__apply_∂r!(output, matrix, input)

Local banded-matrix apply for derivative operators.

_compute_gradient_fallback(l, m, cache)

Fallback gradient computation using finite differences on the Y_l^m grid.

_gauss_legendre_point(n, i)

Compute the i-th Gauss-Legendre quadrature point and weight for n points.

_get_gauss_legendre_from_shtns(sht_config, nlat)

Extract Gauss-Legendre quadrature points and weights. SHTnsKit uses Gauss-Legendre points internally.

_wigner_3j(l1, l2, l3, m1, m2, m3)

Compute Wigner 3j symbol using Racah formula.

apply_all_topography_corrections!(fields, topography, config)

Apply all enabled topography corrections to the simulation fields.

This is the main entry point for applying topography coupling during the simulation timestep.

Arguments

• fields: Named tuple or struct containing velocity, magnetic, temperature fields
• topography: TopographyData containing ICB and CMB topography
• config: TopographyCouplingConfig (optional, uses global if not provided)

apply_clapeyron_correction!(state::StefanState{T}, pressure_field) where T

Apply Clapeyron slope correction to the melting temperature at ICB.

The melting temperature varies with pressure: Tm(P) = Tm0 + (dT_m/dP) ΔP

This modifies the effective temperature boundary condition at the ICB.

apply_composition_topography_correction!(composition_field, topography::TopographyData,
                                         config::TopographyCouplingConfig)

Apply topography corrections to compositional boundary conditions.

Composition follows the same mathematical structure as temperature, so we can reuse the thermal correction functions.

apply_magnetic_correction_at_boundary!(poloidal, toroidal, topo_field,
                                       gaunt, ε, config, location, bc_type)

Apply magnetic topography corrections at a specific boundary.

Arguments

• bc_type: :insulatinginner, :insulatingouter, or :conducting_inner

apply_magnetic_topography_correction!(magnetic_field, topography::TopographyData,
                                      config::TopographyCouplingConfig)

Apply topography corrections to magnetic boundary conditions.

This modifies the boundary values of the magnetic field's poloidal and toroidal components to account for topographic coupling.

Arguments

• magnetic_field: SHTnsMagneticFields or similar structure with toroidal/poloidal fields
• topography: TopographyData containing ICB and/or CMB topography
• config: TopographyCouplingConfig with coupling parameters

apply_thermal_correction_at_boundary!(spectral, cache, boundary_values,
                                      bc_inner, bc_outer, topo_field,
                                      gaunt, ε, config, location, T_cond)

Apply thermal topography corrections at a specific boundary.

apply_thermal_topography_correction!(temperature_field, topography::TopographyData,
                                     config::TopographyCouplingConfig;
                                     T_cond::Union{Function, Nothing}=nothing)

Apply topography corrections to thermal boundary conditions.

Arguments

• temperature_field: SHTnsTemperatureField or similar scalar field structure
• topography: TopographyData containing ICB and/or CMB topography
• config: TopographyCouplingConfig with coupling parameters
• T_cond: Optional conductive temperature profile function T_cond(r)

apply_velocity_correction_at_boundary!(poloidal, toroidal, topo_field,
                                       gaunt, ε, config, location)

Apply velocity topography corrections at a specific boundary.

apply_velocity_topography_correction!(velocity_field, topography::TopographyData,
                                      config::TopographyCouplingConfig)

Apply topography corrections to velocity boundary conditions.

This modifies the boundary values of the velocity field's poloidal and toroidal components to account for topographic coupling.

Arguments

• velocity_field: SHTnsVelocityFields or similar structure with toroidal/poloidal fields
• topography: TopographyData containing ICB and/or CMB topography
• config: TopographyCouplingConfig with coupling parameters

assemble_magnetic_boundary_operator(l::Int, topo::TopographyField{T},
                                    gaunt::GauntTensorCache{T},
                                    config::TopographyCouplingConfig,
                                    location::BoundaryLocation) where T

Assemble the magnetic boundary condition operator matrix row for mode l.

Returns a sparse representation of coupling coefficients.

Arguments

• l: Target spherical harmonic degree
• topo: Topography field at boundary
• gaunt: Gaunt tensor cache
• config: Coupling configuration
• location: INNERBOUNDARY (ICB) or OUTERBOUNDARY (CMB)

assemble_thermal_boundary_operator(l::Int, topo::TopographyField{T},
                                   gaunt::GauntTensorCache{T},
                                   config::TopographyCouplingConfig,
                                   bc_type::Symbol) where T

Assemble the thermal boundary condition operator matrix row for mode l.

Arguments

• l: Target spherical harmonic degree
• topo: Topography field at boundary
• gaunt: Gaunt tensor cache
• config: Coupling configuration
• bc_type: :dirichlet or :neumann

assemble_velocity_boundary_operator(l::Int, topo::TopographyField{T},
                                    gaunt::GauntTensorCache{T},
                                    config::TopographyCouplingConfig,
                                    bc_type::Symbol) where T

Assemble the boundary condition operator matrix row for mode l.

Returns a sparse representation of coupling coefficients to other modes.

Arguments

• l: Target spherical harmonic degree
• topo: Topography field at boundary
• gaunt: Gaunt tensor cache
• config: Coupling configuration
• bc_type: :impermeability, :noslip, or :stressfree

compute_associated_legendre(l::Int, m::Int, x::T) where T

Compute the associated Legendre polynomial P_l^m(x). Uses the standard normalization (not Schmidt or fully normalized).

compute_boundary_derivative_cache(field, ∂r, ∂²r, domain) -> BoundaryDerivativeCache

Compute boundary values and radial derivatives for all modes using the full radial profile (MPI-safe). This is expensive but avoids per-mode Allreduce inside tight coupling loops.

compute_boundary_heat_flux(temperature_field, topography::TopographyData,
                           config::TopographyCouplingConfig,
                           location::BoundaryLocation;
                           k::Float64=1.0) -> Matrix

Compute the heat flux at a topographic boundary.

The normal heat flux on the true surface (to first order) is: qn = -k [∂r T - ε ∇H h · ∇H T + εh ∂_rr T]

Returns heat flux on the (θ, φ) grid.

compute_boundary_heat_flux_spectral(temperature_field, r::T, side::Symbol) where T

Compute spectral coefficients of heat flux at a boundary.

Arguments

• temperature_field: Temperature field
• r: Boundary radius
• side: :inner or :outer

Returns vector of spectral coefficients for ∂_r T at the boundary.

compute_cmb_insulating_correction(l, m, poloidal, toroidal, topo,
                                  gaunt, ro, config) -> (P_corr, T_corr)

Compute topography correction to CMB insulating boundary condition.

Flat-sphere conditions:

• ∂r P{b,ℓm} + (ℓ+1)/ro P{b,ℓm} = 0
• T_{b,ℓm} = 0

With topography, the poloidal condition becomes: [∂r P + (ℓ+1)/ro P]{ℓm} + ε Σ h^o{LM} (α^o ∂_r P + β^o T + γ^o P) = 0

and toroidal: T + εho ∂r T = 0

compute_cross_gaunt_tensor(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int,
                           cache::GauntTensorCache{T}) where T

Compute the cross Gaunt integral: G^{(×)}{l1,m1,l2,m2,L,M} = ∫ Y{l1}^{m1*} r̂ · (∇H Y{l2}^{m2} × ∇H YL^M) dΩ

This integral couples toroidal and poloidal modes through topography. Uses SHTnsKit for gradient evaluation when available.

compute_dirichlet_thermal_correction(l, m, spectral, topo, gaunt, rb, config, dTcond_dr)

Compute topography correction to Dirichlet thermal BC for mode (l, m).

From PDF equation 38: Θ{ℓm} + ε Σ h^b{LM} G{ℓm,ℓ'm',LM} ∂r Θ{ℓ'm'} = [Tb - Tcond(rb)]{ℓm} - ε Σ h^b{LM} G{ℓm,00,LM} ∂r T_cond

The correction to the RHS involves:

1. Shift term: h · ∂_r Θ (couples modes)
2. Conductive correction: h · ∂r Tcond (modifies boundary value)

compute_gaunt_from_wigner3j(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int)

Compute Gaunt coefficient using Wigner 3j symbols (analytic formula).

For the Gaunt integral with conjugate: G = ∫ Y{l1}^{m1*} Y{l2}^{m2} Y_L^M dΩ

Using Yl^{m*} = (-1)^m Yl^{-m}, the formula is: G = (-1)^{m1} * sqrt((2l1+1)(2l2+1)(2L+1)/(4π)) * (l1 l2 L; 0 0 0) * (l1 l2 L; -m1 m2 M)

This is faster than numerical integration for production use.

compute_gaunt_tensor(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int,
                     cache::GauntTensorCache{T}) where T

Compute the basic Gaunt integral using SHTnsKit for spherical harmonic evaluation: G{l1,m1,l2,m2,L,M} = ∫ Y{l1}^{m1*} Y{l2}^{m2} YL^M dΩ

Uses numerical quadrature on Gauss-Legendre × uniform φ grid.

compute_gradient_gaunt_tensor(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int,
                              cache::GauntTensorCache{T}) where T

Compute the gradient Gaunt integral using the analytic identity: G^{(∇)}{l1,m1,l2,m2,L,M} = (1/2)[ℓ2(ℓ2+1) + L(L+1) - ℓ1(ℓ1+1)] G{l1,m1,l2,m2,L,M}

This is more efficient and accurate than numerical gradient computation.

compute_icb_insulating_correction(l, m, poloidal, toroidal, topo,
                                  gaunt, ri, config) -> (P_corr, T_corr)

Compute topography correction to ICB insulating boundary condition.

Flat-sphere conditions:

• ∂r P{b,ℓm} - ℓ/ri P{b,ℓm} = 0
• T_{b,ℓm} = 0

With topography, the poloidal condition becomes: [∂r P - ℓ/ri P]{ℓm} + ε Σ h^i{LM} (α^i ∂_r P + β^i T + γ^i P) = 0

compute_impermeability_correction(l, m, poloidal, toroidal, topo,
                                  gaunt, rb, config) -> Complex

Compute the topography correction to the impermeability condition for mode (l, m).

The flat-sphere condition is: uᵣ = ℓ(ℓ+1)/r² P_{ℓm} = 0

With topography: uᵣ - ε(∇H h)·uₜ + εh ∂r uᵣ = 0

In spectral space (schematic): ℓ(ℓ+1)/r² P{ℓm} + ε Σ{LM,ℓ'm'} h{LM} [ G · ∂r(ℓ'(ℓ'+1)/r² P{ℓ'm'}) # shift term - G^{(∇)} · ∂r P{ℓ'm'} # slope term (poloidal) - G^{(×)} · T{ℓ'm'} # slope term (toroidal) ] = 0

compute_neumann_thermal_correction(l, m, spectral, topo, gaunt, rb, config,
                                   dTcond_dr, d2Tcond_dr2)

Compute topography correction to Neumann thermal BC for mode (l, m).

From PDF equation 39: ∂r Θ{ℓm} - ε Σ h^b{LM} G^{(∇)}{ℓm,ℓ'm',LM} Θ{ℓ'm'} + ε Σ h^b{LM} G{ℓm,ℓ'm',LM} ∂rr Θ{ℓ'm'} = -q{b,ℓm}/k - ∂r Tcond δ{ℓ0}δ{m0} - ε Σ h^b{LM} G{ℓm,00,LM} ∂rr Tcond

The correction involves:

1. Slope term: -∇H h · ∇H Θ (gradient coupling)
2. Shift term: h · ∂_rr Θ (second derivative coupling)
3. Conductive correction: h · ∂rr Tcond

compute_normal_velocity_spectral(velocity_field, r::T,
                                 location::BoundaryLocation=OUTER_BOUNDARY) where T

Compute spectral coefficients of normal (radial) velocity at a boundary.

uₙ = uᵣ = ℓ(ℓ+1)/r² P at the boundary

Returns vector of spectral coefficients for uᵣ.

compute_noslip_correction(l, m, poloidal, toroidal, topo,
                          gaunt, rb, config) -> (Complex, Complex)

Compute topography correction to no-slip condition for mode (l, m).

The flat-sphere condition is: uₜ = 0

With topography: uₜ + εh ∂r uₜ = U{b,t}

Returns (poloidalcorrection, toroidalcorrection).

compute_potential_matching_coefficients(l::Int, topo::TopographyField{T},
                                       gaunt::GauntTensorCache{T},
                                       location::BoundaryLocation) where T

Compute the coefficients for matching the core magnetic field to an external potential field across a topographic boundary.

For an insulating exterior, B = -∇Φ where:

• Φext = Σ a{ℓm} r^{-(ℓ+1)} Yℓ^m (exterior, r > ro)
• Φint = Σ b{ℓm} r^ℓ Yℓ^m (interior, r < ri)

Continuity of B_t at the true surface introduces topography coupling.

compute_stefan_flux(state::StefanState{T}) -> Vector{Complex{T}}

Compute the Stefan flux contribution F_{ℓm} for topography evolution.

F{ℓm} = kic ∂n Θ^{ic}{ℓm} - k ∂n Θ{ℓm}

Returns spectral coefficients of the net heat flux imbalance.

compute_stefan_flux_with_topography(state::StefanState{T}, temperature_ic,
                                    temperature_oc, topo_data::TopographyData,
                                    gaunt::GauntTensorCache{T},
                                    config::TopographyCouplingConfig) where T

Compute Stefan flux including topography corrections to the normal derivative.

The linearized normal derivative at r = ri is: ∂n ≈ ∂r - ε ∇H hi · ∇H + εhi ∂rr

This adds coupling between modes through Gaunt tensors.

compute_stressfree_correction(l, m, poloidal, toroidal, topo,
                              gaunt, rb, config) -> Complex

Compute topography correction to stress-free condition for mode (l, m).

The flat-sphere condition is: ∂_r(uₜ/r) = 0

With topography: ∂r(uₜ/r) = (ε/rb)(∇H h) ∂r uᵣ

The RHS involves the slope of topography coupling with radial velocity gradient.

create_random_topography(spectrum::Function, radius::Float64,
                         location::BoundaryLocation; lmax::Int=32,
                         seed::Int=0) -> TopographyField

Create random topography with prescribed power spectrum.

Arguments

• spectrum: Function l -> power at degree l (e.g., l -> l^(-2) for red spectrum)
• radius: Boundary radius
• location: INNERBOUNDARY or OUTERBOUNDARY
• lmax: Maximum degree
• seed: Random seed (0 for no seed)

Example

# Create topography with l^(-2) spectrum
topo = create_random_topography(l -> 1.0/max(l,1)^2, 1.0, OUTER_BOUNDARY)

create_spherical_harmonic_topography(l::Int, m::Int, amplitude::Float64,
                                     radius::Float64, location::BoundaryLocation;
                                     lmax::Int=32) -> TopographyField

Create single spherical harmonic topography h = amplitude * Y_l^m(θ, φ). For m >= 0, sets the real (cosine) part. For m < 0, sets the imaginary (sine) part via conjugate symmetry, since only |m| entries are stored.

create_topography_data(; icb_coeffs=nothing, cmb_coeffs=nothing,
                        icb_radius=nothing, cmb_radius=nothing,
                        lmax=32, mmax=32, epsilon=0.01) -> TopographyData

Create TopographyData from spectral coefficients.

Arguments

• icb_coeffs: Complex vector of ICB topography coefficients h_{LM}
• cmb_coeffs: Complex vector of CMB topography coefficients h_{LM}
• icb_radius: Inner boundary radius (default: from parameters)
• cmb_radius: Outer boundary radius (default: from parameters)
• lmax, mmax: Maximum spherical harmonic degree/order
• epsilon: Topography amplitude parameter

Example

# Create CMB topography with degree-2 pattern
cmb_h = zeros(ComplexF64, 6)
cmb_h[lm_to_index(2, 0, 32)] = 0.1  # h_{2,0}
topo = create_topography_data(cmb_coeffs=cmb_h, cmb_radius=1.0, lmax=32)

create_uniform_topography(amplitude::Float64, radius::Float64,
                          location::BoundaryLocation; lmax::Int=32) -> TopographyField

Create uniform (degree-0) topography h = amplitude everywhere.

disable_topography!()

Disable all topography coupling.

enable_topography!(; kwargs...)

Enable topography coupling with optional configuration.

Arguments

• epsilon::Float64=0.01: Topography amplitude parameter
• velocity::Bool=true: Enable velocity coupling
• magnetic::Bool=true: Enable magnetic coupling
• thermal::Bool=true: Enable thermal coupling
• stefan::Bool=false: Enable Stefan condition
• slope_terms::Bool=true: Include slope terms
• shift_terms::Bool=true: Include shift terms
• lmax_topo::Int=-1: Maximum topography degree (-1 for auto)

Example

enable_topography!(epsilon=0.05, stefan=true)

evaluate_spherical_harmonic_gradient_grid(l::Int, m::Int, cache::GauntTensorCache{T}) where T

Evaluate ∇H Yl^m on the Gauss-Legendre × uniform φ grid using SHTnsKit.

Returns (∇θ, ∇φ) matrices of complex values.

evaluate_spherical_harmonics_grid(l::Int, m::Int, cache::GauntTensorCache{T}) where T

Evaluate Y_l^m on the Gauss-Legendre × uniform φ grid using SHTnsKit.

Returns a (nlat, nlon) matrix of complex values.

evaluate_topography(field::TopographyField{T}, config) where T -> Matrix{T}

Evaluate topography in physical space using SHTnsKit synthesis.

Arguments

• field: TopographyField containing spectral coefficients
• config: SHTnsKit configuration for the transform

Returns h(θ, φ) as a nlat × nlon matrix on the Gauss-Legendre grid.

evaluate_topography(field::TopographyField{T}, theta::Vector{T},
                    phi::Vector{T}) where T -> Matrix{T}

Evaluate topography in physical space on an arbitrary (θ, φ) grid. This uses direct summation and is slower than the SHTnsKit version.

Returns h(θi, φj) as a nlat × nlon matrix.

evaluate_topography_fallback(field::TopographyField{T}, config) where T

Fallback evaluation when SHTnsKit synthesis fails.

gaunt_on_the_fly(l1, m1, l2, m2, L, M; use_wigner=true)

Compute a single Gaunt coefficient on-the-fly. Useful when only a few coefficients are needed.

gauss_legendre_nodes(n::Int) -> (nodes, weights)

Compute Gauss-Legendre quadrature nodes and weights. Returns nodes in [-1, 1] and corresponding weights.

get_coefficient(field::TopographyField{T}, l::Int, m::Int) where T

Get the h_{l,m} coefficient from the topography field.

get_cross_gaunt(cache::GauntTensorCache{T}, l1, m1, l2, m2, L, M) where T

Get the cross Gaunt tensor G^{(×)}_{l1,m1,l2,m2,L,M} from cache.

get_gaunt_tensor(cache::GauntTensorCache{T}, l1, m1, l2, m2, L, M) where T

Get the basic Gaunt tensor G_{l1,m1,l2,m2,L,M} from cache. Returns 0 if not found (implies zero due to selection rules).

get_gradient_gaunt(cache::GauntTensorCache{T}, l1, m1, l2, m2, L, M) where T

Get the gradient Gaunt tensor G^{(∇)}_{l1,m1,l2,m2,L,M} from cache.

get_spectral_boundary_value(field, l::Int, m::Int, location::BoundaryLocation=OUTER_BOUNDARY)

Get the boundary value of spectral coefficient (l, m).

get_spectral_coefficient(field, l::Int, m::Int)

Get spectral coefficient (l, m) from a field.

get_spectral_radial_derivative(field, l::Int, m::Int, r,
                               location::BoundaryLocation=OUTER_BOUNDARY;
                               ∂r, domain)

Compute the radial derivative of spectral coefficient (l, m) at a boundary. Requires access to the radial derivative matrix and domain.