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 ODINNWe fetch the paths with the files for the available glaciers on disk
rgi_paths = get_rgi_paths()Dict{String, String} with 56 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.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"
"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"
"RGI60-11.01238" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01238"
"RGI60-11.03005" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03005"
"RGI60-08.00087" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00087"
"RGI60-11.00787" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00787"
"RGI60-08.00203" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00203"
"RGI60-07.01193" => "per_glacier/RGI60-07/RGI60-07.01/RGI60-07.01193"
⋮ => ⋮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/120.08333333333333333We 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
)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.0e-21, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(false, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, false, 4, "", 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), 4), Hyperparameters{Float64, Int64}(1, 1, Float64[], Any[Optimisers.Adam(eta=0.02, beta=(0.9, 0.999), epsilon=1.0e-8), Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}(10, LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.BackTracking{Float64, Int64}(0.0001, 0.5, 0.1, 5, 3, Inf, nothing), nothing, Returns{Nothing}(nothing), Optim.Flat(), true)], 0.0, [2, 2], 1), 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), SciMLBase.NoAD(), 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", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :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: 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 ****
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 => 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.0We 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.glaciers1-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 ****
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{Sleipnir.Model{SIA2Dmodel{Float64, Law{MatrixCacheGlacierId, Sleipnir.GenInputsAndApply{@NamedTuple{}, ODINN.var"#252#262"{Float64, Float64}}, Sleipnir.GenInputsAndApply{@NamedTuple{}, typeof(Sleipnir.emptyVJPWithInputs)}, Sleipnir.GenInputsAndApply{@NamedTuple{}, typeof(Sleipnir.emptyVJPWithInputs)}, ODINN.var"#254#264", Int64, Sleipnir.GenInputsAndApply{@NamedTuple{}, ODINN.var"#255#265"}, CustomVJP}, 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}, ODINN.TrainableComponents{GriddedInv{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(θ = ViewAxis(1:1088, Axis(var"1" = ViewAxis(1:1088, ShapedAxis((34, 32))),)),)}}}}, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyIC, SIA2D_A_target, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(A = ViewAxis(1:1088, Axis(var"1" = ViewAxis(1:1088, ShapedAxis((34, 32))),)),)}}}}}, Sleipnir.ModelCache{SIA2DCache{Float64, Int64, MatrixCacheGlacierId, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, Array{Float64, 0}, Array{Float64, 0}, ScalarCacheNoVJP, ScalarCacheNoVJP}, Nothing}, AbstractGlacier, Results{Sleipnir.Results{Float64, Int64}, TrainingStats{Float64, Int64}}}(Sleipnir.Model{SIA2Dmodel{Float64, Law{MatrixCacheGlacierId, Sleipnir.GenInputsAndApply{@NamedTuple{}, ODINN.var"#252#262"{Float64, Float64}}, Sleipnir.GenInputsAndApply{@NamedTuple{}, typeof(Sleipnir.emptyVJPWithInputs)}, Sleipnir.GenInputsAndApply{@NamedTuple{}, typeof(Sleipnir.emptyVJPWithInputs)}, ODINN.var"#254#264", Int64, Sleipnir.GenInputsAndApply{@NamedTuple{}, ODINN.var"#255#265"}, CustomVJP}, 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}, ODINN.TrainableComponents{GriddedInv{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(θ = ViewAxis(1:1088, Axis(var"1" = ViewAxis(1:1088, ShapedAxis((34, 32))),)),)}}}}, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyIC, SIA2D_A_target, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(A = ViewAxis(1:1088, Axis(var"1" = ViewAxis(1:1088, ShapedAxis((34, 32))),)),)}}}}}(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, A: --- Param to invert ---
Matrix per glacier
θ: ComponentVector of length 1088
), nothing, 1-element Vector{AbstractGlacier} distributed over regions 11 (x1)
RGI60-11.03638
, 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.0e-21, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(false, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, false, 4, "", 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), 4), Hyperparameters{Float64, Int64}(1, 1, Float64[], Any[Optimisers.Adam(eta=0.02, beta=(0.9, 0.999), epsilon=1.0e-8), Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}(10, LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.BackTracking{Float64, Int64}(0.0001, 0.5, 0.1, 5, 3, Inf, nothing), nothing, Returns{Nothing}(nothing), Optim.Flat(), true)], 0.0, [2, 2], 1), 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), SciMLBase.NoAD(), 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", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :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()))), Results{Sleipnir.Results{Float64, Int64}, TrainingStats{Float64, Int64}}(Sleipnir.Results{Float64, Int64}[], TrainingStats{Float64, Int64}(nothing, Float64[], 0, nothing, ComponentArrays.ComponentVector[], ComponentArrays.ComponentVector[], nothing, Dates.DateTime("0000-01-01T00:00:00"))))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 BFGS
[ 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.