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...

```julia 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?