Classical inversion tutorial

This tutorial provides a simple example on how to perform a classical gridded inversion in ODINN.jl. For this, we generate a synthetic dataset using a forward simulation, and then we use this dataset to perform the classical inversion. The goal of this classical inversion is to retrieve the matrix of A values associated to the Glen coefficient/rigidity in Glen's Law that was used to generate the results of a forward simulation.

Step 1: Parameter and glacier initialization

using ODINN

We fetch the paths with the files for the available glaciers on disk

rgi_paths = get_rgi_paths()

Dict{String, String} with 83 entries:
  "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.02787" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02787"
  "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646"
  "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270"
  "RGI60-11.02773" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02773"
  "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-11.00597" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00597"
  "RGI60-11.02801" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02801"
  "RGI60-11.01876" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01876"
  "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207"
  "RGI60-11.02739" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02739"
  "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"
  "RGI60-07.00065" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00065"
  "RGI60-11.02346" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02346"
  "RGI60-01.00570" => "per_glacier/RGI60-01/RGI60-01.00/RGI60-01.00570"
  ⋮                => ⋮

Define which glacier RGI IDs we want to work with

rgi_ids = ["RGI60-11.03638"]

1-element Vector{String}:
 "RGI60-11.03638"

Define the time step for the simulation output, in this case, a month.

δt = 1/12

0.08333333333333333

We now define the parameters used for the simulation

params = Parameters(
    simulation = SimulationParameters(
        use_MB = false,
        tspan = (2010.0, 2015.0),
        test_mode = false,
        multiprocessing = false, # We are processing only one glacier
        rgi_paths = rgi_paths,
        gridScalingFactor = 4), # Downscale the glacier grid to speed-up this example for the GitHub servers
    hyper = Hyperparameters(
        batch_size = length(rgi_ids), # Set batch size equals size of the dataset
        epochs = [2, 2], # [35,30]
        optimizer = [
            ODINN.Adam(0.02),
            ODINN.LBFGS(
                linesearch = ODINN.LineSearches.BackTracking(iterations = 5)
            )
        ]),
    physical = PhysicalParameters(
        minA = 8e-21,
        maxA = 8e-17),
    UDE = UDEparameters(
        optim_autoAD = ODINN.NoAD(),
        empirical_loss_function = LossH() # Loss function based on ice thickness
    ),
    solver = Huginn.SolverParameters(step = δt) # Save simulation every one month
)

Parameters
  physical    ρ = 900.0 · A ∈ [8.0e-21, 8.0e-17] · C ∈ [8.5e-20, 8.0e-17]
  simulation  tspan = (2010.0, 2015.0) · ✓ iceflow ✗ MB ✓ velocities
  solver      RDPK3Sp35 · reltol = 1.0e-12 · maxiters = 100000
  hyper       epochs = [2, 2] · batch_size = 1 · optimizer = [Adam, LBFGS]
  UDE         target = :A · "AD+AD" · ContinuousAdjoint · loss = LossH

Step 2: Generate synthetic ground truth data with a forward simulation

We define a synthetic law to generate the synthetic dataset. For this, we use some tabular data from Cuffey and Paterson (2010) [2] in a law relating ice temperature with the coefficient A. This law is already available in ODINN.jl, which we specify to be in a gridded format which means that it varies spatially (i.e. non-scalar).

A_law = CuffeyPaterson(scalar = false)

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

**** Model ****

------- Detailed info about the Shallow Ice Approximation (SIA) -------
The SIA is a 2D low order model that approximates the glacier dynamics.
In the ODINN ecosystem this can be parameterized through different variables, and with the current status of the model:
 - H̄ is the ice thickness on the dual grid
 - ∇S is the ice surface slope, given by ∇S=∇(B+H) with B the bedrock elevation
 - D is the diffusivity
 - U is the flux
 - A is the rheology of the ice in Glen's flow law
 - n is the exponent in Glen's flow law
 - ρ and g are resp. the ice density and the gravity constant
 - C is the sliding coefficient in the Weertman sliding law
 - p is the stress exponent in the Weertman sliding law
 - q is the effective-pressure exponent in the Weertman sliding law
> These details will be shown only once. To show them in future prints set model.show_status .= 1;

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: (:T,) -> Matrix{Float64}   (↧@start  )
      C: ConstantLaw -> Array{Float64, 0}
      n: ConstantLaw -> Array{Float64, 0}
      p: ConstantLaw -> Array{Float64, 0}
      q: ConstantLaw -> Array{Float64, 0}
  where
      T => averaged_gridded_long_term_temperature

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

We initialize the glaciers with all the necessary data:

glaciers = initialize_glaciers(rgi_ids, params)

1-element Vector{Glacier2D} distributed over regions 11 (x1)
RGI60-11.03638

Time snapshots where to store data for the inversion:

tstops = collect(2010:δt:2015)

61-element Vector{Float64}:
 2010.0
 2010.0833333333333
 2010.1666666666667
 2010.25
 2010.3333333333333
 2010.4166666666667
 2010.5
 2010.5833333333333
 2010.6666666666667
 2010.75
    ⋮
 2014.3333333333333
 2014.4166666666667
 2014.5
 2014.5833333333333
 2014.6666666666667
 2014.75
 2014.8333333333333
 2014.9166666666667
 2015.0

We generate the synthetic dataset using the forward simulation. This will generate a dataset with the ice thickness and surface velocities for each glacier at each time step. The dataset will be used to make the inversion hereafter.

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

glaciers = prediction.glaciers

1-element Vector{AbstractGlacier} distributed over regions 11 (x1)
RGI60-11.03638

Now we compute the spatially varying A to have a ground truth for the comparison at the end of this tutorial.

A_ground_truth = zeros(size(prediction.glaciers[1].H₀))
A_ground_truth[1:(end - 1), 1:(end - 1)] .= eval_law(
    prediction.model.iceflow.A, prediction, 1,
    (; T = get_input(iAvgGriddedTemp(), prediction, 1, tstops[1])), nothing)
A_ground_truth[prediction.glaciers[1].H₀ .== 0] .= NaN;

Step 3: Model specification to perform a classical inversion

After this forward simulation, we restart the iceflow model to be ready for the inversions

trainable_model = GriddedInv(params, glaciers, :A)
A_law = LawA(params; scalar = false)
model = Model(
    iceflow = SIA2Dmodel(params; A = A_law),
    mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0),
    regressors = (; A = trainable_model)
)

**** Model ****

------- Detailed info about the Shallow Ice Approximation (SIA) -------
The SIA is a 2D low order model that approximates the glacier dynamics.
In the ODINN ecosystem this can be parameterized through different variables, and with the current status of the model:
 - H̄ is the ice thickness on the dual grid
 - ∇S is the ice surface slope, given by ∇S=∇(B+H) with B the bedrock elevation
 - D is the diffusivity
 - U is the flux
 - A is the rheology of the ice in Glen's flow law
 - n is the exponent in Glen's flow law
 - ρ and g are resp. the ice density and the gravity constant
 - C is the sliding coefficient in the Weertman sliding law
 - p is the stress exponent in the Weertman sliding law
 - q is the effective-pressure exponent in the Weertman sliding law
> These details will be shown only once. To show them in future prints set model.show_status .= 1;

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: () -> Matrix{Float64}   (↧@start  custom VJP  ✅ precomputed)
      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
Learnable components
  A: --- Param to invert ---
    Matrix per glacier
    θ: ComponentVector of length 1088

***************

Step 4: Perform the inversion by optimizing the model

We specify the type of simulation we want to perform

inversion = Inversion(model, glaciers, params)

Inversion
  glaciers    1 glacier
  model       iceflow = SIA2Dmodel · mass_balance = TImodel1 · learnable =
      A: --- Param to invert ---
        Matrix per glacier
        θ: ComponentVector of length 1088
  parameters  
      physical    ρ = 900.0 · A ∈ [8.0e-21, 8.0e-17] · C ∈ [8.5e-20, 8.0e-17]
      simulation  tspan = (2010.0, 2015.0) · ✓ iceflow ✗ MB ✓ velocities
      solver      RDPK3Sp35 · reltol = 1.0e-12 · maxiters = 100000
      hyper       epochs = [2, 2] · batch_size = 1 · optimizer = [Adam, LBFGS]
      UDE         target = :A · "AD+AD" · ContinuousAdjoint · loss = LossH
  cache       (nothing)
  results     ✗  not yet run

And finally, we just run the simulation

run!(inversion)

[ Info: Optimizing with ADAM
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [    1 /     4]     Loss:2.29993e+00
Iteration: [    2 /     4]     Loss:2.10338e+00     Improvement: -8.55 %
Iteration: [    3 /     4]     Loss:2.10338e+00     Improvement: 0.00 %
[ Info: Optimizing with LBFGS
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [    4 /     4]     Loss:2.10338e+00     Improvement: 0.00 %
Iteration: [    5 /     4]     Loss:1.07828e+00     Improvement: -48.74 %
Iteration: [    6 /     4]     Loss:4.56952e-01     Improvement: -57.62 %

Now that the model has been optimized, we retrieve the inverted parameters. These parameters do not correspond directly to the values of A. What it defines instead is a parameterization of A to ensure positiveness through a tanh function.

θ = inversion.results.stats.θ

ComponentVector{Float64}(A = (1 = [-0.00010001000133356149 -0.00010001000133356149 … -0.00010001000133356149 -0.00010001000133356149; -0.00010001000133356149 -0.00010001000133356149 … -0.00010001000133356149 -0.00010001000133356149; … ; -0.00010001000133356149 -0.00010001000133356149 … -0.00010001000133356149 -0.00010001000133356149; -0.00010001000133356149 -0.00010001000133356149 … -0.00010001000133356149 -0.00010001000133356149]))

We map the parameters to the values of A by evaluating the law.

A = zeros(size(inversion.glaciers[1].H₀))
inn1(A) .= eval_law(inversion.model.iceflow.A, inversion, 1, (;), θ)
A[inversion.glaciers[1].H₀ .== 0] .= NaN;

Step 5: Compare the inverted parameter to the synthetic ground truth

Finally we visualize the inverted A.

plot_gridded_data(A, inversion.results.simulation[1]; colormap = :YlGnBu, logPlot = true)

We can compare it to the ground truth A values:

plot_gridded_data(A_ground_truth, inversion.results.simulation[1]; colormap = :YlGnBu, logPlot = true)

Unsurprisingly the inverted A is noisy in comparison to the ground truth. This is because the inversion requires regularization. The fact that we are using only the ice thickness might also not help. Adding a second target observation in the loss function, i.e. the ice surface velocities, would help. For more information on how to define regularizations, see the Optimization section.

This page was generated using Literate.jl.