Make Predictive Glacier Forward Model Compatible with Dual Numbers

Specific Domains Modelling & Simulations
1

Hi Julia Community,

I have an in place predictive model for glaciers, which is well tested and works perfectly fine with any <:AbstractFloat. Now what I would like to do is make this model compatible with ForwardDiff.Dual, with the purpose of making it usable for a classic inversion for some parameters A and n. What is the best way to make sure that this forward model can be used with DualNumbers?

I have been working for days on this and tried several things such as using a customtype const DualFloat = Union{AbstractFloat, ForwardDiff.Dual} or using PreallocationTools.jl to setup caches for both dual numbers and Floats. However I just can’t seem to make sure that the forward model is compatible with DualNumbers.

I am going to provide you with the workflow:

params = Parameters(OGGM = OGGMparameters(working_dir=working_dir,
                                                multiprocessing=false,
                                                workers=1,
                                                ice_thickness_source = "Farinotti19"),
                            simulation = SimulationParameters(use_MB=true,
                                                            use_iceflow = true,
                                                            velocities = true,
                                                            use_glathida_data = true,
                                                            tspan=(2010.0, 2015.0),
                                                            working_dir = working_dir,
                                                            multiprocessing=false,
                                                            workers=1),
                            solver = SolverParameters(reltol=1e-12))

glacier = initialize_glaciers(rgi_ids, params)[1]

model = Model(
            iceflow = SIA2Dmodel(params), 
            mass_balance = TImodel1(params),
    )

    
# Initialize glacier ice flow model
initialize_iceflow_model!(model.iceflow, glacier_idx, glacier, params)
params.solver.tstops = define_callback_steps(params.simulation.tspan, params.solver.step)
stop_condition(u,t,integrator) = Sleipnir.stop_condition_tstops(u,t,integrator, params.solver.tstops) 

#closure
function action!(integrator)
        if params.simulation.use_MB 
            # Compute mass balance
            MB_timestep!(model, glacier, params.solver.step, integrator.t)
            apply_MB_mask!(integrator.u, glacier, model.iceflow)
        end
end
cb_MB = DiscreteCallback(stop_condition, action!)

# Run iceflow PDE for this glacier
du = params.simulation.use_iceflow ? SIA2D! : noSIA2D!
results = simulate_iceflow_PDE!(simulation, model, params, cb_MB; du = du)


function simulate_iceflow_PDE!(
    simulation::SIM, 
    model::Sleipnir.Model, 
    params::Sleipnir.Parameters, 
    cb::DiscreteCallback; 
    du = SIA2D!) where {SIM <: Simulation}

    # Define problem to be solved
    iceflow_prob = ODEProblem{true,SciMLBase.FullSpecialize}(du, model.iceflow.H, params.simulation.tspan, tstops=params.solver.tstops, simulation)
    iceflow_sol = solve(iceflow_prob, 
                        params.solver.solver, 
                        callback=cb, 
                        tstops=params.solver.tstops, 
                        reltol=params.solver.reltol, 
                        save_everystep=params.solver.save_everystep, 
                        progress=params.solver.progress, 
                        progress_steps=params.solver.progress_steps)
    # @show iceflow_sol.destats
    # Compute average ice surface velocities for the simulated period
    model.iceflow.H .= iceflow_sol.u[end]
    map!(x -> ifelse(x>0.0,x,0.0), model.iceflow.H, model.iceflow.H)
    

    glacier_idx = simulation.model.iceflow.glacier_idx
    glacier::Sleipnir.Glacier2D = simulation.glaciers[glacier_idx[]]

    # Surface topography
    model.iceflow.S .= glacier.B .+ model.iceflow.H 

    # Update simulation results
    results = Sleipnir.create_results(simulation, glacier_idx[], iceflow_sol; light=true)  #not so important at the moment

    return results
end

To keep this post as concise as possible I refer you to GitHub, for some important functions:

Please let me know if you have any clues on how to tackle this problem; every bit of help is greatly appreciated!

2

What exactly is the error you are getting, i.e., method missing, conversion error etc?
Just skimming the SIAD code, you often restrict the types using a single type parameters, which forces all appearance to be the same type. This might be too restrictive for ForwardDiff which passes duals only to some arguments.
Further, dual numbers are not AbstractFloats, but rather Reals which is the supertype of AbstractFloat:

julia> using ForwardDiff

julia> typeof(ForwardDiff.Dual(1.2, 0)) |> supertype
Real

julia> subtypes(Real)
5-element Vector{Any}:
 AbstractFloat
 AbstractIrrational
 ForwardDiff.Dual
 Integer
 Rational

# Define some function similarly restricted to your examples
julia> f(x::T, y::T) where {T<:Real} = let z::T = x + y; x * z end
f (generic function with 1 method)

# More liberal version
julia> g(x::Real, y::Real) = let z = x + y; x * z end
g (generic function with 1 method)

julia> ForwardDiff.derivative(x -> f(x, 1.2), 0.4)
ERROR: MethodError: no method matching f(::ForwardDiff.Dual{ForwardDiff.Tag{var"#1#2", Float64}, Float64, 1}, ::Float64)

julia> ForwardDiff.derivative(x -> g(x, 1.2), 0.4)
2.0
3

Instead of real I defined a new constant which is the union of AbstractFloat and Dual Numbers. The reason I did this is because I generally want to avoid the model to be compatible with integers.

Most common errors I am facing are method errors corresponding to wanting to converse DualNumbers to Float64. I managed to fix some of those (for example instead by using S .= B .+H use S = B .+ H which fixes allocation) , however specifically the method error that happen when dealing with reference values are the ones that break up.

The example you provided seems like something i can utilise! However when working with matrices I should rather say for example

g(x::Matrix{<:Real}, y::Matrix{<:Real}) = let z = x + y; x * z end

Right? Otherwise

4

Why would you wanna do that? Integers are promoted to a more general type anyways, as soon as they hit any operation involving non-integer types. Maybe I’m missing something here, but don’t see a reason why passing integers should not be allowed?

Indeed, can be difficult to predict/compute the required type of containers in generic code. As a start, it’s often easier to not use in-place mutations … first get it to work, then worry about allocations and make it fast again. Maybe PreallocationTools.jl can help here (have not used it myself though).

Yes, that’s the idea. You probably also want to drop the type annotations on any of the local variables, e.g., B = glacier.B instead of B::Matrix{F} = glacier.B in your SIA2D! function.

1 Like