Hi Julia Community,
I’ve been trying to make a glacier SIA2D forward model compatible with dual numbers using PreallocationTools.jl, however I run into the following warning using the solver RDPK3Sp35:
┌ Warning: dt(2.2737367544323206e-13) <= dtmin(2.2737367544323206e-13) at t=2010.097388613536, and step error estimate = 52082.68508677296. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/8XHkk/src/integrator_interface.jl:599
I am using RDPK3Sp35 because with just float numbers this is benchmarked to be the best one for this problem. However this issue also arises when using other solvers.
Specifically when applying Rosenbrock23() as a solver method I get the following:
ERROR: First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:
1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
Rosenbrock23(autodiff=false)). More details can befound at
https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
2. Improving the compatibility of `f` with ForwardDiff.jl automatic
differentiation (using tools like PreallocationTools.jl). More details
can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
3. Defining analytical Jacobians. More details can be
found at https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction
Note: turning off automatic differentiation tends to have a very minimal
performance impact (for this use case, because it's forward mode for a
square Jacobian. This is different from optimization gradient scenarios).
However, one should be careful as some methods are more sensitive to
accurate gradients than others. Specifically, Rodas methods like `Rodas4`
and `Rodas5P` require accurate Jacobians in order to have good convergence,
while many other methods like BDF (`QNDF`, `FBDF`), SDIRK (`KenCarp4`),
and Rosenbrock-W (`Rosenbrock23`) do not. Thus if using an algorithm which
is sensitive to autodiff and solving at a low tolerance, please change the
algorithm as well.
And the following warning:
┌ Warning: The supplied DiffCache was too small and was enlarged. This incurs allocations
│ on the first call to `get_tmp`. If few calls to `get_tmp` occur and optimal performance is essential,
│ consider changing 'N'/chunk size of this DiffCache to 25.
└ @ PreallocationTools ~/.julia/packages/PreallocationTools/mJSsc/src/PreallocationTools.jl:151
This leads me to believe that my implementation is not correct. I did the following, to my Forward model:
iceflow_caches = (
H̄ = DiffCache(model.iceflow.H̄),
S = DiffCache(model.iceflow.S),
D = DiffCache(model.iceflow.D),
Dx = DiffCache(model.iceflow.Dx),
Dy = DiffCache(model.iceflow.Dy),
∇S = DiffCache(model.iceflow.∇S),
∇Sx = DiffCache(model.iceflow.∇Sx),
∇Sy = DiffCache(model.iceflow.∇Sy),
Fx = DiffCache(model.iceflow.Fx),
Fy = DiffCache(model.iceflow.Fy),
Fxx = DiffCache(model.iceflow.Fxx),
Fyy = DiffCache(model.iceflow.Fyy),
dSdx = DiffCache(model.iceflow.dSdx),
dSdy = DiffCache(model.iceflow.dSdy),
dSdx_edges = DiffCache(model.iceflow.dSdx_edges),
dSdy_edges = DiffCache(model.iceflow.dSdy_edges)
)
p = (simulation,iceflow_caches)
# Define problem to be solved
iceflow_prob = ODEProblem{true,SciMLBase.FullSpecialize}(du, model.iceflow.H, params.simulation.tspan, tstops=params.solver.tstops, p)
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)
function SIA2D!(dH::Matrix{<:Real}, H::Matrix{<:Real}, p, t::F) where {F <: Real}
simulation = p[1]
caches = p[2]
# Retrieve parameters
SIA2D_model::SIA2Dmodel = simulation.model.iceflow
glacier::Sleipnir.Glacier2D = simulation.glaciers[simulation.model.iceflow.glacier_idx[]]
params::Sleipnir.Parameters = simulation.parameters
int_type = simulation.parameters.simulation.int_type
H̄ = get_tmp(caches.H̄,H)
A = SIA2D_model.A
n = SIA2D_model.n
B = glacier.B
S = get_tmp(caches.S,H)
dSdx = get_tmp(caches.dSdx,H)
dSdy = get_tmp(caches.dSdy,H)
D = get_tmp(caches.D,H)
Dx = get_tmp(caches.Dx,H)
Dy = get_tmp(caches.Dy,H)
dSdx_edges = get_tmp(caches.dSdx_edges,H)
dSdy_edges = get_tmp(caches.dSdy_edges,H)
∇S = get_tmp(caches.∇S,H)
∇Sx = get_tmp(caches.∇Sx,H)
∇Sy = get_tmp(caches.∇Sy,H)
Fx = get_tmp(caches.Fx,H)
Fy = get_tmp(caches.Fy,H)
Fxx = get_tmp(caches.Fxx,H)
Fyy = get_tmp(caches.Fyy,H)
Δx::F = glacier.Δx
Δy::F = glacier.Δy
Γ = SIA2D_model.Γ
ρ::F = simulation.parameters.physical.ρ
g::F = simulation.parameters.physical.g
# First, enforce values to be positive
map!(x -> ifelse(x>0.0,x,0.0), H, H)
# Update glacier surface altimetry
S .= B .+ H
# All grid variables computed in a staggered grid
# Compute surface gradients on edges
diff_x!(dSdx, S, Δx)
diff_y!(dSdy, S, Δy)
avg_y!(∇Sx, dSdx)
avg_x!(∇Sy, dSdy)
∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n[] - 1)/2)
# @infiltrate
avg!(H̄, H)
Γ[] = 2.0 * A[] * (ρ * g)^n[] / (n[]+2) # 1 / m^3 s
D .= Γ[] .* H̄.^(n[] + 2) .* ∇S
# Compute flux components
@views diff_x!(dSdx_edges, S[:,2:end - 1], Δx)
@views diff_y!(dSdy_edges, S[2:end - 1,:], Δy)
# Cap surface elevaton differences with the upstream ice thickness to
# imporse boundary condition of the SIA equation
η₀ = params.physical.η₀
dSdx_edges .= @views @. min(dSdx_edges, η₀ * H[2:end, 2:end-1]/Δx)
dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
dSdy_edges .= @views @. min(dSdy_edges, η₀ * H[2:end-1, 2:end]/Δy)
dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)
avg_y!(Dx, D)
avg_x!(Dy, D)
Fx .= .-Dx .* dSdx_edges
Fy .= .-Dy .* dSdy_edges
# Flux divergence
diff_x!(Fxx, Fx, Δx)
diff_y!(Fyy, Fy, Δy)
inn(dH) .= .-(Fxx .+ Fyy)
end
Any idea what I could be doing wrong?
For completeness some of the extra functions shown.
# Helper functions for the staggered grid
#### Non-allocating functions ####
diff_x!(O, I, Δx) = @views @. O = (I[begin+1:end,:] - I[1:end-1,:]) / Δx
diff_y!(O, I, Δy) = @views @. O = (I[:,begin+1:end] - I[:,1:end - 1]) / Δy
avg!(O, I) = @views @. O = (I[1:end-1,1:end-1] + I[2:end,1:end-1] + I[1:end-1,2:end] + I[2:end,2:end]) * 0.25
avg_x!(O, I) = @views @. O = (I[1:end-1,:] + I[2:end,:]) * 0.5
avg_y!(O, I) = @views @. O = (I[:,1:end-1] + I[:,2:end]) * 0.5
#### Allocating functions ####
"""
avg(A)
4-point average of a matrix
"""
@views avg(A) = 0.25 .* ( A[1:end-1,1:end-1] .+ A[2:end,1:end-1] .+ A[1:end-1,2:end] .+ A[2:end,2:end] )
"""
avg_x(A)
2-point average of a matrix's X axis
"""
@views avg_x(A) = 0.5 .* ( A[1:end-1,:] .+ A[2:end,:] )
"""
avg_y(A)
2-point average of a matrix's Y axis
"""
@views avg_y(A) = 0.5 .* ( A[:,1:end-1] .+ A[:,2:end] )
"""
diff_x(A)
2-point differential of a matrix's X axis
"""
@views diff_x(A) = (A[begin + 1:end, :] .- A[1:end - 1, :])
"""
diff_y(A)
2-point differential of a matrix's Y axis
"""
@views diff_y(A) = (A[:, begin + 1:end] .- A[:, 1:end - 1])
"""
inn(A)
Access the inner part of the matrix (-2,-2)
"""
@views inn(A) = A[2:end-1,2:end-1]
"""
inn1(A)
Access the inner part of the matrix (-1,-1)
"""
@views inn1(A) = A[1:end-1,1:end-1]