Improving DAE performance

Specific Domains Modelling & Simulations
1

I’m looking to solve two implicitly-defined nonlinear ODEs. The system diverges toward the final integration z=H, but in a way that I expect, T\sim1/(H-z) and u\sim\log(H-z), plus constants.

What I’m hoping to do is integrate the system to something close to 1 and subtract off this diverging behaviour from the solution to find the non-singular constants that remain.

The first equation depends only on dT/dz and z, while the second equation defines du/dz as a function of dTdz and z.

This is a MWE:

using DifferentialEquations, Sundials, ComponentArrays

function ddz!(res,dy,y,p,z)
    Q = p.Q
    T_z,U_z = dy
    
    res[1] = Q*abs(Q) + z^2*(p.H-z) * T_z^2 * abs(1 - Q^2/((p.H-z)^2*T_z))^p.α
    res[2] = U_z + (p.H-z)*T_z / Q

end

p = ComponentArray(Q=-0.5, α = 1.8, H = 1.)

dy0 = [3., 3.]
y0 = [1., 1.]
zspan = [0.01, p.H - 0.001]
prob1 = DAEProblem(ddz!, dy0, y0, zspan, p, differential_vars = [true,true])

solve(prob1,IDA(linear_solver=:GMRES),abstol=1e-9,reltol=1e-9)

I’m having two main problems:

  1. Perhaps most problematic, my solution changes when I make (seemingly) innocuous changes to ddz!. For example, the above example gives u(1 - 0.001)=8.37. However, if I multiply the last line of the function by a negative:
    res[2] = - U_z - (p.H-z)*T_z / Q,
    I get u(1 - 0.001)=0.094. The results continue to change if I try to multiply the second equation by Q or manipulate the first equation in any way.
    This seems to stop happening for different linear_solver options, but I don’t know how to pick one.

  2. I think I need quite low tolerances. Increasing them to, say, 1e-4 changes the final u. However, for different parameters in the full problem, I sometimes get convergence errors even for relatively high tolerances.

I was thinking this would all be improved with a manually defined jacobian, but I can’t get it to run.

function ddz_jac(J,du,u,p,gamma,z)
    Q = p.Q
    T_z,U_z = du
    
    J[1,1] = gamma*z^2*(p.H-z) * (2*T_z + Q^2/(p.H-z)^2 * 2 * T_z / (T_z - Q^2/(p.H-z)^2)) * abs(1-Q^2/((p.H-z)^2*T_z))^2
    J[1,2] = 0.
    J[2,1] = gamma * (p.H-z) / Q
    J[2,2] = gamma
    nothing
end

DAE_F = DAEFunction(ddz!,jac = ddz_jac)
prob2 = DAEProblem(DAE_F, dy0, y0, zspan, p, differential_vars = [true,true])

Running solve(prob2,IDA()) gives me an error:
[IDAS ERROR] IDACalcIC Newton/Linesearch algorithm failed to converge.

And running solve(prob2,DFBDF()) gives:
type Nothing has no field α, which I don’t understand in this context. I tried manually setting my exponent in the function and jacobian to a fixed number, and I still get this error.

Any help would be greatly appreciated!

2

Have you tried defining this with ModelingToolkit and seeing what structural_simplify and a mass matrix method would do? Generally that does a lot better than the implicit DAEProblem.

2 Likes
3

I wasn’t aware of structural_simplify, but this greatly helps convergence. It also gives me many more solver options to play with to find the one that works best for my problem.

Thank you!

4

No problem.

And yes the benchmarks very consistently show that the implicit DAE form of IDA is just slower than the transformed version.

https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/TransistorAmplifier/

https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/ChemicalAkzoNobel/

https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/OregoDAE/

https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/ROBERDAE/

There are some clear numerical/mathematical arguments as to why that’s the case as well, but we’re starting to have enough benchmarks to make it clear without requiring any background in how the methods work. So yeah… we’re going to be generally trying to push more things down this path over time, and this is why we haven’t put as much work into DFBDF. IDA is fine, but there’s a fundamental limit to its performance by keeping the f definition implicit.

2 Likes
5

Totally makes sense.

I had even looked at those benchmarks, but just thought I was stuck with IDA.

It’s totally possible I’m alone in my misunderstanding as I lack a rigorous background in DEs. But if you think others might fall victim to the same problem, a sentence in the DiffEqDocs faq might be helpful (not to suggest more work for you. Or perhaps later as you continue to push down this path). It’s always my go-to when I get stuck.

6

If you missed it than others would. I’ll add something to the FAQ.

1 Like