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.

[caojialin000]

and as last, it will always export NaN data, which will cause model crash.

[glwagner]

I don't think that output is from the Kelvin-Helmholtz instability example. I'm not sure what code you're running

[caojialin000]

well, this is not KH instability case exactly. This model is setup by myself to caculate some specific question about symmetric instability. Here, I just want to borrow functions used to caculate growth rate from KH case.

[caojialin000]

`"""
grow_instability!(simulation, energy)

Grow an instability by running simulation.

Estimates the growth rate σ of the instability
using the fractional change in volume-mean kinetic energy,
over the course of the simulation

energy(t₀ + Δτ) / energy(t₀) ≈ exp(2 σ Δτ)

where t₀ is the starting time of the simulation and t₀ + Δτ
the ending time of the simulation. We thus find that the growth rate
is measured by

σ = log(energy(t₀ + Δτ) / energy(t₀)) / (2 Δτ) .
"""
function grow_instability!(simulation, energy)
# Initialize
# σ = []
simulation.model.clock.iteration = 0
t₀ = simulation.model.clock.time = 0
compute!(energy)
energy₀ = energy[1, 1, 1]

# Grow
run!(simulation)  #,pickup = "./model_checkpoint_iteration8451.jld2"

# Analyze
compute!(energy)
energy₁ = energy[1, 1, 1]
Δτ = simulation.model.clock.time - t₀

# ½(u² + v²) ~ exp(2 σ Δτ)
σ = growth_rate = log(energy₁ / energy₀) / 2Δτ
# push!(σ, grow_instability!(simulation, energy))
return growth_rate

end
nothing # hide`

this is the key function caculating growth_rate. Unfortunately, I came across some bug. So I want to export velocity field first, and then use the formula to caculate directly. Could you please offer some advice?
Thanks~

[glwagner]

The velocity field components are accessible from the simulation:

u, v, w = simulation.model.velocities

is that what you're looking for? (I'm not sure what you mean by "export").

[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...

[glwagner]

I'm going to convert this to a discussion. The example linked above illustrates how to do linear stability calculations with Oceananigans. Any other tips people have to add are welcome!

[caojialin000]

My field of research is the study of instability mechanisms, and I am a huge fan of Oceananigans.jl. This really helps me a lot, especially KH intability you put on Oceananigans.example. I thank you for that. However, I do have some issues need handling.
(1) when I put the code on oceananigans KH instability on GPU, here comes so many bug. Could you please make a GPU version for KH instability? especially, the functions used to caculate growth-rate and rescale.
(2) it would be wonderful if you guys can put on some ohter type of instability. Well, I think the key is how you analyze the instability parameters like energy and so on.
hope for your answer, please.

[navidcy]

Hi @caojialin000,

1. Thanks for pointing out that the KH example does not work on GPU. It'd be great if you can help us fix it. There are a few tips at https://clima.github.io/OceananigansDocumentation/stable/simulation_tips/ that could be helpful. Usually it's a few minor changes that are needed, e.g. some constant appended with const. We'll very much appreciate a PR if you figure it out!

2. Unfortunately, we don't have any plans at the moment for putting more stability problems together. But we would very much welcome a contribution by you if you'd like! If you have any idea of a stability calculation then it'd be nice if you first put it forward here so we can all have a discussion before you go ahead putting work to it towards a PR.

[caojialin000]

thanks. Oceananigans is a great platform for us all. I will do my best to make some contribution. 

…

---- Replied Message ---- From Navid C. ***@***.***> Date 9/10/2023 16:27 To ***@***.***> Cc jialin ***@***.***> , ***@***.***> Subject Re: [CliMA/Oceananigans.jl] linear stability calculator? (Discussion #3000) Hi @caojialin000, Thanks for pointing out that the KH example does not work on GPU. It'd be great if you can help us fix it. There are a few times at https://clima.github.io/OceananigansDocumentation/stable/simulation_tips/ that could be helpful. Usually it's a few minor changes that are needed, e.g. some constant appended with const. We'll very much appreciate a PR if you figure it out! Unfortunately, we don't have any plans at the moment for putting more stability problems together. But we would very much welcome a contribution by you if you'd like! If you have any idea of a stability calculation then it'd be nice if you first put it forward here so we can all have a discussion before you go ahead putting work to it towards a PR. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[francispoulin]

I will add that the Shallow Water model example looks at a barotropic instability, in case that is of interest to you.

[caojialin000]

thanks~look forward to your new work, it really helps me a lot. 

…

---- Replied Message ---- From Francis J. ***@***.***> Date 9/11/2023 01:16 To ***@***.***> Cc jialin ***@***.***> , ***@***.***> Subject Re: [CliMA/Oceananigans.jl] linear stability calculator? (Discussion #3000) I will add that the Shallow Water model example looks at a barotropic instability, in case that is of interest to you. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>