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

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 for of the forward PDE.

ContinuousVJP

DiscreteAdjoint{VJP <: 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.

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

DiscreteVJP

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

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

Constructor for FunctionalInversion 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

• FunctionalInversion: A new instance of the FunctionalInversion struct.

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

An object representing a functional inversion simulation (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 functional inversion and of the forward simulations.

FunctionalInversionBinder{FI <: FunctionalInversion, CA <: ComponentArray}

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

Fields

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

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.

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

Create an Inversion object using the provided model, glaciers, and parameters.

Arguments

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

Returns

• inversion: An Inversion object initialized with the provided model, glaciers, and parameters.

Inversion <: Simulation

A mutable struct that represents an inversion simulation.

Fields

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

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.

mutable struct InversionResults{F <: Real}

A structure to store the results of an inversion simulation.

Fields

• rgi_id::Union{String, Nothing}: The RGI identifier, which can be a string or nothing.
• A::F: Area parameter.
• n::F: Exponent parameter.
• C::Matrix{F}: Coefficient matrix.
• H_pred::Matrix{F}: Predicted height matrix.
• H_obs::Matrix{F}: Observed height matrix.
• H_diff::Matrix{F}: Difference between predicted and observed height matrices.
• V_pred::Matrix{F}: Predicted volume matrix.
• V_obs::Matrix{F}: Observed volume matrix.
• V_diff::Matrix{F}: Difference between predicted and observed volume matrices.
• MSE::F: Mean squared error.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• lon::F: Longitude.
• lat::F: Latitude.

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.

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

Feed-forward neural network.

Fields

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

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

Mutable struct containing the results of a functional 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 a functional 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.

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.

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

SciMLSensitivityAdjoint

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 SciMLSensitivityAdjoint().
• 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.

Base.copyto!(
    dest::FunctionalInversionBinder{FI, CA},
    src::FunctionalInversionBinder{FI, CA},
) where {FI <: FunctionalInversion, 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::FunctionalInversionBinder{FI, CA},
) where {FI <: FunctionalInversion, CA <: ComponentArray}

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

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(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters,
)

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.

Returns

• A_law: A Law{Array{Float64, 0}} 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 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 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.
• The init_cache function initializes the cache with a scalar zero.

Example

```julia nnmodel = NeuralNetwork(params) Alaw = LawA(nn_model, params)

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.

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

```julia nnmodel = NeuralNetwork(params) boundsH = (0.0, 300.0) bounds∇S = (0.0, 0.5) Ulaw = LawU(nnmodel, params; maxNN = 50.0, prescalebounds = [boundsH, 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

```julia nnmodel = NeuralNetwork(params) boundsT = (-25.0, 0.0) boundsH = (0.0, 500.0) Ylaw = LawY(nnmodel, params; prescalebounds = [boundsT, boundsH])

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

Model(;
    iceflow::Union{IFM, Nothing},
    mass_balance::Union{MBM, Nothing},
    regressors::NamedTuple = (;),
    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.
• machine_learning::Union{MLM, Nothing}: The machine learning model to be used. Can be a single model or nothing.

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

```julia T, A = TAAlaw(simulation, glacier_idx, θ, 2010.0)

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

Compute an out-of-place adjoint step of the Shallow Ice Approximation PDE. Given an output gradient, it backpropagates the gradient to the inputs H and A. To some extent, this function is equivalent to VJPλ∂SIA∂Hcontinuous and VJPλ∂SIA∂θcontinuous.

Arguments:

• ∂dH::Matrix{R}: Output gradient to backpropagate.
• 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:

• ∂H::Matrix{R}: Input gradient wrt H.
• ∂A::F: Input gradient wrt A.

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

Implementation of the continuous adjoint 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 adjoint 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:

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

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.

_batch_iceflow_UDE(
    container::FunctionalInversionBinder,
    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

```julia mymodel = StatefulLuxLayer{true}(archi, nothing, st) y = predNN([1.0, 2.0], mymodel, θ, prescalefn, postscalefn)

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

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

Arguments:

• container::FunctionalInversionBinder: 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

function create_interpolation(A::Matrix; n_interp_half::Int)

Function to create an intepolation for AD computation combining uniform and quantiles.

create_results(θ, simulation::FunctionalInversion, 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::FunctionalInversion: Simulation structure that contains all the required information about the functional inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

define_iceflow_prob(
    simulation::FunctionalInversion,
    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.
• inputs::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

```julia result = evallaw(simulation.model.iceflow.A, simulation, glacieridx, (; T=273.15), θ)

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

```julia 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::FunctionalInversion)

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

Arguments

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

Returns

• A vector of FunctionalInversion 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 FunctionalInversion 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.
• If the simulation results are empty, the function creates batches with empty results objects.

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::FunctionalInversion, 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::FunctionalInversion: Simulation structure that contains all the required information about the functional inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

grad_loss_iceflow!(θ, simulation::FunctionalInversion, 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::FunctionalInversion, 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=1.,
) 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_iceflow_transient(θ, simulation::FunctionalInversion, 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::FunctionalInversion: Simulation structure that contains all the required information about the functional 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::FunctionalInversion)

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.

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

```julia 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.

run!(simulation::FunctionalInversion)

Run the training process for a given FunctionalInversion simulation.

Arguments

• simulation::FunctionalInversion: The simulation object containing the parameters and settings for the functional 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.

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::FunctionalInversionBinder,
    cb::SciMLBase.DECallback,
    iceflow_prob::ODEProblem,
)

Make a forward simulation of the iceflow UDE.

BFGS Training

trainUDE!(simulation::FunctionalInversion; saveevery_iter::Bool=false, logger::Union{<: TBLogger, Nothing}=nothing)

Trains UDE based on the current FunctionalInversion.

ADAM Training

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.

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.

InpH̄ <: 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.

InpTemp <: AbstractInput

Input that represents the long term air temperature of a glacier. It is computed using the OGGM data over a period predefined in Gungnir.

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 = (avgy(diffx(S) / Δx).^2 .+ avgx(diffy(S) / Δy).^2).^(1/2)

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.

struct SIA2DCache{R <: Real, I <: Integer, A_CACHE, C_CACHE, n_CACHE, Y_CACHE, U_CACHE, ∂A∂θ_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.
• ∂A∂θ_CACHE: Type used for caching ∂A∂θ in the VJP computation which is a scalar.
• ∂Y∂θ_CACHE: Type used for caching ∂Y∂θ in the VJP computation which is a scalar.
• ∂U∂θ_CACHE: Type used for caching ∂U∂θ in the VJP computation which is a scalar.

Fields

• A::A_CACHE: Flow rate factor.
• n::n_CACHE: Flow 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.
• C::C_CACHE: Sliding coefficient.
• Y::Y_CACHE: Hybrid diffusivity.
• U::U_CACHE: Diffusive velocity.
• ∂A∂H::A_CACHE: Buffer for VJP computation.
• ∂A∂θ::∂A∂θ_CACHE: Buffer for VJP computation.
• ∂Y∂H::Y_CACHE: Buffer for VJP computation.
• ∂Y∂θ::∂Y∂θ_CACHE: Buffer for VJP computation.
• ∂U∂H::U_CACHE: Buffer for VJP computation.
• ∂U∂θ::∂U∂θ_CACHE: Buffer for VJP computation.
• 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.

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.
• 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_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::OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm = RDPK3Sp35(),
                  reltol::F = 1e-12,
                  step::F = 1.0/12.0,
                  tstops::Union{Nothing,Vector{F}} = nothing,
                  save_everystep = false,
                  progress::Bool = true,
                  progress_steps::I = 10,
                  maxiters::I = Int(1e5),
                ) where {F <: AbstractFloat, I <: Integer}

Arguments

• solver::OrdinaryDiffEq.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 for the callbacks. These are mainly used to run the surface mass balance model. Defaults to 1.0/12.0 (i.e. a month).
• tstops::Union{Nothing, Vector{F}}: Optional vector of time points where the solver should stop. Defaults to nothing.
• save_everystep::Bool: Whether to save the solution at every step. 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::OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm: The algorithm used for solving differential equations.
• reltol::F: The relative tolerance for the solver.
• step::F: The step size for the solver.
• tstops::Union{Nothing, Vector{F}}: Optional vector of time points where the solver should stop for the callbacks.
• save_everystep::Bool: Flag indicating whether to save the solution at every step.
• 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.

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

Law that represents a constant A in the SIA.

Arguments:

• A::F: Rheology factor A.

CuffeyPaterson()

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:

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

Model(;
    iceflow::Union{IFM, Nothing},
    mass_balance::Union{MBM, Nothing}
) where {IFM <: IceflowModel, MBM <: MBmodel}

Initialize Model at Huginn level (no machine learning model).

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

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.

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, θ) -> 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. If such a frequency is set (via callback_freq), the update is done using a PeriodicCallback. Otherwise, no callback is used for that component.

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{Float64, Float64}: A tuple representing the start and end times.
• step::Float64: The step size for generating the callback steps.

Returns

• Vector{Float64}: 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},
) 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 field.

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

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.

Example

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

glaciers = generate_ground_truth(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.

run!(simulation::Prediction)

In-place run of the model.

function simulate_iceflow_PDE!(
    simulation::SIM,
    cb::DiscreteCallback,
    du
) where {SIM <: Simulation}

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

Description

This function iterates over the glaciers in the Prediction object and generates the simulated thickness data (H) and corresponding time steps (t). I then computes the surface ice velocity data. 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.

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::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::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.
• batch_id::Union{Nothing, I}: Optional batch identifier for simulations using Reverse Diff. Defaults to nothing.

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.

If batch_id is provided, the function updates the mass balance for the specified batch; otherwise, it updates the mass balance for the entire model.

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)

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

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

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.

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{T}(init_cache)

Creates a constant law of type Law{T} 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

• T: The type of the internal state. 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)

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 is not 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; smoothing=false, test=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.
• 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,
    slope::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    dist_border::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    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.
• slope::Matrix{F}: Slope matrix.
• dist_border::Matrix{F}: Distance to border matrix.
• 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.
• slope::Matrix{F}: Surface slope matrix.
• dist_border::Matrix{F}: Distance to the glacier border matrix.
• 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}(; f!, init_cache, callback_freq=nothing, inputs=nothing)

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, AbstractFloat}: Optional. If provided, the law is treated as a callback law and is only applied every callback_freq time units.
• inputs::Union{Nothing, Tuple{<:AbstractInput}}: Optional. Provides automatically generated inputs passed to f! at runtime.

Type Parameters

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

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,
)

# TODO: create a meaningful example of Law with glacier-shaped 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.
• thresDate::Float: Threshold in days to accept or reject a surface velocity snapshot. This allows selecting surface velocity data that is close enough to the current time step in the loss function.

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

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

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.
• machine_learning::MLM: Represents the machine learning component, which is an instance of MLM.

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.
• MLM: A subtype of AbstractEmptyModel representing the type of the machine learning model.

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,
    maxTlaw::Float64 = 1.0,
    minTlaw::Float64 = -25.0,
    noise_A_magnitude::Float64 = 5e-18
    )

Keyword arguments

- `ρ`: Density of ice.
- `g`: Gravitational acceleration.
- `ϵ`: Regularization used in the square root of norms for AD numerical stability.
- `η₀`: Initial viscosity.
- `maxA`: Maximum value for `A` (Glen's coefficient).
- `minA`: Minimum value for `A` (Glen's coefficient).
- `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 value for A (Glen's coefficient).
• minA::F: Minimum value for A (Glen's coefficient).
• maxTlaw::F: Maximum value of Temperature used in simulations on fake law.
• minTlaw::F: Minimum value of Temperature used in simulations on fake 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),
    θ::Union{Nothing,ComponentArray{F}} = nothing,
    loss::Union{Nothing,Vector{F}} = nothing
) 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.
• θ::Union{Nothing, ComponentArray{F}}: Machine learning model parameters.
• loss::Union{Nothing, Vector{F}} Vector with evolution of loss function.

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,
    overwrite_climate::Bool = false,
    use_glathida_data::Bool = false,
    tspan::Tuple{F, F} = (2010.0,2015.0),
    step::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).
• 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::F: Time step for the 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.
• 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::F: Time step for the 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, isGridGlacierAligned::Bool = false, ) where {F <: AbstractFloat}

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

Important remarks:

• 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 component of surface velocity.
• vy::Union{Vector{Matrix{F}}, Nothing}: y component of surface velocity.
• 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.
• 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

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

create_results(
    simulation::SIM,
    glacier_idx::I,
    solution,
    loss = nothing;
    light = false,
    tstops::Union{Vector, Nothing} = nothing,
    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.
• loss=nothing: The loss value, default is nothing.
• light=false: A boolean flag to indicate if only the first and last steps of the solution should be used.
• 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

• 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

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.

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

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

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.
   • Triggers garbage collection to free up memory.

get_cumulative_climate!(climate, period; gradient_bounds=[-0.009, -0.003])

Calculate and update the cumulative climate data for a given period.

Keyword arguments

• climate::Climate: The climate object containing raw climate data.
• period::Period: The period for which to calculate the cumulative climate data.
• gradient_bounds::Vector{Float64}: Optional. The bounds within which to clamp the gradient values. Default is [-0.009, -0.003].

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.
• 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.
• test::Bool: An optional boolean flag indicating whether to run in test mode. Default is false.
• 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
) 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.

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.

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 to find a pixel with value zero.

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

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,
    arrowsize::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.
• arrowsize::Float64: Arrow size 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,
    H::Vector{Matrix{Float64}},
    glacier::Glacier2D,
    simuparams::SimulationParameters,
    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.
• H: A vector of matrices containing the ice thickness over time. This should be replaced by a Results instance in the future once Results no longer depends on an iceflow model.
• glacier: A glacier instance.
• simuparams: The simulation parameters.
• 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)".

ratio_max(v, vabs)

Compute the maximum ratio between v and vabs at points where the value of vabs is not a NaN.

reverseForHeatmap(
    inp::Matrix{F},
    x::Vector{F},
    y::Vector{F}
) where {F <: AbstractFloat}

Out-of-place reverse of a matrix based on the values of the x and y axes. This function corrects the orientation so that the heatmap is displayed correctly.

Arguments

• inp::Matrix{F}: The matrix to reverse.
• x::Vector{F}: Values of the x axis.
• y::Vector{F}: Values of the y axis.

Returns

• Out-of-place copy of inp that has been reversed if needed.

save_results_file!(results_list::Vector{Results{F, I}}, simulation::SIM; path::Union{String,Nothing}=nothing) where {F <: AbstractFloat, I <: Integer, SIM <: Simulation}

Save the results of a simulation to a file.

Arguments

• results_list::Vector{Results{F, I}}: A vector containing the results of the simulation.
• simulation::SIM: The simulation object containing the parameters and results.
• path::Union{String,Nothing}: Optional. The path where the results file will be saved. If not provided, a default path will be used.

Description

This function saves the results of a simulation 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.

smooth!(A)

Smooths the interior of a 2D array A using a simple averaging method. The function modifies the array A in place.

Arguments

• A::AbstractMatrix: A 2D array to be smoothed.

Details

The function updates the interior elements of A (excluding the boundary elements) by adding a weighted average of the second differences along both dimensions. The boundary elements are then set to the values of their nearest interior neighbors to maintain the boundary conditions.

stop_condition_tstops(u, t, integrator, tstops)

Check if the current time t is in the list of stop times tstops.

Arguments

• u: The current state of the system (not used in this function).
• t: The current time.
• integrator: The integrator object (not used in this function).
• tstops: A collection of times at which the integration should stop.

Returns

• Bool: true if t is in tstops, otherwise false.

trim_period(period, climate)

Adjusts the given period to fit within the bounds of the climate data, ensuring it aligns with hydrological years.

Arguments

• period::UnitRange{Date}: The initial date range to be trimmed.
• climate::AbstractArray: The climate data array, which should have a time dimension Ti.

Returns

• UnitRange{Date}: The adjusted date range that fits within the climate data's time bounds.

Details

• If the start of the climate data is later than the start of the period, the period is adjusted to start from October 1st of the year of the climate data's start.
• If the end of the climate data is earlier than the end of the period, the period is adjusted to end on September 30th of the year of the climate data's end.