Change to DifferentialEquations

[ChrisRackauckas]

This changes the internal time stepping solvers to use the DifferentialEquations.jl interface. This makes it compatible with all of the solvers on DifferentialEquations, including OrdinaryDiffEq.jl, Sundials.jl, ODEInterface.jl, ODE.jl, etc. The OrdinaryDiffEq.jl requirement is only needed because of the default algorithm setting Tsit5(). The formulation was also changed to in-place updates. Additionally, the passing of keyword args makes the full interface available.

This should not only make you immune to further deprecation of ODE.jl, but also should make this much more efficient.

@smldis

[codecov-io]

Codecov Report

Merging #2 into master will increase coverage by 0.79%. The diff coverage is 50%.

@@            Coverage Diff             @@
##           master       #2      +/-   ##
==========================================
+ Coverage   76.79%   77.58%   +0.79%     
==========================================
  Files          13       13              
  Lines         655      647       -8     
==========================================
- Hits          503      502       -1     
+ Misses        152      145       -7

Impacted Files	Coverage Δ
src/PDEs.jl	50% <50%> (+4.28%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 538e450...5ad733d. Read the comment docs.

[smldis]

just in case this gets merged, this is the patch for the readme file: https://github.com/smldis/GalerkinSparseGrids.jl/tree/patch-1

[ABAtanasov]

This looks like a great addition, and I hope to be able to add it very soon. The one issue that I've been having is with Tsit5 is as follows: One of the most recent tests that I've added (commented out on lines 320-347 in the branch you have) concerns the actual time evolution of the wave equations. One of the requirements I impose is that the total energy of the solution (L2 norm of f + L2 norm of \dot f) does not increase as the wave equation evolves in time. The hope is, of course, that this stays as close to constant as possible, but since approximation error prevents this ideal scenario, the stability condition requires no new energy entering into the system. From my test, this looks like its violated with Tsit5, where in fact the energy increases very slightly, and this seems to indicate instability. I'm looking into fixing this, and I have a hunch that there's a relatively simple solution. @eschnett perhaps you may have an idea about this.

Thanks for this contribution, and looking forward to merging it soon. ~Alex

[ChrisRackauckas]

One of the most recent tests that I've added (commented out on lines 320-347 in the branch you have) concerns the actual time evolution of the wave equations. One of the requirements I impose is that the total energy of the solution (L2 norm of f + L2 norm of \dot f) does not increase as the wave equation evolves in time. The hope is, of course, that this stays as close to constant as possible, but since approximation error prevents this ideal scenario, the stability condition requires no new energy entering into the system. From my test, this looks like its violated with Tsit5, where in fact the energy increases very slightly, and this seems to indicate instability. I'm looking into fixing this, and I have a hunch that there's a relatively simple solution. @eschnett perhaps you may have an idea about this.

That has nothing to do with instability. That's actually expected. Only symplectic integrators will keep 1st integral properties like energy constant. Any other method, which includes the Dormand-Prince method, will not.

What you're seeing is actually mostly a difference in default tolerances. DiffEq defaults to reltol=1e-3 which is the standard default for most packages you've ever used (R/Python/MATLAB/C/Fortran). ODE.jl has a (pretty absurd) default of reltol=1e-7. But that default means that, at the level you were testing at you didn't see an energy change. You would see no change in the energy if the Tsit5() solution was made more exact by dropping the tolerance. This method is much more efficient, so it should get below your previous test at something like reltol=1e-5 (I'd have to check on this problem). But the user can pass in keyword arguments to handle that as they please.

The other option is to setup the problem as partitioned and for use with a symplectic integrator. We don't have very many symplectic integrators right now, but that is an option and will be very useful as we get more.

Another option is to use a post-step manifold projection.

http://docs.juliadiffeq.org/latest/features/callback_library.html#Manifold-Conservation-and-Projection-1

That's a pretty cheap way to conserve both order and integration properties. This is a way to ensure that any method conserves energy "exactly", whether it's symplectic or not, regardless of the tolerance.

Obviously there's a tradeoff here between computational cost and conservation properties, so it's moreso where you want to put the default. I left it at the common interface defaults so the users could override as they see fit, but you could easily inject your own defaults in there.

[ABAtanasov]

Yes, I see what you mean, though to clarify the energy isn't constant even with the high tolerance of ODE.jl, but rather decreases very slightly (on the order of 1e-10 per crossing time) and monotonically in time. Regardless, I'll be testing Tsit5/7 before merging around this next weekend. It'll be a welcome change, since ODE.jl has giving me a lot of grief in the package tests.

[ChrisRackauckas]

Yes, I see what you mean, though to clarify the energy isn't constant even with the high tolerance of ODE.jl, but rather decreases very slightly (on the order of 1e-10 per crossing time) and monotonically in time.

Yes, that is to be expected since it's using similar methods to DP5()/Tsit5() and DP8(), which are RK methods and are not symplectic.

[eschnett]

I don't know exactly which system is tested here, but: the sparse grid representation is not actually conservative, in the sense that the derivative operator is losing certain modes during evolution. This means that we expect the energy to drop, but only very slightly, and be conserved in the continuum limit. If using a particular integrator has the energy increase, then this means that the time integration error is larger than the spatial discretization error, which is undesirable here, as this is a test case for the spatial discretization.

The statement above (regarding energy non-increase for stability) assumes a perfect time integrator without any temporal discretization error. In practice, we just use a very high accuracy for the time integrator for these tests.

[ChrisRackauckas]

In practice, we just use a very high accuracy for the time integrator for these tests.

You're using a high order integrator, but not necessarily a "high accuracy" integrator. The accuracy is set by the tolerance, and high order just means more efficient at lower tolerances. But the tolerance isn't actually set. Before it was using the ODE.jl default (reltol=1e-7) while now it's using the DiffEq default (reltol=1e-3). If you want high accuracy, you just change the tolerance.

[miguelraz]

Hello all! I liked this project, so I gave it a spin and brought it up to speed in for Julia 1.0 over here, except for the ODE.jl upgrades that @ChrisRackauckas made. I have all tests passing locally on the v1fixes branch, but there is a lot of work to be done that I am glad to take on. I am unsure how to rebase on top of this branch, so I would appreciate any advise on how to merge the changes into this repo.