Boundary Topography

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

get_topography_config() -> TopographyCouplingConfig

Get the current topography coupling configuration.

set_topography_config!(config::TopographyCouplingConfig)

Set the topography coupling configuration.

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)

disable_topography!()

Disable all topography coupling.

is_topography_enabled() -> Bool

Check if topography coupling is enabled.

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!

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

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

precompute_gaunt_tensors!(cache::GauntTensorCache{T};
                          verbose::Bool=true, use_wigner::Bool=true) where T

Precompute all non-zero Gaunt tensors up to lmax.

Arguments

• verbose::Bool=true: Print progress information
• use_wigner::Bool=true: Use analytic Wigner 3j formula (faster and more accurate)

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!

initialize_stefan_state!(state::StefanState{T}, temperature_ic, temperature_oc,
                         velocity_field; gaunt_cache=nothing) where T

Initialize the Stefan state from current field values.

Arguments

• state: StefanState to initialize
• temperature_ic: Inner core temperature field (or spectral coefficients)
• temperature_oc: Outer core temperature field
• velocity_field: Velocity field for computing uₙ
• gaunt_cache: Optional pre-computed Gaunt tensors

update_icb_topography!(state::StefanState{T}, dt::T, velocity_field,
                       temperature_ic, temperature_oc;
                       topo_data::Union{TopographyData, Nothing}=nothing,
                       gaunt::Union{GauntTensorCache{T}, Nothing}=nothing,
                       config::Union{TopographyCouplingConfig, Nothing}=nothing) where T

Update the ICB topography based on the Stefan condition.

From Eq. 22/40: ε ∂t hi = uₙ + (kic ∂n Tic - k ∂n T) / (ρ L)

or equivalently with Stefan number: ε ∂t hi = uₙ + (1/St) (λ ∂n Θic - ∂_n Θ)

Arguments

• state: StefanState to update
• dt: Time step
• velocity_field: Current velocity field
• temperature_ic: Inner core temperature
• temperature_oc: Outer core temperature
• topo_data: Optional TopographyData for coupled evolution
• gaunt: Optional Gaunt cache for topography coupling
• config: Optional coupling configuration

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.

config is retained so reads index modes in the SAME canonical m-major order the cache was filled in (via local_spectral_storage_slot / get_mode_index).

__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 (B_r = λP/r² ⇒ exterior vacuum P ∝ r^{-l}):

• ∂r P{b,lm} + l/ro P{b,lm} = 0
• T_{b,lm} = 0

With topography, the poloidal condition becomes: [∂r P + l/ro P]{lm} + ε Σ 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: Θ{lm} + ε Σ h^b{LM} G{lm,l'm',LM} ∂r Θ{l'm'} = [Tb - Tcond(rb)]{lm} - ε Σ h^b{LM} G{lm,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)[l2(l2+1) + L(L+1) - l1(l1+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 (B_r = λP/r² ⇒ interior vacuum P ∝ r^{l+1}):

• ∂r P{b,lm} - (l+1)/ri P{b,lm} = 0
• T_{b,lm} = 0

With topography, the poloidal condition becomes: [∂r P - (l+1)/ri P]{lm} + ε Σ 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ᵣ = l(l+1)/r² P_{lm} = 0

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

In spectral space (schematic): l(l+1)/r² P{lm} + ε Σ{LM,l'm'} h{LM} [ G · ∂r(l'(l'+1)/r² P{l'm'}) # shift term - G^{(∇)} · ∂r P{l'm'} # slope term (poloidal) - G^{(×)} · T{l'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 Θ{lm} - ε Σ h^b{LM} G^{(∇)}{lm,l'm',LM} Θ{l'm'} + ε Σ h^b{LM} G{lm,l'm',LM} ∂rr Θ{l'm'} = -q{b,lm}/k - ∂r Tcond δ{l0}δ{m0} - ε Σ h^b{LM} G{lm,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ᵣ = l(l+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{lm} r^{-(l+1)} Yl^m (exterior, r > ro)
• Φint = Σ b{lm} r^l Yl^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_{lm} for topography evolution.

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

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.

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.

get_stefan_diagnostics(state::StefanState{T}) where T

Get diagnostic information about the Stefan condition state.

get_topography_coefficients(field::TopographyField{T}) -> Vector{Complex{T}}

Get complex spectral coefficients from topography field.

gradient_gaunt_from_basic(l1::Int, l2::Int, L::Int, G_basic::T) where T

Compute gradient Gaunt using the identity: G^{(∇)} = (1/2)[l2(l2+1) + L(L+1) - l1(l1+1)] G

This avoids computing gradients explicitly.

index_to_lm(idx::Int, lmax::Int) -> (l, m)

Convert linear index back to (l, m).

initialize_gaunt_cache!(topo::TopographyData{T}, lmax_field::Int;
                        precompute::Bool=true) where T

Initialize and optionally precompute Gaunt tensor cache for topography coupling.

is_no_slip_boundary(field, location::BoundaryLocation) -> Bool

Check if the field uses Dirichlet tangential conditions at the boundary.

is_stress_free_boundary(field, location::BoundaryLocation) -> Bool

Check if the field has stress-free boundary conditions at the given location.

lm_to_index(l::Int, m::Int, lmax::Int) -> Int

Convert (l, m) to linear index for spectral coefficient storage. Uses the convention: index = l(l+1)/2 + m + 1 for m >= 0

lm_to_spectral_index(l::Int, m::Int, config) -> Int

Convert (l, m) to spectral array index using SHTnsKit convention.

load_topography_from_array(h_physical::Matrix{T}, radius::Float64,
                           location::BoundaryLocation, config;
                           lmax::Int=-1) where T

Load topography from a physical space array h(θ, φ) by transforming to spectral space.

Arguments

• h_physical: 2D array of topography values on (θ, φ) grid
• radius: Boundary radius
• location: INNERBOUNDARY or OUTERBOUNDARY
• config: SHTnsKit configuration for transform
• lmax: Maximum degree to retain (-1 for config.lmax)

load_topography_from_file(filename::String, location::BoundaryLocation;
                          radius::Union{Float64, Nothing}=nothing) -> TopographyField

Load topography from a NetCDF file.

Expected NetCDF structure:

• Dimensions: l, m (or lm for 1D storage)
• Variables: hreal, himag (or h for complex)
• Attributes: lmax, radius (optional)

print_stefan_summary(state::StefanState)

Print a summary of the Stefan condition state.

print_topography_summary(topography::TopographyData)

Print a summary of the current topography configuration.

save_topography_to_file(field::TopographyField, filename::String)

Save topography field to NetCDF file.

set_topography_coefficients!(field::TopographyField{T}, coeffs::Vector) where T

Set spectral coefficients for a topography field.

Accepts either:

• Complex vector: splits into real and imaginary parts
• Real vector: sets real parts only (imaginary = 0)

update_icb_topography_semiimplicit!(state::StefanState{T}, dt::T, velocity_field,
                                    temperature_ic, temperature_oc;
                                    theta::T=0.5) where T

Update ICB topography using a semi-implicit (θ-scheme) time integration.

h^{n+1} = h^n + dt [θ RHS^{n+1} + (1-θ) RHS^n]

where RHS = (1/ε) [uₙ + F/(ρL)]

Arguments

• theta::T=0.5: Implicitness parameter (0.5 = Crank-Nicolson, 1.0 = backward Euler)

update_topography_statistics!(field::TopographyField{T}) where T

Update RMS and max amplitude statistics for the topography field.