Computing parameter sensitivities with autodiff

This example demonstrates how to parameterize and differentiate a simulation with two parameters (summer and winter n-factors) using forward-mode automatic simulation.

TODO: add more detail/background

Set up forcings and boundary conditions similarly to other examples:

using CryoGrid
forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
grid = CryoGrid.Presets.DefaultGrid_5cm
initT = initializer(:T, tempprofile)
tile = CryoGrid.Presets.SoilHeatTile(
    :T,
    TemperatureBC(forcings.Tair, NFactor(nf=Param(0.5), nt=Param(0.9))),
    GeothermalHeatFlux(0.053u"W/m^2"),
    soilprofile,
    initT;
    freezecurve=PainterKarra(),
    grid=grid
)
tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
u0, du0 = @time initialcondition!(tile, tspan);

 15.345525 seconds (16.12 M allocations: 1.295 GiB, 99.99% compilation time)

Collect model parameters

p = CryoGrid.parameters(tile)

CryoGridParams{Float64} with 2 parameters
Parameters:
┌───────────┬───────────┬─────────┬────────┬───────┬────────┐
│ component │ fieldname │     val │  layer │   idx │   name │
│  UnionAll │    Symbol │ Float64 │ Symbol │ Int64 │ Symbol │
├───────────┼───────────┼─────────┼────────┼───────┼────────┤
│   NFactor │        nf │     0.5 │    top │     1 │     nf │
│   NFactor │        nt │     0.9 │    top │     2 │     nt │
└───────────┴───────────┴─────────┴────────┴───────┴────────┘

Create the CryoGridProblem.

prob = CryoGridProblem(tile, u0, tspan, p, saveat=24*3600.0);

Solve the forward problem with default parameter settings:

sol = @time solve(prob)
out = CryoGridOutput(sol)

CryoGridOutput with 2921 time steps (2010-10-30T00:00:00 to 2011-10-30T00:00:00) and 1 variable:
    T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (218, 2921)

Define a "loss" function; here we'll just take the mean over the final temperature field.

using Statistics
function loss(p)
    local u0, _ = initialcondition!(tile, tspan, p)
    local prob = CryoGridProblem(tile, u0, tspan, p, saveat=24*3600.0)
    local sol = solve(prob, CGEuler());
    local out = CryoGridOutput(sol)
    return mean(ustrip.(out.T[:,end]))
end

loss (generic function with 1 method)

ForwardDiff provides tools for forward-mode automatic differentiation.

using ForwardDiff
pvec = vec(p)

Compute gradient with forward diff:

grad = @time ForwardDiff.gradient(loss, pvec)

This page was generated using Literate.jl.