Huginn.jl API

HalfarParameters(; λ=0.0, H₀=3600.0, R₀=750000.0, n=3.0, A=1e-16, f=0.0, ρ=900.0, g=9.81)

Holds parameters for the Halfar similarity solution of the shallow ice approximation (SIA).

Parameters

• λ::AbstractFloat=0.0: Mass balance coefficient (used to model accumulation/ablation).
• H₀::AbstractFloat=3600.0: Dome height at initial time t₀ [m].
• R₀::AbstractFloat=750000.0: Ice sheet margin radius at t₀ [m].
• n::AbstractFloat=3.0: Glen flow law exponent.
• A::AbstractFloat=1e-16: Flow rate factor in Glen's law [Pa⁻ⁿ yr⁻¹].
• f::AbstractFloat=0.0: Fraction of isostatic bed depression (0 for fully grounded ice).
• ρ::AbstractFloat=900.0: Ice density in kg/m³ (same default as Sleipnir.PhysicalParameters).
• g::AbstractFloat=9.81: Gravitational acceleration [m/s²].

Notes

Default parameters set as in Bueler (2005) "Exact solutions and verification of numerical models for isothermalice sheets", experiment B.

IceflowModel

An abstract type representing the base model for ice flow simulations. All specific ice flow models should subtype this abstract type.

Prediction{CACHE} <: Simulation

A mutable struct that represents a prediction simulation.

Fields

• model::Sleipnir.Model: The model used for the prediction.
• glaciers::Vector{Sleipnir.AbstractGlacier}: A vector of glaciers involved in the prediction.
• parameters::Sleipnir.Parameters: The parameters used for the prediction.
• results::Vector{Results}: A vector of results obtained from the prediction.

Prediction(model::Sleipnir.Model, glaciers::Vector{G}, parameters::Sleipnir.Parameters) where {G <: Sleipnir.AbstractGlacier}

Create a Prediction object using the given model, glaciers, and parameters.

Arguments

• model::Sleipnir.Model: The model used for prediction.
• glaciers::Vector{G}: A vector of glacier objects, where each glacier is a subtype of Sleipnir.AbstractGlacier.
• parameters::Sleipnir.Parameters: The parameters used for the prediction.

Returns

• Prediction: A Prediction object based on the input values.

SIA2DCache{
    R <: Real,
    I <: Integer,
    A_CACHE,
    C_CACHE,
    n_CACHE,
    p_CACHE,
    q_CACHE,
    n_H_CACHE,
    n_∇S_CACHE,
    Y_CACHE,
    U_CACHE
} <: SIAmodel

Store and preallocated all variables needed for running the 2D Shallow Ice Approximation (SIA) model efficiently.

Type Parameters

• R: Real number type used for physical fields.
• I: Integer type used for indexing glaciers.
• A_CACHE, C_CACHE, n_CACHE: Types used for caching A, C, and n, which can be scalars, vectors, or matrices.
• Y_CACHE: Type used for caching Y which is a matrix.
• U_CACHE: Type used for caching U which is a matrix.

Fields

• A::A_CACHE: Flow rate factor.
• C::C_CACHE: Sliding coefficient.
• n::n_CACHE: Flow law exponent.
• p::n_CACHE: Sliding law exponent.
• q::n_CACHE: Sliding law exponent.
• n_H::n_CACHE: Exponent used for the power of H when using the Y law.
• n_∇S::n_CACHE: Exponent used for the power of ∇S when using the Y law.
• Y::Y_CACHE: Hybrid diffusivity.
• U::U_CACHE: Diffusive velocity.
• H₀::Matrix{R}: Initial ice thickness.
• H::Matrix{R}: Ice thickness.
• H̄::Matrix{R}: Averaged ice thickness.
• S::Matrix{R}: Surface elevation.
• dSdx::Matrix{R}: Surface slope in the x-direction.
• dSdy::Matrix{R}: Surface slope in the y-direction.
• D::Matrix{R}: Diffusivity.
• Dx::Matrix{R}: Diffusivity in the x-direction.
• Dy::Matrix{R}: Diffusivity in the y-direction.
• dSdx_edges::Matrix{R}: Surface slope at edges in the x-direction.
• dSdy_edges::Matrix{R}: Surface slope at edges in the y-direction.
• ∇S::Matrix{R}: Norm of the surface gradient.
• ∇Sy::Matrix{R}: Surface gradient component in the y-direction.
• ∇Sx::Matrix{R}: Surface gradient component in the x-direction.
• Fx::Matrix{R}: Flux in the x-direction.
• Fy::Matrix{R}: Flux in the y-direction.
• Fxx::Matrix{R}: Second derivative of flux in the x-direction.
• Fyy::Matrix{R}: Second derivative of flux in the y-direction.
• V::Matrix{R}: Velocity magnitude.
• Vx::Matrix{R}: Velocity in the x-direction.
• Vy::Matrix{R}: Velocity in the y-direction.
• Γ::A_CACHE: Basal shear stress.
• MB::Matrix{R}: Mass balance.
• MB_history::Vector{Matrix{R}}: Mass balance history.
• MB_times::Vector{R}: Mass balance time steps.
• MB_mask::BitMatrix: Boolean mask for applying the mass balance.
• MB_total::Matrix{R}: Total mass balance field.
• glacier_idx::I: Index of the glacier for use in simulations with multiple glaciers.
• A_prep_vjps, C_prep_vjps, n_prep_vjps, Y_prep_vjps and U_prep_vjps: Structs that contain the prepared VJP functions for the adjoint computation and for the different laws. Useful mainly when the user does not provide the VJPs and they are automatically inferred using DifferentiationInterface.jl which requires to store precompiled functions. When no gradient is computed, these structs are nothing.

SIA2Dmodel(A, C, n, Y, U, n_H, n_∇S)
SIA2Dmodel(params; A, C, n, Y, U, n_H, n_∇S)

Create a SIA2Dmodel, representing a two-dimensional Shallow Ice Approximation (SIA) model.

The SIA model describes glacier flow under the assumption that deformation and basal sliding dominate the ice dynamics. It relies on:

• Glen's flow law for internal deformation, with flow rate factor A and exponent n,

• A sliding law governed by coefficient C,

• Optionally the user can provide either:

  • A specific diffusive velocity U such that D = U * H
  • A modified creep coefficient Y that takes into account the ice thickness such that D = (C + Y * 2/(n+2)) * (ρ*g)^n * H^(n_H+1) * |∇S|^(n_∇S-1) where n_H and n_∇S are optional parameters that control if the SIA should use the n law or not. This formulation is denoted as the hybrid diffusivity in the code.

This struct stores the laws used to compute these three parameters during a simulation. If not provided, default constant laws are used based on glacier-specific values.

Arguments

• A: Law for the flow rate factor. Defaults to a constant value from the glacier.
• C: Law for the sliding coefficient. Defaults similarly.
• n: Law for the flow law exponent. Defaults similarly.
• p: Law for the sliding law exponent (basal drag). Defaults similarly.
• q: Law for the sliding law exponent (normal stress). Defaults similarly.
• Y: Law for the hybrid diffusivity. Providing a law for Y discards the laws of A, C and n.
• U: Law for the diffusive velocity. Defaults behavior is to disable it and in such a case it is computed from A, C and n. Providing a law for U discards the laws of A, C, n and Y.
• n_H::F: The exponent to use for H in the SIA equation when using the Y law (hybrid diffusivity). It should be nothing when this law is not used.
• n_∇S::F: The exponent to use for ∇S in the SIA equation when using the Y law (hybrid diffusivity). It should be nothing when this law is not used.
• Y_is_provided::Bool: Whether the diffusivity is provided by the user through the hybrid diffusivity Y or it has to be computed from the SIA formula from A, C and n.
• U_is_provided::Bool: Whether the diffusivity is provided by the user through the diffusive velocity U or it has to be computed from the SIA formula from A, C and n.
• n_H_is_provided::Bool: Whether the H exponent is prescribed by the user, or the one of the n law has to be used. This flag is used only when a law for Y is used.
• n_∇S_is_provided::Bool: Whether the ∇S exponent is prescribed by the user, or the one of the n law has to be used. This flag is used only when a law for Y is used.
• apply_A_in_SIA::Bool: Whether the value of the A law should be computed each time the SIA is evaluated.
• apply_C_in_SIA::Bool: Whether the value of the C law should be computed each time the SIA is evaluated.
• apply_n_in_SIA::Bool: Whether the value of the n law should be computed each time the SIA is evaluated.
• apply_p_in_SIA::Bool: Whether the value of the p law should be computed each time the SIA is evaluated.
• apply_q_in_SIA::Bool: Whether the value of the q law should be computed each time the SIA is evaluated.
• apply_Y_in_SIA::Bool: Whether the value of the Y law should be computed each time the SIA is evaluated.
• apply_U_in_SIA::Bool: Whether the value of the U law should be computed each time the SIA is evaluated.

SIAmodel

An abstract type representing the Shallow Ice Approximation (SIA) models. This type is a subtype of IceflowModel and serves as a base for all SIA-specific models.

A mutable struct that holds parameters for the solver.

SolverParameters{F <: AbstractFloat, I <: Integer, ST <: OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm}

Fields

• solver::ST: The algorithm used for solving differential equations.
• supply_eigen_est::Bool: Whether to give ROCK2 and ROCK4 our own estimate of the spectral radius instead of letting them compute it (see with_eigen_est). Other solvers ignore it. Only useful for finite difference gradient checks.
• reltol::F: The relative tolerance for the solver.
• abstol::F: The absolute tolerance for the solver, in metres of ice. This is the tolerance that matters in practice, since reltol * H is much smaller over most of the glacier. It is tightened for runs longer than 5 years unless scale_abstol is false (see effective_abstol).
• scale_abstol::Bool: Whether to tighten abstol for long runs (see effective_abstol). Set to false to use abstol as it is.
• adaptive::Bool: Whether the solver chooses its own step size. Adaptive stepping makes the solution discontinuous in the parameters — an arbitrarily small change can flip which steps are accepted — which is harmless for a forward run but makes finite differences meaningless. Gradient checks want adaptive = false; forward runs don't.
• dt::F: Fixed step size used when adaptive is false, in years. Ignored otherwise.
• step::F: The step size that controls at which frequency the results must be saved.
• tstops::Vector{F}: Optional vector of time points where the solver should stop to store the results.
• save_everystep::Bool: Flag indicating whether to save the solution at every step computed by the solver.
• progress::Bool: Flag indicating whether to show progress during the solving process.
• progress_steps::I: The number of steps between progress updates.
• maxiters::I: Maximum number of iterations to perform in the iceflow solver.

Constructs a SolverParameters object with the specified parameters or using default values.

SolverParameters(;
    solver::OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm = RDPK3Sp35(),
    supply_eigen_est::Bool = false,
    reltol::F = 1e-12,
    abstol::F = 1e-3,
    scale_abstol::Bool = true,
    adaptive::Bool = true,
    dt::F = 1.0/120.0,
    step::F = 1.0/12.0,
    tstops::Vector{Sleipnir.Float} = Vector{Sleipnir.Float}(),
    save_everystep = false,
    progress::Bool = true,
    progress_steps::I = 10,
    maxiters::I = Int(1e5),
) where {F <: AbstractFloat, I <: Integer}

Arguments

• solver::OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm: The ODE solver algorithm to use. Defaults to RDPK3Sp35().
• supply_eigen_est::Bool: Whether to give ROCK2 and ROCK4 our own estimate of the spectral radius (see with_eigen_est). Only useful for finite difference gradient checks. Defaults to false.
• reltol::F: The relative tolerance for the solver. Defaults to 1e-12.
• abstol::F: The absolute tolerance for the solver, in metres of ice. Defaults to 1e-3, for a run of 5 years (see effective_abstol).
• scale_abstol::Bool: Whether to tighten abstol for long runs. Defaults to true.
• adaptive::Bool: Whether the solver picks its own step size. Defaults to true; gradient checks against finite differences need false.
• dt::F: Fixed step in years, used only when adaptive is false. Defaults to 1.0/120.0.
• step::F: The step size that controls at which frequency the solution should be computed and returned in the results. Defaults to 1.0/12.0 (i.e. a month).
• tstops::Vector{Sleipnir.Float}: Optional vector of time points where the solver should stop. Defaults to an empty vector.
• save_everystep::Bool: Whether to save the solution at every step computed by the solver. Defaults to false.
• progress::Bool: Whether to show progress during the solving process. Defaults to true.
• progress_steps::I: The number of steps between progress updates. Defaults to 10.
• maxiters::I: Maximum number of iterations to perform in the iceflow solver. Defaults to 1e5.

Returns

• solver_parameters: A SolverParameters object constructed with the specified parameters.

iAvgGriddedTemp <: AbstractInput

Input that represents the long term air temperature over the glacier grid. It is computed using the OGGM climate data over a period predefined in Gungnir (i.e. around 30 years).

iAvgScalarTemp <: AbstractInput

Input that represents the long term air temperature over the whole glacier. It is computed using the OGGM climate data over a period predefined in Gungnir (i.e. around 30 years).

iCPDD <: AbstractInput

Input that represents the cumulative positive degree days (PDD) over the last time window window. It is computed by summing the daily PDD values from t - window to t using the glacier's climate data.

iH̄ <: AbstractInput

Input that represents the ice thickness in the SIA. It is the averaged ice thickness computed on the dual grid, that is H̄ = avg(H) which is different from the ice thickness solution H.

iTopoRough{F<:AbstractFloat} <: AbstractInput

Input that represents the topographic roughness of the glacier. It is computed as the curvature of the glacier bed (or surface) over a specified window size. The curvature can be calculated in different directions (flow, cross-flow, or both) and using different curvature types (scalar or variability).

i∇S <: AbstractInput

Input that represents the surface slope in the SIA. It is computed using the bedrock elevation and the ice thickness solution H. The spatial differences are averaged over the opposite axis:

S = B + H
∇S = (avg_y(diff_x(S) / Δx) .^ 2 .+ avg_x(diff_y(S) / Δy) .^ 2) .^ (1/2)

ConstantA(A::F) where {F <: AbstractFloat}

Law that represents a constant A in the SIA.

Arguments:

• A::F: Rheology factor A.

CuffeyPaterson(; scalar::Bool = true)

Create a rheology law for the flow rate factor A. The created law maps the long term air temperature to A using the values from Cuffey & Peterson through polyA_PatersonCuffey that returns a polynomial which is then evaluated at a given temperature in the law.

Halfar(halfar_params::HalfarParameters) -> (_halfar::Function, t₀_years::Float64)

Constructs the Halfar similarity solution to the SIA for a radially symmetric ice sheet dome following Bueler (2005) "Exact solutions and verification of numerical models for isothermalice sheets"

Arguments

• halfar_params::HalfarParameters: A struct containing physical and geometric parameters for the Halfar solution, including dome height, margin radius, Glen exponent, and other constants.

Returns

• _halfar::Function: A function (x, y, t) -> H that evaluates the ice thickness H at position (x, y) and time t (in years).
• t₀_years::Float64: The characteristic timescale t₀ (in years) of the solution, based on the specified parameters.

Description

The solution has the form:

\[H(r, t) = H₀ (t / t₀)^(-α) [1 - ((t / t₀)^(-β) (r / R₀))^((n+1)/n)]^{n / (2n + 1)}\]

Halfar_velocity(halfar_params::HalfarParameters)

Same as Halfar(halfar_params), but instead of returning a function that gives the ice thickness as a function of space and time, this returns the ice surface velocity according to the Shallow Ice Approximation.

MB_diagnostics(simulation, iceflow_sol)

Mass balance accumulated over each mass balance window, and the times it's reported at.

Populates the MB and t_MB fields of Results. With mass balance in the ice flow right hand side, no snapshot is recorded during the solve, so it's rebuilt from the states saved on window edges with a midpoint rule,

\[MB_k ≈ ṁ\left(\frac{H_{k-1} + H_k}{2}, \frac{t_{k-1} + t_k}{2}\right) (t_k - t_{k-1})\]

on the window grid, not the result grid — MB and t_MB don't depend on step. Without mass balance in the right hand side (use_MB = false) nothing is recorded, so both are empty.

MB_solver_stops(simulation, tstops::Vector{F}) where {F <: AbstractFloat}

Times the solver must stop at, and times it must save at, for a given result grid tstops.

No callback saves the solution on the result grid, so saving has to be asked for explicitly: create_results reads the result grid out of the solution and fails on states that are missing. This applies to every run, with or without mass balance.

With mass balance in the right hand side, the solver must also stop at every window edge (ṁ is piecewise constant in time, so the right hand side jumps there), and those edge states are saved for MB_diagnostics to read.

SIA2D!(
    dH::Matrix{R},
    H::Matrix{R},
    simulation::SIM,
    t::Real,
    θ,
) where {R <:Real, SIM <: Simulation}

Simulates the evolution of ice thickness in a 2D shallow ice approximation (SIA) model. Works in-place.

Arguments

• dH::Matrix{R}: Matrix to store the rate of change of ice thickness.
• H::Matrix{R}: Matrix representing the ice thickness.
• simulation::SIM: Simulation object containing model parameters and state.
• t::Real: Current simulation time. Deliberately not tied to the element type of the state, so that the state can be dual or tracked while time stays a plain float.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Details

This function updates the ice thickness H and computes the rate of change dH using the shallow ice approximation in 2D. It retrieves necessary parameters from the simulation object, enforces positive ice thickness values, updates glacier surface altimetry and computes surface gradients. It then applies the necessary laws that are not updated via callbacks (A, C, n or U depending on the use-case) and computes the flux components, and flux divergence. Finally it adds the surface mass balance source term (see add_MB!).

Notes

• The function operates on a staggered grid for computing gradients and fluxes.
• Surface elevation differences are capped using upstream ice thickness to impose boundary conditions.
• The function modifies the input matrices dH and H in-place.

See also SIA2D

SyntheticC(params::Sleipnir.Parameters)

Creates a synthetic law for calculating the parameter C using a nonlinear sigmoid transformation based on the ratio of CPDD (cumulative positive degree days) to topo_roughness (topographic roughness). The law is parameterized by minimum and maximum values (Cmin, Cmax) from params.physical, and applies a sigmoid scaling to smoothly interpolate between these bounds.

TemperateA()

Law that represents a constant A in the SIA for temperate ice (0°C).

The value of A is set to 2.4e-24 s⁻¹Pa⁻³ which come from Cuffey & Peterson.

V_from_H(
    simulation::SIM,
    H::Matrix{F},
    t::Real,
    θ,
) where {F <: AbstractFloat, SIM <: Simulation}

Compute surface velocity from ice thickness using the SIA model. It relies on surface_V to compute Vx and Vy and it additionally computes the magnitude of the velocity V.

Arguments:

• simulation::SIM: The simulation structure used to retrieve the physical parameters.
• H::Matrix{F}: The ice thickness matrix.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Returns:

• Vx: x axis component of the surface velocity.
• Vy: y axis component of the surface velocity.
• V: Magnitude velocity.

add_MB!(dH, H, simulation, t::Real)

Add the surface mass balance rate ṁ(H, t) to the ice flow tendency dH, in-place.

Surface mass balance is a source term of the ice flow equation,

\[∂H/∂t = -∇·(D∇S) + ṁ(H, t)\]

so the solver integrates it together with the flux divergence and controls its error with the same tolerances. Because ṁ reads the surface elevation S = B + H, the elevation feedback ∂ṁ/∂H is part of the Jacobian and is picked up by the adjoint.

Ice flow right hand sides call this unconditionally; it returns immediately when the mass balance cache is inactive.

Arguments

• dH: Rate of change of ice thickness, updated in-place.
• H: Current ice thickness.
• simulation: Simulation object holding the model, the caches and the glaciers.
• t::Real: Current simulation time.

apply_MB_mask!(H, ifm::SIA2DCache)

Apply the mass balance (MB) mask to the iceflow model in-place. This function ensures that no MB is applied on the borders of the glacier to prevent overflow.

Arguments:

• H: Ice thickness.
• ifm::SIA2DCache: Iceflow cache of the SIA2D that provides the mass balance information and that is modified in-place.

averageV(
    θ,
    simulation,
    tspan_velocity::Tuple{F,F},
    step::F,
    ts::Vector{F},
    Hs::Vector{Matrix{F}},
) where {F <: AbstractFloat}

Compute average surface velocity from ice thickness over a given time window and using a constant time stepping. This function uses the SIA model to map the ice thickness to instant surface velocity.

Arguments:

• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.
• simulation::SIM: The simulation structure used to retrieve the physical parameters.
• tspan_velocity::Tuple{F,F}: Time window over which the average surface velocity is computed. It can be different than the simulation period.
• step::F: Constant time step used to compute the average surface velocity.
• ts::Vector{F}: Vector representing time associated to each ice thickness matrix.
• Hs::Vector{Matrix{F}}: Vector ice thickness matrices.
• H::Matrix{F}: The ice thickness matrix.
• t::R: Current simulation time.

Returns:

• avg_Vx: x axis component of the average surface velocity.
• avg_Vy: y axis component of the average surface velocity.
• avg_V: Magnitude velocity.

avg!(O, I)

Compute the average of adjacent elements in the input array I and store the result in the output array O.

Arguments

• O: Output array where the averaged values will be stored.
• I: Input array from which the adjacent elements will be averaged.

Details

This function uses the @views macro to avoid creating temporary arrays and the @. macro to broadcast the operations. The averaging is performed by taking the mean of each 2x2 block of elements in I and storing the result in the corresponding element in O.

avg(A::AbstractArray)

Compute the average of each 2x2 block in the input array A. The result is an array where each element is the average of the corresponding 2x2 block in A.

Arguments

• A::AbstractArray: A 2D array of numerical values.

Returns

• A 2D array of the same type as A, where each element is the average of a 2x2 block from A.

avg_surface_V!(simulation::SIM, t::R, θ) where {SIM <: Simulation, R <: Real}

Calculate the average surface velocity for a given simulation.

Arguments

• simulation::SIM: A simulation object of type SIM which is a subtype of Simulation.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Description

This function computes the average surface velocity components (Vx and Vy) and the resultant velocity (V) for the ice flow model within the given simulation. It first calculates the surface velocities at the initial and current states, then averages these velocities and updates the ice flow model's velocity fields.

Notes

• The function currently uses a simple averaging method and may need more datapoints for better interpolation. # TODO: Add more datapoints to better interpolate this

avg_x!(O, I)

Compute the average of adjacent elements along the first dimension of array I and store the result in array O.

Arguments

• O: Output array where the averaged values will be stored.
• I: Input array from which adjacent elements will be averaged.

avg_x(A::AbstractArray)

Compute the average of adjacent elements along the first dimension of the array A.

Arguments

• A::AbstractArray: Input array.

Returns

• An array of the same type as A with one less element along the first dimension, containing the averages of adjacent elements.

avg_y!(O, I)

Compute the average of adjacent elements along the second dimension of array I and store the result in array O.

Arguments

• O: Output array where the averaged values will be stored.
• I: Input array from which the adjacent elements will be averaged.

avg_y(A::AbstractArray)

Compute the average of adjacent elements along the second dimension of the input array A.

Arguments

• A::AbstractArray: An array of numeric values.

Returns

• An array of the same type as A containing the averages of adjacent elements along the second dimension.

batch_iceflow_PDE!(glacier_idx::I, simulation::Prediction) where {I <: Integer}

Solve the Shallow Ice Approximation iceflow PDE in-place for a given set of laws prescribed in the simulation object. It creates the iceflow problem, the necessary callbacks and solve the PDE.

Arguments:

• glacier_idx::I: Integer ID of the glacier.
• simulation::Prediction: Simulation object that contains all the necessary information to solve the iceflow.

Returns

• A Results instance that stores the iceflow solution.

build_callback(model::SIA2Dmodel, cache::SIA2DCache, glacier_idx::Real, θ, tspan) -> CallbackSet

Return a CallbackSet that updates the cached values of A, C, n and U at provided time intervals.

Each law can optionally specify a callback frequency via callback_freq.

• If callback_freq > 0, a PeriodicCallback is used to update the corresponding component at regular intervals.
• If callback_freq == 0, a PresetTimeCallback is used to trigger the update only at the initial time (taken from tspan[1]).
• If no callback is specified for a component, a dummy CallbackSet is returned.

Arguments:

• model::SIA2Dmodel: The ice flow model definition.
• cache::SIA2DCache: Model cache for efficient component access and updates.
• glacier_idx::Real: Index of the glacier in the simulation.
• θ: Optional parameter for law evaluation.
• tspan: Tuple or floats specifying the simulation time span. Used to determine initial callback time when freq == 0.

Returns:

• A CallbackSet containing all the callbacks for periodic or preset updates of model components.

d2dx(f::Matrix{T}, i::Int, j::Int, Δx::Float64) where T <: Real

Compute the second central difference in the x-direction at (i,j).

d2dxy(f::Matrix{T}, i::Int, j::Int, Δx::Float64, Δy::Float64) where T <: Real

Compute the mixed second central difference (∂²f/∂x∂y) at (i,j).

d2dy(f::Matrix{T}, i::Int, j::Int, Δy::Float64) where T <: Real

Compute the second central difference in the y-direction at (i,j).

define_callback_steps(tspan::Tuple{F, F}, step::F) where {F <: AbstractFloat}

Defines the times to stop for the DiscreteCallback given a step and a timespan.

Arguments

• tspan::Tuple{F, F}: A tuple representing the start and end times.
• step::F: The step size for generating the callback steps.

Returns

• Vector{F}: A vector of callback steps within the specified time span.

diff_x!(O, I, Δx)

Compute the finite difference of array I along the first dimension and store the result in array O. The difference is computed using the spacing Δx.

Arguments

• O: Output array to store the finite differences.
• I: Input array from which finite differences are computed.
• Δx: Spacing between points in the first dimension.

Notes

• The function uses @views to avoid copying data when slicing arrays.
• The operation is performed in-place, modifying the contents of O.

diff_x(A::AbstractArray)

Compute the difference along the first dimension of the array A.

Arguments

• A::AbstractArray: Input array.

Returns

• An array of the same type as A containing the differences along the first dimension.

diff_y!(O, I, Δy)

Compute the finite difference along the y-axis and store the result in O.

Arguments

• O: Output array where the result will be stored.
• I: Input array from which the finite difference is computed.
• Δy: The spacing between points in the y-direction.

Description

This function calculates the finite difference along the y-axis for the input array I and stores the result in the output array O. The calculation is performed using the formula:

O = (I[:,begin+1:end] - I[:,1:end - 1]) / Δy

The @views macro is used to avoid copying data when slicing the array.

diff_y(A::AbstractArray)

Compute the difference between adjacent elements along the second dimension (columns) of the input array A.

Arguments

• A::AbstractArray: An array of numeric values.

Returns

• An array of the same type as A containing the differences between adjacent elements along the second dimension.

effective_abstol(abstol, tspan; verbose = true)

Tighten abstol for long runs. abstol is given for a run of ABSTOL_REFERENCE_YEARS (5 years), and it is divided by how many times longer the run is: a 30 year run uses abstol / 6. Shorter runs keep abstol as it is. We tried loosening it for them too, but that made short inversions and gradient checks noticeably less accurate.

We do this because the solver error grows with the length of the run. On RGI60-11.03638, with abstol = 1e-3, the RMS error in H was 0.016, 0.129 and 0.223 m after 5, 15 and 30 years, compared to a converged reference (RDPK3Sp35, abstol = 1e-10, fixed dt = 1/2400). Between 5 and 30 years this is roughly t^1.4. Tightening abstol helps less than one would expect, since the error only went down as abstol^0.4. This is probably because the mass balance is not smooth in H (linear interpolation in the elevation lookup table and accumulation ramp). With this rescaling, the error grows roughly linearly with the length of the run. Keeping it constant would need a tolerance hundreds of times smaller for 30 years, which would be too expensive.

These numbers come from a single glacier. If you need a given accuracy, set abstol yourself and scale_abstol = false in SolverParameters.

The adjustment is logged, use verbose = false to silence it.

generate_ground_truth(
    glaciers::Vector{G},
    params::Sleipnir.Parameters,
    model::Sleipnir.Model,
    tstops::Vector{F};
    store::Tuple=(:H, :V, :dhdt),
) where {G <: Sleipnir.AbstractGlacier, F <: AbstractFloat}

Generate ground truth data for a glacier simulation by using the laws specified in the model and running a forward model. It returns a new vector of glaciers with updated thicknessData, velocityData and dhdtData fields based on the store argument.

Arguments

• glaciers::Vector{G}: A vector of glacier objects of type G, where G is a subtype of Sleipnir.AbstractGlacier.
• params::Sleipnir.Parameters: Simulation parameters.
• model::Sleipnir.Model: The model to use for the simulation.
• tstops::Vector{F}: A vector of time steps at which the simulation will be evaluated.
• store::Tuple: Which generated simulation products to store. It can include :H, :V, :avgV and/or :dhdt.

Description

1. Runs a forward model simulation for the glaciers using the provided laws, parameters, model, and time steps.
2. Build a new vector of glaciers and store the simulation results as ground truth in the glaciers struct. For each glacier it populates
   • thicknessData field if store contains :H,
   • velocityData if store contains :V or :avgV,
   • dhdtData if store contains :dhdt.

Example

glaciers = [glacier1, glacier2] # dummy example
params = Huginn.Parameters() # to be filled
model = Huginn.Model() # to be filled
tstops = 0.0:1.0:10.0

glaciers = generate_ground_truth(glaciers, params, model, tstops)

generate_ground_truth_prediction(
    glaciers::Vector{G},
    params::Sleipnir.Parameters,
    model::Sleipnir.Model,
    tstops::Vector{F},
) where {G <: Sleipnir.AbstractGlacier, F <: AbstractFloat}

Wrapper for generate_ground_truth that also updates the glaciers field of the Prediction object.

Arguments

• glaciers::Vector{G}: A vector of glacier objects of type G, where G is a subtype of Sleipnir.AbstractGlacier.
• params::Sleipnir.Parameters: Simulation parameters.
• model::Sleipnir.Model: The model to use for the simulation.
• tstops::Vector{F}: A vector of time steps at which the simulation will be evaluated.

Description

This function calls generate_ground_truth to generate ground truth data for the glaciers using the provided laws, parameters, model, and time steps. In addition, it updates the glaciers field of the Prediction object with the newly generated glaciers containing the ground truth data.

Example

glaciers = [glacier1, glacier2] # dummy example
params = Huginn.Parameters() # to be filled
model = Huginn.Model() # to be filled
tstops = 0.0:1.0:10.0

prediction = generate_ground_truth_prediction(glaciers, params, model, tstops)

generate_result(placeholder_sim::SIM, A, n) where {SIM <: Simulation}

Generate the result of a simulation by initializing the model with the specified parameters and running the simulation.

Arguments

• simulation::SIM: An instance of a type that is a subtype of Simulation.
• A: The parameter to set for simulation.model.iceflow.A.
• n: The parameter to set for simulation.model.iceflow.n.

Returns

• result: The first result from the simulation's results.

inn(A::AbstractArray)

Extracts the inner part of a 2D array A, excluding the first and last rows and columns.

Arguments

• A::AbstractArray: A 2D array from which the inner part will be extracted.

Returns

• A subarray of A containing all elements except the first and last rows and columns.

inn1(A::AbstractArray)

Returns a view of the input array A excluding the last row and the last column.

Arguments

• A::AbstractArray: The input array from which a subarray view is created.

Returns

• A view of the input array A that includes all elements except the last row and the last column.

polyA_PatersonCuffey()

Returns a function of the coefficient A as a polynomial of the temperature. The values used to fit the polynomial come from Cuffey & Peterson.

precompute_all_VJPs_laws!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation::Prediction,
    glacier_idx::Integer,
    t::Real,
    θ,
)

Function that does nothing and its existence is just to support multiple dispatch. The implementation that is useful is available in ODINN when simulation is a Inversion object.

project_curvatures(H, eₚ, eₛ)

Computes the scalar second derivative of the surface in a specific direction.

Arguments

• H: Hessian matrix (2x2).
• eₚ: Principal direction vector 1 (2x1).
• eₛ: Principal direction vector 2 (2x1).

Returns

• Kₚ: Curvature in the direction of eₚ.
• Kₛ: Curvature in the direction of eₛ.

run!(simulation::Prediction)

In-place run of the model.

simulate_iceflow_PDE!(
    simulation::SIM,
    cb::SciMLBase.DECallback,
    du,
    tstops::Vector{F},
) where {SIM <: Simulation, F <: AbstractFloat}

Make forward simulation of the iceflow PDE determined in du in-place and create the results.

Results are sampled on tstops. When mass balance runs in the right hand side, the solver is also forced to stop on the mass balance window edges — discontinuities of the right hand side, not result time steps.

spectral_radius(simulation)

Estimate of the spectral radius of the ice flow right hand side, in yr⁻¹. It is not a guaranteed upper bound, see the warning in with_eigen_est.

Two terms contribute: the flux divergence, a diffusion estimated as 4 D (1/Δx² + 1/Δy²) with the largest D, and the mass balance source's diagonal Jacobian ∂ṁ/∂H.

The mass balance term isn't a correction — where ice is thin the ramp makes ∂ṁ/∂H large, and where A is small D is negligible, so mass balance can be the only term. Dropping it underestimates the spectral radius and the stabilised solver picks too few stages.

The mass balance term is a bound from the lookup table, not the exact maximum at the current state: this runs inside the solver, where the adjoint hands back an augmented [H; θ] vector that a state-dependent version couldn't read anyway.

surface_V!(H::Matrix{<:Real}, simulation::SIM, t::R, θ) where {SIM <: Simulation, R <: Real}

Compute the surface velocities of a glacier using the Shallow Ice Approximation (SIA) in 2D.

Arguments

• H::Matrix{<:Real}: The ice thickness matrix.
• simulation::SIM: The simulation object containing parameters and model information.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Returns

• Vx: The x-component of the surface velocity.
• Vy: The y-component of the surface velocity.

Description

This function updates the glacier surface altimetry and computes the surface gradients on edges using a staggered grid. It then calculates the surface velocities based on the Shallow Ice Approximation (SIA) model.

Details

• params: The simulation parameters.
• iceflow_model: The ice flow model from the simulation.
• glacier: The glacier object from the simulation.
• B: The bedrock elevation matrix.
• H̄: The average ice thickness matrix.
• dSdx, dSdy: The surface gradient matrices in x and y directions.
• ∇S, ∇Sx, ∇Sy: The gradient magnitude and its components.
• Γꜛ: The surface stress.
• D: The diffusivity matrix.
• A: The flow rate factor.
• n: The flow law exponent.
• C: The sliding coefficient.
• Δx, Δy: The grid spacing in x and y directions.
• ρ: The ice density.
• g: The gravitational acceleration.

The function computes the surface gradients, averages the ice thickness, and calculates the surface stress and diffusivity. Finally, it computes the surface velocities Vx and Vy based on the gradients and diffusivity.

surface_V(
    H::Matrix{R},
    simulation::SIM,
    t::Real,
    θ,
) where {R <: Real, SIM <: Simulation}

Compute the surface velocities of a glacier using the Shallow Ice Approximation (SIA) in 2D.

Arguments

• H::Matrix{R}: Ice thickness matrix.
• simulation::SIM: Simulation object containing parameters and model information.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Returns

• Vx: Matrix of surface velocities in the x-direction.
• Vy: Matrix of surface velocities in the y-direction.

Details

This function computes the surface velocities of a glacier by updating the glacier surface altimetry and calculating the surface gradients on the edges. It uses a staggered grid approach to compute the gradients and velocities.

Notes

• The function assumes that the simulation object contains the necessary parameters and model information.

thickness_velocity_data(
    prediction::Prediction,
    tstops::Vector{F};
    store::Tuple=(:H, :V, :dhdt),
) where {F <: AbstractFloat}

Return a new vector of glaciers with the simulated thickness ice velocity and dhdt data for each of the glaciers.

Arguments

• prediction::Prediction: A Prediction object containing the simulation results and associated glaciers.
• tstops::Vector{F}: A vector of time steps (of type F <: AbstractFloat) at which the simulation was evaluated.
• store::Tuple: Which generated simulation products to store. It can include :H, :V, :avgV and/or :dhdt.

Description

This function iterates over the glaciers in the Prediction object and generates the simulated data based on the store argument at corresponding time steps (t). If store includes :H, then the ice thickness is stored. If store includes :V or :avgV, then it computes the surface ice velocity data and stores it. These two options are mutually exclusive. When :avgV is provided, only one snapshot of ice surface velocity is computed and it corresponds to the ice surface velocity at the closest time to (tspan[1]+tspan[2])/2. If store includes :dhdt, then it computes the mean surface elevation change and stores it. A new vector of glaciers is created and each glacier is a copy with an updated thicknessData, velocityData and dhdtData fields.

Notes

• The function asserts that the time steps (ts) in the simulation results match the provided tstops. If they do not match, an error is raised.

Returns

A new vector of glaciers where each glacier is a copy of the original one with the updated thicknessData, velocityData and dhdtData based on the values provided in store.

with_eigen_est(alg, simulation)

Give ROCK2 and ROCK4 our own estimate of the spectral radius, computed with spectral_radius, instead of letting them compute it. Other solvers are returned unchanged. It is only used when supply_eigen_est = true in SolverParameters.

By default, ROCK estimates the spectral radius with a power iteration. The result of this iteration changes slightly between two close parameter values, and so can the number of stages used by the solver. This doesn't matter for a forward run, but it adds noise to the loss, which breaks finite difference gradient checks. Our estimate doesn't have this problem.

This is not a speed-up. On RGI60-11.03638 over 2010–2015, the solver needed 1130 RHS evaluations with our estimate and 1157 without it. Multiplying our estimate by 2 or 5 to be on the safe side costs 25% and 83% more evaluations.

∇slope(S::Matrix{T}, Δx::T, Δy::T) where T <: Real

Compute the magnitude of the surface slope ∇S for a scalar field S defined on a rectilinear grid.

apply_all_callback_laws!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation,
    glacier_idx::Integer,
    t::Real,
    θ,
)

Applies the different laws required by the SIA2D glacier model for a given glacier and simulation state. If U_is_provided is true in SIA2D_model and U is a callback law, it applies the law for U only. Otherwise if Y_is_provided and Y is a callback law, it applies the law for Y only. Finally, if U_is_provided and Y_is_provided are false, the function checks and applies the laws for A, C, and n, if they are defined as "callback" laws (i.e., handled as callbacks by the ODE solver). Results are written in-place to the cache for subsequent use in the simulation step.

Arguments

• SIA2D_model: The model object containing the laws (A, C, n, Y and U).
• SIA2D_cache: A cache object to store the evaluated values of the laws (A, C, n, Y and U) for the current step.
• simulation: The simulation object.
• glacier_idx::Integer: Index of the glacier being simulated, used to select data for multi-glacier simulations.
• t::Real: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Notes

• The function mutates the contents of SIA2D_cache.
• Only "callback" laws are applied.
• This function is typically used in the manual adjoint and in the tests where only portions of the code are applied and we need to apply all the laws used in the iceflow model.

function apply_all_non_callback_laws!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation,
    glacier_idx::Integer,
    t::Real,
    θ,
)

Applies the different laws required by the SIA2D glacier model for a given glacier and simulation state. If U_is_provided is true in SIA2D_model and U is not a callback law, it applies the law for U only. Otherwise if Y_is_provided and Y is not a callback law, it applies the law for Y only. Finally, if U_is_provided and Y_is_provided are false, the function checks and applies the laws for A, C, and n, unless they are defined as "callback" laws (i.e., handled as callbacks by the ODE solver). Results are written in-place to the cache for subsequent use in the simulation step.

Arguments

• SIA2D_model: The model object containing the laws (A, C, n, Y and U).
• SIA2D_cache: A cache object to store the evaluated values of the laws (A, C, n, Y and U) for the current step.
• simulation: The simulation object.
• glacier_idx::Integer: Index of the glacier being simulated, used to select data for multi-glacier simulations.
• t::Real: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Notes

• The function mutates the contents of SIA2D_cache.
• "Callback" laws are skipped, as they are expected to be handled outside this function.
• This function is typically called at each simulation time step for each glacier.

function initcache( iceflowmodel::SIA2Dmodel, glacier::AbstractGlacier, glacier_idx::I, θ, ) where {IF <: IceflowModel, I <: Integer}

Initialize iceflow model data structures to enable in-place mutation.

Keyword arguments

• iceflow_model: Iceflow model used for simulation.
• glacier_idx: Index of glacier.
• glacier: Glacier to provide basic initial state of the ice flow model.
• θ: Optional parameters of the laws.