OrdinaryDiffEq.jl v6 LinearSolve.jl

[ChrisRackauckas]

I know CLIMA is one of the few projects making use of the OrdinaryDiffEq.jl linear solver interface. The latest breaking release's breaking part is that the old linear solver part was ripped out and replaced with a full library https://github.com/SciML/LinearSolve.jl . The nice thing of course is that if you define LinearSolve.jl solvers you can now just TRBDF2(linsolve=PETSCLinsolve()) and it'll use that internally. Let me know if you need some help updating your linear solvers, but also we should now thing about more generally pulling out generically good linear solvers to there, along with the expanded preconditioners.

(Note better preconditioning is then coming to OrdinaryDiffEq.jl. Now that this is done it's not hard and I hope to do it by the end of the month)

[simonbyrne]

Thanks, we will check it out. Have you kept the Wfact interface?

[ChrisRackauckas]

That wasn't touched. What are you using that for? I can't think of a legitimate use since if you just use a DiffEqOperator you will get the W in its lazy form in the linear solver.

[simonbyrne]

We use it for solving the Schur form: basically we have an efficient way of solving (I - γ * J) Y = Y₀ by defining an extra variable (p for pressure in our case), which satisfies a wave equation (I - γ^2 G) p = γ * h(p₀, ...), and then the remainder of the variables can be computed as a linear update (Y = Y₀ + γ K p). So we define a Wfact should be an object which forms the tridiagonal matrix of the operator (I - γ^2 G), and then a special linear solver that defines an ldiv! method which solves for p, then computes Y from p.

I'd have to look to see if we can use DiffEqOperator (perhaps if J is constant, not sure generally).

See https://github.com/CliMA/ClimaCore.jl/issues/41

[ChrisRackauckas]

Yeah if you make your J a linear operator then the linear solver gives you a WOperator which is a lazy representation of I - γ * J to work with. The most common use case of this is that WOperator has mul! defined and this makes Jacobian-free Newton Krylov automatic, but you should be able to use that to instead target the linear solve directly. That would be more standard. Though I don't like code churn so I wouldn't plan on actually removing Wfact for a few years, even though it is undocumented.