Huginn.jl API

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

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

Parameters

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

Notes

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

IceflowModel

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

Prediction{CACHE} <: Simulation

A mutable struct that represents a prediction simulation.

Fields

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

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

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

Arguments

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

Returns

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

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

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

Type Parameters

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

Fields

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

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

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

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

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

• A sliding law governed by coefficient C,

• Optionally the user can provide either:

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

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

Arguments

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

SIAmodel

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

A mutable struct that holds parameters for the solver.

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

Fields

• solver::ST: The algorithm used for solving differential equations.
• reltol::F: The relative tolerance for the solver.
• step::F: The step size that controls at which frequency the results must be saved.
• tstops::Vector{F}: Optional vector of time points where the solver should stop to store the results.
• save_everystep::Bool: Flag indicating whether to save the solution at every step computed by the solver.
• progress::Bool: Flag indicating whether to show progress during the solving process.
• progress_steps::I: The number of steps between progress updates.
• maxiters::I: Maximum number of iterations to perform in the iceflow solver.

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

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

Arguments

• solver::OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm: The ODE solver algorithm to use. Defaults to RDPK3Sp35().
• reltol::F: The relative tolerance for the solver. Defaults to 1e-12.
• step::F: The step size that controls at which frequency the solution should be computed and returned in the results. Defaults to 1.0/12.0 (i.e. a month).
• tstops::Vector{Sleipnir.Float}: Optional vector of time points where the solver should stop. Defaults to an empty vector.
• save_everystep::Bool: Whether to save the solution at every step computed by the solver. Defaults to false.
• progress::Bool: Whether to show progress during the solving process. Defaults to true.
• progress_steps::I: The number of steps between progress updates. Defaults to 10.
• maxiters::I: Maximum number of iterations to perform in the iceflow solver. Defaults to 1e5.

Returns

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

iAvgGriddedTemp <: AbstractInput

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

iAvgScalarTemp <: AbstractInput

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

iCPDD <: AbstractInput

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

iH̄ <: AbstractInput

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

iTopoRough{F<:AbstractFloat} <: AbstractInput

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

i∇S <: AbstractInput

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

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

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

Law that represents a constant A in the SIA.

Arguments:

• A::F: Rheology factor A.

CuffeyPaterson(; scalar::Bool = true)

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

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

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

Arguments

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

Returns

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

Description

The solution has the form:

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

Halfar_velocity(halfar_params::HalfarParameters)

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

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

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

Arguments

• dH::Matrix{R}: Matrix to store the rate of change of ice thickness.
• H::Matrix{R}: Matrix representing the ice thickness.
• simulation::SIM: Simulation object containing model parameters and state.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Details

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

Notes

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

See also SIA2D

SyntheticC(params::Sleipnir.Parameters)

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

TemperateA()

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

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

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

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

Arguments:

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

Returns:

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

apply_MB_mask!(H, ifm::SIA2DCache)

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

Arguments:

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

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

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

Arguments:

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

Returns:

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

avg!(O, I)

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

Arguments

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

Details

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

avg(A::AbstractArray)

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

Arguments

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

Returns

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

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

Calculate the average surface velocity for a given simulation.

Arguments

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

Description

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

Notes

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

avg_x!(O, I)

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

Arguments

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

avg_x(A::AbstractArray)

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

Arguments

• A::AbstractArray: Input array.

Returns

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

avg_y!(O, I)

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

Arguments

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

avg_y(A::AbstractArray)

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

Arguments

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

Returns

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

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

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

Arguments:

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

Returns

• A Results instance that stores the iceflow solution.

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

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

Each law can optionally specify a callback frequency via callback_freq.

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

Arguments:

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

Returns:

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

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

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

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

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

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

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

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

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

Arguments

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

Returns

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

diff_x!(O, I, Δx)

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

Arguments

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

Notes

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

diff_x(A::AbstractArray)

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

Arguments

• A::AbstractArray: Input array.

Returns

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

diff_y!(O, I, Δy)

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

Arguments

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

Description

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

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

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

diff_y(A::AbstractArray)

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

Arguments

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

Returns

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

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

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

Arguments

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

Description

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

Example

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

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

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

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

Arguments

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

Description

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

Example

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

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

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

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

Arguments

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

Returns

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

inn(A::AbstractArray)

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

Arguments

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

Returns

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

inn1(A::AbstractArray)

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

Arguments

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

Returns

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

polyA_PatersonCuffey()

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

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

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

project_curvatures(H, eₚ, eₛ)

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

Arguments

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

Returns

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

run!(simulation::Prediction)

In-place run of the model.

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

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

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

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

Arguments

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

Returns

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

Description

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

Details

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

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

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

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

Arguments

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

Returns

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

Details

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

Notes

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

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

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

Arguments

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

Description

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

Notes

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

Returns

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

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

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

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

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

Arguments

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

Notes

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

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

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

Arguments

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

Notes

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

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

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

Keyword arguments

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