Zero Order Extinction Reaction

Is there any way to achieve the positive_domain effect with the Catalyst.jl syntax?

In my example (decomposition of CaC2 yielding Calcium gas and solid Carbon):
CaC2 → Ca(g) + 2C(s)
I am looking for a way to “remove” Ca from the reacting system. However, I get negative mass:

  • if I override the mass kinetic rate law (using => ) to set a zero-order extinction
    dif = @reaction Dif, Ca=>0
  • or set a negative reaction rate
    dif = @reaction -Dif, 0 --> Ca

plot_1

using a positive_domain callback, the code works as expected:

plot_2

here is the full example:

using DataFrames, Plots
using Catalyst, DifferentialEquations, IfElse
using StatsBase, Optimization, OptimizationOptimJL, OptimizationOptimisers
using OptimizationBBO, OptimizationNLopt
using LsqFit

function positive_domain()
    condition(u,t,integrator) = minimum(u) < 0
    affect!(integrator) = integrator.u .= integrator.uprev
    return DiscreteCallback(condition,affect!,save_positions = (false,false)) 
end

R = 8.314 #J/mol - Universal Gas Constant
MW = [64., 40., 12.]
solids = [1, 1, 1]
@variables t 
D = Differential(t)
@parameters A E Dif
hr = 20
T_iso = 1850. 
time_step = (T_iso-1000)/hr #minus the starting time
T(t) = IfElse.ifelse(t<time_step, hr*t+1273, T_iso+273) 
k(t) = (A*1E4) * exp(-(E*1E3)/R/(T(t)))

rxn = @reaction	k(t), CaC2 --> Ca + 2C
dif = @reaction -Dif, 0 --> Ca #Ca => 0.
@named rs = ReactionSystem([rxn, dif], t)
#osys = convert(ODESystem, rs)

p = (A => 5.9E2, E => 300, Dif => 2E-2)
tspan = (0., 120.)
u0 = [:CaC2 => 1.0, :Ca => 0.0, :C => 0.0]

oprob = ODEProblem(rs, u0, tspan, p)
osol = solve(oprob, AutoTsit5(Rosenbrock23()), saveat = 1., reltol=1E-12, abstol=1E-12)#, callback=positive_domain())

plot(osol, width=2, title="Reaction System", xlabel="Time (min)", ylabel="Species (molar)")
plot!(twinx(), osol.t, T.(osol.t) .-273, legend=:bottom, label="T [C]")

Yes this is something that will be very numerically sensitive and using some kind of callback is generally a good idea. The isoutofdomain solver option and PositiveDomain callback in the DiffEqCallbacks.jl library were designed for this kind of application. The one that you have there introduces error so I would recommend the two standard approaches over that.

Thanks @ChrisRackauckas ! I opted for the PositiveDomain option and it seems ok.

osol = solve(oprob, Tsit5(), saveat = 1., reltol=1E-12, abstol=1E-12, callback = PositiveDomain([1.0]; save = false, abstol = 1E-9))

I have not done any Benchmark but the Optimization block I am using to fit the parameters appears to work smoother.