Dead species coming back alive #1651
Description
This was noted in the MWE https://github.com/HelgavonLichtenstein/Food-Webs by @HelgavonLichtenstein.
I created a debugging script https://github.com/HelgavonLichtenstein/Food-Webs/pull/2 to see what's going on. The end result is captured in:
(A\b)[101] # -2.6218322568017613e-23
(sparse(A)\b)[101] # -0.0 Basically, the linear solver introduces a numerical error which then revives a dead quantity. It's dead because:
function condition(out, bioS, t, integrator)
# extinction threshold for species
out[1:(num_plant+num_animal)] .= bioS[(num_nutrient+1):(index_vessel-1)] .- 1e-6
end
function affect!(integrator, event_idx)
integrator.p[1][event_idx+num_nutrient] = 0.0
integrator.u[event_idx+num_nutrient] = 0.0
@show event_idx[], integrator.t
endenforces u == 0 and du == 0 for that component, yet it can still revive because:
if repeat_step
linres = dolinsolve(integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
du = cache.fsalfirst, u = u, p = p, t = t, weight = weight, solverdata = (;gamma = dtgamma))
else
linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
du = cache.fsalfirst, u = u, p = p, t = t, weight = weight, solverdata = (;gamma = dtgamma))
end
@inbounds @simd ivdep for i in eachindex(u)
k1[i] = -linres.u[i]
end
integrator.destats.nsolve += 1
@inbounds @simd ivdep for i in eachindex(u)
u[i] = uprev[i] + a21*k1[i]
end
@show u[101], uprev[101], k1[101]in the numerical error of k1 introduced by the factorization. Sparse matrices can impose this as a true zero, though that can lose performance for smaller sizes. The right answer is probably just doing u[101] = 0 in the ODE, at least for performance.