# Make Predictive Glacier Forward Model Compatible with Dual Numbers

**URL:** <https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977>\
**Category:** Modelling & Simulations\
**Created:** [January 18, 2024, 5:47pm UTC](https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977 "2024-01-18T17:47:08Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![vivekag7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivekag7/32/202869_2.png) [@vivekag7](https://discourse.julialang.org/u/vivekag7)\
**Post date:** [January 18, 2024, 5:47pm UTC](https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977/1 "2024-01-18T17:47:08Z")

</div>

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:

```julia
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:

- [initialize\_iceflow\_model! and model.iceflow = SIA2Dmodel](https://github.com/vivekag7/Huginn.jl/blob/main/src/models/iceflow/SIA2D/SIA2D.jl)
- [SIA2D predictive forward model](https://github.com/vivekag7/Huginn.jl/blob/main/src/models/iceflow/SIA2D/SIA2D_utils.jl)

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

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 18, 2024, 7:46pm UTC](https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977/2 "2024-01-18T19:46:47Z")

</div>

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 `AbstractFloat`s, but rather `Real`s which is the supertype of `AbstractFloat`:

```julia-repl
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

```

---

<div class="post-metadata">

**Author:** ![vivekag7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivekag7/32/202869_2.png) [@vivekag7](https://discourse.julialang.org/u/vivekag7)\
**Post date:** [January 18, 2024, 10:03pm UTC](https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977/3 "2024-01-18T22:03:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 18, 2024, 10:49pm UTC](https://discourse.julialang.org/t/make-predictive-glacier-forward-model-compatible-with-dual-numbers/108977/4 "2024-01-18T22:49:55Z")

</div>

> [@vivekag7](#):
>
> 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.

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?

> [@vivekag7](#):
>
> 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.

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](https://github.com/SciML/PreallocationTools.jl) can help here (have not used it myself though).

> [@vivekag7](#):
>
> 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

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.