Added preconditioner based on Taylor expansion #158
Conversation
|
Please add a small example or a test of this new preconditioner. Looks interesting, but I'd like to see this at work (and so does my CI monster). |
some of them available as preconditioners.
|
😭 seemed like a good idea ... |
|
Oh don't worry! I am quite eager to work on this, but it's gonna be a whole thing! I realised that I was constructing linear multistep methods as preconditioners and these are not necessarily stable. For the Piline equation, I could get really good results (fourth order local error in one sweep with three nodes) but for advection and heat equation, the schemes were unstable. I just decided to close this PR so you don't get an email every time I make a commit :) |
root condition for convergence
|
So apparently there is a theorem that there are no A-stable LMMs with order higher than two, which makes these preconditioners somewhat useless... But they are convergent at least. |
|
How does the second order preconditionner looks like ? I've been messing a bit with high order sweeps during my PhD to increase explicit SDC stability, and found out that getting higher order was not helping much ... but I found an explicit second order sweep that was as expensive as Forward Euler in term of computation cost, and wasn't that bad at the end. But I never got to write it using the QDelta matrix formulation (that was waaaay too obscure for me at that time 😅) |
|
The first order LMM that is derived in this way is backward Euler and the second order LMM that is derived here is the trapezoidal rule, both of which are A-stable, so that is good. Yes, the QDelta is confusing. I just don't think about it too much and just construct time marching schemes and plug them in there how I learned from other preconditioners. But why this works out is unclear to me, since you subtract in on the rhs. and have the full Q. But it looks like it's working.. :D |
|
We don't have to merge this in the master branch by the way. I will maybe look more into the regions of stability (which depend on the nodes that you are using, both kind and number) but this may not be useful enough to put it in the master.. |
|
Actually, I don't know if there is already a functionality in pySDC that automatically computes and plots the regions of stability for given nodes and QDelta ... that would quite useful to compare preconditionners |
|
My plan was to make script that executed dahlquist problems for some patch in the complex plane and shows if the scheme is stable. Do you think it's possible to do this analytically in an automated way? I haven't spend to many thoughts on this. I am guessing with the boundary locus method it must be possible somehow... |
|
Don't know what the boundary locus method is ... The way I imagine it : you discretize the complex plane into lam_{i,j} = x_i + y_j i, then put all those lam_{i,j} values in a diagonal matrix A, then use SDC to solve dU/dt = AU with u_0 = 1, on one time step delta_t = 1, to get a numerical solution u_{i,j} for each lam_{i,j}. Then you look at the absolute values of u_{i,j} and plot it in the complex plane. Unitary isocontours represent the stability contour, and regions with |u_{i,j}| <= 1 the numerically stable regions. Here are some examples of stability contours I got with this approach ... using FE sweep
and trying some explicit multi-steps preconditionners (SO = second order, TO = third order)
PS : and I forgot, but the stability regions also depend on the number of iterations ... ;) |
|
Ok that was what I was about to do, so that's good! The many-tests-simulateneously Dahlquist equation is already implemented in pySDC. I saw that the stability region depends on the number of iterations in my tests, but I thought that was some bogus numerics... Even if nothing comes of this project, I have definitely learned a lot already, and will learn lots more, so it wasn't a waste of time ;) |
|
The plots are here: https://github.com/Parallel-in-Time/pySDC/blob/35535b6b1bcdd4038d77047897f7ea8ff5b2f682/pySDC/playgrounds/Preconditioners/LMM_preconditioner.ipynb Doesn't look too good... |
This is correct but there are A(alpha)-stable LMMs of order more than two. So unless you are targeting stiff oscillatory problems like we did in the FWSW-SDC paper, you might get away with A(alpha) stable methods. |
Do you know how to construct particularly stable LMMs? So far I did not come up with something better than trial and error.. |
Probably to link with what is currently investigated in https://github.com/tlunet/collgen |
This is a preconditioner built in a similar way to the extrapolation based error estimate. I was just playing around with diagonal preconditioners and somehow arrived at this one, which is not at all diagonal, but it's higher order than one (kind of). At least I got a fourth order method with two sweeps in the Pi-line problem, but then fifth order with three sweeps.
Was this done before? Does any sort of analysis exist for this? Is it worth pursuing? I have no idea, I was just messing around a bit.