Growing cell population using callback

I think the insert function had a small bug too, e.g. the second cell would be sometimes set to zero instead of 50% of the parent.

This seems to work:

using OrdinaryDiffEq

function f(du,u,p,t)
    for i in 1:length(u)
        du[i] = p.α*u[i]
    end
end

function condition(u,t,integrator)
    return 1.0 - maximum(u) # has any cell hit 1? 
end

function affect!(integrator)
    
    u = integrator.u
    
    idx = findall(y-> y >= 1.0, u) # find all cells which have hit concentration == 1
    
    resize!(integrator,length(u)+length(idx)) # resize solution to accomodate all new divisions

    for id in reverse(idx) # for each cell hitting 1, divide:
        u[id+2:end] .= @view u[id+1:end-1]
        u[id] *= 0.5
        u[id+1] = u[id]
    end
    
    nothing
end

callback = ContinuousCallback(condition,affect!, rootfind=SciMLBase.RightRootFind)
u0 = [0.2]
tspan = (0.0,10.0)
p = (; α = 0.3)
prob = ODEProblem(f,u0,tspan, p)
sol = solve(prob,Tsit5(),callback=callback);
3 Likes