search

SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.84k stars 224 forks source link

Compatibility with ApproxFun #847

Closed mjyshin closed 2 years ago

mjyshin commented 2 years ago

In the documentation, it says that there is full compatibility with ApproxFun types, but I am having trouble timestepping (nonlinear) PDEs forward as ODEs of ApproxFun Funs. Using the periodic example in the now-deprecated DiffEqApproxFun.jl's readme, I would imagine it might be possible to do something like

using ApproxFun, DifferentialEquations

S = Fourier()
u0=Fun(θ->cos(cos(θ-0.1))-cos(cos(0-0.1)),S)

cosx = Fun(x->cos(x),S)
model(u,θ,t) = u'' + (cosx+1)*u'

tspan = (0.0,1.0)

prob = ODEProblem(model,u0,tspan)

which builds the problem type just fine, but then be able to solve using

u = solve(prob,SomeStiffSolver())

which currently doesn't work. Also, I liked the general idea of putting nontrivial boundary conditions in the callback argument (I think this is different from how ode15s is used in MATLAB's Chebfun because there the BCs are baked into the nonlinear operator object), although it might be nice to have something in the same vein as BVProblem like

function bc!(residual,u,θ,t)
    residual[1] = u(0) - l # left Dirichlet condition
    residual[2] = u'(2π) - r # right Neumann condition
end

(no doubt much easier said than done). I am wondering why the DiffEqApproxFun.jl project (and SpectralTimeStepping.jl) has been discontinued? The readme says that there are "better ways to do this"... What are those ways?

ChrisRackauckas commented 2 years ago

Yeah this is a complex one. DiffEqApproxFun was a fun project but ultimately it wasn't the most efficient way to do things. There's a few things going on that made it a bit awkward to connect. Simply putting Fun objects into the differential equation solver works (or at least used to work), and so we had quite a few nice demos showcasing it.

That said, mixing it with the stiff ODE solvers was always precarious because of how the adaptivity works. ApproxFun.jl always tries to expand the function in the basis to floating point eps accuracy. The stiff ODE solvers are trying to solve a linear solve (Krylov etc. if you choose it) to some tolerance, some nonlinear solve to a tolerance, and then checking that the error of a step is within tolerance. This tolerance of course is usually much higher than 1e-16, so ApproxFun ends up doing a ton of work to get something that is more accuracy than necessary. On top of that, this iterative process only ever can increase the number of Fourier modes. So while it was working, we saw that it could easily overflow the memory, even when using L-stable solvers (our hypothesis at first was that L-stability was required in order to achieve stability with this adaptivity. While that may be the case, this numerical issue is moreso simply a result of reasonable adaptivity rules clashing). So in the end, one could use automated code generation to get something to work here, but it wasn't the fastest implementation and it only really worked well with non-stiff ODE solvers, where we really wanted the implicit methods and such.

The better ways to do this are to do things like, generating code for a pseudospectral discretization and only adaptive between steps, not operations or stages. This can be done by pairing a pseudospectral discretization with a callback that resize!s the integrator based on the integrator.EEst estimates and such. And since the time of ApproxFun, we've also built out an entirely new set of automation on PDEs. Specifically, the ModelingToolkit.jl PDESystem interface gives someone a generic symbolic way of defining a PDE that can be consumed by any discretizer, of which there are currently MethodOfLines.jl for finite difference and NeuralPDE.jl for physics-informed neural network PDE solvers. We plan to then extend this implementation to include NeuralOperators.jl for Fourier Neural Operator and DeepONet solvers, a wrapper to Gridap.jl/FEniCS.jl, and, you guessed it, a pseudospectral discretizer (most likely built on ApproxFun.jl). Not only would this make the syntax a lot cleaner (it would allow a user to simply describe the equations symbolically), it would allow for more optimal code, allow for putting this all in an interface where the solution handling is the same between all PDE solvers, and allow for symbolic manipulations to improve the PDE's formulation before the numerical solve occurs. Because of how this is strictly better, we have halted most of the tests with ApproxFun.jl for now, are finishing MethodOfLines.jl, and then will continue that project to the spectral parts.

In the meantime, there are some documentation examples for how to use ApproxFun-generated operators for things like transforms and application of derivatives in a given basis. You can find examples of doing this in the SciMLBenchmarks, such as

https://benchmarks.sciml.ai/html/MOLPDE/burgers_spectral_wpd.html

And that's the kind of workflow that we plan to generate from the symbolic PDE descriptions.