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 = Huginn.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
Huginn.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
pdiff = 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 Prediction from the previous example (i.e. a forward run). This first tutorial keeps things simple, and since we are not using machine learning models, we will only use the Model type to specify the iceflow and mass balance models. These functionalities are mainly covered by 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 come with a default constructor, which will provide default values in case you don't want to tune those. The main types of parameters are:

  • Simulation parameters: SimulationParameters includes all the parameters related to the ODINN.jl simulation, including the number of workers, the timespan of the simulation or the working directory.
  • Hyperparameters: Hyperparameters includes all the necessary hyperparameters for a machine learning model.
  • UDEparameters: 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 need to specify a list of RGI IDs of the glacier we want to work with. Specifying an RGI region is also possible. 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{ODINN.SciMLSensitivityAdjoint}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(true, true, true, true, 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, 0.1643835616438356), 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}
  alpha: Float64 1.0
  scaled: Bool false
, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}
  delta: Float64 0.1
  sigma: Float64 0.9
  alphamax: Float64 Inf
  rho: Float64 5.0
  epsilon: Float64 1.0e-6
  gamma: Float64 0.66
  linesearchmax: Int64 50
  psi3: Float64 0.1
  display: Int64 0
  mayterminate: Base.RefValue{Bool}
  cache: Nothing nothing
  check_flatness: Bool 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, nothing, false, true, 10, 100000), UDEparameters{ODINN.SciMLSensitivityAdjoint}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), ADTypes.AutoEnzyme(), ODINN.SciMLSensitivityAdjoint(), "AD+AD", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :A), 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}
  alpha: Float64 1.0
  scaled: Bool false
, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}
  delta: Float64 0.1
  sigma: Float64 0.9
  alphamax: Float64 Inf
  rho: Float64 5.0
  epsilon: Float64 1.0e-6
  gamma: Float64 0.66
  linesearchmax: Int64 50
  psi3: Float64 0.1
  display: Int64 0
  mayterminate: Base.RefValue{Bool}
  cache: Nothing nothing
  check_flatness: Bool 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 three 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.
  • 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.
  • Machine learning model: MLmodel is the machine learning model (e.g. a neural network) which will be used as part of a hybrid model based on a Universal Differential Equation.

The model is initialized using the Model constructor:

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

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 OGGM, mostly providing a Julia interface to automatize this. 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)

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 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, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}, Array{Float64, 0}}, Nothing}}(Sleipnir.Model{SIA2Dmodel{Float64, ConstantLaw{Array{Float64, 0}, Huginn.var"#7#8"}, ConstantLaw{Array{Float64, 0}, Huginn.var"#9#10"}, ConstantLaw{Array{Float64, 0}, Huginn.var"#11#12"}, NullLaw, NullLaw}, TImodel1{Float64}, Nothing}(, , nothing), nothing, , Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ODINN.SciMLSensitivityAdjoint}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(true, true, true, true, 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, 0.1643835616438356), 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}
  alpha: Float64 1.0
  scaled: Bool false
, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}
  delta: Float64 0.1
  sigma: Float64 0.9
  alphamax: Float64 Inf
  rho: Float64 5.0
  epsilon: Float64 1.0e-6
  gamma: Float64 0.66
  linesearchmax: Int64 50
  psi3: Float64 0.1
  display: Int64 0
  mayterminate: Base.RefValue{Bool}
  cache: Nothing nothing
  check_flatness: Bool 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, nothing, false, true, 10, 100000), UDEparameters{ODINN.SciMLSensitivityAdjoint}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), ADTypes.AutoEnzyme(), ODINN.SciMLSensitivityAdjoint(), "AD+AD", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :A), 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}
  alpha: Float64 1.0
  scaled: Bool false
, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}
  delta: Float64 0.1
  sigma: Float64 0.9
  alphamax: Float64 Inf
  rho: Float64 5.0
  epsilon: Float64 1.0e-6
  gamma: Float64 0.66
  linesearchmax: Int64 50
  psi3: Float64 0.1
  display: Int64 0
  mayterminate: Base.RefValue{Bool}
  cache: Nothing nothing
  check_flatness: Bool false
, nothing, nothing, Optim.Flat()))), Sleipnir.Results[])

And once we have the Prediction object, we can run it using Huginn.run!:

Huginn.run!(prediction)
[ Info: Running forward in-place PDE ice flow model
      From worker 2:	Processing glacier RGI60-11.03638 for PDE forward simulation
      From worker 4:	Processing glacier RGI60-11.01450 for PDE forward simulation
      From worker 3:	Processing glacier RGI60-11.02346 for PDE forward simulation
      From worker 3:	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:02

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

This page was generated using Literate.jl.