You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
There is a dimension mis-match between the (modified) state vector and the ipiv pivot vector after a ContinuousCallback in which an element of the state vector is removed via deleteat! which only happens above a certain initial particle population (state vector) size.
Expected behavior
I expect this callback should work without issue for any (reasonable) particle population size. I have a feeling that somewhere before appleaccelerate.jl, A.ipiv is not being updated after the callback, because A.factors shows the right dimensionality for the new state vector, while A.ipiv does not.
Minimal Reproducible Example 👇
using DifferentialEquations
using Random
using Distributions
using Plots
Random.seed!(1010)
# ODE Function
function const_nanoparticle_radius_change_dynamics(du, u, p, t)
rate_constants = p
R = u # radii of particles
for i in eachindex(R)
du[i] = rate_constants[i]
end
return nothing
end
# Callback condition for removing particles when radius is at/close to 0
function condition(u, t, integrator)
# Trigger the event if the order of magnitude of the radius goes below the threshold
minimum(u) < 1e-4 ? 0 : 1
end
# Callback to modify the state vector if above condition is met
function affect!(integrator)
original_size = length(integrator.u)
idxs = findall(r -> r <= 1e-4, integrator.u)
new_size = original_size - length(idxs)
# Remove the identified radii that have effectively gone to zero
deleteat!(integrator.u, idxs)
deleteat_non_user_cache!(integrator, idxs)
resize!(integrator, new_size)
resize_non_user_cache!(integrator, new_size)
nothing
end
# Number of nanoparticles, and the mean size and standard deviation of particle sizes
num_particles = 11
mean_size = 12
std = 0.1*mean_size
# Generate random initial particle size distribution
initial_radii = rand(Normal(mean_size, std), num_particles)
# Generate random rate rate constants for each particle (growing or shrinking)
random_rates = rand(Uniform(-0.5,0.5),num_particles)
# Make input vectors
u0 = initial_radii
p = random_rates
tspan = (0, 70)
# Make ODE problem
prob = ODEProblem(const_nanoparticle_radius_change_dynamics, u0, tspan, p)
# Callback to remove a particle from the simulation if it effectively "disappears" due to shrinking
disappearing_callback = ContinuousCallback(condition, affect!)
# Solve
@time sol = solve(prob, Rosenbrock23(), maxiters=1e6, abstol=1e-5, reltol=1e-5, isoutofdomain=(u,p,t)->any(x->x<0.0, u), callback=disappearing_callback)
# Plot number of particles decrease over time
plot(sol.t, map((x) -> length(x), sol[:]), lw = 3,
ylabel = "Number of Nanoparticles", xlabel = "Time")
Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin22.4.0)
CPU: 12 × Apple M3 Pro
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)
Additional context
This minimal example is for a simulation of a kinetic model for nanoparticle sintering (shrinking or growing over time) with constant rates. The process has particles that "vanish" during the simulation because they effectively merge with other, bigger particles. So the callback is meant to remove particles (the respective particle radius) from the simulation after they reach a minimum threshold. Interestingly, for a num_particles < 11, this error doesn't show up. But for num_particles >= 11, I get this (what I think is a linear algebra method) error.
The text was updated successfully, but these errors were encountered:
An easy workaround is Rosenbrock23(linsolve = LUFactorization()) with using LinearSolve. The reason is because RecursiveFactorization.jl needs to store the ipiv but other LU factorizations like that one do not. But this means what you're really needing is to support resizing in LinearSolve.jl and then the algorithms would need to resize individually, where many algorithms wouldn't need to do anything but some algorithms would need to resize the cache.
Describe the bug 🐞
There is a dimension mis-match between the (modified) state vector and the
ipiv
pivot vector after a ContinuousCallback in which an element of the state vector is removed viadeleteat!
which only happens above a certain initial particle population (state vector) size.Expected behavior
I expect this callback should work without issue for any (reasonable) particle population size. I have a feeling that somewhere before appleaccelerate.jl,
A.ipiv
is not being updated after the callback, becauseA.factors
shows the right dimensionality for the new state vector, whileA.ipiv
does not.Minimal Reproducible Example 👇
Error & Stacktrace⚠️
Environment (please complete the following information):
using Pkg; Pkg.status()
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
versioninfo()
Additional context
This minimal example is for a simulation of a kinetic model for nanoparticle sintering (shrinking or growing over time) with constant rates. The process has particles that "vanish" during the simulation because they effectively merge with other, bigger particles. So the callback is meant to remove particles (the respective particle radius) from the simulation after they reach a minimum threshold. Interestingly, for a num_particles < 11, this error doesn't show up. But for num_particles >= 11, I get this (what I think is a linear algebra method) error.
The text was updated successfully, but these errors were encountered: