Oceananigans.jl: linear stability calculator?

[francispoulin]

Inspired by @glwagner's script from this morning with a beautiful movie of Rayeigh-Taylor instability, I thought it would be great to set up code to solve the linear stability problem for the example (and for others). I'm not sure if others are interested so I will work to do this sometime, and this issue is a way for me not to forget about it.

[glwagner]

@navidcy put together an example that uses the "power method" for this purpose for Kelvin-Helmholz instability:

https://clima.github.io/OceananigansDocumentation/stable/generated/kelvin_helmholtz_instability/

@kburns thinks we can generalize this method to find all eigenvalues / eigenfunctions with something called "Arnoldi iteration".

Not sure this is what you had in mind, but its a start perhaps. I'm not sure how hard it would be to design a general abstraction for this purpose. If its not too much code I think it could be appropriate to go in the src. Another option is to create an independent package that uses Oceananigans (with the possibility of using other models too) for this purpose.

[francispoulin]

Yes, that is pretty much what I had in mind but for a different physical problem.

There are a variety of ways one can do this. Below are a few that come to mind.

1. Power method: worked well for Kelvin-Helmholtz instability and should work well for other problems. Two nice features are that it's pretty simple, and it's already coded up for one example so should be easy to adapt. One not so nice feature is that it only tells us the most unstable mode (spatial structure, growth rates and phase speed) but it doesn't tell us about any other unstable modes.

2. Arnoldi Method: A quick google search came across this library. I have not tried it but could try adapting it to the above mentioned KH problem.

3. Standard eigenvals: If we assume periodicity in two directions it boils down to a 1D eigenvalue problem that should be easy to solve using eigenvals. There should be a direct solver that forms a full matrix and gives you all the eigenfunctions but the indirect uses Arnoldi and gives you some.

As for where they appera I don't kown but at some point, after I get a ShallowWaterModel off the ground, I will play around with different solvers and see what looks promising. This won't happen this week I'm pretty sure.

[kburns]

The issue with the power method is that its convergence depends on the _difference_ between the two largest eigenvalues, and this can be really slow. Arnoldi iteration does a lot better, and gives you multiple eigenmodes. I think it should also be easy to add since most implementations (including ArnoldiMethod.jl, it looks like) are matrix-free. Directly calculating the full spectrum might be harder -- Dedalus can do this because it forms sparse matrices for the global linear operator, but I don't think oceananigans does this?

[francispoulin]

Thanks @kburns for the suggestion. I will look at Dedalus to see what they do and learn from them. I don't think we would want to find the full spectrum unless the problem is one-dimensional, but I guess it depends on what we want it to do. Lots of potential here I'm sure.

[kburns]

@glwagner or @ali-ramadhan would know better -- I don't think oceananigans forms a global LHS matrix for timestepping, but it might be feasible to build it just for eigenvalue problems. But on the other hand, the oceananigans forward evaluation is so fast that I think Arnoldi would probably really rip if you can do it all on the gpu!

[glwagner]

There might be fancy ways to build the matrix using existing abstractions, but I'm not sure due to the pressure solver + fractional step stuff. With immersed boundaries I think we are out of luck for getting a matrix too, since that again makes non-trivial modifications to the discretization and time-stepping algorithm. A matrix free method that uses the user API might be easier to maintain as well. Arnoldi sounds like a pretty sweet idea...