As an alternative to DAE definition or SavingCallbacks, I’m using a FunctionCallingCallback with DifferentialEquations.jl to solve algebraic states for an ODE. The callback would need to occur after the integration timestep but before the integration function. Maybe I’m missing it but it appears that this callback type only evaluates post integration function. Is there a way to call the function pre integration function?
What do you mean by the “integration function”?
In the MWE below I’m thinking euler! is my “integration function” and what I referred to as “integration timestep” is what I would consider “solver stuff.”
If what I’m doing now is:
solver stuff → euler! → euler_cb! → solver stuff
I’m hoping to do:
solver stuff → euler_cb! → euler! → solver stuff
Looking at the solution below, we can see that at index 5, (t = 2), the body now has ω due to torque starting at t = 1, but H is still zeros until the next full timestep, which the euler equation is dependent on. My goal is for H to be calculated in euler_cb! ahead of euler!.
Understanding that I could integrate T directly to get H and preserve ODE form, I’m just using this as an example for how I would do this more generally in other applications.
using LinearAlgebra,DifferentialEquations,ComponentArrays,UnPack
x = ComponentArray(
θ = zeros(3), #angles
ω = zeros(3), #rates
T = zeros(3), #torque
H = zeros(3), #ang momentum
)
p = (
J = 1000*I(3), #inertia tensor
)
# integration function?
function euler!(d,x,p,t)
@unpack ω,H,T = x
@unpack J = p
d.θ = ω
d.ω = inv(J)*(T - ω×H)
nothing
end
# FunctionCallback function
function euler_cb!(u,t,s)
s.u.H = s.p.J * s.u.ω
nothing
end
cb = FunctionCallingCallback(euler_cb!)
# start torquing to have something to look at
function torque_cb!(s)
s.u.T[1] = 10
end
torque_cb = PresetTimeCallback(1,torque_cb!)
# sim
prob = ODEProblem(euler!, x, (0, 5), p)
sol = solve(prob,callback = CallbackSet(cb,torque_cb), adaptive = false, dt = 1.)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 8-element Vector{Float64}:
0.0
1.0
1.0
1.0
2.0
3.0
4.0
5.0
u: 8-element Vector{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(θ = 1:3, ω = 4:6, T = 7:9, H = 10:12)}}}}:
ComponentVector{Float64}(θ = [0.0, 0.0, 0.0], ω = [0.0, 0.0, 0.0], T = [0.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.0, 0.0, 0.0], ω = [0.0, 0.0, 0.0], T = [0.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.0, 0.0, 0.0], ω = [0.0, 0.0, 0.0], T = [0.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.0, 0.0, 0.0], ω = [0.0, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.004999999999999974, 0.0, 0.0], ω = [0.009999999999999998, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.019999999999999934, 0.0, 0.0], ω = [0.019999999999999997, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [9.999999999999998, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.044999999999999915, 0.0, 0.0], ω = [0.029999999999999995, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [19.999999999999996, 0.0, 0.0])
ComponentVector{Float64}(θ = [0.07999999999999988, 0.0, 0.0], ω = [0.039999999999999994, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [29.999999999999996, 0.0, 0.0])
sol.u[5]
ComponentVector{Float64}(θ = [0.004999999999999974, 0.0, 0.0], ω = [0.009999999999999998, 0.0, 0.0], T = [10.0, 0.0, 0.0], H = [0.0, 0.0, 0.0])
What you’re looking for are stage_limiters. These are functions called between stages. They are documented in the solver docs:
More and more solvers are getting support for that over time.