More DAE benchmarks

[ChrisRackauckas]

https://archimede.dm.uniba.it/~testset/testsetivpsolvers/?page_id=26#DAE

Edit: URL changed to https://archimede.uniba.it/~testset/testsetivpsolvers/?page_id=26#DAE

• [x] Chemical Akzo Nobel
• [ ] Andrews' squeezing mechanism problem
• [x] Transistor amplifier problem
• [ ] Charge-pump problem
• [ ] Two bit adding unit problem
• [ ] Car axis problem
• [ ] Fekete problem
• [x] Slider crank problem
• [ ] Water tube problem
• [ ] NAND gate problem
• [ ] Wheelset problem

[ChrisRackauckas]

If anyone wants to help, what needs to be done is a benchmark file similar to the ROBER DAE one needs to be created:

https://github.com/SciML/SciMLBenchmarks.jl/blob/master/benchmarks/DAE/ROBERDAE.jmd

To do this, you can go to the site and find the benchmark problem description and translate it to the ModelingToolkit.jl form:

https://archimede.dm.uniba.it/~testset/src/problems/andrews.f

If you don't want to setup the solvers, posting the MTK code for the benchmark would be amazing.

[caseykneale]

I got a start on the first one, little stuck figuring out the DAE problem, can't find the constructor for that, (yet) Incase I run out of time(looks like I might) I'll dump the code I have here to maybe save someone some effort:

long code block...

using ModelingToolkit
using OrdinaryDiffEq
using DiffEqDevTools, ODEInterfaceDiffEq
using LinearAlgebra
using DASSL, DASKR
using Sundials

ModelingToolkit.@parameters begin 
      k₁=18.7 
      k₂=0.58 
      k₃=0.09 
      k₄=0.42 
      kbig=34.4 
      kla=3.3 
      ks=115.83 
      po2=0.9 
      hen=737 
end

@variables begin 
      t 
      y₁(t)  = 1.0 
      y₂(t)  = 1.0 
      y₃(t)  = 1.0
      y₄(t)  = 1.0
      y₅(t)  = 1.0
      y₆(t)  = 0.0
end

D = Differential(t)
eqs = [
      r₁  ~ k₁ * (y₁^4.)*sqrt(y₂)
      r₂  ~ k₂ * y₃ * y₄
      r₃  ~ k₂/kbig * y₁ * y₅
      r₄  ~ k₃*y₁*(y₄^2)
      r₅  ~ k₄*(y₆^2)*sqrt(y₂)
      fin ~ kla*(po2/hen-y₂)
      D(y₁) ~ -2. * r₁ + r₂ - r₃ - r₄
      D(y₂) ~ -0.5 * r₁ - r₄ - 0.5*r₅ + fin
      D(y₃) ~ r₁ - r₂ + r₃
      D(y₄) ~ -r₂ + r₃ - 2. * r₄
      D(y₅) ~ r₂ - r₃ + r₅
      D(y₆) ~ ks * y₁ * y₄ - y₆
]

ModelingToolkit.@named sys = ModelingToolkit.ODESystem(eqs)
sys = ode_order_lowering(sys)
u0 = [ 
      y₁  => 0.444,
      y₂  => 0.00123,
      y₃  => 0,
      y₄  => 0.007,
      y₅  => 0,
      y₆  => ks*y₁*y₄#115.83*0.444*0.007#,
]

tspan = (0.0, 180.0)
simpsys = structural_simplify(sys)
mmprob = ODEProblem(simpsys, u0, tspan)
sol = solve(mmprob, Tsit5(),abstol=1/10^14,reltol=1/10^14)
sol(180.0)
#compare to actual solution
# y1_ref = 0.1150794920661702
# y2_ref = 0.1203831471567715 * 10^-2
# y3_ref = 0.1611562887407974
# y4_ref = 0.3656156421249283 * 10^-3
# y5_ref = 0.1708010885264404 * 10^-1
# y6_ref = 0.4873531310307455 * 10^-2

odaeprob = ODAEProblem(simpsys,u0,tspan)
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);
ode_ref_sol[end]

#broken???
daeprob = DAEProblem(sys,u0,[],tspan)
ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);

using Plots
plot(sol,vars=(y₁));
plot!(sol,vars=(y₂));
plot!(sol,vars=(y₃));
plot!(sol,vars=(y₄));
plot!(sol,vars=(y₅));
plot!(sol,vars=(y₆))

[ChrisRackauckas]

You just need to use the ODE form to predict the initial conditions of the DAE. Here's a full example:

using ModelingToolkit
using OrdinaryDiffEq
using DiffEqDevTools, ODEInterfaceDiffEq
using LinearAlgebra
using DASSL, DASKR
using Sundials

ModelingToolkit.@parameters begin
      k₁=18.7
      k₂=0.58
      k₃=0.09
      k₄=0.42
      kbig=34.4
      kla=3.3
      ks=115.83
      po2=0.9
      hen=737
end

@variables begin
      t
      y₁(t)  = 0.444
      y₂(t)  = 0.00123
      y₃(t)  = 0.0
      y₄(t)  = 0.007
      y₅(t)  = 1.0
      y₆(t)  = 115.83*0.444*0.007 # ks*y₁*y₄
end

D = Differential(t)

r₁  = k₁ * (y₁^4.)*sqrt(y₂)
r₂  = k₂ * y₃ * y₄
r₃  = k₂/kbig * y₁ * y₅
r₄  = k₃*y₁*(y₄^2)
r₅  = k₄*(y₆^2)*sqrt(y₂)
fin = kla*(po2/hen-y₂)

eqs = [

      D(y₁) ~ -2. * r₁ + r₂ - r₃ - r₄
      D(y₂) ~ -0.5 * r₁ - r₄ - 0.5*r₅ + fin
      D(y₃) ~ r₁ - r₂ + r₃
      D(y₄) ~ -r₂ + r₃ - 2. * r₄
      D(y₅) ~ r₂ - r₃ + r₅
      D(y₆) ~ ks * y₁ * y₄ - y₆
]

ModelingToolkit.@named sys = ModelingToolkit.ODESystem(eqs)

tspan = (0.0, 180.0)
simpsys = structural_simplify(sys)
mmprob = ODEProblem(simpsys, u0, tspan)
sol = solve(mmprob, Rodas4(),abstol=1/10^14,reltol=1/10^14)
sol(180.0)
#compare to actual solution
# y1_ref = 0.1150794920661702
# y2_ref = 0.1203831471567715 * 10^-2
# y3_ref = 0.1611562887407974
# y4_ref = 0.3656156421249283 * 10^-3
# y5_ref = 0.1708010885264404 * 10^-1
# y6_ref = 0.4873531310307455 * 10^-2

odaeprob = ODAEProblem(simpsys,u0,tspan)
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);
ode_ref_sol[end]

du = mmprob.f(mmprob.u0,mmprob.p,0.0)
du0 = D.(states(sys)) .=> du
daeprob = DAEProblem(sys,du0,u0,tspan)
ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);

using Plots
plot(sol,vars=(y₁));
plot!(sol,vars=(y₂));
plot!(sol,vars=(y₃));
plot!(sol,vars=(y₄));
plot!(sol,vars=(y₅));
plot!(sol,vars=(y₆))

[AayushSabharwal]

I might have misunderstood something here, but the matrix M is given as:

So shouldn't this:

D(y₆) ~ ks * y₁ * y₄ - y₆

Be this:

0. ~ ks * y₁ * y₄ - y₆

?

[ChrisRackauckas]

yes it should. Can you PR a correction?