Implemented framework for spectral methods and Rayleigh Benard Convection with ultraspherical method #476
Conversation
Convection with ultraspherical method
|
I don't have anything useful to say except I am very happy about having an ultra-spherical method in pySDC. So ultra. Can we invent an ultra-SDC method, please? |
|
Can we also have ultra-spectral deferred corrections?? |
|
Maybe this is a stupid question, but why do you call the good ol' man "Chebychov"? I fail to find a reasonable reference for this way of writing his name. Google wants to refer me to "Chebychev", which to me is the far more common name? |
Of course we can. I have been thinking about this only superficially and I don't think it's actually useful. Basically, it would be a different quadrature matrix and you would work in spectral space in time. But the advantage of the ultraspherial method is that it is sparse and hence more easy to solve directly. So maybe we should do ultraspherical Runge-Kutta instead. |
Wikipedia convinced me that the "o" is more in line with the original Russian pronunciation. Notice also that the German spelling is "Tschebyschow". But I can change it to the more common spelling if you want ;) |
Well, I would prefer Tchebychev 😉 (according to Wikipedia ...) |
|
Since we're using English here, I'd prefer the common "Chebyshev". As Wikipedia says: "Currently, the English transliteration Chebyshev has gained widespread acceptance". This is also about findability and consistency with the literature (not the one from Russia, that is). |
Funny how Wikipedia has become a reference now, while 15 years ago it wasn't always considered as trustworthy. Cannot wait a few years when we will use "As ChatGPT says : ..." 😅 |
|
It turns out that it's actually very easy to cache LU decompositions using both scipy and cupy. This is particularly good on GPUs because once solved systems can be efficiently solved again with no detour to CPU needed. The LU decomposition is computed on CPU also when using GPUs. I will work out some adaptive step size selection controller that tries to minimise LU decompositions. Notice that SDC is ideally suited to this because the number of stages is often lower for a given order than in DIRK. That means we can get high order with fewer expensive LU decompositions. Although SDIRK methods might be most efficient still ;) |
|
|
||
| Note that in implicit Euler, the right hand side will be composed of the initial conditions. We don't want that in lines that don't depend on time. Therefore, we multiply the right hand side by the mass matrix. This means you can only do algebraic conditions that add up to zero. But you can easily overload this function with something more generic if needed. | ||
|
|
||
| We use a tau method to enforce boundary conditions in Chebychov methods. This means we replace a line in the system matrix by the polynomials evaluated at a boundary and put the value we want there in the rhs at the respective position. Since we have to do that in spectral space along only the axis we want to apply the boundary condition to, we transform back to real space after applying the mass matrix, and then transform only along one axis, apply the boundary conditions and transform back. Then we transform along all dimensions again. If you desire speed, you may wish to overload this function with something less generic that avoids a few transformations. |
There was a problem hiding this comment.
I removed the need for the comment as well as the comment by implementing more efficient boundary conditions here.
|
Are we good to go now? |
Well, I would like to have a more closer look, but I first have to finish correcting my 210 exams ... can you let me until the end of the week ? |
There was a problem hiding this comment.
I'm sincerely impressed ! Lot of comments are requests of a clearer documentation, and a few comments deals with code design problematics, but nothing that seems very bad in general.
My advice : maybe read first all the comments before modifying anything, some are related ...
|
Thanks a lot for taking the time to look through all this, @tlunet! I agree with most of your comments, but before I start implementing them, I want to address some big picture things. Personally, I am not a big fan of properties because I don't think the individual helpers need the same interface and I want to keep it as modular as possible by use of Finally, I am not happy with the implementation of the transforms, as are you. The 1D bases have a reasonably clean implementation of non-distributed 1D transforms, but the |
Well, you don't need And if you do need |
Seems to be a TODO that can be put somewhere then (describing in details the current issues with the code), and solved later when we have more feedback |
Properties are the python getters ... and they can be extended by setters for more functionalities : for instance, if you want to change the grid at some point (for, e.g, adaptive meshing), then you can simply do |
|
Sorry for the delay on that, I guess I can be happy with all the changes done, thanks @brownbaerchen ! I'll merge it then ... 😉 |
This is quite a large PR and I do apologise. But I think it's finally time to merge some stuff back to the main repository.
I implemented a new spectral helper, which allows to build 2D rectangular geometries with spectral discretizations via direct products of 1D spectral bases. Supported 1D bases are FFTs, Chebychov and Ultraspherical.
Then I implemented a few problems with this, namely heat equation, Burger's equation and Rayleigh Benard Convection. The goal was to make it simple to implement new problems with this. It's a lot more manual than Dedalus. But with a rough understanding of the discretizations it should not be too hard to implement new problems.
This PR does not include any projects done with any of the new problems. I want to merge before because this is already a lot of changes and my branches keep diverging. Also, while I did make an effort to write clean code, I don't have the time to go back and refactor everything too rigorously. Especially the transforms are kind of a mess because I use FFTs to do DCTs and I have not found a clean way to do it, especially with padding. And honestly, I don't know if I will find a cleaner way.. If you spot any large BS or have some hints, I will try, though ;)
The main goal was to do Rayleigh Benard Convection. I will run some more testing, but I was able to confirm the numerical order of SDC with fixed step size when turbulence develops (see here) and I copied the spacial part from Dedalus rather closely.
Speaking of Dedalus, I cannot compete in terms of runtime. My code is slower by a factor of maybe 10. It's hard to compare exactly because of SDC and RK and so on. I will probably shave off a few transforms, but Dedalus is optimised in a way that I will not be able to compete with anytime soon. They cache matrix decompositions and things like that to make solving more efficient. Initial tests with porting to GPUs also didn't look promising as GPUs don't excel at solving sparse linear systems. I haven't given up yet, but no promises.
On the other hand, with this implementation it's much easier to see what's actually going on. The optimisation and generality of Dedalus makes some details hard to follow. So maybe this is still of interest to @chelseajohn. Let me know what you think.
If you are interested in the details of Chebychov or Ultraspherical methods, I recommend this paper, which is a rather concise summary.