API

AbstractAdjointMethod

Abstract type representing the flavor of AD and adjoint to be used to compute the gradient of the cost function. There are two parts where one can play with how the gradient is propagated: the iceflow model VJP and the adjoint of the ODE solver. The VJP of the iceflow model can be computed using either AD (Zygote or Enzyme), the discrete, or the continuous adjoint of the iceflow model. As for the computation of the adjoint of the ODE solution, it can be handled by SciMLSensitivity, or computed using the adjoint implemented in ODINN.

AbstractVJPMethod

Abstract type representing the flavor of AD to be used to compute the VJP inside the gradient of the cost function.

ContinuousAdjoint{
    F <: AbstractFloat,
    I <: Integer,
    VJP <: AbstractVJPMethod,
    MBVJP <: AbstractVJPMethod,
} <: AbstractAdjointMethod

Continuous adjoint of SIA2D with manual implementation of the backward in the ODE scheme.

Fields

• VJP_method::VJP: Type of AbstractVJPMethod used to compute VJPs inside adjoint calculation.
• solver::Any: The solver to be used for adjoint.
• reltol::F: Relative tolerance to be used in the ODE solver of the adjoint.
• abstol::F: Absolute tolerance to be used in the ODE solver of the adjoint.
• dtmax::F: Maximum time step to be used in the ODE solver of the adjoint.
• interpolation: Interpolation method to be used to interpolate the variables in the computation of the adjoint. Currently only :Linear is supported.
• n_quadrature::I: Number of nodes used in the Gauss quadrature for the numerical integration of the loss function.
• MB_VJP::MBVJP: Type of AbstractVJPMethod used to compute the MB VJP inside adjoint calculation.

ContinuousVJP{ADTYPE <: DI.AbstractADType} <: AbstractVJPMethod

Continuous manual implementation of the VJP required inside the adjoint calculation. It relies in the continuous expresion for the adjoint operation based on the functional formula of the forward PDE.

Fields

• regressorADBackend::ADTYPE: Specifies the AD backend used for the laws when their associated VJPs functions are not provided. The type parameter ADTYPE must be a subtype of DI.AbstractADType.

Constructor

• The default constructor allows specifying the backend via the regressorADBackend keyword argument, defaulting to DI.AutoMooncake().

DiffusivityRegularization(; reg = TikhonovRegularization())

Regularization for diffusivity fields using a specified spatial operator.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator applied to diffusivity.

DiscreteAdjoint{
    VJP <: AbstractVJPMethod,
    MBVJP <: AbstractVJPMethod,
} <: AbstractAdjointMethod

Discrete adjoint of SIA2D with manual implementation of the backward in the ODE scheme.

Fields

• VJP_method: Type of AbstractVJPMethod used to compute VJPs inside adjoint calculation.
• MB_VJP::MBVJP: Type of AbstractVJPMethod used to compute the MB VJP inside adjoint calculation.

DiscreteVJP{ADTYPE <: DI.AbstractADType} <: AbstractVJPMethod

Discrete manual implementation of the VJP required inside the adjoint calculation. This implements the pullback function for the function to differentiate.

Fields

• regressorADBackend::ADTYPE: Specifies the AD backend used for the laws when their associated VJPs functions are not provided. The type parameter ADTYPE must be a subtype of DI.AbstractADType.

Constructor

• The default constructor allows specifying the backend via the regressorADBackend keyword argument, defaulting to DI.AutoMooncake().

Struct to provide a dummy gradient. It does not have to be the true gradient. Mainly used to test that the optimization pileline works independenly of the gradient calculation.

DummyAdjoint

Fields:

• grad::Function: In-place function f(du, u; kwargs) that fills the first argument du with the gradient values.

Enzyme implementation of VJP used inside the adjoint calculation.

EnzymeVJP

FunctionalModel <: TrainableModel

Abstract type representing functional learnable components of the model. This is a subtype of TrainableModel. Typically used for functional inversions.

GlacierWideInv{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier invertible parameter container. GlacierWideInv wraps a ComponentVector (θ) that stores one scalar parameter per glacier and implements the PerGlacierModel interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier parameter vector (one scalar value per glacier).

Constructor

GlacierWideInv(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    var::Symbol,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• var::Symbol: Symbol naming the field on each glacier to use as the initial value.

Example

GlacierWideInv(params, glaciers, :A)

GriddedInv{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier invertible parameter container. GriddedInv wraps a ComponentVector (θ) that stores one matrix parameter per glacier and implements the PerGlacierModel interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier parameter vector (one matrix per glacier).

Constructor

GriddedInv(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    var::Symbol,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• var::Symbol: Symbol naming the field on each glacier to use as the initial value.

Example

GriddedInv(params, glaciers, :A)

Hyperparameters(; current_epoch::Int64 = 1, current_minibatch::Int64 = 1, loss_history::Vector{Float64} = Vector{Float64}(), optimizer::Union{Optim.FirstOrderOptimizer, Flux.Optimise.AbstractOptimiser, Optimisers.AbstractRule} = BFGS(initial_stepnorm=0.001), loss_epoch::Float64 = 0.0, epochs::Int64 = 50, batch_size::Int64 = 15)

Constructs a Hyperparameters object with the specified parameters.

Arguments

• current_epoch::Int64: The current epoch number. Defaults to 1.
• current_minibatch::Int64: The current minibatch number. Defaults to 1.
• loss_history::Vector{Float64}: A vector to store the history of loss values. Defaults to an empty vector.
• optimizer::Union{Optim.FirstOrderOptimizer, Flux.Optimise.AbstractOptimiser, Optimisers.AbstractRule}: The optimizer to be used. Defaults to BFGS(initial_stepnorm=0.001).
• loss_epoch::Float64: The loss value for the current epoch. Defaults to 0.0.
• epochs::Int64: The total number of epochs. Defaults to 50.
• batch_size::Int64: The size of each minibatch. Defaults to 15.

Returns

• A Hyperparameters object initialized with the provided values.

mutable struct Hyperparameters{F <: AbstractFloat, I <: Integer} <: AbstractParameters

A mutable struct that holds hyperparameters for training a machine learning model.

Keyword arguments

• current_epoch::I: The current epoch number.
• current_minibatch::I: The current minibatch number.
• loss_history::Vector{F}: A vector storing the history of loss values.
• optimizer::Union{Optim.FirstOrderOptimizer, Flux.Optimise.AbstractOptimiser, Optimisers.AbstractRule}: The optimizer used for training.
• loss_epoch::F: The loss value for the current epoch.
• epochs::I: The total number of epochs for training.
• batch_size::I: The size of each minibatch.

InitialCondition{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier initial condition container. InitialCondition wraps a ComponentVector (θ) that stores one matrix per glacier and implements the InitialCondition interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier matrix of initial condition.

Constructor

InitialCondition(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    initialization::Symbol = :Farinotti2019,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• initialization::Symbol: Symbol providing the way the initial condition should be initialized.

Example

InitialCondition(params, glaciers, :Farinotti2019)

InitialThicknessRegularization(; reg = TikhonovRegularization(), t₀ = 1994.0)

A composite regularization type designed for initial ice thickness. It combines a simple spatial regularization (e.g., TikhonovRegularization) with a reference initial time.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): The spatial regularization operator applied to the initial field. By default, a Tikhonov (Laplacian-based) regularization is used.
• t₀::AbstractFloat = 1994.0: The reference initial time (e.g., year) at which the regularization applies.

function Inversion(
    model::M,
    glaciers::Vector{G},
    parameters::P
) where {G <: Sleipnir.AbstractGlacier, M <: Sleipnir.Model, P <: Sleipnir.Parameters}

Constructor for Inversion struct with glacier model information, glaciers, and parameters.

Arguments

• model::Sleipnir.Model: The model used for the simulation.
• glaciers::Vector{G}: A vector of glaciers involved in the simulation.
• parameters::Sleipnir.Parameters: The parameters used for the simulation.

Returns

• Inversion: A new instance of the Inversion struct.

mutable struct Inversion{MODEL, CACHE, GLACIER, RES} <: Simulation

An object representing an inversion simulation. It can involve at the same time a classical inversion and a functional inversion (i.e. the inversion of a function using some data-driven regressor).

Fields

• model::Sleipnir.Model: The model used for the simulation.
• glaciers::Vector{Sleipnir.AbstractGlacier}: A vector of glaciers involved in the simulation.
• parameters::Sleipnir.Parameters: The parameters used for the simulation.
• results::ODINN.Results: A ODINN.Results instance to store the results of the inversion and of the forward simulations.

InversionBinder{FI <: Inversion, CA <: ComponentArray}

Struct used for the binding with SciMLSensitivity. It is defined as a SciMLStructure and it contains the inversion structure and the vector of parameters to differentiate.

Fields

• simulation::FI: Inversion instance.
• θ::CA: ComponentArray that contains the parameters used to differentiate the iceflow.

InversionParameters{F<:AbstractFloat}(;
    initial_conditions::Vector{F} = [1.0],
    lower_bound::Vector{F} = [0.0],
    upper_bound::Vector{F} = [Inf],
    regions_split::Vector{Int} = [1, 1],
    x_tol::F = 1.0e-3,
    f_tol::F = 1.0e-3,
    solver = BFGS()
)

Initialize the parameters for the inversion process.

Arguments

• initial_conditions::Vector{F}: Starting point for optimization.
• lower_bound::Vector{F}: Lower bounds for optimization variables.
• upper_bound::Vector{F}: Upper bounds for optimization variables.
• regions_split::Vector{Int}: Defines the amount of region split based on altitude and distance to border for the inversion process.
• x_tol::F: Tolerance for variables convergence.
• f_tol::F: Tolerance for function value convergence.
• solver: Optimization solver to be used.

InversionParameters{F<:AbstractFloat}

A mutable struct that holds parameters for inversion processes. This struct is a subtype of AbstractParameters.

Fields

• initial_conditions::Vector{F}: A vector of initial conditions.
• lower_bound::Vector{F}: A vector specifying the lower bounds for the parameters.
• upper_bound::Vector{F}: A vector specifying the upper bounds for the parameters.
• regions_split::Vector{Int}: A vector indicating how the regions are split.
• x_tol::F: The tolerance for the solution's x-values.
• f_tol::F: The tolerance for the function values.
• solver::Any: The solver to be used for the inversion process.

L2Sum{I <: Integer} <: AbstractSimpleLoss

Struct that defines an L2 sum loss. The sum is defined only on pixels inside the glacier. The parameter distance controls the pixels that should be used to compute the sum. In order for a pixel to be used, it should be at least at distance from the glacier border. The mask defining the glacier borders are computed using the ground truth ice thickness.

$loss(a,b) = sum_{i\in\text{mask}} (a[i]-b[i])^2 / normalization$

Fields

• distance::I: Distance to border.

LogSum{I <: Integer, F <: AbstractFloat} <: AbstractSimpleLoss

Struct that defines a Logarithmic sum loss.

$loss(a,b) = log^2( (a + ϵ) / (b + ϵ) ) / normalization$

Fields

• distance::I: Distance to border.
• ϵ::F: Epsilon used inside the loss function to handle division by zero and log(0). It somehow represents the minimum value the loss function should be sensible to.

LossH{L <: AbstractSimpleLoss} <: AbstractLoss

Struct that defines the ice thickness loss.

Fields

• loss::L: Type of loss to use for the ice thickness. Default is L2Sum().

LossHV{
    F <: AbstractFloat,
    LH <: AbstractLoss,
    LV <: AbstractLoss,
} <: AbstractLoss

Struct that defines the ice thickness and ice velocity loss. It consists of two fields that define the ice thickness and ice velocity loss. It also has a scaling coefficient that balances the ice velocity term in the loss.

$loss(\hat H,H) = loss_H(\hat H,H) + scaling * loss_V(\hat V,V)$

with $\hat V$ computed from $\hat H$ for the SIA.

Fields

• hLoss::LH: Type of loss to use for the ice thickness. Default is LossH().
• vLoss::LV: Type of loss to use for the ice velocity. Default is LossV().
• scaling::F: Scaling of the ice velocity term. Default is 1.

LossV{L <: AbstractSimpleLoss} <: AbstractLoss

Struct that defines the ice velocity loss.

Fields

• loss::L: Type of loss to use for the ice velocity. Default is L2Sum().
• component::Symbol: Component of the velocity field used in the loss. Options include :xy for both x and y component, and :abs for the norm/magnitude of the velocity.
• scale_loss::Bool: Whether to scale the loss function with the reference ice velocity magnitude.

MatrixCacheGlacierId <: Cache

A cache structure for storing a matrix value (Float64 2D array) along with their associated vector-Jacobian products (VJP). It also stores the glacier ID as an integer. This is typically used to invert a spatially varying field per glacier. Fields:

• value::Array{Float64, 2}: The cached matrix value.
• vjp_inp::Array{Float64, 2}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.
• glacier_id::Int64: Glacier ID in the list of simulated glaciers.

MatrixCacheInterp(nodes_H, nodes_∇S, interp_θ)

A mutable cache structure for storing interpolation data on a 2D grid, used to efficiently evaluate and reuse interpolated matrices and their gradients. This interpolation makes complex inversions feasible since it allows the precomputation of all gradients before the solving the reverse PDE associated to the adjoint variable.

Fields

• value::Array{Float64, 2}: Matrix to store value of the diffusivity.
• nodes_H::Vector{Float64}: Grid nodes corresponding to the first interpolation dimension, typically representing values of ice thickness H.
• nodes_∇S::Vector{Float64}: Grid nodes corresponding to the second interpolation dimension, typically representing absolute values of slope ∇S.
• interp_θ::Interpolations.GriddedInterpolation{Vector{Float64}, 2, Matrix{Vector{Float64}}, Interpolations.Gridded{InterPolations.Linear{InterPolations.Throw{OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}: A gridded linear interpolation object mapping (nodes_H, nodes_∇S) to parameter vectors θ. Used to compute interpolated parameters and enable fast evaluation during repeated model calls.

MultiLoss(; losses = (L2Sum(),), λs = (1.0,))

Combines multiple loss functions into a single weighted objective.

MultiLoss enables composing several individual loss terms—each possibly representing a different physical constraint, data fidelity term, or regularization penalty—into a single differentiable loss function.

Keyword Arguments (Constructor)

• losses::Tuple = (L2Sum(),): A tuple of loss objects (each subtype of AbstractLoss) to be combined.
• λs::Tuple = (1.0,): A tuple of scalar weights or hyperparameters corresponding to each loss term.

Fields (Struct)

• losses::TL: Tuple of loss functions.
• λs::TS: Tuple of weighting coefficients.

NeuralNetwork{
    ChainType <: Lux.Chain,
    ComponentVectorType <: ComponentVector,
    NamedTupleType <: NamedTuple,
} <: FunctionalModel

Feed-forward neural network.

Fields

• architecture::ChainType: Flux.Chain neural network architecture
• θ::ComponentVectorType: Neural network parameters
• st::NamedTupleType: Neural network status

No VJP flavor when the contribution of a given term should not be computed inside the adjoint calculation (e.g. MB).

NoVJP

PerGlacierModel <: TrainableModel

Abstract type representing per glacier optimizable components of the model. This is a subtype of TrainableModel. Typically used for classical inversions.

mutable struct Results{RES <: Sleipnir.Results, STAT <: TrainingStats}

Mutable struct containing the results of an inversion. This object stores both the results of the optimization through TrainingStats and the simulation results of the forward evaluations using the optimized variables through Sleipnir.Results. It expands the functionalities offered by Sleipnir.Results to allow saving the results of an inversion. Multiple dispatch is used to select either Sleipnir.Results or ODINN.Results.

Fields

• simulation::Vector{RES}: Vector of Sleipnir.Results representing the results of a forward evaluation for each glacier.

• stats::STAT: Training statistics gathered during the optimization.

  function Results( simulation::Vector{<: Sleipnir.Results}, stats::TrainingStats, )

Constructor for the Results object used to track statistics during training and the results of the forward evaluations simulated with the optimized variables.

Arguments

• simulation::Vector{<: Sleipnir.Results}: Vector of Sleipnir.Results associated to the forward simulation of each glacier.
• stats::TrainingStats: Training statistics.

RheologyRegularization(; reg = TikhonovRegularization())

Regularization of the gridded rheology A in the context of classical inversions. It can include a spatial smoothing operator through the field reg.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator.

SIA2D_D_target(; interpolation=:None, n_interp_half=20,
                 prescale=nothing, postscale=nothing)

Inversion of general diffusivity as a function of physical parameters.

D(H, ∇S, θ) = H * NN(H, ∇S; θ)

So now we are learning the velocity field given by D * ∇S. This inversion is similar to learnign the velocity field assuming that this is parallel to the gradient in surface ∇S.

Arguments

• interpolation::Symbol = :None: Specifies the interpolation method. Options include :Linear, :None.
• n_interp_half::Int = 20: Half-width of the interpolation stencil. Determines resolution of interpolation.
• prescale::Union{Fin, Nothing} = nothing: Optional prescaling function or factor applied before parametrization. Must be of type Fin or nothing.
• postscale::Union{Fout, Nothing} = nothing: Optional postscaling function or factor applied after parametrization. Must be of type Fout or nothing.

Type Parameters

• Fin: Type of the prescale function or operator.
• Fout: Type of the postscale function or operator.

Supertype

• AbstractSIA2DTarget: Inherits from the abstract target type for 2D SIA modeling.

Returns

• An instance of SIA2D_D_target configured with optional scaling and interpolation parameters.

ScalarCacheGlacierId <: Cache

A cache structure for storing a scalar value as a zero-dimensional array of Float64 along with their associated vector-Jacobian products (VJP). It also stores the glacier ID as an integer. This is typically used to invert a single scalar per glacier. Fields:

• value::Array{Float64, 0}: The cached scalar value.
• vjp_inp::Array{Float64, 0}: VJP with respect to inputs. Must be defined but never used in practice since this cache is used for classical inversions and the law does not have inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.
• glacier_id::Int64: Glacier ID in the list of simulated glaciers.

Struct that defines the SciMLSensitivity adjoint flavor. This is the default behavior in ODINN.

SciMLSensitivityAdjoint

TikhonovRegularization(; operator = :laplacian, distance = 3)

A simple regularization type implementing Tikhonov regularization (also known as ridge regularization) for inverse problems.

This struct includes both the forward and reverse (adjoint) operators, which are required for the computation of the gradients with respect to the model parameters.

Keyword Arguments (Constructor)

• operator::Symbol = :laplacian: The regularization operator to use. Currently, only :laplacian is implemented, which penalizes large gradients by applying the Laplacian operator.
• distance::Integer = 3: A width parameter to determine how far from the margin evaluate the loss.

Fields (Struct)

• operator_forward::Function: The forward regularization operator (e.g., ∇²).
• operator_reverse::Function: The reverse-mode (VJP) of the operator (e.g., VJP_λ_∂∇²a_∂a).
• distance::Integer: The distance parameter controlling the extent of regularization.

TrainingStats(;
    retcode::Union{String, Nothing} = nothing,
    losses::Vector{F} = Float64[],
    niter::I = 0,
    θ::Union{ComponentVector, Nothing} = nothing,
    θ_hist::Union{Vector{ComponentVector}, Nothing} = ComponentVector[],
    ∇θ_hist::Union{Vector{ComponentVector}, Nothing} = ComponentVector[]
) where {F <: AbstractFloat, I <: Integer}

Constructor for TrainingStats object used to store important information during training.

Arguments

• retcode: Report code of the optimization.
• losses: Vector storing the value of the loss function at each iteration.
• niter: Total number of iterations/epochs.
• θ: Parameters of neural network after training
• θ_hist: History of parameters of neural network during training
• ∇θ_hist: History of gradients training

mutable struct TrainingStats{F <: AbstractFloat, I <: Integer}

An object with the information of the training.

Fields

• retcode::Union{String, Nothing}: Report code of the optimization.
• losses::Vector{F}: Vector storing the value of the loss function at each iteration.
• niter::I: Total number of iterations/epochs.
• θ::Union{<: ComponentVector, Nothing}: Parameters of neural network after training
• θ_hist::Vector{<: ComponentVector}: History of parameters of neural network during training
• ∇θ_hist::Vector{<: ComponentVector}: History of gradients training
• lastCall::DateTime: Last time the callback diagnosis was called. This is used to compute the time per iteration.

UDEparameters(; sensealg, optim_autoAD, grad, optimization_method, empirical_loss_function, target) where {ADJ <: AbstractAdjointMethod}

Create a UDEparameters object for configuring the sensitivity analysis and optimization of a Universal Differential Equation (UDE).

Keyword Arguments

• sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm: The sensitivity algorithm to use for adjoint calculations. Defaults to GaussAdjoint(autojacvec=SciMLSensitivity.EnzymeVJP()).
• optim_autoAD::AbstractADType: The automatic differentiation type for optimization. Defaults to Optimization.AutoEnzyme().
• grad::ADJ: The adjoint gradient computation method. Defaults to ContinuousAdjoint().
• optimization_method::String: The optimization method to use. Must be either "AD+AD" or "AD+Diff". Defaults to "AD+AD".
• empirical_loss_function::AbstractLoss: The loss function to use for optimization. Defaults to LossH().
• target::Union{Symbol, Nothing}: The target variable for optimization. Defaults to :A.

Returns

• A UDEparameters object configured with the specified sensitivity, optimization, and loss settings.

Description

This function creates a UDEparameters object that encapsulates the configuration for sensitivity analysis, optimization, and loss computation in a Universal Differential Equation (UDE) framework. It verifies that the provided optimization_method is valid and constructs the solver parameters accordingly.

Notes

• The optimization_method must be either "AD+AD" (automatic differentiation for both forward and backward passes) or "AD+Diff" (automatic differentiation combined with finite differences).
• The empirical_loss_function determines how the loss is computed during optimization.

A mutable struct that holds parameters for a UDE (Universal Differential Equation).

UDEparameters{ADJ <: AbstractAdjointMethod} <: AbstractParameters

Fields

• sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm: The sensitivity algorithm used for adjoint sensitivity analysis.
• optimization_method::String: The optimization method to be used.
• target::Symbol: The target variable for the optimization.

struct VJPsPrepLaw <: AbstractPrepVJP

A container struct that holds all objects needed to compute vector-Jacobian products (VJPs) for a law using DifferentiationInterface.

Fields:

• f_θ_first: Function to evaluate the law with parameters θ as the first argument.
• f_inp_first: Function to evaluate the law with inputs as the first argument.
• prep_θ: Precomputed gradient preparation for parameters θ.
• prep_inp: Precomputed gradient preparation for inputs.

This struct is used to prepare the VJP computation with DifferentiationInterface (DI). Depending on the AD backend, DI might require to precompile code and this struct stores the results. This allows each VJP call to be fast in the adjoint PDE by reusing the preparation results.

VelocityRegularization(; reg = TikhonovRegularization(), components = :abs, distance = 3)

Regularization for velocity fields, combining a spatial smoothing operator with optional component control.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator.
• components::Symbol = :abs: Determines which velocity components to regularize (e.g. :abs, :x, :y).
• distance::Integer = 3: Distance to glacier margin.

Base.copyto!(
    dest::InversionBinder{FI, CA},
    src::InversionBinder{FI, CA},
) where {FI <: Inversion, CA <: ComponentArray}

Overload Base.copyto! as we need a way to copy the SciMLStructure. It is used in SciMLSensitivity to differentiate the callbacks.

Base.zero(
    p::InversionBinder{FI, CA},
) where {FI <: Inversion, CA <: ComponentArray}

Overload Base.zero as we need a way to copy the SciMLStructure. It is used in SciMLSensitivity to differentiate the callbacks.

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

Precomputes the vector-Jacobian products (VJPs) for all laws used in the SIA2D ice flow model for a given glacier, time, and model parameters.

Depending on which target (U, Y, or neither) is provided in SIA2D_model, this function checks if the corresponding law supports VJP precomputation and, if so, triggers the appropriate precompute routine for that law. If neither U nor Y is provided, precomputes VJPs for the A, C, and n laws.

Arguments

• SIA2D_model::SIA2Dmodel: The model containing the configuration and laws used for SIA2D ice flow.
• SIA2D_cache::SIA2DCache: A cache object holding intermediate values and storage relevant for precomputations.
• simulation::Inversion: Simulation object containing global simulation parameters.
• glacier_idx::Integer: Index of the glacier being simulated.
• t::Real: Current time in the simulation.
• θ: Model parameters or state variables for the simulation step.

Notes

• This routine is intended as a preparatory step for manual adjoint.
• Only laws supporting VJP precomputation are processed.

run!(simulation::Inversion)

Run the training process for a given Inversion simulation.

Arguments

• simulation::Inversion: The simulation object containing the parameters and settings for the inversion process.

Description

This function initiates the training of a Universal Differential Equation (UDE) for the provided simulation. It prints a message indicating the start of the training process, calls the train_UDE! function to perform the training, and collects the results in results_list. The results are intended to be saved using Sleipnir.save_results_file!, but this step is currently commented out and will be enabled once the optimization is working. Finally, the garbage collector is triggered to free up memory.

Notes

• The Sleipnir.save_results_file! function call is currently commented out and should be enabled once the optimization process is confirmed to be working.
• The garbage collector is explicitly run using GC.gc() to manage memory usage.

CallbackOptimizationSet(θ, l; callbacks)

Helper to combine callbacks for Optimization function. This executes the action of each callback. (equivalent to CallbackSet for DifferentialEquations.jl)

ComponentVector2Vector(cv::ComponentVector)

Transform a ComponentVector into a Vector of same length. This function creates a new Vector and does not mutate the original ComponentVector.

Arguments:

• cv::ComponentVector: Input ComponentVector.

Diffusivity(target::SIA2D_D_target; H, ∇S, θ, iceflow_model, glacier, params)

Compute the effective diffusivity field for a 2D shallow ice model using the diagnostic target and a predicted velocity matrix U.

This function uses a learned or specified model to estimate the velocity matrix U, then calculates the diffusivity as either H .* U (if dimensions match) or the averaged H times U if dimensions differ by one grid cell (staggered grid). Errors if dimensions are incompatible.

Arguments

• target::SIA2D_D_target: Diagnostic target object defining interpolation and scaling rules.

Keyword Arguments

• H: Ice thickness.
• ∇S: Ice surface slope.
• θ: Parameters of the model.
• iceflow_model: Iceflow model used for simulation.
• glacier: Glacier data.
• params: Model parameters.

Returns

• A matrix of diffusivity values with the same shape as H or staggered by one cell, depending on U.

Throws

• An error if the dimensions of U and H are not compatible for diffusivity calculation.

Notes

Supports both grid-matched and staggered configurations by averaging H where necessary.

Gauss Quadratrue for numerical integration

LawA(params::Sleipnir.Parameters; scalar::Bool=true)

Construct a law that defines an ice rheology A per glacier to invert. This can be either a spatially varying A or a scalar value per glacier based on the value of scalar.

Arguments

• params::Sleipnir.Parameters: Parameters struct used to retrieve the minimum and maximum values of A for scaling the parameter to invert.
• scalar::Bool: Whether the ice rheology to invert is a scalar per glacier, or a spatially varying A per glacier (matrix to invert).

function LawA(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    precompute_VJPs::Bool = true,
)

Constructs a law object for the creep coefficient A in the SIA based on a neural network that takes as input the long term air temperature. The creep coefficient A with this law is a scalar. See also SIA2D_A_target.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct used to retrieve the minimum and maximum values of A for scaling of the neural network output.
• precompute_VJPs::Bool: If true, enables precomputation of vector-Jacobian products before solving the adjoint PDE for efficient autodiff.

Returns

• A_law: A Law{ScalarCache} instance that computes the creep coefficient A based on an input temperature using the neural network. The law scales the output to the physical bounds defined by params.

Notes

• The VJP is computed automatically using DifferentiationInterface.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to be between params.physical.minA and params.physical.maxA.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote when using DifferentiationInterface.
• The init_cache function initializes the cache with a scalar zero for the forward placeholder, and with a vector of zeros for the VJP placeholder.

Example

nn_model = NeuralNetwork(params)
A_law = LawA(nn_model, params; precompute_VJPs = false)

LawU(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    max_NN::Union{F, Nothing} = 50.0,
    prescale_bounds::Union{Vector{Tuple{F,F}}, Nothing} = [(0.0, 300.0), (0.0, 0.5)],
) where {F <: AbstractFloat}

Constructs a law object for the diffusive velocity U in the SIA based on a neural network that takes as input the ice thickness H̄ and the surface slope ∇S. The diffusive velocity U with this law is a matrix and the diffusivity in the SIA is obtained through D = U * H̄. See also SIA2D_D_target.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct. Not used for the moment but kept as an argument to keep consistency with other equivalent functions LawA and LawY.
• max_NN::Union{F, Nothing}: Expected maximum value of the neural network output. If set to nothing, no postscaling is applied.
• prescale_bounds::Union{Vector{Tuple{F,F}}, Nothing}: Vector of tuples where each tuple defines the lower and upper bounds of the input for scaling. If set to nothing, no prescaling is applied.
• precompute_interpolation::Bool: Determines which cache to use depending if interpolation is used or not for the evaluation of gradients.
• precompute_VJPs::Bool: Determines is VJPs are stored in the cache during the reverse step.

Returns

• U_law: A Law{Array{Float64, 2}} instance that computes the diffusive velocity U based on the ice thickness H̄ and the surface slope ∇S using the neural network. The law scales the output using the max_NN argument.

Notes

• The computation is compatible with Zygote for automatic differentiation.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to match max_NN.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote.
• The init_cache function initializes the cache with a zero matrix.

Example

nn_model = NeuralNetwork(params)
bounds_H = (0.0, 300.0)
bounds_∇S = (0.0, 0.5)
U_law = LawU(nn_model, params; max_NN = 50.0, prescale_bounds = [bounds_H, bounds_∇S])

LawY(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    max_NN::Union{F, Nothing} = nothing,
    prescale_bounds::Vector{Tuple{F,F}} = [(-25.0, 0.0), (0.0, 500.0)],
) where {F <: AbstractFloat}

Constructs a law object for the hybrid diffusivity Y in the SIA based on a neural network that takes as input the long term air temperature and the ice thickness H̄. The hybrid diffusivity Y with this law is a matrix as it depends on the ice thickness. This law is used in an hybrid setting where the n exponent in the mathematical expression of the diffusivity is different from the one used to generate the ground truth. The goal of this law is to retrieve the missing part of the diffusivity. Please refer to SIA2D_D_hybrid_target for a mathematical definition.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct used to retrieve the maximum value of A for scaling of the neural network output.
• max_NN::Union{F, Nothing}: Expected maximum value of the neural network output. If not specified, the law takes as an expected maximum value params.physical.maxA.
• prescale_bounds::Vector{Tuple{F,F}}: Vector of tuples where each tuple defines the lower and upper bounds of the input for scaling.

Returns

• Y_law: A Law{Array{Float64, 2}} instance that computes the hybrid diffusivity Y based on an input temperature and ice thickness using the neural network. The law scales the output to the physical bounds defined by params.

Notes

• The computation is compatible with Zygote for automatic differentiation.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to match the maximum value which is either max_NN or params.physical.maxA.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote.
• The init_cache function initializes the cache with a zero matrix.

Example

nn_model = NeuralNetwork(params)
bounds_T = (-25.0, 0.0)
bounds_H = (0.0, 500.0)
Y_law = LawY(nn_model, params; prescale_bounds = [bounds_T, bounds_H])

This function allows to extend the Wrapper layers define in Lux to matrices operations.

Model(;
    iceflow::Union{IFM, Nothing} = nothing,
    mass_balance::Union{MBM, Nothing} = nothing,
    regressors::Union{NamedTuple, Nothing} = nothing,
    target::Union{TAR, Nothing} = nothing,
) where {IFM <: IceflowModel, MBM <: MBmodel, TAR <: AbstractTarget}

Creates a new model instance using the provided iceflow, mass balance, and machine learning components.

Arguments

• iceflow::Union{IFM, Nothing}: The iceflow model to be used. Can be a single model or nothing.
• mass_balance::Union{MBM, Nothing}: The mass balance model to be used. Can be a single model or nothing.
• regressors::Union{NamedTuple, Nothing}: The regressors to be used in the laws.

Returns

• model: A new instance of Sleipnir.Model initialized with the provided components.

Constructor for the Parameters type. Since some of the subtypes of parameters are defined in different packages of the ODINN ecosystem, this constructor will call the constructors of the different subtypes and return a Parameters object with the corresponding subtypes. The Parameters mutable struct is defined in Sleipnir.jl using abstract types, which are later on defined in the different packages of the ODINN ecosystem.

Parameters(;
        physical::PhysicalParameters = PhysicalParameters(),
        simulation::SimulationParameters = SimulationParameters(),
        solver::SolverParameters = SolverParameters(),
        hyper::Hyperparameters = Hyperparameters(),
        UDE::UDEparameters = UDEparameters()
        inversion::InversionParameters = InversionParameters()
        )

Keyword arguments

• physical::PhysicalParameters: Physical parameters for the simulation.
• simulation::SimulationParameters: Parameters related to the simulation setup.
• solver::SolverParameters: Parameters for the solver configuration.
• hyper::Hyperparameters: Hyperparameters for the model.
• UDE::UDEparameters: Parameters specific to the UDE (Universal Differential Equation).
• inversion::InversionParameters: Parameters for inversion processes.

Currently just used for Enzyme

Inverse with batch

Compute gradient glacier per glacier

T_A_Alaw(simulation::Simulation, glacier_idx::Integer, θ, t::AbstractFloat)

Evaluate the A law when it defines a mapping between the long term air temperature and the creep coefficient A and return both the input temperature T and the computed creep coefficient A.

Arguments

• simulation::Simulation: The simulation object containing model data and parameters.
• glacier_idx::Integer: Index specifying which glacier to evaluate.
• θ: Model parameters to be used in the law.
• t::AbstractFloat: The time at which to evaluate the law. For this law it is useless but in the general setting, a law needs a time t in order to retrieve the inputs. For the sake of consistency, this input was kept.

Returns

• (T, A): A tuple containing:

  • T: The input long term air temperature for the specified glacier.
  • A: The evaluated creep coefficient for the specified glacier.

Details

• The function checks that the inputs of the A law are exactly as expected (long term air temperature only).
• Retrieves the long term air temperature for the specific glacier.
• Evaluates the creep coefficient using the law.
• Returns both the temperature and creep coefficient as a tuple. Since the cache of A is a zero dimensional array, it is converted to float before returning the value.

Example

T, A = T_A_Alaw(simulation, glacier_idx, θ, 2010.0)

VJP_λ_∂SIA∂H_continuous(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the continuous VJP of the SIA2D equation with respect to H. Given λ and H, it returns the VJP of λ^T * ∂(SIA2D)/∂H (H).

Arguments:

• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• dλ::Matrix{R}: Jacobian vector product, also called input gradient in reverse-mode AD.

VJP_λ_∂SIA∂H_discrete(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <:Real, SIM <: Simulation}

Implementation of the discrete VJP of the SIA2D equation with respect to H. Given λ and H, it returns the VJP of λ^T * ∂(SIA2D)/∂H (H).

Arguments:

• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• dλ::Matrix{R}: Jacobian vector product, also called input gradient in reverse-mode AD.

VJP_λ_∂SIA∂θ_continuous(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the continuous VJP of the SIA2D equation with respect to θ. Given λ, H and θ, it returns the VJP of λ^T * ∂(SIA2D)/∂θ (θ).

Arguments:

• θ: Vector of parameters
• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• ∂θ: Jacobian vector product with respect to θ, also called input gradient in reverse-mode AD. It has the same type as θ.

VJP_λ_∂SIA∂θ_discrete(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the discrete VJP of the SIA2D equation with respect to θ. Given λ, H and θ, it returns the VJP of λ^T * ∂(SIA2D)/∂θ (θ).

Arguments:

• θ: Vector of parameters
• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• ∂θ: Jacobian vector product with respect to θ, also called input gradient in reverse-mode AD. It has the same type as θ.

VJP_λ_∂∇²a_∂a(λ::Matrix{R}, a::Matrix{R}, Δx::R, Δy::R) where {R<:Real}

Computes the vector-Jacobian product (VJP) of the Laplacian operator ∇² with respect to its input field a. This function effectively propagates sensitivities (adjoints) λ backward through the Laplacian, as required in adjoint or reverse-mode differentiation.

Arguments

• λ::Matrix{R}: Adjoint field associated with the Laplacian output.
• a::Matrix{R}: Input scalar field to the Laplacian operator.
• Δx::R: Grid spacing in the x-direction.
• Δy::R: Grid spacing in the y-direction.

Returns

• Matrix{R}: The adjoint (VJP) with respect to a, i.e. ∂⟨λ, ∇²a⟩/∂a.

Vector2ComponentVector(v::Vector, cv_template::ComponentVector)

Transform a vector v to a ComponentVector that has the same structure as cv_template. This function creates a new ComponentVector and copies the values of v explicitly. The arguments v and cv_template must be of the same length.

Arguments:

• v::Vector: Vector whose values are copied.
• cv_template::ComponentVector: ComponentVector whose structure is used to create a new ComponentVector.

Function to evaluate derivatives of ice surface velocity in D inversion.

TODO: This functions right now just make a call to the regular functions used for the calculation of the adjoint. This is not correct, but we keep it as this for now until we figure out how to do this in the case of the D inversion.

_batch_iceflow_UDE(
    container::InversionBinder,
    glacier_idx::Integer,
    iceflow_prob::ODEProblem,
)

Define the callbacks to be called by the ODE solver, solve the ODE and create the results.

_ml_model_postscale(
    Y::Vector,
    max_NN,
)

Applies an exponential transformation to each element in Y, then rescales the result by multiplying with max_NN. For each element, the transformation is: max_NN * exp((Y - 1.0) / Y)

Arguments

• Y::Vector: Values to be post-processed.
• max_NN: Scalar representing the maximum value for rescaling.

Returns

• The rescaled values after applying the exponential transformation.

_ml_model_prescale(
    X::Vector,
    prescale_bounds::Vector{Tuple{F, F}},
) where {F <: AbstractFloat}

Scales each element of the input vector X using the corresponding bounds from prescale_bounds. For each index i, X[i] is normalized based on the interval specified in prescale_bounds[i] using the normalize function. This function is typically used to ensure that the scales of the inputs of a neural network are comparable to each other.

Arguments

• X::Vector: A vector of input values to be normalized.
• prescale_bounds::Vector{Tuple{F, F}}: A vector of tuples specifying the lower and upper bounds for normalization of each corresponding element in X.

Returns

• A vector where each element is the normalized value of the corresponding input, using the specified bounds.

Notes

• The length of X and prescale_bounds must be equal.

_pred_NN(inp::Vector{F}, smodel, θ, prescale, postscale) where {F <: AbstractFloat}

Compute the output of a neural network model on the input vector inp.

Arguments

• inp::Vector{F}: Input vector of floats.
• smodel: The neural network model.
• θ: Parameters for the neural network model.
• prescale: Function to scale the input vector before passing it to the model.
• postscale: Function to scale the model output.

Returns

• The single (scalar) output value from the neural network after applying prescale to the input, evaluating the model, and then applying postscale. The result is extracted via only.

Notes

• The function assumes that the neural network, when evaluated, returns an iterable with exactly one element.
• Using only will throw an error if the output is not exactly one element.

Example

mymodel = StatefulLuxLayer{true}(archi, nothing, st)
y = _pred_NN([1.0, 2.0], mymodel, θ, prescale_fn, postscale_fn)

aggregate∇θ(∇θ::Vector{<: ComponentArray}, θ, submodels::TrainableComponents)

Aggregate the vector of gradients ∇θ as a single ComponentArray. The argument ∇θ is the vector of all the gradients computed for each glacier. This function aggregates them based on the optimizable components of submodels.

backward_loss(
    lossType::MultiLoss,
    H_pred::Matrix{F},
    H_ref,
    V_ref, Vx_ref, Vy_ref,
    t::F,
    glacier_idx::Integer,
    θ,
    simulation,
    normalization::F,
    Δt,
) where {F <: AbstractFloat}

Computes the gradient of a composite loss defined by a MultiLoss object with respect to both the predicted field H_pred and model parameters θ.

Each sub-loss's backward gradient is weighted by its corresponding coefficient in lossType.λs and summed to form the total gradient.

Arguments

• lossType::MultiLoss: Composite loss object containing individual losses and weights.
• H_pred::Matrix{F}: Predicted ice thickness.
• H_ref::Matrix{F}: Reference ice thickness.
• t::F: Current time or simulation step.
• glacier_idx::Integer: Glacier id in the list of glaciers in simulation.
• θ: Model parameters used in the simulation.
• simulation: Simulation object providing necessary context for gradient computation.
• normalization::F: Normalization factor applied within each individual loss.
• Δt: Named tuple containing the time step to use for the approximation of continuous in time loss terms. For example if LossH is used, there must be a term Δt.H containing the time step since the last computation of the ice thickness loss term. If the current time t where the loss is evaluated does not correspond to a time step of the LossH term, then the value of Δt.H has no impact.

Returns

• (∂L∂H, ∂L∂θ): Tuple containing:

  • ∂L∂H::Matrix{F}: Gradient of the composite loss with respect to H_pred.
  • ∂L∂θ: Gradient of the composite loss with respect to model parameters θ.

batch_loss_iceflow_transient(
    container::InversionBinder,
    glacier_idx::Integer,
    iceflow_prob::ODEProblem,
)

Solve the ODE, retrieve the results and compute the loss.

Arguments:

• container::InversionBinder: SciMLStruct that contains the simulation structure and the vector of parameters to optimize.
• glacier_idx::Integer: Index of the glacier.
• iceflow_prob::ODEProblem: Iceflow problem defined as an ODE with respect to time.

callback_diagnosis(θ, l, simulation; save::Bool = false, tbLogger::Union{<: TBLogger, Nothing} = nothing)

Callback function to track and diagose training. It includes print and updates in simulation::Simulation. It also logs training statistics with tensorboard if tbLogger is provided.

callback_plots_A(θ, l, simulation)

Callback function to generate plots during training.

Normalization of D to cap at a maximum physical value

create_interpolation(A::Matrix; n_interp_half::Int) -> Vector{Float64}

Construct a one-dimensional interpolation grid from the elements of a matrix A by flattening it and delegating to create_interpolation(::Vector). This is a convenience method that allows users to pass a 2D array to the function create_interpolation(A::Vector) directly without manually reshaping it.

create_interpolation(
    A::Vector;
    n_interp_half::Int,
    dilation_factor = 1.0,
    minA_unif::Union{F, Nothing} = nothing,
    minA_quantile::Union{F, Nothing} = nothing,
    maxA_unif::Union{F, Nothing} = nothing,
    maxA_quantile::Union{F, Nothing} = nothing
    ) where {F <: AbstractFloat}

Construct a one-dimensional interpolation grid from the data in A, combining uniformly spaced and quantile-based sampling points.

This hybrid interpolation grid provides both coverage of the entire range of values and higher resolution in regions where A has dense data, making it useful for interpolation or machine learning applications that need balanced sampling.

Arguments

• A::Vector: Input data vector (typically containing positive values).
• n_interp_half::Int: Number of points used for both the uniform and quantile-based subsets of the interpolation grid.
• dilation_factor::Real = 1.0: Optional multiplier applied to maximum(A) to slightly extend the grid beyond the data range (useful to avoid extrapolation issues).
• minA_unif::Union{F, Nothing} = nothing: Minimum value used for the uniform interpolation
• minA_quantile::Union{F, Nothing} = nothing: Maximum value used for the uniform interpolation
• maxA_unif::Union{F, Nothing} = nothing: Minimum value used for the quantile interpolation
• maxA_quantile::Union{F, Nothing} = nothing: Maximum value used for the quantile interpolation

Returns

A sorted, unique vector of interpolation nodes combining:

• n_interp_half uniformly spaced values between 0 and dilation_factor * maximum(A)
• n_interp_half quantile-based values computed from the positive entries of A

create_results(θ, simulation::Inversion, mappingFct)

Given the parameters θ, solve the iceflow problem for all the glaciers and aggregate the results for all of them. This function is typically used at the end of a training once θ has been optimized and one wants to run one last forward simulation in order to retrieve statistics about each of the iceflow problems.

Arguments:

• θ: Parameters to use for the forward simulation.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

define_iceflow_prob(
    simulation::Inversion,
    glacier_idx::Integer,
)

Given a simulation struct and a glacier_idx, build the iceflow problem that has to be solved in the ODE solver. In practice, the returned iceflow problem is used inside simulate_iceflow_UDE! through remake. The definition of the iceflow problem has to be done outside of the gradient computation, otherwise Zygote fails at differentiating it.

enable_multiprocessing(params::Sleipnir.Parameters) -> Int

Configures and enables multiprocessing based on the provided simulation parameters.

Arguments

• params::Sleipnir.Parameters: A parameter object containing simulation settings, including the number of workers (params.simulation.workers) and whether multiprocessing is enabled (params.simulation.multiprocessing).

Behavior

• If multiprocessing is enabled (params.simulation.multiprocessing == true) and the specified number of workers (params.simulation.workers) is greater than 0:

  • Adds the required number of worker processes if the current number of processes (nprocs()) is less than the specified number of workers.
  • Suppresses precompilation output on the worker processes and ensures the ODINN module is loaded on all workers.
  • If the specified number of workers is 1, removes all worker processes.

Returns

• The number of worker processes (nworkers()) after configuration.

Notes

• This function uses @eval to dynamically add or remove worker processes.
• Precompilation output is suppressed on workers to reduce noise in the console.

eval_law(law::AbstractLaw, simulation::Simulation, glacier_idx::Integer, inputs::NamedTuple, θ)

Evaluates a law on the specified glacier within a simulation context and for a user defined input.

Arguments

• law::AbstractLaw: The law object to be evaluated. Must provide a function f and an init_cache method.
• simulation::Simulation: The simulation context, containing model parameters and machine learning components.
• glacier_idx::Integer: Index identifying which glacier in the simulation to evaluate the law for.
• input_values::NamedTuple: Input data required by the law and provided by the user.
• θ: Weights used in the law to make inference. This can be nothing when the law has no parameter.

Returns

• The updated cache after evaluating the law. The cache contains the result of the law's computation for the specified glacier and inputs.

Details

• The function initializes a cache for the law using init_cache.
• If the simulation has a machine learning model, the model's parameters (θ) are updated in-place with the provided θ.
• The law's function is then called with the cache, inputs, and parameters. The result is stored in the cache and the cache is returned.
• In future versions, the design may change so that only inputs and θ are needed, with the cache handled separately so that no simulation is required.

Example

result = eval_law(simulation.model.iceflow.A, simulation, glacier_idx, (; T = 273.15), θ)    # Initialize the cache to be able to make an inference of the law

evaluate_H₀(
    θ::ComponentArray,
    glacier::Glacier2D,
    filter::Symbol,
    glacier_id::Integer,
)

Evaluate the initial ice thickness H₀ for a given glacier, optionally applying a smooth thresholding function.

Arguments

• θ::ComponentArray: A ComponentArray containing glacier parameters.

• glacier::Glacier2D: Glacier for which to evaluate H₀.

• filter::Symbol: Specifies the smoothing function to apply to the raw initial condition:

  • :identity: applies the identity function (no change).
  • :softplus: applies the softplus function log(1 + exp(x)) to ensure positivity.
  • :Zang1980: applies the σ_zang function (Zang 1980) as a smooth positivity threshold.

• glacier_id::Integer: Index of the glacier in order to retrieve the parameters of the IC in θ.

Returns

• A numeric value or array representing the filtered initial ice thickness for the specified glacier.

evaluate_∂H₀(
    θ::ComponentArray,
    glacier::Glacier2D,
    filter::Symbol,
    glacier_id::Integer,
)

Evaluate the derivative of the initial ice thickness H₀ for a given glacier, optionally applying a smooth thresholding function.

Arguments

• θ::ComponentArray: A ComponentArray containing glacier parameters.

• glacier::Glacier2D: Glacier for which to evaluate ∂H₀.

• filter::Symbol: Specifies the smoothing function to apply to the raw initial condition:

  • :identity: applies the identity function (no change).
  • :softplus — applies the softplus function log(1 + exp(x)) to ensure positivity.
  • :Zang1980 — applies the σ_zang function (Zang 1980) as a smooth positivity threshold.

• glacier_id::Integer: Index of the glacier in order to retrieve the parameters of the IC in θ.

Returns

• A numeric value or array representing the filtered initial ice thickness for the specified glacier.

feed_input_cache!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation,
    glacier_idx::Integer,
    θ,
    result
)

Populate the input cache of an SIA2DCache instance with interpolation nodes for ice thickness (H) and surface slope (∇S), based on the results of a previous forward simulation. This function is required just when results of the forward pass are required to evaluate the elements of the cache in the reverse step.

This function prepares the interpolation knots used later by the reverse evaluation of the adjoint SIA2D model.

Right now, this function is just required for the inversion w.r.t to D, which is indicated by the boolean variable SIA2D_model.U_is_provided. Other inversions may not required the definition of this function.

Arguments

• SIA2D_model::SIA2Dmodel: The 2D shallow-ice approximation model instance.
• SIA2D_cache::SIA2DCache: The cache object that stores precomputed interpolation nodes.
• simulation: The simulation object containing glacier configurations and model settings.
• glacier_idx::Integer: Index of the glacier within simulation.glaciers for which the cache is being populated.
• θ: Model parameters (not directly used in this function but included for interface consistency).
• result: Output of a previous forward run containing ice thickness fields H.

fourier_feature(v, n::Integer=10, random=false, σ=5.0)

Generates a Fourier feature embedding of a vector v, optionally using randomized projections.

Arguments

• v: Input vector to be transformed (typically a coordinate or feature vector).
• n::Integer=10: Number of Fourier features to generate (default is 10).
• random::Bool=false: Whether to use random Fourier features (default is false).
• σ::Float64=5.0: Standard deviation of the normal distribution used for random feature projection (only used if random=true).

Returns

• A 2n-dimensional vector consisting of sine and cosine features of the transformed input vector.

Notes

Fourier features help to overcome spectral bias in neural networks and can further help to learn higher frequncy components of the function faster. For more information, see Tancik et. al (2020), "Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional Domains".

Example

v = [0.5, 1.0]
features = fourier_feature(v, n = 4, random = true, σ = 2.0)

generate_batches(simulation::S; shuffle=false) where {S <: Simulation}

Generate a data loader for batching simulations, optionally shuffling the batches.

Arguments

• simulation::S: A Simulation object (or subtype of Simulation) containing the data to be batched.
• shuffle::Bool=false: A flag indicating whether to shuffle the batches. Defaults to false.

Returns

• A DataLoader object that provides batched access to the simulation data.

Description

This function creates a DataLoader for batching the provided simulation object. The DataLoader allows for efficient iteration over the simulation data in batches. The batch size is set to 1 by default, and the shuffle flag determines whether the batches are shuffled. If shuffle is enabled, a warning is logged to indicate that the batches used for parallelization are being shuffled.

Notes

• The batch size is fixed at 1 in this implementation. To modify the batch size, you may need to adjust the DataLoader initialization.
• Shuffling the batches may affect reproducibility and parallelization behavior.

generate_simulation_batches(simulation::Inversion)

Generate batches of simulations from a Inversion object for parallel or batched processing.

Arguments

• simulation::Inversion: A Inversion object containing the model, glaciers, parameters, results, and statistics for the simulation.

Returns

• A vector of Inversion objects, each representing a batch of simulations. Each batch contains a subset of glaciers, models, and results from the original simulation.

Description

This function splits the glaciers and associated data in the simulation object into smaller batches for processing. Each batch is represented as a new Inversion object. The number of batches is determined by the nbatches variable (currently set to 1). If the simulation results are empty, the function creates batches with empty results. Otherwise, it includes the corresponding results for each glacier in the batches.

Notes

• The number of glaciers (ninstances) must be divisible by the number of batches (nbatches). An assertion is used to enforce this condition.
• The function currently defaults to nbatches = 1, meaning no actual batching is performed. This can be updated to use simulation.parameters.hyper.batchsize for dynamic batching. # TODO: we need to change how this is done here, with a manual =1 batch size
• If the simulation results are empty, the function creates batches with empty results objects. #simulation.parameters.hyper.batchsize

get_NN(θ_trained, ft; lightNN=false)

Generates a neural network.

Arguments

• θ_trained: Pre-trained neural network parameters (optional).
• ft: Float type used.
• lightNN: Boolean that determines if a light architecture is returned or not.

Returns

• UA: Lux.Chain neural network architecture.
• θ: Neural network parameters.
• st: Lux state.

grad_loss_iceflow!(dθ, θ, simulation::Inversion, mappingFct)

Compute the gradient with respect to θ for all the glaciers and assign the result in-place to dθ.

Arguments:

• dθ: Gradient of the parameters where the computed gradient should be stored.
• θ: Parameters to differentiate.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

grad_loss_iceflow!(θ, simulation::Inversion, mappingFct)

Compute the gradient with respect to θ for all the glaciers and return the result out-of-place. See the in-place implementation for more information.

grad_parallel_loss_iceflow!(θ, simulation::Inversion, glacier_idx::Integer)

Compute the gradient with respect to θ for a particular glacier and return the computed gradient. This function defines the iceflow problem and then calls Zygote to differentiate batch_loss_iceflow_transient with respect to θ. It uses the SciMLSensitivity implementation under the hood to compute the adjoint of the ODE.

function loss(
    lossType::LogSum,
    a::Matrix{F},
    b::Matrix{F},
    mask::BitMatrix,
    normalization::F,
) where {F <: AbstractFloat}

Compute logarithmic loss function for ice velocity fields following Morlighem, M. et al., "Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica". Geophys. Res. Lett. 37, (2010). Given a minimum velocity ϵ the absolute velocity given by a and b, it computes the sum of

log^2( (a + ϵ) / (b + ϵ) )

It has been shown that this loss function enables robust estimation of drag coefficient.

loss(
    lossType::MultiLoss,
    H_pred::Matrix{F},
    H_ref,
    V_ref, Vx_ref, Vy_ref,
    t::F,
    glacier_idx::Integer,
    θ,
    simulation,
    normalization::F,
    Δt,
) where {F <: AbstractFloat}

Computes the weighted composite loss for a prediction H_pred against a reference H_ref using a MultiLoss object.

Each individual loss in lossType.losses is evaluated and multiplied by its corresponding weight in lossType.λs. The final loss is the sum of these weighted contributions.

Arguments

• lossType::MultiLoss: Composite loss object containing individual losses and weights.
• H_pred::Matrix{F}: Predicted ice thickness.
• H_ref::Matrix{F}: Reference ice thickness.
• t::F: Current time or simulation step.
• glacier_idx::Integer: Glacier id in the list of glaciers in simulation.
• θ: Model parameters used in the simulation.
• simulation: Simulation object providing necessary context for loss evaluation.
• normalization::F: Normalization factor applied within each individual loss.
• Δt: Named tuple containing the time step to use for the approximation of continuous in time loss terms. For example if LossH is used, there must be a term Δt.H containing the time step since the last computation of the ice thickness loss term. If the current time t where the loss is evaluated does not correspond to a time step of the LossH term, then the value of Δt.H has no impact.

Returns

• F: The total scalar loss, computed as the sum of weighted individual losses.

loss_iceflow_transient(θ, simulation::Inversion, mappingFct)

Given the parameters θ, this function:

1. Solves the iceflow problem for all the glaciers.
2. Computes the loss function defined as the sum of the loss functions for each of the glaciers. The loss function of each glacier depends on the type of loss. Refer to empirical_loss_function in the UDE parameters for more information. The loss function is transient meaning that the state of the glacier is compared to a reference at different time steps over the simulated period.
3. Return the value of the loss function.

Arguments:

• θ: Parameters to use for the forward simulation.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

merge_batches(results::Vector)

Merge simulation results from multiple batches into a single collection.

Arguments

• results::Vector: A vector where each element is a collection of results (e.g., arrays or vectors) from a batch.

Returns

• A single collection containing all the merged results from the input batches.

Description

This function takes a vector of results from multiple simulation batches and merges them into a single collection using vertical concatenation (vcat). It is useful for combining results that were processed in parallel or in separate batches.

normalize(X; lims::Tuple{F, F}; method = :shift) where {F <: AbstractFloat}

Normalize a variable by using an affine transformation defined by some input lower and upper bounds (m, M) and transforming to O(1) scale.

Arguments

• X: Input value.
• lims::Tuple{F, F}: Lower and upper bounds to use in the affine transformation.
• method::Symbol: Method to scale data.

Returns

• The input variable scaled by the affine transformation.

parallel_loss_iceflow_transient(θ, simulation::Inversion)

Loop over a list of glaciers to process. When multiprocessing is enabled, each call of this function has a dedicated process. This function calls batch_loss_iceflow_transient which returns both the loss and the result structure. The function keeps only the loss.

plot_law(law::AbstractLaw, simulation::Simulation, inputs::NamedTuple, θ; glacier_idx=1, idx_fixed_input=0)

Plot a law function with one or two input variables.

Arguments

• law::AbstractLaw: The law to plot (e.g., a sliding law or creep law).
• simulation::Simulation: The simulation containing glaciers and parameters.
• inputs::NamedTuple: Named tuple of input variables for the law (e.g., (T=iTemp(), H̄=iH̄())).
• θ: Parameters for the law (can be nothing for laws without parameters).

Keyword Arguments

• glacier_idx::Integer=1: Index of the glacier to use for extracting input values (for 2D inputs).
• idx_fixed_input::Integer=0: For two-input laws, index (1 or 2) of the input to fix at its mean value. If 0, plots a 3D surface.
• plot_full_input_range::Bool=false: If true and plotting a 1D law, uses the full physical range for that input.
• ground_truth_law::Union{AbstractLaw, Nothing}=nothing: Optional ground truth law to overlay on the plot for comparison.

Returns

• A plot figure (1D line plot, 2D scatter, or 3D surface) and saves it to the simulation's working directory.

Examples

# 1D plot (temperature input for Cuffey-Paterson law)
plot_law(A_law, simulation, (T = iTemp(),), nothing)

pretraining(architecture::Lux.Chain;
            X::Matrix,
            Y::Matrix,
            nepochs::Int=3000,
            lossfn::GenericLossFunction=MSLELoss(; agg=mean, epsilon=1e-10),
            rng::AbstractRNG=Random.default_rng())

Pretrains a neural network model using a input and output.

Arguments

• architecture::Lux.Chain: The neural network architecture to be trained.
• X::Matrix: Input feature matrix where each column is a feature vector.
• Y::Matrix: Target output matrix corresponding to the inputs in X.
• nepochs::Int=3000: Number of training epochs (default is 3000).
• lossfn::GenericLossFunction=MSLELoss(...): Loss function used for training. Defaults to Mean Squared Logarithmic Error.
• rng::AbstractRNG=Random.default_rng(): Random number generator used for parameter initialization.

Returns

• architecture: The trained neural network architecture.
• θ_pretrain: Trained parameters of the neural network.
• st_pretrain: Internal states of the trained model.

Notes

Pretrainign helps to reduce the number of total epochs required to train the UDE by selecting a physical meaningful initialization for the model. The function initializes the model parameters and states using Lux.setup, then performs training using a custom train_model! function with ADAM optimizer. Loss values are printed every 100 epochs during training.

Example

using Lux, Random
arch = Chain(Dense(10 => 20, relu), Dense(20 => 1))
X = rand(10, 100)
Y = rand(1, 100)
model, params, state = pretraining(arch; X = X, Y = Y)

printProgressLoss(iter, total_iters, loss, improvement)

Print function to track training.

random_matrix(mean, std, corr_length)

Generate a random matrix with entries drawn from a multivariate normal distribution whose mean is given by mean and whose covariance decays exponentially with grid distance.

Arguments

• mean::AbstractMatrix{<:Real}: Matrix specifying the spatial mean values at each grid point. Entries equal to 0.0 are treated as inactive and skipped in sampling.
• std::Real: Standard deviation scaling factor for the covariance kernel.
• corr_length::Real: Correlation length parameter controlling how fast correlations decay with Euclidean distance between grid points.

Returns

• H_sample::Matrix{Float64}: A random realization of the same size as mean, with correlated entries drawn from MvNormal(mean, Σ), where Σ[i,j] = std * exp(-‖coords[i] - coords[j]‖ / corr_length).

safe_slice(obj, ind::Integer)

Return a sliced object obj if ind > 0, otherwise return 0.0.

This function saves the results of an inversion to a file in JLD2 format. If the path argument is not provided, the function will create a default path based on the current project directory. The results are saved in a file named prediction_<nglaciers>glaciers_<tspan>.jld2, where <nglaciers> is the number of glaciers in the simulation and <tspan> is the simulation time span.

simulate_iceflow_UDE!(
    container::InversionBinder,
    cb::SciMLBase.DECallback,
    iceflow_prob::ODEProblem,
    tstops,
)

Make a forward simulation of the iceflow UDE.

splitθ(θ, glacier_idx::Integer, optimizableComponent::TrainableModel)

Given a ComponentVector θ, a glacier_idx and an optimizableComponent, extract the content of θ relevant for the given optimizableComponent and glacier ID glacier_idx.

BFGS optim

train_UDE!(
    simulation::Inversion;
    save_every_iter::Bool = false,
    logger::Union{<: TBLogger, Nothing} = nothing
)

Trains UDE based on the current Inversion.

ADAM optim

update_training_state!(simulation::S, l) where {S <: Simulation}

Update the training state to determine if the training has completed an epoch. If an epoch is completed, reset the minibatches, update the history loss, and increment the epoch count.

Arguments

• simulation: The current state of the simulation or training process.
• l: The current loss value or other relevant metric.

Returns

• None. This function updates the state in-place.

σ_zang(x; β = 2.0)

Smooth thresholding function for enforcing non-negativity and zero values for negative values following I. Zang, "A smoothing-out technique for min—max optimization" (1980).

Arguments

• x::Real: Input value to be thresholded.
• β::Real: (optional) Parameter controlling the transition zone width. Default is 2.0.

∂σ_zang(x; β = 2.0)

Derivative of the smooth thresholding function σ_zang.

Arguments

• x::Real: Input value to be thresholded.
• β::Real: (optional) Parameter controlling the transition zone width. Default is 2.0.

∇²(a::Matrix{F}, Δx::F, Δy::F) where {F<:AbstractFloat}

Computes the 2D Laplacian operator of a scalar field a on a regular grid using finite differences and staggered (dual–primal) averaging.

Arguments

• a::Matrix{F}: 2D scalar field to differentiate.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.

Returns

• Matrix{F}: Approximation of the Laplacian ∇²a with boundary values set to 0.0.

HalfarParameters(; λ=0.0, H₀=3600.0, R₀=750000.0, n=3.0, A=1e-16, f=0.0, ρ=910.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=910.0: Ice density [kg/m³].
• 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(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.

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.

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_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.

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

SolverParameters(;
    solver::OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm = RDPK3Sp35(),
    reltol::F = 1e-12,
    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().
• reltol::F: The relative tolerance for the solver. Defaults to 1e-12.
• 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.

A mutable struct that holds parameters for the solver.

SolverParameters{F <: AbstractFloat, I <: Integer}

Fields

• solver::OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm: The algorithm used for solving differential equations.
• reltol::F: The relative tolerance for the solver.
• 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.

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.

SIA2D!(
    dH::Matrix{R},
    H::Matrix{R},
    simulation::SIM,
    t::R,
    θ,
) 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::R: Current simulation time.
• θ: 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.

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.

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.

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.

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.

generate_ground_truth(
    glaciers::Vector{G},
    params::Sleipnir.Parameters,
    model::Sleipnir.Model,
    tstops::Vector{F};
    store::Tuple=(:H, :V),
) 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 and velocityData 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 and/or :V.

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 the thicknessData field if store contains :H and it populates velocityData if store contains :V.

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.

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),
) where {F <: AbstractFloat}

Return a new vector of glaciers with the simulated thickness and ice velocity 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 and/or :V.

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, then it computes the surface ice velocity data and store it. A new vector of glaciers is created and each glacier is a copy with an updated thicknessData and velocityData 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 and velocityData based on the values provided in store.

∇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.

MBmodel <: AbstractModel

An abstract type representing a mass balance model in the Muninn package. This serves as a base type for all specific mass balance models, ensuring they adhere to a common interface and can be used interchangeably within the ODINN framework.

TImodel <: MBmodel

An abstract type representing a temperature index mass balance models within the ODINN framework. This type serves as a parent type for more specialized mass balance models, ensuring they adhere to a common interface defined by the MBmodel abstract type.

TImodel1(params::Sleipnir.Parameters; DDF::F = 7.0/1000.0, acc_factor::F = 1.0/1000.0) where {F <: AbstractFloat}

Create a temperature index model with one degree-day factor (DDF) with the given parameters.

Arguments

• params::Sleipnir.Parameters: The simulation parameters.
• DDF::F: Degree-day factor (default is 7.0/1000.0).
• acc_factor::F: Accumulation factor (default is 1.0/1000.0).

Returns

• TI1_model: An instance of TImodel1 with the specified parameters.

TImodel1{F <: AbstractFloat}

A structure representing a temperature index model with degree-day factor and accumulation factor.

Keyword arguments

• DDF::F: Degree-day factor, which is a coefficient used to convert temperature into melt.
• acc_factor::F: Accumulation factor, which is a coefficient used to adjust the accumulation of mass.

Type Parameters

• F: A subtype of AbstractFloat representing the type of the factors.

TImodel2(params::Sleipnir.Parameters; DDF_snow::F = 3.0/1000.0, DDF_ice::F = 6.0/1000.0, acc_factor::F = 1.0/1000.0) where {F <: AbstractFloat}

Create a temperature-index model with two degree-day factors (TImodel2) for mass balance calculations.

Arguments

• params::Sleipnir.Parameters: The parameters object containing simulation settings.
• DDF_snow::F: Degree-day factor for snow (default: 3.0/1000.0).
• DDF_ice::F: Degree-day factor for ice (default: 6.0/1000.0).
• acc_factor::F: Accumulation factor (default: 1.0/1000.0).

Returns

• TI2_model: An instance of the TImodel2 with the specified parameters.

TImodel2{F <: AbstractFloat}

A type representing a temperature-index model with parameters for snow and ice degree-day factors, and an accumulation factor.

Keyword arguments

• DDF_snow::F: Degree-day factor for snow, which determines the melt rate of snow per degree above the melting point.
• DDF_ice::F: Degree-day factor for ice, which determines the melt rate of ice per degree above the melting point.
• acc_factor::F: Accumulation factor, which scales the accumulation of snow.

Type Parameters

• F: A subtype of AbstractFloat, representing the numeric type used for the model parameters.

MB_timestep!(cache, model::Model, glacier::G, step::F, t) where {F <: AbstractFloat, G <: AbstractGlacier}

Simulates a mass balance timestep for a given glacier model.

Arguments

• cache: The model cache to update.
• model::Model: The glacier model containing iceflow and mass balance information.
• glacier::G: The glacier object containing climate and DEM data.
• step::F: The timestep duration.
• t: The current time.

Description

This function performs the following steps:

1. Computes the period from the previous timestep to the current time.
2. Retrieves cumulative climate data for the specified period.
3. Downscales the climate dataset to a 2D grid based on the glacier's DEM.
4. Computes the mass balance for the glacier and updates the model's iceflow mass balance.

MB_timestep(model::Model, glacier::G, step::F, t::F) where {F <: AbstractFloat, G <: AbstractGlacier}

Calculate the mass balance (MB) for a glacier over a given timestep.

Keyword arguments

• model::Model: The model containing mass balance parameters.
• glacier::G: The glacier object containing climate and DEM data.
• step::F: The timestep duration.
• t::F: The current time.

Returns

• MB::Matrix{F}: The computed mass balance matrix for the given timestep.

Details

1. Computes the period between the current time t and the previous step t - step.
2. Retrieves cumulative climate data for the specified period.
3. Downscales the climate data to a 2D grid based on the glacier's DEM.
4. Computes the mass balance using the downscaled climate data.

compute_MB(
    mb_model::TImodel1,
    climate_2D_period::Climate2Dstep,
    step::AbstractFloat,
)

Compute the mass balance (MB) for a given mass balance model and climate period.

Arguments

• mb_model::TImodel1: The mass balance model containing parameters such as accumulation factor (acc_factor) and degree-day factor (DDF).
• climate_2D_period::Climate2Dstep: The climate data for a specific period, including snow accumulation (snow) and positive degree days (PDD).
• step::AbstractFloat: The step used to update MB. This scales the MB so that the accumulation and degree-day factors are scaled monthly.

Returns

• A numerical array representing the computed mass balance, calculated as the difference between the product of the accumulation factor and snow, and the product of the degree-day factor and positive degree days.

AbstractData

Abstract type that represents data. Used to implement ThicknessData and SurfaceVelocityData.

AbstractGlacier

An abstract type representing a glacier. This serves as a base type for different glacier implementations in the Sleipnir package.

AbstractInput

Abstract type representing an input source for a Law. Concrete subtypes must implement:

• default_name(::ConcreteInput): returns the default field name (as a Symbol) under which this input will be accessed in the law's input NamedTuple.
• get_input(::ConcreteInput, simulation, glacier_idx, t): retrieves the actual input value given the simulation context, glacier index, and time t.

AbstractLaw

Abstract type representing a synthetic law. Currently it's only used for testing by making easier to create dumb laws, but in the future it may be cleaner to use different concrete type of laws (for example CallbackLaw, ContinuousLaw, or LearnableLaw)

Concrete subtypes must implement:

• apply_law!(::ConcreteLaw, state, simulation, glacier_idx, t, θ)
• init_cache(::ConcreteLaw, glacier, glacier_idx)
• law_VJP_input(::ConcreteLaw, cache, simulation, glacier_idx, t, θ)
• law_VJP_θ(::ConcreteLaw, cache, simulation, glacier_idx, t, θ)
• precompute_law_VJP(::ConcreteLaw, cache, vjpsPrepLaw, simulation, glacier_idx, t, θ)
• cache_type(::ConcreteLaw)
• is_callback_law(::ConcreteLaw)
• is_precomputable_law_VJP(::ConcreteLaw)
• callback_freq(::ConcreteLaw)
• inputs(::ConcreteLaw)
• inputs_defined(::ConcreteLaw)
• apply_law_in_model(::ConcreteLaw)

AbstractModel

An abstract type that serves as a base for all model types in ODINN.

AbstractParameters

An abstract type that serves as a base for all parameter-related types in the ODINN ecosystem.

AbstractPrepVJP

Abstract type representing the preparation of Vector-Jacobian Product (VJP) computations for the laws. Subtypes of AbstractPrepVJP are used to handle any precomputations or setup required before evaluating VJPs, such as configuring automatic differentiation backends or precompiling code. This type provides a flexible interface for implementing custom VJP preparation strategies for different laws. It is for internal use only and not exposed to the user.

Cache

Abstract type for defining a cache struct of a law.

Mandatory field:

• value: Store the result of a forward evaluation.

Optional fields:

• vjp_inp: Store the result of the evaluation of the vector-Jacobian product (VJP) with respect to the inputs.
• vjp_θ: Store the result of the evaluation of the vector-Jacobian product (VJP) with respect to the parameters θ.

Notes:

• If a concrete subtype does not implement vjp_inp and vjp_θ, the law should not be used for gradient computation, and therefore for inversions.

A mutable struct representing a 2D climate for a glacier with various buffers and datasets.

Climate2D{CLIMRAW <: RasterStack, CLIMRAWSTEP <: RasterStack, CLIMSTEP <: ClimateStep, CLIM2DSTEP <: Climate2Dstep, F <: AbstractFloat}

Fields

• raw_climate::CLIMRAW: Raw climate dataset for the whole simulation.

• climate_raw_step::CLIMRAWSTEP: Raw climate trimmed for the current step to avoid memory allocations.

• climate_step::ClimateStep: Climate data for the current step.

• climate_2D_step::Climate2Dstep: 2D climate data for the current step to feed to the mass balance (MB) model.

• longterm_temps::Vector{F}: Long-term temperatures for the ice rheology.

• avg_temps::F: Intermediate buffer for computing average temperatures.

• avg_gradients::F: Intermediate buffer for computing average gradients.

• ref_hgt::F: Reference elevation of the raw climate data.

  Climate2D( rgi_id, params::Parameters, S::Matrix{<: AbstractFloat}, Coords::Dict, )

Initialize the climate data given a RGI ID, a matrix of surface elevation and glacier coordinates.

Arguments

• rgi_id: The glacier RGI ID.
• params::Parameters: The parameters containing simulation settings and paths.
• S::Matrix{<: AbstractFloat}: Matrix of surface elevation used to initialize the downscaled climate data.
• Coords::Dict: Coordinates of the glacier.

Description

This function initializes the climate data for a glacier by:

1. Creating a dummy period based on the simulation time span and step.

2. Loading the raw climate data from a NetCDF file.

3. Calculating the cumulative climate data for the dummy period.

4. Downscaling the cumulative climate data to a 2D grid.

5. Retrieving long-term temperature data for the glacier.

6. Returning the climate data, including raw climate data, cumulative climate data, downscaled 2D climate data, long-term temperatures, average temperatures, and average gradients.

   Climate2D( rawclimate::RasterStack, climaterawstep::RasterStack, climatestep::ClimateStep, climate2Dstep::Climate2Dstep, longtermtemps::Vector{<: AbstractFloat}, avgtemps::AbstractFloat, avggradients::AbstractFloat, refhgt::AbstractFloat, )

Initialize the climate data with the fields provided as arguments. Refer to the list of fields for a complete description of the arguments.

Climate2Dstep{F <: AbstractFloat}

A mutable struct representing a 2D climate time step with various climate-related parameters.

Keyword arguments

• temp::Matrix{F}: Temperature matrix.
• PDD::Matrix{F}: Positive Degree Days matrix.
• snow::Matrix{F}: Snowfall matrix.
• rain::Matrix{F}: Rainfall matrix.
• gradient::F: Gradient value.
• avg_gradient::F: Average gradient value.
• x::Vector{F}: X-coordinates vector.
• y::Vector{F}: Y-coordinates vector.
• ref_hgt::F: Reference height.

ClimateStep{F <: AbstractFloat}

Mutable struct that represents a climate step before downscaling.

Keyword arguments

• prcp::F: Cumulative precipitation for the given period.
• temp::F: Cumulative temperature at the reference elevation for the given period.
• gradient::F: Cumulative temperature gradient for the given period.
• avg_temp::F: Average temperature over the time step.
• avg_gradient::F: Average temperature gradient over the time step.
• ref_hgt::F: Reference elevation of the raw climate data.

ConstantLaw{CACHE_TYPE}(init_cache)

Creates a constant law of type ConstantLaw{CACHE_TYPE} that holds a fixed value for the entire simulation.

This is useful to inject glacier-specific or global constants into the simulation without modifying them over time. The update function is a no-op, and only the init_cache function matters.

Arguments

• init_cache::Function: A function init_cache(simulation, glacier_idx, θ)::T that provides the constant value.

Type Parameters

• CACHE_TYPE: The type of the cache. Must be specified manually and should match the return type of init_cache.

Examples

# Same value for all glaciers
n_law = ConstantLaw{Float64}(Returns(4.))

# Value depending on the glacier
n_law = ConstantLaw{Float64}((sim, i, θ) -> sim.glaciers[i].n)

# Learned value
n_law = ConstantLaw{Float64}((sim, i, θ) -> θ.n)

Container

Abstract type that defines a container to be used in the PDE solver. It is useful to retrieve the simulation object when applying callback laws.

CustomVJP <: VJPType

Indicates that a law uses a custom-defined function for VJP computation. This is used when the VJP is provided manually rather than computed automatically.

DIVJP <: VJPType

Indicates that a law uses the default VJP computation provided by DifferentiationInterface.jl. This is used when no custom VJP function is provided, and the Vector-Jacobian Product (VJP) is computed automatically using DifferentiationInterface.jl.

GenInputsAndApply{IN, F}

Given a tuple of AbstractInputs and a function f, returns a callable struct. This struct with_input_f can be evaluated as a function that generates the inputs and then applies the function's law with_input_f.f. It is for internal use only and it isn't exposed to the user.

function Glacier1D(; rgiid::String = "", climate::Union{Climate1D, Nothing} = nothing, H₀::Vector{F} = Vector{Sleipnir.Float}([]), S::Vector{F} = Vector{Sleipnir.Float}([]), B::Vector{F} = Vector{Sleipnir.Float}([]), V::Vector{F}= Vector{Sleipnir.Float}([]), A::Union{F, Nothing} = nothing, C::Union{F, Nothing} = nothing, n::Union{F, Nothing} = nothing, w₀::Union{Vector{F}, Nothing} = nothing, λ::Union{Vector{F}, Nothing} = nothing, slope::Vector{F} = Vector{Sleipnir.Float}([]), distborder::Vector{F} = Vector{Sleipnir.Float}([]), Coords::Dict{String, Vector{Float64}} = Dict{String, Vector{Float64}}("lon" => [], "lat" => []), Δx::F = 0, Δy::F = 0, nx::I = 0, ny::I = 0, ) where {F <: AbstractFloat, I <: Integer}

Constructor for empty 2D Glacier object.

Glacier2D(
    glacier::Glacier2D;
    thicknessData::Union{<: ThicknessData, Nothing} = nothing,
    velocityData::Union{<: SurfaceVelocityData, Nothing} = nothing,
)

Copies a Glacier2D object and updates the thickness and/or surface velocity data.

Arguments

• glacier::Glacier2D: The original glacier struct.
• thicknessData::Union{<: ThicknessData, Nothing}: Thickness data structure that is used to store the reference values. Default is nothing which keeps the existing thickness data.
• velocityData::Union{<: SurfaceVelocityData, Nothing}: Surface velocity data structure that is used to store the reference values. Default is nothing which keeps the existing surface velocity data.

Returns

• A Glacier2D object that is a copy of the original one with the thickness and/or surface velocity data updated.

Glacier2D(
    rgi_id::String,
    params::Parameters;
    masking::Union{Int, Nothing, BitMatrix} = 2,
    smoothing=false
)

Build glacier object for a given RGI ID and parameters.

Arguments

• rgi_id::String: The RGI ID of the glacier.

• params::Parameters: A Parameters object containing simulation parameters.

• masking::Union{Int, Nothing, BitMatrix}: Type of mask applied to the glacier to determine regions with no ice.

  • When masking is an Int, the mask is based on the initial ice thickness H₀ and it is set to true for pixels outside at a distance of the glacier borders greater than the value of masking.
  • When masking is set to nothing, the mask is set to a BitMatrix full of falses.
  • When masking is a BitMatrix, this matrix is used for the mask. Defaults to 2.

• smoothing::Bool=false: Optional; whether to apply smoothing to the initial ice thickness. Default is false.

• test::Bool=false: Optional; test flag. Default is false.

Returns

• glacier::Glacier2D: A Glacier2D object initialized with the glacier data.

Description

This function loads and initializes the glacier data for a given RGI ID. It retrieves the initial ice thickness conditions based on the specified source in the parameters, applies optional smoothing, and initializes the glacier's topographical and velocity data. The function also handles Mercator projection for the glacier coordinates and filters glacier borders in high elevations to avoid overflow problems.

Notes

• The function reverses the matrices for ice thickness, bedrock, and other data to match the required orientation.
• If the Mercator projection includes latitudes larger than 80°, a warning is issued.
• If the glacier data is missing, the function updates a list of missing glaciers and issues a warning.

Glacier2D(;
    rgi_id::String = "",
    name::String = "",
    climate::Climate2D = nothing,
    H₀::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    H_glathida::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    S::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    B::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    V::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    Vx::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    Vy::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    A::F = 0.0,
    C::F = 0.0,
    n::F = 0.0,
    p::F = 0.0,
    q::F = 0.0,
    slope::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    dist_border::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    mask::BitMatrix = BitMatrix([;;]),
    Coords::Dict{String, Vector{Float64}} = Dict{String, Vector{Float64}}("lon" => [], "lat" => []),
    Δx::F = 0.0,
    Δy::F = 0.0,
    nx::I = 0,
    ny::I = 0,
    cenlon::F = NaN,
    cenlat::F = NaN,
    params_projection::Dict{String, Float64} = Dict{String, Float64}(),
    thicknessData::THICKDATA = nothing,
    velocityData::SURFVELDATA = nothing,
) where {
    F <: AbstractFloat,
    I <: Integer,
    THICKDATA <: Union{<: ThicknessData, Nothing},
    SURFVELDATA <: Union{<: SurfaceVelocityData, Nothing},
}

Constructs a Glacier2D object with the given parameters, including default ones.

Arguments

• rgi_id::String: The RGI identifier for the glacier.
• name::String: The name of the glacier if available.
• climate::Climate2D: The climate data associated with the glacier.
• H₀::Matrix{F}: Initial ice thickness matrix.
• H_glathida::Matrix{F}: Ice thickness matrix from GLATHIDA.
• S::Matrix{F}: Surface elevation matrix.
• B::Matrix{F}: Bed elevation matrix.
• V::Matrix{F}: Ice velocity magnitude matrix.
• Vx::Matrix{F}: Ice velocity in the x-direction matrix.
• Vy::Matrix{F}: Ice velocity in the y-direction matrix.
• A::F: Flow law parameter.
• C::F: Sliding law parameter.
• n::F: Flow law exponent.
• p::F: Power law exponent associated to Weertman sliding law (Power associated to basal drag).
• q::F: Power law exponent associated to Weertman sliding law (Power associated to normal pressure).
• slope::Matrix{F}: Slope matrix.
• dist_border::Matrix{F}: Distance to border matrix.
• mask::BitMatrix: Boolean matrix representing the glacier mask, where true values indicate regions constrained by the mask (i.e., no-ice zones)
• Coords::Dict{String, Vector{Float64}}: Coordinates dictionary with keys "lon" and "lat".
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• cenlon::F: Central longitude of the glacier.
• cenlat::F: Central latitude of the glacier.
• params_projection::Dict{String, Float64}: Projection parameters that allows mapping the regional grid to global WGS84 coordinates.
• thicknessData::THICKDATA: Thickness data structure that is used to store the reference values.
• velocityData::SURFVELDATA: Surface velocity data structure that is used to store the reference values.

Returns

• A Glacier2D object with the specified parameters.

A mutable struct representing a 2D glacier. Notice that all fields can be empty by providing nothing as the default value.

/!\ WARNING /!\ Glacier objects should not be constructed manually, but rather through the initialize_glaciers function.

Glacier2D{F <: AbstractFloat, I <: Integer, CLIM <: Climate2D, THICKDATA <: Union{<: ThicknessData, Nothing}, SURFVELDATA <: Union{<: SurfaceVelocityData, Nothing}}

Fields

• rgi_id::String: The RGI (Randolph Glacier Inventory) identifier for the glacier.
• name::String: The name of the glacier if available.
• climate::CLIM: The climate data associated with the glacier.
• H₀::Matrix{F}: Initial ice thickness matrix.
• H_glathida::Matrix{F}: Ice thickness matrix from the GLATHIDA dataset.
• S::Matrix{F}: Surface elevation matrix.
• B::Matrix{F}: Bedrock elevation matrix.
• V::Matrix{F}: Ice velocity magnitude matrix.
• Vx::Matrix{F}: Ice velocity in the x-direction matrix.
• Vy::Matrix{F}: Ice velocity in the y-direction matrix.
• A::F: Flow law parameter.
• C::F: Sliding law parameter.
• n::F: Flow law exponent.
• p::F: Power law exponent associated to Weertman sliding law (Power associated to basal drag).
• q::F: Power law exponent associated to Weertman sliding law (Power associated to normal pressure).
• slope::Matrix{F}: Surface slope matrix.
• dist_border::Matrix{F}: Distance to the glacier border matrix.
• mask::BitMatrix: Boolean matrix representing the glacier mask, where true values indicate regions constrained by the mask (i.e., no-ice zones)
• Coords::Dict{String, Vector{Float64}}: Coordinates dictionary with keys as coordinate names and values as vectors of coordinates.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• cenlon::F: Longitude of the glacier center.
• cenlat::F: Latitude of the glacier center.
• params_projection::Dict{String, Float64}: Projection parameters that allows mapping the regional grid to global WGS84 coordinates.
• thicknessData::THICKDATA: Thickness data structure that is used to store the reference values.
• velocityData::SURFVELDATA: Surface velocity data structure that is used to store the reference values.

IntegratedTrajectoryMapping <: VelocityMapping

Integrated trajectory mapping. This mapping is closer to reality as it consists in integrating over time the instantaneous ice surface velocities along ice flow trajectories in a Lagrangian way. This integrated velocity is then compared to the velocity of the datacube. It has not been implemented yet but its computational cost will likely be expensive.

Fields

• spatialInterp::Symbol: The spatial interpolation to use to map the ice surface velocity grid to the glacier grid. For the moment only :nearest is supported.

Law{T}(;
    inputs = nothing,
    f!, f_VJP_input! = nothing, f_VJP_θ! = nothing,
    init_cache, callback_freq = nothing,
    p_VJP! = nothing,
    max_value = NaN, min_value = NaN, name = :unknown,
) where{T}

Defines a physical or empirical law applied to a glacier model that mutates an internal state T at each simulation time step.

# ❌ Will not work: Float64 is immutable, so cache .= ... has no effect
Law{Float64}(;
    f! = (cache, _, _, t, θ) -> cache = θ.scale * sin(2π * t + θ.shift),
    init_cache = (_, _, _) -> 0.0,
)

# ✅ Correct: using a 0-dimensional array allows in-place mutation
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
)

Arguments

• f!::Function: A function with signature f!(cache::T, simulation, glacier_idx, t, θ) that updates the internal state. If inputs are provided, the function instead takes the form f!(cache::T, inputs, θ).
• init_cache::Function: A function init_cache(simulation, glacier_idx, θ)::T that initializes the internal state for a given glacier.
• callback_freq::Union{Nothing, Real, Period, Month}: Optional. If provided, the law is treated as a callback law and is only applied every callback_freq time units. If callback_freq is set to zero, then the law is applied only once at the beginning of the simulation. If callback_freq is set to nothing (default), then the law is applied at every iteration. If callback_freq is provided as a Period or Month, it is converted to a float value in a yearly basis.
• f_VJP_input!: A function with signature (cache::T, simulation, glacier_idx, t, θ) that updates cache.vjp_inp which is the VJP with respect to the inputs.
• f_VJP_θ!: A function with signature (cache::T, simulation, glacier_idx, t, θ) that updates cache.vjp_θ which is the VJP with respect to the parameters θ.
• p_VJP!: A function with signature (cache::T, vjpsPrepLaw, simulation, glacier_idx, t, θ) that performs the precomputation of the VJPs.
• inputs::Union{Nothing, Tuple{<:AbstractInput}}: Optional. Provides automatically generated inputs passed to f! at runtime.
• max_value::Float64: Optional. The maximum value that the law can take, used for plotting and capping the function output.
• min_value::Float64: Optional. The minimum value that the law can take, used for plotting and capping the function output.
• name::Symbol: A name for the law, used for identification and plotting.

Type Parameters

• T: The type of the internal state. Must be specified manually and should match the return type of init_cache.

Notes

• Refer to the tutorials in the documentation for a complete description of the VJP options.

Examples

# A law applied at every timestep, storing a scalar value
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
)

# A callback law applied once per month (assuming time in years)
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
    callback_freq = 1 / 12,
)

MatrixCache <: Cache

A cache structure for storing a two-dimensional array of Float64 values along with their associated vector-Jacobian products (VJP). This is typically used for spatially varying laws. Fields:

• value::Array{Float64, 2}: The cached matrix.
• vjp_inp::Array{Float64, 2}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.

MatrixCacheNoVJP <: Cache

A mutable cache structure for storing a two-dimensional array of Float64 values. This is typically used for spatially varying laws. This struct is intended for use cases where the law is not differentiated, and hence the vector-Jacobian products (VJP) are not required. Fields:

• value::Array{Float64, 2}: The cached matrix.

MeanDateVelocityMapping <: VelocityMapping

Mean date velocity mapping. It is the most simple mapping one can build and it consists in taking the 2D vector field of ice velocity associated to a given mean date and compare it to the instantaneous ice surface velocity obtained from the ice flow model. It is valid only for ice surface velocities estimated from short time windows since the velocity can vary within this time window.

Fields

• spatialInterp::Symbol: The spatial interpolation to use to map the ice surface velocity grid to the glacier grid. For the moment only :nearest is supported.

Model{IFM <: AbstractEmptyModel, MBM <: AbstractEmptyModel, TC <: AbstractEmptyModel}

A mutable struct that represents a model with three components: iceflow, mass balance, and machine learning.

Model(
    iceflow::IFM,
    mass_balance::MBM,
    trainable_components::TC,
) where {IFM <: AbstractEmptyModel, MBM <: AbstractEmptyModel, TC <: AbstractEmptyModel}

Model(;iceflow, mass_balance) = Model(iceflow, mass_balance, nothing)

Initialize Model (no machine learning model).

Keyword arguments

• iceflow::IFM}: Represents the iceflow component, which is an instance of IFM.
• mass_balance::Union{MBM, Vector{MBM}}: Represents the mass balance component, which is an instance of MBM.
• trainable_components::TC: Represents the trainable components, which is an instance of TC.

Type Parameters

• IFM: A subtype of AbstractEmptyModel representing the type of the iceflow model.
• MBM: A subtype of AbstractEmptyModel representing the type of the mass balance model.
• TC: A subtype of AbstractEmptyModel representing the type of the trainable components.

ModelCache{IFC, MBC}

Cache struct that holds the internal state or memory buffers for the components of a Model.

Typically used to store per-glacier preallocated buffers or intermediate results that persist across time steps during simulation.

Fields

• iceflow::IFC: Cache associated with the iceflow model.
• mass_balance::MBC: Cache associated with the mass balance model.

Type Parameters

• IFC: Cache type for the iceflow model.
• MBC: Cache type for the mass balance model.

NullLaw <: AbstractLaw

This struct represents a law that is not used in the iceflow model.

Parameters(; physical::PhysicalParameters = PhysicalParameters(), simulation::SimulationParameters = SimulationParameters())

Constructs a Parameters object with the given physical and simulation parameters.

Arguments

• physical::PhysicalParameters: An instance of PhysicalParameters (default: PhysicalParameters()).
• simulation::SimulationParameters: An instance of SimulationParameters (default: SimulationParameters()).

Returns

• A Parameters object initialized with the provided physical and simulation parameters.

Notes

• If simulation.multiprocessing is enabled, multiprocessing is configured with the specified number of workers.

mutable struct Parameters{PPHY <: AbstractEmptyParams, PSIM <: AbstractEmptyParams, PHY <: AbstractEmptyParams,
    PSOL <: AbstractEmptyParams, PUDE <: AbstractEmptyParams, PINV <: AbstractEmptyParams}

A mutable struct that holds various parameter sets for different aspects of a simulation or model.

Fields

• physical::PPHY: Physical parameters.
• simulation::PSIM: Simulation parameters.
• hyper::PHY: Hyperparameters.
• solver::PSOL: Solver parameters.
• UDE::PUDE: Universal Differential Equation (UDE) parameters.
• inversion::PINV: Inversion parameters.

Type Parameters

• PPHY: Type of the physical parameters, must be a subtype of AbstractEmptyParams.
• PSIM: Type of the simulation parameters, must be a subtype of AbstractEmptyParams.
• PHY: Type of the hyperparameters, must be a subtype of AbstractEmptyParams.
• PSOL: Type of the solver parameters, must be a subtype of AbstractEmptyParams.
• PUDE: Type of the UDE parameters, must be a subtype of AbstractEmptyParams.
• PINV: Type of the inversion parameters, must be a subtype of AbstractEmptyParams.

Initialize the physical parameters of a model.

PhysicalParameters(;
    ρ::Float64 = 900.0,
    g::Float64 = 9.81,
    ϵ::Float64 = 1e-10,
    η₀::F = 1.0,
    maxA::Float64 = 8e-17,
    minA::Float64 = 8.5e-20,
    maxC::Float64 = 8e-17, # TODO: to be revised
    minC::Float64 = 8.5e-20,
    maxTlaw::Float64 = 1.0,
    minTlaw::Float64 = -25.0,
    noise_A_magnitude::Float64 = 5e-18
    )

Keyword arguments

- `ρ`: Ice density
- `g`: Gravitational acceleration.
- `ϵ`: Regularization used in the square root of norms for AD numerical stability.
- `η₀`: Factor to cap surface elevation differences with the upstream ice thickness to impose boundary condition in the iceflow equation
- `maxA`: Maximum value for `A` (Glen's coefficient)
- `minA`: Minimum value for `A` (Glen's coefficient)
- `maxC`: Maximum value of sliding coefficient `C`
- `minC`: Minimum value of sliding coefficient `C`
- `maxTlaw`: Maximum value of Temperature used in simulations on fake law
- `minTlaw`: Minimum value of Temperature used in simulations on fake law
- `noise_A_magnitude`: Magnitude of noise added to A

A structure representing physical parameters used in simulations.

PhysicalParameters{F <: AbstractFloat}

Fields

• ρ::F: Density of ice.
• g::F: Gravitational acceleration.
• ϵ::F: Regularization used in the square root of norms for AD numerical stability.
• η₀::F: Initial viscosity.
• maxA::F: Maximum A.
• minA::F: Minimum A.
• maxC::F: Maximum C.
• minC::F: Minimum C.
• maxTlaw::F: Maximum temperature according to some law.
• minTlaw::F: Minimum temperature according to some law.
• noise_A_magnitude::F: Magnitude of noise in A.

Results(glacier::G, ifm::IF;
    rgi_id::String = glacier.rgi_id,
    H::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    H_glathida::Matrix{F} = glacier.H_glathida,
    H_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    S::Matrix{F} = zeros(Sleipnir.Float, size(ifm.S)),
    B::Matrix{F} = zeros(Sleipnir.Float, size(glacier.B)),
    V::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vx::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vy::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    V_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vx_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vy_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    date_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    date1_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    date2_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    Δx::F = glacier.Δx,
    Δy::F = glacier.Δy,
    lon::F = glacier.cenlon,
    lat::F = glacier.cenlat,
    nx::I = glacier.nx,
    ny::I = glacier.ny,
    t::Vector{F} = Vector{Sleipnir.Float}([]),
    tspan::Tuple{F, F} = (NaN, NaN),
) where {G <: AbstractGlacier, F <: AbstractFloat, IF <: AbstractModel, I <: Integer}

Construct a Results object for a glacier simulation.

Arguments

• glacier::G: The glacier object, subtype of AbstractGlacier.
• ifm::IF: The model object, subtype of AbstractModel.
• rgi_id::String: The RGI identifier for the glacier. Defaults to glacier.rgi_id.
• H::Vector{Matrix{F}}: Ice thickness matrices. Defaults to an empty vector.
• H_glathida::Matrix{F}: Ice thickness from GlaThiDa. Defaults to glacier.H_glathida.
• H_ref::Vector{Matrix{F}}: Reference ice thickness. Defaults to an empty vector.
• S::Matrix{F}: Surface elevation matrix. Defaults to a zero matrix of the same size as ifm.S.
• B::Matrix{F}: Bed elevation matrix. Defaults to a zero matrix of the same size as glacier.B.
• V::Vector{Matrix{F}}: Velocity magnitude matrix. Defaults to an empty vector.
• Vx::Vector{Matrix{F}}: Velocity in the x-direction matrix. Defaults to an empty vector.
• Vy::Vector{Matrix{F}}: Velocity in the y-direction matrix. Defaults to an empty vector.
• V_ref::Vector{Matrix{F}}: Reference velocity magnitude matrix. Defaults to an empty vector.
• Vx_ref::Vector{Matrix{F}}: Reference velocity in the x-direction matrix. Defaults to an empty vector.
• Vy_ref::Vector{Matrix{F}}: Reference velocity in the y-direction matrix. Defaults to an empty vector.
• date_Vref::Vector{F}: Date of velocity observation (mean of date1 and date2). Defaults to an empty vector.
• date1_Vref::Vector{F}: First date of velocity acquisition. Defaults to an empty vector.
• date2_Vref::Vector{F}: Second date of velocity acquisition. Defaults to an empty vector.
• Δx::F: Grid spacing in the x-direction. Defaults to glacier.Δx.
• Δy::F: Grid spacing in the y-direction. Defaults to glacier.Δy.
• lon::F: Longitude of the glacier grid center. Defaults to glacier.cenlon.
• lat::F: Latitude of the glacier grid center. Defaults to glacier.cenlat.
• nx::I: Number of grid points in the x-direction. Defaults to glacier.nx.
• ny::I: Number of grid points in the y-direction. Defaults to glacier.ny.
• tspan::Tuple(F, F): Timespan of the simulation.
• θ::Union{Nothing, ComponentArray{F}}: Model parameters. Defaults to nothing.
• loss::Union{Nothing, Vector{F}}: Loss values. Defaults to nothing.

Returns

• results::Results: A Results object containing the simulation results.

mutable struct Results{F <: AbstractFloat, I <: Integer}

A mutable struct to store the results of simulations.

Fields

• rgi_id::String: Identifier for the RGI (Randolph Glacier Inventory).
• H::Vector{Matrix{F}}: Vector of matrices representing glacier ice thickness H over time.
• H_glathida::Matrix{F}: Optional matrix for Glathida ice thicknesses.
• H_ref::Vector{Matrix{F}}: Reference data for ice thickness.
• S::Matrix{F}: Glacier surface altimetry.
• B::Matrix{F}: Glacier bedrock.
• V::Matrix{F}: Glacier ice surface velocities.
• Vx::Matrix{F}: x-component of the glacier ice surface velocity V.
• Vy::Matrix{F}: y-component of the glacier ice surface velocity V.
• V_ref::Matrix{F}: Reference data for glacier ice surface velocities V.
• Vx_ref::Matrix{F}: Reference data for the x-component of the glacier ice surface velocity Vx.
• Vy_ref::Matrix{F}: Reference data for the y-component of the glacier ice surface velocity Vy.
• date_Vref::Vector{F}: Date of velocity observation (mean of date1 and date2).
• date1_Vref::Vector{F}: First date of velocity acquisition.
• date2_Vref::Vector{F}: Second date of velocity acquisition.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• lon::F: Longitude of the glacier grid center.
• lat::F: Latitude of the glacier grid center.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• tspan::Vector{F}: Time span of the simulation.

ScalarCache <: Cache

A cache structure for storing a scalar value as a zero-dimensional array of Float64 along with their associated vector-Jacobian products (VJP). This is typically used for constant per glacier laws. Fields:

• value::Array{Float64, 0}: The cached scalar value.
• vjp_inp::Array{Float64, 0}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.

ScalarCacheNoVJP <: Cache

A mutable cache structure for storing a scalar value as a zero-dimensional array of Float64. This is typically used for constant per glacier laws. This struct is intended for use cases where the law is not differentiated, and hence the vector-Jacobian products (VJP) are not required. Fields:

• value::Array{Float64, 0}: The cached scalar value.

Simulation

An abstract type representing a generic simulation. This type is intended to be subclassed by specific simulation types to provide a common interface and shared functionality for all simulations.

Constructor for SimulationParameters type, including default values.

SimulationParameters(;
    use_MB::Bool = true,
    use_iceflow::Bool = true,
    plots::Bool = true,
    use_velocities::Bool = true,
    f_surface_velocity_factor::F = 1.0,
    overwrite_climate::Bool = false,
    use_glathida_data::Bool = false,
    tspan::Tuple{F, F} = (2010.0, 2015.0),
    step_MB::F = 1/12,
    multiprocessing::Bool = true,
    workers::I = 4,
    working_dir::String = "",
    test_mode::Bool = false,
    rgi_paths::Dict{String, String} = Dict{String, String}(),
    ice_thickness_source::String = "Farinotti19",
    mapping::VM = MeanDateVelocityMapping(),
    gridScalingFactor::I = 1,
) where {I <: Integer, F <: AbstractFloat, VM <: VelocityMapping}

Keyword arguments

• use_MB::Bool: Whether to use mass balance (default: true).
• use_iceflow::Bool: Whether to use ice flow (default: true).
• plots::Bool: Whether to generate plots (default: true).
• use_velocities::Bool: Whether to calculate velocities (default: true).
• f_surface_velocity_factor::F: Numerical factor representing the ratio between depth integrated ice velocity and surface velocity (default: 1.0).
• overwrite_climate::Bool: Whether to overwrite climate data (default: false).
• use_glathida_data::Bool: Whether to use GLATHIDA data (default: false).
• float_type::DataType: Data type for floating point numbers (default: Float64).
• int_type::DataType: Data type for integers (default: Int64).
• tspan::Tuple{F, F}: Time span for the simulation (default: (2010.0, 2015.0)).
• step_MB::F: Time step for the MB simulation (default: 1/12).
• multiprocessing::Bool: Whether to use multiprocessing (default: true).
• workers::I: Number of workers for multiprocessing (default: 4).
• working_dir::String: Working directory for the simulation (default: "").
• test_mode::Bool: Whether to run in test mode (default: false).
• rgi_paths::Dict{String, String}: Dictionary of RGI paths (default: Dict{String, String}()).
• ice_thickness_source::String: Source of ice thickness data, either "Millan22" or "Farinotti19" (default: "Farinotti19").
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• gridScalingFactor::I: Grid downscaling factor, used to speed-up the tests. Default value is 1 which means no downscaling is applied.

Returns

• simulation_parameters: A new SimulationParameters object.

Throws

• AssertionError: If ice_thickness_source is not "Millan22" or "Farinotti19".

Notes

• If the global variable ODINN_OVERWRITE_MULTI is set to true, multiprocessing is is enabled in any case and the number of workers specified in the simulation parameters must correspond to the number of processes with which Julia has been started. This is to allow the documentation to build successfully in ODINN as we cannot change the number of process in the CI.

A structure to hold simulation parameters for a simulation in ODINN.

struct SimulationParameters{I <: Integer, F <: AbstractFloat, VM <: VelocityMapping} <: AbstractParameters

Fields

• use_MB::Bool: Flag to indicate whether mass balance should be used.
• use_iceflow::Bool: Flag to indicate whether ice flow should be used.
• plots::Bool: Flag to indicate whether plots should be generated.
• use_velocities::Bool: Flag to indicate whether velocities should be calculated.
• f_surface_velocity_factor::F: Numerical factor representing the ratio between depth integrated ice velocity and surface velocity.
• overwrite_climate::Bool: Flag to indicate whether to overwrite climate data.
• use_glathida_data::Bool: Flag to indicate whether to use GLATHIDA data.
• tspan::Tuple{F, F}: Time span for the simulation.
• step_MB::F: Time step for the MB simulation.
• multiprocessing::Bool: Flag to indicate whether multiprocessing should be used.
• workers::I: Number of workers for multiprocessing.
• working_dir::String: Directory for working files.
• test_mode::Bool: Flag to indicate whether to run in test mode.
• rgi_paths::Dict{String, String}: Dictionary of RGI paths.
• ice_thickness_source::String: Source of ice thickness data.
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• gridScalingFactor::I: Grid downscaling factor, used to speed-up the tests. Default value is 1 which means no downscaling is applied.

Constructs SurfaceVelocityData using data from Rabatel et. al (2023) with the given parameters, including default ones.

function SurfaceVelocityData(; x::Union{Vector{F}, Nothing} = nothing, y::Union{Vector{F}, Nothing} = nothing, lat::Union{Vector{F}, Nothing} = nothing, lon::Union{Vector{F}, Nothing} = nothing, vx::Union{Vector{Matrix{F}}, Nothing} = nothing, vy::Union{Vector{Matrix{F}}, Nothing} = nothing, vabs::Union{Vector{Matrix{F}}, Nothing} = nothing, vxerror::Union{Vector{F}, Nothing} = nothing, vyerror::Union{Vector{F}, Nothing} = nothing, vabserror::Union{Vector{F}, Nothing} = nothing, date::Union{Vector{DateTime}, Nothing} = nothing, date1::Union{Vector{DateTime}, Nothing} = nothing, date2::Union{Vector{DateTime}, Nothing} = nothing, dateerror::Union{Vector{Day}, Vector{Millisecond}, Nothing} = nothing, flag::Union{BitMatrix, Nothing} = nothing, isGridGlacierAligned::Bool = false, ) where {F <: AbstractFloat}

Constructor for ice surface velocity data based on Rabatel et. al (2023).

Important remarks:

• Velocities values are reported in m/yr. Positive velocities correspond to the direction of increasing index of the glacier. When the glacier is oriented in east-west and south-north (see latitude and coordinate ordering), positive velocities of the ice surface velocity correspond to positive east-west and south-north velocity component.
• The error in velocity is unique per timestamp, rather than being pixel distributed.
• The error in the absolute velocities vabs_error is overestimated.

References:

• Rabatel, A., Ducasse, E., Millan, R. & Mouginot, J. Satellite-Derived Annual Glacier Surface Flow Velocity Products for the European Alps, 2015–2021. Data 8, 66 (2023).

A mutable struct representing a surface velocity data. Notice that all fields can be empty by providing nothing as the default value.

SurfaceVelocityData{F <: AbstractFloat} <: AbstractData

Fields

• x::Union{Vector{F}, Nothing}: Easting of observation.
• y::Union{Vector{F}, Nothing}: Northing of observation.
• lat::Union{Vector{F}, Nothing}: Latitude of observation.
• lon::Union{Vector{F}, Nothing}: Longitude of observation.
• vx::Union{Vector{Matrix{F}}, Nothing}: x / longitudinal component of surface velocity. Positive velocities correspond to the direction of increasing index of the glacier x-coordinate.
• vy::Union{Vector{Matrix{F}}, Nothing}: y / latitudinal component of surface velocity. Positive velocities correspond to the direction of increasing index of the glacier y-coordinate.
• vabs::Union{Vector{Matrix{F}}, Nothing}: Absolute ice surface velocity.
• vx_error::Union{Vector{F}, Nothing}: Error in vx
• vy_error::Union{Vector{F}, Nothing}: Error in vy
• vabs_error::Union{Vector{F}, Nothing}: Error in vabs.
• date::Union{Vector{DateTime}, Nothing}: Date of observation (mean of date1 and date2)
• date1::Union{Vector{DateTime}, Nothing}: First date of acquisition.
• date2::Union{Vector{DateTime}, Nothing}: Second date of acquisition.
• date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}: Error in date.
• flag::Union{BitMatrix, Nothing}: Flag indicating whether a pixel is considered as reliable or not.
• isGridGlacierAligned::Bool: Whether the data have been gridded to the glacier grid or not.

Simple time series of ice thickness data to test transient inversion

VJPType

Abstract type representing the mode of Vector-Jacobian Product (VJP) computation used in a law. Subtypes of VJPType define how the VJP is evaluated (e.g., custom or using DifferentiationInterface.jl).

VelocityMapping

Abstract type representing the mapping to use in order to map the ice velocity products onto the glacier grid. It contains all needed information to build both the spatial projection, and how to interpolate the data in time.

DummyClimate2D(;
    longterm_temps::Vector{F} = []
) where {F <: AbstractFloat}

Dummy climate initialization for very specific use cases where we don't have climate data and we need to build a minimalistic climate with only a few data. For the moment it supports only the initialization of the long term temperatures. It returns a minimalistic Climate2D instance.

Arguments:

• longterm_temps::Vector{F}: Long term temperatures.

ReverseUTMercator(x::F, y::F; k=0.9996, cenlon=0.0, cenlat=0.0, x0=0.0, y0=0.0, zone::Union{Nothing, Int}=nothing, hemisphere=nothing) where {F <: AbstractFloat}

Transverse Mercator Projection. This function reprojects latitude/longitude into northing/easting coordinates.

Keyword arguments

- `k`: scale factor of the projection
- `cenlon`: Central longitude used in the projection
- `cenlat`: Central latitude used in the projection
- `x0`: Shift in easting
- `y0`: Shift in northing
- `zone` : Zone of the projection
- `hemisphere`: Either :north or :south

UTMercator(x::F, y::F; k=0.9996, cenlon=0.0, cenlat=0.0, x0=0.0, y0=0.0, zone::Union{Nothing, Int}=nothing, hemisphere=nothing) where {F <: AbstractFloat}

Transverse Mercator Projection. This function reprojects northing/easting coordinates into latitude/longitude.

Keyword arguments

- `k`: scale factor of the projection
- `cenlon`: Central longitude used in the projection
- `cenlat`: Central latitude used in the projection
- `x0`: Shift in easting
- `y0`: Shift in northing
- `zone` : Zone of the projection
- `hemisphere`: Either :north or :south

apply_all_callback_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)

This function is a placeholder and must be implemented for your custom model type.

It is intended to apply all callback laws in the simulation to update the model cache for a given glacier at time t with parameters θ. By default, calling this function will throw an error to indicate that the user should provide their own implementation tailored to their model.

Arguments

• model::AbstractModel: The model instance.
• cache: The cache object storing state variables.
• simulation: The simulation context.
• glacier_idx: Index identifying the glacier.
• t: The current simulation time.
• θ: The parameter vector.

Throws

• Always throws an error: "This function should not be called. Implement apply_all_callback_laws! for your own model."

apply_all_non_callback_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)

This function is a placeholder and must be implemented for your custom model type.

It is intended to apply all non-callback laws in the simulation to update the model cache for a given glacier at time t with parameters θ. By default, calling this function will throw an error to indicate that the user should provide their own implementation tailored to their model.

Arguments

• model::AbstractModel: The model instance.
• cache: The cache object storing state variables.
• simulation: The simulation context.
• glacier_idx: Index identifying the glacier.
• t: The current simulation time.
• θ: The parameter vector.

Throws

• Always throws an error: "This function should not be called. Implement apply_all_non_callback_laws! for your own model."

apply_t_cumul_grad!(climate_2D_step::Climate2Dstep, S::Matrix{F}) where {F <: AbstractFloat}

Apply temperature and precipitation gradients based on the positive degree day (PDD) and on the elevation matrix S to the climate data in climate_2D_step.

Arguments

• climate_2D_step::Climate2Dstep: The climate data structure containing temperature, PDD, gradients, and reference height.
• S::Matrix{F}: A matrix of elevations.

Description

This function updates the temperature and PDD fields in climate_2D_step by applying the respective gradients based on the difference between the elevation matrix S and the reference height. Negative PDD values are cropped to zero. Additionally, the function adjusts the rain and snow fractions based on the updated temperature values.

apply_t_grad!(climate::RasterStack, dem::Raster)

Apply temperature gradients to the climate data based on the digital elevation model (DEM).

Arguments

• climate::RasterStack: A RasterStack object containing climate data, including temperature and gradient information.
• dem::Raster: A Raster object representing the digital elevation model (DEM) data.

Description

This function adjusts the temperature data in the climate object by applying the temperature gradients. The adjustment is based on the difference between the mean elevation from the DEM data and a reference height specified in the metadata of the climate object.

block_average(mat::Matrix{F}, n::Int) where {F <: AbstractFloat}

Downsamples a matrix by averaging non-overlapping n x n blocks. Returns a matrix of the block-averaged values with size (div(X, n), div(Y, n)) where (X, Y) = size(mat).

Arguments

• mat::Matrix{F}: Input 2D matrix.
• n::Int: Block size for downsampling. Both matrix dimensions must be divisible by n.

block_average_pad_edge(mat::Matrix{F}, n::Int) where {F <: AbstractFloat}

Downsamples a matrix by averaging n x n blocks, using edge-replication padding when the matrix dimensions are not divisible by n. Edge padding replicates the last row/column values to expand the matrix so that both dimensions are divisible by n. Returns a matrix of averaged values with size (ceil(Int, X/n), ceil(Int, Y/n)).

Arguments

• mat::Matrix{F}: Input 2D matrix.
• n::Int: Block size for downsampling.

build_affect(law::AbstractLaw, cache, glacier_idx, θ)

Return a !-style function suitable for use in a callback, which applies the given law to update the cache for a specific glacier and parameters θ, using the simulation time.

combine_velocity_data(refVelocities; merge=false)

Combine multiple ice surface velocity datasets into a single SurfaceVelocityData object.

Arguments

• refVelocities::Vector{SurfaceVelocityData}: A vector of ice surface velocity datasets to combine. Each element must have the same grid alignment (isGridGlacierAligned must be true for all).
• merge::Bool=false: If true, velocities with the same date are averaged, and corresponding date ranges (date1, date2) are reduced to their min/max. If false, data are simply concatenated.

Returns

• SurfaceVelocityData: A single object containing the combined velocity data, including vx, vy, vabs and their associated errors, as well as coordinate (x, y, lat, lon) and date information. The isGridGlacierAligned field reflects whether all input datasets were aligned.

Notes

• The function asserts that all input datasets are aligned on the same grid.
• When merge=true, velocities and errors are averaged over datasets sharing the same date.
• Uses nanmean for averaging to handle missing data.
• date_error is set to nothing when merging. # Check all surfaces are on the same grid as the glacier

create_results(
    simulation::SIM,
    glacier_idx::I,
    solution,
    tstops::Vector{F};
    processVelocity::Union{Nothing, Function} = nothing,
) where {SIM <: Simulation, I <: Integer}

Create a Results object from a given simulation and solution.

Arguments

• simulation::SIM: The simulation object of type Simulation.
• glacier_idx::I: The index of the glacier within the simulation.
• solution: The solution object containing all the steps including intermediate ones.
• tstops::Vector{F}: The list of time steps to use to construct the results.
• processVelocity::Union{Nothing, Function}=nothing: Post processing function to map the ice thickness to the surface velocity. It is called before creating the results. It takes as inputs simulation, ice thickness (matrix) and the associated time and returns 3 variables Vx, Vy, V which are all matrix. Defaults is nothing which means no post processing is applied.

Returns

• results: A Results object containing the processed simulation data.

Details

The function processes the solution to select the last value for each time step. It then constructs a Results object containing various attributes from the simulation and the iceflow model.

downscale_2D_climate!(glacier::Glacier2D)

Update the 2D climate structure for a given glacier by downscaling climate data.

Arguments

• glacier::Glacier2D: The glacier object containing the climate data to be downscaled.

Description

This function updates the 2D climate structure of the given glacier by:

1. Updating the temperature, PDD (Positive Degree Days), snow, and rain fields in the 2D climate step with the corresponding values from the climate step.
2. Updating the gradients and average gradients in the 2D climate step.
3. Applying temperature gradients and computing the snow/rain fraction for the selected period by reprojecting the current S with the RasterStack structure.

Notes

# Update 2D climate structure

• The function modifies the glacier object in place.

downscale_2D_climate(climate_step::ClimateStep, S::Matrix{<: AbstractFloat}, Coords::Dict)

Downscales climate data to a 2D grid based on the provided matrix of surface elevation and coordinates.

Arguments

• climate_step::ClimateStep: A struct containing climate data for a specific time step. Expected fields are:

  • "avg_temp": Average temperature.
  • "temp": Temperature.
  • "prcp": Precipitation.
  • "gradient": Temperature gradient.
  • "avg_gradient": Average temperature gradient.
  • "ref_hgt": Reference height.

• S::Matrix{<: AbstractFloat}: Surface elevation data.

• Coords::Dict: A dictionary with keys "lon" and "lat" for longitude and latitude coordinates.

Returns

• Climate2Dstep: A Climate2Dstep object containing the downscaled climate data with fields:

  • temp: 2D array of temperature.
  • PDD: 2D array of positive degree days.
  • snow: 2D array of snow precipitation.
  • rain: 2D array of rain precipitation.
  • gradient: Temperature gradient.
  • avg_gradient: Average temperature gradient.
  • x: Longitude coordinates.
  • y: Latitude coordinates.
  • ref_hgt: Reference height.

Description # Create dummy 2D arrays to have a base to apply gradients afterwards

This function creates dummy 2D arrays based on the provided surface elevation data and applies the climate step data to these arrays. It then constructs a Climate2Dstep object with the downscaled climate data and applies temperature gradients to compute the snow/rain fraction for the selected period.

emptyPrepVJP(cache, vjpPrep, simulation, glacier_idx, t, θ)
emptyPrepVJPWithInputs(cache, vjpPrep, inputs, θ)

Function that defines an empty !-style function for the preparation of the VJPs. The two methods define the two possible signatures for the function that updates the cache in-place. Trying to apply this function will yield an error. It is for internal use only and it isn't exposed to the user.

emptyVJP(cache, simulation, glacier_idx, t, θ)
emptyVJPWithInputs(cache, inputs, θ)

Function that defines an empty !-style function for the VJPs. The two methods define the two possible signatures for the function that updates the cache in-place. It is for internal use only and it isn't exposed to the user.

fake_interpolated_datacube()

Create a fake datacube of ice surface velocity time series. It corresponds to the interpolated data.

fake_multi_datacube()

Create a fake datacube of ice surface velocity time series. It corresponds to the filtered multi source data.

fillNaN!(A::AbstractArray, fill::Number=zero(eltype(A)))

Replace all NaN values in the array A with the specified fill value.

Arguments

• A::AbstractArray: The array in which NaN values will be replaced.
• fill::Number: The value to replace NaN with. Defaults to zero(eltype(A)).

fillNaN(A::AbstractArray, fill::Number=zero(eltype(A)))

Replace all NaN values in the array A with the specified fill value. If no fill value is provided, it defaults to the zero value of the element type of A.

Arguments

• A::AbstractArray: The input array that may contain NaN values.
• fill::Number: The value to replace NaNs with. Defaults to zero(eltype(A)).

Returns

• An array of the same type and shape as A, with all NaN values replaced by fill.

fillZeros!(A::AbstractArray, fill::Number=NaN)

Replace all zero elements in the array A with the specified fill value.

Arguments

• A::AbstractArray: The array in which to replace zero elements.
• fill::Number: The value to replace zero elements with. Defaults to NaN.

fillZeros(A::AbstractArray, fill::Number=NaN) -> AbstractArray

Replace all zero elements in the array A with the specified fill value.

Arguments

• A::AbstractArray: The input array in which zero elements are to be replaced.
• fill::Number: The value to replace zero elements with. Defaults to NaN.

Returns

• AbstractArray: A new array with zero elements replaced by the fill value.

filter_missing_glaciers!(rgi_ids::Vector{String}, params::Parameters)

Filter out glaciers that cannot be processed from the given list of RGI IDs.

Arguments

• rgi_ids::Vector{String}: A vector of RGI IDs representing glaciers.
• params::Parameters: A Parameters object containing simulation parameters.

Description

This function filters out glaciers from the provided rgi_ids list based on two criteria:

1. Glaciers that are marked as level 2 in the RGI statistics CSV file.
2. Glaciers listed in the missing_glaciers.jld2 file located in the params.simulation.working_dir directory.

Notes

TODO: see if this is necessary, otherwise remove

• The RGI statistics CSV file is downloaded from a remote server.
• If the missing_glaciers.jld2 file is not available, a warning is logged and the function skips this filtering step. # Check which glaciers we can actually process # TODO: see if this is necessary, otherwise remove

generate_raw_climate_files(rgi_id::String, simparams::SimulationParameters)

Generate raw climate files for a given RGI (Randolph Glacier Inventory) ID and simulation parameters.

Arguments

• rgi_id::String: The RGI ID for which to generate raw climate files.
• simparams::SimulationParameters: The simulation parameters containing the time span and RGI paths.

Description

This function generates raw climate files for a specified RGI ID if they do not already exist. It retrieves raw climate data, ensures the desired period is covered, crops the data to the desired time period, and saves the raw climate data to disk.

Details

1. Constructs the path to the RGI directory using the provided rgi_id and simparams.

2. Checks if the raw climate file for the specified time span already exists.

3. If the file does not exist:

   • Retrieves the raw climate data.
   • Ensures the desired period is covered by the climate data.
   • Crops the climate data to the desired time period.
   • Saves the cropped climate data to disk. # Initialize RGI path to be accessible outside the try block
   • Triggers garbage collection to free up memory.

get_cumulative_climate!(climate, t::AbstractFloat, step::AbstractFloat, gradient_bounds=[-0.009, -0.003])
get_cumulative_climate!(climate, period::StepRange{Date, Day}, gradient_bounds=[-0.009, -0.003])

Calculate and update the cumulative climate data for a given period. The user can choose between providing a specific time t and a time step step, or a time period defined by period.

Keyword arguments

• climate::Climate: The climate object containing raw climate data.
• gradient_bounds::Vector{Float64}: Optional. The bounds within which to clamp the gradient values. Default is [-0.009, -0.003]. Optional parameters to specify the time period:
• t::AbstractFloat: Time at which the cumulative climate data should be computed.
• step::AbstractFloat: Time step used to compute the cumulative climate data. Together with t they define a time period. or
• period::StepRange{Date, Day}: The time period for which to compute the cumulative climate data.

Updates

• climate.climate_raw_step: The raw climate data for the given period.
• climate.avg_temps: The average temperature for the given period.
• climate.avg_gradients: The average gradient for the given period.
• climate.climate_step.prcp: The cumulative precipitation for the given period.
• climate.climate_step.temp: The cumulative temperature for the given period.
• climate.climate_step.gradient: The cumulative gradient for the given period.
• climate.climate_step.avg_temp: The average temperature for the given period.
• climate.climate_step.avg_gradient: The average gradient for the given period.
• climate.climate_step.ref_hgt: The reference height from the raw climate data.

get_cumulative_climate(
    climate::RasterStack,
    gradient_bounds::Vector{Float64}=[-0.009, -0.003],
)

Calculate cumulative climate statistics from the given climate data.

Keyword arguments

• climate::RasterStack: A RasterStack object containing temperature, precipitation, and gradient data.
• gradient_bounds::Vector{Float64}: A two-element vector specifying the lower and upper bounds for the gradient values. Defaults to [-0.009, -0.003].

Returns

• climate_sum::ClimateStep: A struct containing the following fields:

  • "temp": The sum of positive degree days (PDDs) from the temperature data.
  • "prcp": The sum of precipitation data.
  • "gradient": The sum of gradient data, clipped within the specified bounds.
  • "avg_temp": The average temperature.
  • "avg_gradient": The average gradient.
  • "ref_hgt": The reference height from the climate metadata.

Notes

• The temperature data is modified to only include positive degree-day values (PDDs).
• The gradient data is clipped within the specified bounds to ensure plausible values.

get_glathida!(glaciers::Vector{G}, params::Parameters; force=false) where {G <: Glacier2D}

Retrieve and process glacier thickness data for a vector of Glacier2D objects.

Arguments

• glaciers::Vector{Glacier2D}: A vector of Glacier2D objects for which the glacier thickness data is to be retrieved.
• params::Parameters: A Parameters object containing simulation parameters.
• force::Bool=false: A boolean flag indicating whether to force the retrieval of glacier thickness data.

Returns

• gtd_grids::Vector: A vector of glacier thickness data grids.
• glaciers::Vector{Glacier2D}: The updated vector of Glacier2D objects after removing glaciers with no data.

Description

This function retrieves glacier thickness data for each glacier in the input vector using parallel processing. It updates a list of missing glaciers if any glacier has all data points equal to zero. The function then removes glaciers with no data from both the gtd_grids and glaciers vectors and returns the updated vectors.

Notes

• The function uses pmap for parallel processing of glaciers.
• The list of missing glaciers is stored in a JLD2 file located at params.simulation.working_dir/data/missing_glaciers.jld2.
• Glaciers with no data are identified and removed based on the condition that all data points in their thickness grid are zero.

get_glathida_glacier(glacier::Glacier2D, params::Parameters, force)

Retrieve or generate the glathida glacier grid for a given glacier.

Arguments

• glacier::Glacier2D: The glacier object for which the glathida grid is to be retrieved or generated.
• params::Parameters: The parameters object containing simulation settings.
• force: A boolean flag indicating whether to force regeneration of the glathida grid even if it already exists.

Returns

• gtd_grid: A 2D array representing the glathida glacier grid.

Description

This function checks if the glathida glacier grid file (glathida.h5) exists in the specified path. If the file exists and force is false, it reads the grid from the file. Otherwise, it reads the glacier thickness data from a CSV file (glathida_data.csv), computes the average thickness for each grid cell, and saves the resulting grid to an HDF5 file (glathida.h5).

get_longterm_temps(rgi_id::String, params::Parameters, climate::RasterStack) -> Array{Float64}

Calculate the long-term average temperatures for a given glacier.

Arguments

• rgi_id::String: The RGI (Randolph Glacier Inventory) identifier for the glacier.
• params::Parameters: A struct containing simulation parameters, including paths to RGI data.
• climate::RasterStack: A RasterStack object containing climate data.

Returns

• Array{Float64}: An array of long-term average temperatures.

Description

This function retrieves the gridded data for the specified glacier using its RGI identifier. It then applies a temperature gradient to the climate data based on the glacier's topography. Finally, it calculates the long-term average temperatures by grouping the temperature data by year and computing the mean for each group.

get_raw_climate_data(rgi_path::String) -> RasterStack

Load raw climate data from a specified path.

Arguments

• rgi_path::String: The file path to the directory containing the climate data file.

Returns

• RasterStack: A RasterStack object containing the climate data from the specified file.

get_result_id_from_rgi(glacier_id::I, simulation::SIM) where {I <: Integer, SIM <: Simulation}

Extract results of specific simulation from the Simulation object.

Arguments

• glacier_id::I: Numerical ID of glacier used to generate simulation.
• simulation::SIM`: The simulation object containing the parameters and results.

glacierName(rgi_id::String)
glacierName(rgi_ids::Vector{String})

Returns the name(s) of one or multiple glaciers based the given RGI ID(s). It uses the rgi62_stats.csv file from OGGM.

grid(
    glacier::G,
    latitudes::Vector{F},
    longitudes::Vector{F},
    vx::Union{FileArray, Array{Union{Missing, F}, 3}},
    vy::Union{FileArray, Array{Union{Missing, F}, 3}},
    mapping::VM
) where {
    G <: AbstractGlacier,
    F <: AbstractFloat,
    VM <: VelocityMapping,
    FileArray <: Rasters.FileArray
}

Grid velocity data onto the glacier grid following the prescribed mapping. This function maps the 3 dimensional surface velocities (x, y and t) to the glacier grid. The provided surface velocities can be a Rasters.FileArray which happens when the RasterStack is instantiated in lazy mode. In this situation, only the smallest cube that contains all the needed data to construct the mapping is read from disk. The returned velocity variables have shape (nTimes, nx, ny) where nTimes is the number of time steps and (nx, ny) is the size of the glacier grid.

Arguments:

• glacier::G: Glacier instance which determines the glacier on which the velocities are projected onto.
• latitudes::Vector{F}: Vector of latitude values of the original surface velocity grid.
• longitudes::Vector{F}: Vector of longitude values of the original surface velocity grid.
• vx::Union{FileArray, Array{Union{Missing, F}, 3}}: X component of the original surface velocities. It can be either a Rasters.FileArray if the datacube is read in lazy mode, or a plain 3 dimensional array.
• vy::Union{FileArray, Array{Union{Missing, F}, 3}}: Y component of the original surface velocities. It can be either a Rasters.FileArray if the datacube is read in lazy mode, or a plain 3 dimensional array.
• mapping::VM: Mapping to use.

Returns:

• xG: A vector that gives the x coordinates of the glacier grid.
• yG: A vector that gives the y coordinates of the glacier grid.
• vxG: A 3 dimensional array of the x component of the velocity gridded onto the glacier grid.
• vyG: A 3 dimensional array of the y component of the velocity gridded onto the glacier grid.

initialize_glacier(rgi_id::String, parameters::Parameters; smoothing=false)

Initialize a glacier with the given RGI ID and parameters.

Arguments

• rgi_id::String: The RGI (Randolph Glacier Inventory) ID of the glacier.
• parameters::Parameters: A struct containing various parameters required for initialization.
• smoothing::Bool: Optional. If true, apply smoothing to the initial topography. Default is false.
• masking::Union{Int, Nothing, Matrix}: Type of mask applied to the glacier to determine regions with no ice.
• velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}: A dictionary that provides for each RGI ID either the path to the datacube or the RasterStack with velocity data.

Returns

• glacier: An initialized glacier object containing the initial topography and climate data.

initialize_glaciers(
    rgi_ids::Vector{String},
    params::Parameters;
    velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}=Dict(),
)

Initialize glaciers based on provided RGI IDs and parameters.

Arguments

• rgi_ids::Vector{String}: A vector of RGI IDs representing the glaciers to be initialized.
• params::Parameters: A Parameters object containing simulation parameters.
• velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}: A dictionary that provides for each RGI ID either the path to the datacube or the RasterStack with velocity data.

Returns

• glaciers::Vector{Glacier2D}: A vector of initialized Glacier2D objects.

Description

This function performs the following steps:

1. Generates a file for missing glaciers if it does not already exist.
2. Filters out missing glaciers from the provided RGI IDs.
3. Generates raw climate data for the glaciers if necessary.
4. Initializes the glaciers using the provided RGI IDs and parameters.
5. If use_glathida_data is enabled in the simulation parameters, assigns GlaThiDa data to the glaciers.

Errors

• Throws an error if none of the provided RGI IDs have GlaThiDa data.

Warnings

• Issues a warning if not all glaciers have GlaThiDa data available.

Example

# We declare a list of glaciers to be initialized with their RGI IDs
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-08.00203"]
# We initialize those glaciers based on the RGI IDs and the parameters we previously specified
glaciers = initialize_glaciers(rgi_ids, params)

initialize_surfacevelocitydata(
    raster::Union{String, <: RasterStack};
    glacier::Union{G, Nothing}=nothing,
    mapping::VM=MeanDateVelocityMapping(),
    compute_vabs_error::Bool=true,
    flag::Union{String, <: RasterStack, Nothing} = nothing
) where {G <: AbstractGlacier, VM <: VelocityMapping}

Initialize SurfaceVelocityData from Rabatel et. al (2023).

Arguments:

• raster::Union{String, RasterStack}: RasterStack or path of the netCDF file with surface velocity data.
• glacier::Union{G, Nothing}: Glacier associated to the surface velocity datacube. When provided, the surface velocity data are gridded on the glacier grid using the mapping.
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• compute_vabs_error::Bool: Whether to compute the absolute error uncertainty.
• flag::Union{String, <: RasterStack, Nothing}: Option to provide a RasterStack containing a fmask raster and which provides an indicator of whether each pixel of the ice surface velocity data is considered as reliable or not.

initialize_surfacevelocitydata_mask(data::RasterStack, flag::Union{String, <:RasterStack, Nothing}=nothing)

Create a mask for a glacier surface velocity dataset by subsetting and aligning a flag raster to the spatial domain of data.

Arguments

• data::RasterStack: The glacier datacube or surface velocity dataset.
• flag::Union{String, RasterStack, Nothing}: Type of flag to be applied

Description

This function extracts the portion of the flag raster that spatially overlaps with the glacier domain defined by data. Because the flag raster is centered on pixel centers, its bounding box is shifted by half a grid step when subsetting.

is_in_glacier(A::Matrix{F}, distance::I) where {I <: Integer, F <: AbstractFloat}

Return a matrix with booleans indicating if a given pixel is at distance at least distance in the set of non zero values of the matrix. This usually allows discarding the border pixels of a glacier. A positive value of distanceindicates a measurement from inside the glacier, while a negativedistance`` indicates one from outside.

Arguments:

• A::Matrix{F}: Matrix from which to compute the matrix of booleans.
• distance::I: Distance to the border, computed as the number of pixels we need to move from within the glacier to find a pixel with value zero.

Compute distance between one point in the format of a TransverseMercator point and a set of points defined through the coordinates x and y in meters.

mapVelocity(
    velocityMapping::MeanDateVelocityMapping,
    velocityData::SurfaceVelocityData,
    t::AbstractFloat,
)

Retrieve the reference ice surface velocity for a given time step. This mapping uses the nearest snap shot available within a time window whose length is controlled by velocityMapping.thresDate. If no snapshot is found in the time window of length 2*velocityMapping.thresDate, the returned ice velocity components are empty matrices and the returned boolean flag useVel is set to false.

Arguments:

• velocityMapping::MeanDateVelocityMapping: Mapping to map the reference ice velocity to a target time step t.
• velocityData::SurfaceVelocityData: Surface velocity data. This is usually an attribute of a glacier.
• t::AbstractFloat: Current time step.

Returns

• Vx_ref: Matrix of the x-component of the ice surface velocity.
• Vy_ref: Matrix of the y-component of the ice surface velocity.
• V_ref: Matrix of the ice surface velocity magnitude.
• useVel: Boolean indicating whether the returned ice surface velocity can be used or not. The value of this boolean depends on the success of the ice surface velocity mapping at the current time step t.

max_or_empty(A::Array)

Return maximum value for non-empty arrays. This is just required to compute the error in the absolute velocity.

parse_proj(proj::String)

Parses the string containing the information of the projection to filter for important information "+proj=tmerc +lat0=0 +lon0=6.985 +k=0.9996 +x0=0 +y0=0 +datum=WGS84 +units=m +no_defs"

partial_year(float::Float64) -> Float64

Calculate the partial year value based on the given floating-point number.

Arguments

• float::Float64: A floating-point number representing the fraction of the year.

Returns

• Float64: The calculated partial year value.

partial_year(period::Type{<:Period}, float)

Calculate a partial year date based on a floating-point year value.

Arguments

• period::Type{<:Period}: The type of period to use (e.g., Month, Day).
• float::Float64: The floating-point year value.

Returns

• Date: The calculated date corresponding to the partial year.

plot_bias(
    results,
    variables;
    treshold = [0, 0],
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the bias of the glacier integrated volume over the specified time span.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: The variables to be plotted.
• title_mapping::Dict{Symbol, String}: A dictionary mapping variable names to their titles.
• tspan::Tuple{Float64, Float64}: A tuple representing the start and end time for the simulation.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier integrated volume bias.

plot_glacier(results::Results, plot_type::String, variables::Vector{Symbol}; kwargs...) -> Figure

Generate various types of plots for glacier data.

Arguments

• results::Results: The results object containing the data to be plotted.

• plot_type::String: Type of plot to generate. Options are:

  • "heatmaps": Heatmaps for glacier variables like :H, :H₀, :S, :B, :V, :Vx, :Vy, :V_ref.
  • "evolution difference": Temporal difference metrics (between start and end) for a variable, with optional metrics like "hist" (histogram) and "difference".
  • "evolution statistics": Temporal statistical metrics for a variable, with optional metrics like "average", "median", "min", "max", and "std".
  • "integrated volume": Temporal evolution of the integrated ice volume for a variable.
  • "bias": Scatter plot to visualize the bias between two variables.

• variables::Vector{Symbol}: Variables to be plotted, e.g., :H.

Optional Keyword Arguments

• tspan: A tuple representing the start and end time for the simulation.
• metrics: Metrics to visualize, e.g., ["average"] for statistics, ["difference"] for difference.
• scale_text_size::Union{Nothing,Float64}: Optional argument to scale the text size for heatmaps.
• threshold::Vector{F}: Threshold values for filtering data in bias plots.
• figsize::Tuple{Int64, Int64}: Size of the figure.

Returns

• A Figure object containing the desired visualization.

Notes

• Ensure the variables and kwargs match the requirements of the specified plot_type.
• The function routes requests to specific plotting functions based on plot_type.

plot_glacier_difference_evolution(
    results::Results,
    variables::Vector{Symbol},
    title_mapping;
    tspan::Tuple{F,F}=results.tspan,
    metrics::Vector{String}="difference",
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
) where {F<:AbstractFloat}

Plot the evolution of the difference in a glacier variable over time.

Arguments

• results::Results: The simulation results object containing the data to be plotted.
• variables::Vector{Symbol}: The variable to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• tspan::Tuple{F,F}: A tuple representing the start and end time for the simulation.
• metrics::Vector{String}: Metrics to visualize, e.g., ["difference"].
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier difference evolution.

plot_glacier_heatmaps(
    results::Results,
    variables::Vector{Symbol},
    title_mapping::Dict;
    scale_text_size::Union{Nothing,Float64}=nothing,
    timeIdx::Union{Nothing,Int64}=nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    plotContour::Bool=false,
) -> Figure

Plot heatmaps for glacier variables.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping::Dict: A dictionary mapping variable names to their titles and colormaps.
• scale_text_size::Union{Nothing,Float64}: Optional argument to scale the text size.
• timeIdx::Union{Nothing,Int64}:: Optional argument to select the index at which data should be plotted when dealing with vector of matrix. Default is nothing which selects the last element available.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.
• plotContour::Bool: Whether to add a contour plot representing the glacier borders at the beginning of the simulation on top of each of the figures. Default is false.

Returns

• A plot of the glacier heatmaps.

plot_glacier_integrated_volume(
    results,
    variables,
    title_mapping;
    tspan,
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the integrated volume of a glacier variable over time.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: The variable to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• tspan: A tuple representing the start and end time for the simulation.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier integrated volume.

plot_glacier_quivers(
    results::Results,
    variables::Vector{Symbol},
    title_mapping::Dict;
    timeIdx::Union{Nothing,Int64} = nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    lengthscale::Float64 = 0.00001,
    tiplength::Float64 = 0.5,
) -> Figure

Plot quivers for glacier variables.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping::Dict: A dictionary mapping variable names to their titles and colormaps.
• timeIdx::Union{Nothing,Int64}:: Optional argument to select the index at which data should be plotted when dealing with vector of matrix. Default is nothing which selects the last element available.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.
• lengthscale::Float64: Lengthscale of the arrows in the quiver plot.
• tiplength::Float64: Length of the arrow in the quiver plot.

Returns

• A plot of the glacier quivers.

plot_glacier_statistics_evolution(
    results::Results,
    variables::Vector{Symbol},
    title_mapping;
    metrics="median",
    tspan,
    threshold=0.5,
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the evolution of statistics for multiple glacier variables over time.

Arguments

• results::Results: The simulation results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• metrics: Metrics to visualize, e.g., "average", "median", "min", "max", and "std". Default is "median".
• tspan: A tuple representing the start and end time for the simulation.
• threshold: A threshold value to filter the data. Default is 0.5.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier statistics evolution.

plot_glacier_vid(
    plot_type::String,
    results::Results,
    glacier::Glacier2D,
    tspan,
    step,
    pathVideo::String;
    framerate::Int=24,
    baseTitle::String=""
)

Generate various types of videos for glacier data. For now only the evolution of the glacier ice thickness is supported. More types of visualizations will be added in the future.

Arguments

• plot_type: Type of plot to generate. Options are:

  • "thickness": Heatmap of the glacier thickness.

• results: A result object containing the simulation results including ice thickness over time.

• glacier: A glacier instance.

• tspan: The simulation time span.

• step: Time step to use to retrieve the results and generate the video.

• pathVideo: Path of the mp4 file to generate.

Optional Keyword Arguments

• framerate: The framerate to use for the video generation.
• baseTitle: The prefix to use in the title of the frames. In each frame it is concatenated with the value of the year in the form " (t=XXXX)".

plot_gridded_data(
    gridded_data::Union{Vector{Matrix{F}}, Matrix{F}},
    results::Results;
    scale_text_size::Union{Nothing,Float64}=nothing,
    timeIdx::Union{Nothing,Int64}=nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    plotContour::Bool=false,
    colormap = :cool,
    logPlot = false,
) where {F <: AbstractFloat}

Plot a gridded matrix (or a time series of matrices) as a heatmap using metadata from results.

Arguments

• gridded_data::Union{Vector{Matrix{F}}, Matrix{F}}: Single snapshot or time series (defaults to last timestep).
• results::Results: Supplies lon, lat, x, y, rgi_id, Δx and H (mask).
• scale_text_size, figsize, colormap: Optional plotting params.
• timeIdx::Union{Nothing,Int64}: Select timestep when gridded_data is a vector.
• plotContour::Bool: overlay glacier-mask contour from results.H.
• logPlot::Bool: Use log10 colorscale (positive non-NaN values determine range).

Behavior

• Masks out cells where results.H[begin] .<= 0 (set to NaN).
• Adds colorbar, central lon/lat tick, and a Δx-based scale bar in km.
• If plotContour, draws mask boundary lines.
• Returns a CairoMakie.Figure.

Errors

• Asserts gridded_data is non-empty and timeIdx (if provided) is in range.

precompute_all_VJPs_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)