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!