Adaptive exponential Rosenbrock integrators #463
Conversation
It's what is normed to get the error estimate. It should be You need to add a convergence test as well. The convergence test automatically sets |
| w2 = phiv(dt, A, F2, 3; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop) | ||
| u = uprev + dt * (w1[:,2] - 2w1[:,4] + 2w2[:,4]) | ||
| if integrator.opts.adaptive | ||
| utilde = dt * w1[:,2] # embedded method |
There was a problem hiding this comment.
In the language of the paper, this is uhat, right? So then in our setup, the error estimate is the first argument of calculate_residuals which should be u - utilde. However, it's more numerically stable to never construct uhat itself. Instead, define utilde = u - uhat = [-2phi_3 phi_3] = 2dt*(-w1[:,4] + w2[:,4]) as the error estimator.
| @@ -101,4 +101,15 @@ end | |||
| @test sol(1.0) ≈ exp_fun2.analytic(u0,nothing,1.0) | |||
| end | |||
| end | |||
|
|
|||
There was a problem hiding this comment.
Needs a convergence test
| sol_ref = solve(prob, Tsit5(); abstol=abstol, reltol=reltol) | ||
|
|
||
| sol = solve(prob, Exprb32(m=20); adaptive=true, abstol=abstol, reltol=reltol) | ||
| @test isapprox(sol(1.0), sol_ref(1.0); rtol=reltol) |
There was a problem hiding this comment.
I am not sure this test will do so well given that the error estimators are local instead of global. This can be off by an order of magnitude or more even in a sane integration
| @test isapprox(sol(1.0), sol_ref(1.0); rtol=reltol) | ||
| sol = solve(prob_ip, Exprb32(m=20); adaptive=true, abstol=abstol, reltol=reltol) | ||
| @test isapprox(sol(1.0), sol_ref(1.0); rtol=reltol) | ||
| end | ||
| end |
There was a problem hiding this comment.
There should be a regression test. I would do some diagnostic printing to make sure things are right: check the error estimator against the analytical solution and it should be in the same ballpark, make sure it goes to zero sanely as the tolerances decrease, etc. Then when you think you have it right, make it solve the problem at some standard tolerance and put a test that its length cannot grow much. So if you measure it at 23 steps, put it at <26 or something like that. Differences between processors can change it by a step or two due to floating point differences, but not by a whole lot. A big indicator for when the error estimator is off is that the step count will grow, so that's a nice regression test for the future to indicate any time that the error estimate breaks.
|
I added convergence tests for both Exprb32 and Exprb43 but found that they seem to only have first and second order respectively. I have double checked the implementation but can't wrap my head around what's wrong. Exprb43 seems most peculiar, because if there is indeed some mistakes on the coefficients the methods should reduce to first order, not 2. In Hochbrock's review (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.296.2179&rep=rep1&type=pdf) it is stated that the stiff order of an exponential integrator can be different from the non-stiff order, and that many of the previous 4th order methods (e.g. ETDRK4) reduces to second order in a worst case scenario, so I can't help wonder if this is somehow connected. But on the other hand the tests for ETDRK4 clearly indicates 4th order convergence, so why would the Exprb methods be any different? |
Codecov Report
@@ Coverage Diff @@
## master #463 +/- ##
==========================================
- Coverage 95.14% 95.09% -0.06%
==========================================
Files 78 78
Lines 8718 8744 +26
==========================================
+ Hits 8295 8315 +20
- Misses 423 429 +6
Continue to review full report at Codecov.
|
|
Wrap this up to something that passes tests with |
|
I couldn't find a fix for Exprb but after going through the paper in detail I seemed to have found some inconsistency in the formulation. The main paper I'm referencing is the Hochbruck review paper (http://na.math.kit.edu/download/papers/exp_int.pdf). The original paper for adaptive Exprb (https://pdfs.semanticscholar.org/e6cc/1023bba0b4670559f1283c50a8ddcb119373.pdf) contains pretty much the same information. The main issue I'm having is what the Butcher tables of Example 2.23 (Exprb32) and Example 2.24 (Exprb43) actually means. I have assumed they are just like the tables for regular ExpRK methods (with an additional embedded method) in my implementation. This is confirmed by the formulations in (2.46a) and (2.46b). However, the authors proposed an alternative formulation in (2.48a) and (2.48b), which is basically what EPIRK uses (Tokman's work is actually mentioned). I have tried rewriting Exprb using the EPIRK formulation with the same coefficients as the Butcher tables, but this too do not seem to work. In actuality, when looking at (2.48a) and (2.48b) you would assume that Leaving this as a note to self. I shall document this in ExponentialIntegrators along with EPIRK5s3. |
|
Close this as development of exponential integrators go to a separate package (https://github.com/MSeeker1340/ExponentialIntegrators). |
|
Wait, why would the exponential integrators be a separate package? The utilities are a separate thing which are useful on their own, but the integrators are highly tied to the ODEIntegrator |
MSeeker1340
changed the title
|
Readded changes to Exprb when developing ExponentialIntegrators. I shall reimplement the other changes (EPIRK, tests, etc) after merging this and finishing ExponentialUtilities. |
|
@ChrisRackauckas OK to merge now? I'll do a PR for the other changes (with ExponentialUtilities) shortly afterwards. |
Implements the adaptive Rosenbrock integrators exprb32 and exprb43 in #394.
The basic structure is still the same as classical ExpRK methods while the adaptive portion is modified from Tsit5. I'm not 100% sure if this is correct so @ChrisRackauckas can you take a look at the code responsible for adaptation? In particular, I'm not sure what
utildeis supposed to be, although by looking at Tsit5 I think this should be the delta u of the embedded method correct?If this is ok, I will add exprb43 (which is more complicated but follows the same structure) and then merge. I should note that the tests are temporary for the moment like the rest of ExpRK, and we should move to more sophisticated convergence tests in the future.