Fourth order EPIRK methods #410
Conversation
| @@ -49,8 +49,8 @@ end | |||
| prob = ODEProblem(f, u0, (0.0, 1.0)) | |||
| prob_ip = ODEProblem{true}(f_ip, u0, (0.0, 1.0)) | |||
|
|
|||
| dt = 0.1; tol=1e-5 | |||
| Algs = [Exp4] | |||
| dt = 0.05; tol=1e-5 | |||
There was a problem hiding this comment.
was there an issue here?
There was a problem hiding this comment.
If you recall our numerical experiments on phiv_timestep!, the provided tolerance is more of a general guideline than a hard limit. This is especially true for small time steps as there's not much room of improvement for the Krylov part and most of the error comes directly from the truncation error of the schemes.
Using the previous dt would fail the test (by a small margin), so I turned it down a bit.
Now that you mentioned it, comparing the results to Tsit5 is not a good test and should probably be rewritten using test_convergence. It's best if we can find a nonlinear system with known solution which we can test on.
There was a problem hiding this comment.
Oh I just meant the dt change.
We don't need analytical solutions for convergence testing, though we should probably switch to using a better analytical solution when SciML/DiffEqProblemLibrary.jl#8 is there. @Vaibhavdixit02 did you ever add the analytical solution? SciML/SciMLSensitivity.jl#14
Codecov Report
@@ Coverage Diff @@
## master #410 +/- ##
==========================================
+ Coverage 83.24% 83.32% +0.08%
==========================================
Files 78 78
Lines 23140 23258 +118
==========================================
+ Hits 19262 19380 +118
Misses 3878 3878
Continue to review full report at Codecov.
|
| @@ -665,6 +666,69 @@ function perform_step!(integrator, cache::Exp4Cache, repeat_step=false) | |||
| # integrator.k is automatically set due to aliasing | |||
| end | |||
|
|
|||
| function perform_step!(integrator, cache::EPIRK4s3AConstantCache, repeat_step=false) | |||
| @unpack t,dt,uprev,f,p = integrator | |||
| A = f.jac(uprev, p, t) | |||
There was a problem hiding this comment.
The EPIRKs should probably use the standard route for grabbing Jacobians, or everything should be updated at the same time.
There was a problem hiding this comment.
It's fine, keep it. SciML/DiffEqBase.jl#113 makes DiffEqFunction-based ODEProblems the default, so this'll be the new style.
|
Is there any specialization that can be done on semilinear problems that would be helpful? |
Definitely. Both the Jacobian and the remainder term can be calculated more conveniently using the knowledge of the linear term. I think we can do this along with the addition of sparse Jacobian support. I should note that in the arxiv paper the author stated that the remainder term can be taken as the nonlinear term for a semilinear problem. However, I believe this is an error on the author's part because EPIRK methods all use Taylor expansion as their basis and this is different from the semilinear partition. The original paper by Tokman, for example, stated clearly that the remainder term should take this particular form. |
|
This should conclude the fourth order EPIRK methods. I shall leave the optimization part as a whole after the fifth order methods are implemented. |
MSeeker1340
changed the title
|
Yes and with the optimization PR we should make the Jacobians use the full setup: https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/derivative_utils.jl#L22-L33 |
The main reference I'm using for the EPIRK integrators is this preprint by Rainwater and Tokman. The two fourth order EPIRK integrators are EPIRK4s3A and EPIRK4s3B (the s3 refers to them having effectively three stages), which I plan to add to OrdinaryDiffEq.
With Exp4 in place the implementation is not very difficult. I made a small improvement to allow
phiv_timestepto accept sub arrays (from@view) for both input and output, so that different calls tophiv_timestep!can share the same cache.It should be noted that EPIRK4s3A has three "flavors": horizontal, vertical and mixed. I choose to use the mixed scheme as it requires one less call of
phiv_timestepand should have better performance (the paper's numerical results also agree with this).