GitHub

Forward simulation tutorial

This tutorial provides a simple example on how to perform a forward simulation using ODINN.jl.

Running the whole code

using ODINN

# Define the working directory
working_dir = joinpath(ODINN.root_dir, "demos")

# Ensure the working directory exists
mkpath(working_dir)

# Define which glacier RGI IDs we want to work with
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-08.00203"]
rgi_paths = get_rgi_paths()

# Create the necessary parameters
params = Parameters(
    simulation = SimulationParameters(
    working_dir = working_dir,
    tspan = (2010.0, 2015.0),
    multiprocessing = true,
    workers = 4,
    rgi_paths = rgi_paths
)
)

# Specify a model based on an iceflow model, a mass balance model,
# and a machine learning model
model = Model(
    iceflow = SIA2Dmodel(params),
    mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0)
)

# We initialize the glaciers with all the necessary data
glaciers = initialize_glaciers(rgi_ids, params)

# We specify the type of simulation we want to perform
prediction = Prediction(model, glaciers, params)

# And finally, we just run the simulation
run!(prediction)

# Then we can visualize the results of the simulation, e.g. the difference in ice thickness
# between 2010 to 2015 for Argentière glacier
plot_glacier(prediction.results[1], "evolution difference", [:H]; metrics = ["difference"])

Step-by-step explanation of the tutorial

Here we will cover in detail each one of the steps that lead us to run the forward models Prediction from the previous example. This first tutorial keeps things simple, and since we are not using machine learning models, we will only use Model to specify the iceflow and mass balance models. These functionalities are mainly included in Huginn.jl.

Step 1: Parameter initialization

The first step is to initialize and specify all the necessary parameters. In ODINN.jl we have many different types of parameters, specifying different aspects of the model. All the parameter types have a default constructor, which provide default values. The main types of parameters are:

  • Simulation parameters: SimulationParameters includes all the parameters related to the simulation, including the number of workers, the timespan of the simulation or the working directory.
  • Hyper parameters: Hyperparameters includes all the necessary hyperparameters to optimize the model.
  • UDE parameters: UDEparameters contains the parameters related to the training of a Universal Differential Equation.

All these sub-types of parameters are held in a Parameters struct, a general parameters structure to be passed to an ODINN simulation.

First, we specify a list of RGI IDs (associated to the Randolph Glacier Inventory) of the glaciers we want to work with. From these RGI IDs, we will look for the necessary files inside the workspace.

rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-08.00203"]
rgi_paths = get_rgi_paths()

params = Parameters(
    simulation = SimulationParameters(
    working_dir = working_dir,
    tspan = (2010.0, 2015.0),
    multiprocessing = true,
    workers = 4,
    rgi_paths = rgi_paths
)
)
Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.5e-20, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(true, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, true, 4, "/home/runner/work/ODINN.jl/ODINN.jl/demos", false, Dict("RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897", "RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213", "RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147", "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270", "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646", "RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232", "RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174", "RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274", "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207", "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"…), "Farinotti19", MeanDateVelocityMapping(:nearest), 1), Hyperparameters{Float64, Int64}(1, 1, Float64[], Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Float64, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, 0.001, Optim.Flat()), 0.0, 50, 15), SolverParameters{Float64, Int64}(OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-12, 0.08333333333333333, Float64[], false, true, 10, 100000), UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), ADTypes.AutoEnzyme(), ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}(DiscreteVJP{ADTypes.AutoMooncake{Nothing}}(ADTypes.AutoMooncake()), OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-8, 0.08333333333333333, :Linear, 200, EnzymeVJP()), "AD+AD", MultiLoss{Tuple{LossH{L2Sum{Int64}}}, Vector{Float64}}((LossH{L2Sum{Int64}}(L2Sum{Int64}(3)),), [1.0]), :A, :identity), InversionParameters{Float64}([1.0], [0.0], [Inf], [1, 1], 0.001, 0.001, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, nothing, Optim.Flat())))

Step 2: Model specification

The next step is to specify which model(s) we want to use for our simulation. In ODINN we have two different types of model, which are encompassed in a Model structure:

  • Iceflow model: IceflowModel is the ice flow dynamics model that will be used to simulate iceflow. It defaults to a 2D Shallow Ice Approximation (SIA). Check out this glaciology notebook for a very good introduction to the Shallow Ice Approximation.
  • Surface mass balance model: MassBalanceModel is the mass balance model that will be used for simulations. Options here include temperature-index models, or machine learning models coming from MassBalanceMachine.

Trainable components can be embedded inside the iceflow model which can be a neural network to learn a parameterization in the context of Universal Differential Equation, or per glacier values in the context of a classical inversion (like inverting the initial conditions). Refer to the functional inversion tutorial for an example of how to incorporate a neural network inside the iceflow model.

The model is initialized using the Model constructor:

model = Model(
    iceflow = SIA2Dmodel(params),
    mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0)
)
**** Model ****

SIA2D iceflow equation  = ∇(D ∇S)  with D = U H̄
  and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
      Γ = 2A (ρg)^n /(n+2)
      A: ConstantLaw -> Array{Float64, 0}
      C: ConstantLaw -> Array{Float64, 0}
      n: ConstantLaw -> Array{Float64, 0}
      p: ConstantLaw -> Array{Float64, 0}
      q: ConstantLaw -> Array{Float64, 0}

Temperature index mass balance model TImodel1
   DDF = 0.006
   acc_factor = 0.0012
No learnable components
***************

Step 3: Glacier initialization

The third step is to fetch and initialize all the necessary data for our glaciers of interest. This is strongly built on top of the Open Global Glacier Model (OGGM) Python package, mostly providing a Julia interface to automatize preprocessing of the glacier data. The package Gungnir is used to fetch the necessary data from the RGI and other sources. The data is then stored in servers and fetched and read using Rasters.jl directly by Sleipnir.jl when needed.

Then, we initialize those glaciers based on those RGI IDs and the parameters we previously specified.

glaciers = initialize_glaciers(rgi_ids, params)
4-element Vector{Glacier2D} distributed over regions 08 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-08.00203

Step 4: Creating and running a simulation

The final step of the pipeline, is to create an ODINN simulation based on all the previous steps, and then to run it. There are different types of simulations that we can do with ODINN:

  • Prediction: This is a forward simulation, where the initial glacier conditions are run forward in time based on specified parameters and climate data.

This is as simple as doing:

prediction = Prediction(model, glaciers, params)
Prediction{Sleipnir.ModelCache{SIA2DCache{Float64, Int64, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, Array{Float64, 0}, Array{Float64, 0}, ScalarCacheNoVJP, ScalarCacheNoVJP}, Nothing}}(Sleipnir.Model{SIA2Dmodel{Float64, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#7#8"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#9#10"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#11#12"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#13#14"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#15#16"}, NullLaw, NullLaw}, TImodel1{Float64}, Nothing}(SIA2D iceflow equation  = ∇(D ∇S)  with D = U H̄
  and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
      Γ = 2A (ρg)^n /(n+2)
      A: ConstantLaw -> Array{Float64, 0}
      C: ConstantLaw -> Array{Float64, 0}
      n: ConstantLaw -> Array{Float64, 0}
      p: ConstantLaw -> Array{Float64, 0}
      q: ConstantLaw -> Array{Float64, 0}
, Temperature index mass balance model TImodel1
   DDF = 0.006
   acc_factor = 0.0012, nothing), nothing, 4-element Vector{AbstractGlacier} distributed over regions 08 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-08.00203
, Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.5e-20, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(true, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, true, 4, "/home/runner/work/ODINN.jl/ODINN.jl/demos", false, Dict("RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897", "RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213", "RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147", "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270", "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646", "RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232", "RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174", "RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274", "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207", "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"…), "Farinotti19", MeanDateVelocityMapping(:nearest), 1), Hyperparameters{Float64, Int64}(1, 1, Float64[], Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Float64, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, 0.001, Optim.Flat()), 0.0, 50, 15), SolverParameters{Float64, Int64}(OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-12, 0.08333333333333333, Float64[], false, true, 10, 100000), UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), ADTypes.AutoEnzyme(), ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}(DiscreteVJP{ADTypes.AutoMooncake{Nothing}}(ADTypes.AutoMooncake()), OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-8, 0.08333333333333333, :Linear, 200, EnzymeVJP()), "AD+AD", MultiLoss{Tuple{LossH{L2Sum{Int64}}}, Vector{Float64}}((LossH{L2Sum{Int64}}(L2Sum{Int64}(3)),), [1.0]), :A, :identity), InversionParameters{Float64}([1.0], [0.0], [Inf], [1, 1], 0.001, 0.001, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, nothing, Optim.Flat()))), Sleipnir.Results[])

And once we have the Prediction object, we can run the simulation using the function run!:

run!(prediction)
[ Info: Running forward in-place PDE ice flow model
      From worker 3:	Processing glacier RGI60-11.03638 for PDE forward simulation
      From worker 2:	Processing glacier RGI60-11.01450 for PDE forward simulation
      From worker 4:	Processing glacier RGI60-11.02346 for PDE forward simulation
      From worker 4:	Processing glacier RGI60-08.00203 for PDE forward simulation

Progress:  50%|████████████████████▌                    |  ETA: 0:00:00
Progress:  75%|██████████████████████████████▊          |  ETA: 0:00:01
Progress: 100%|█████████████████████████████████████████| Time: 0:00:03

There we go, we have successfully simulated the evolution of 3 glaciers for 5 years in around 1-2 seconds!

Step 5: Visualizing the results

Finally, we can use the plotting functions of ODINN.jl to visualize the results of the simulation. Like the glacier ice thickness evolution:

plot_glacier(prediction.results[1], "evolution difference", [:H]; metrics = ["difference"])

Or the initial glacier ice thickness and the resulting ice surface velocities:

plot_glacier(prediction.results[1], "heatmaps", [:H, :V])

We can also visualize the results for other glaciers, like Aletsch:

plot_glacier(prediction.results[2], "evolution difference", [:H]; metrics = ["difference"])

This page was generated using Literate.jl.