Dear Lawrence, David and Colin,
Thanks a lot for the helpful information. I changed the subject since I thought it might be helpful to have it match the topic of discussion.
Now that I have an example on how to do this in firedrake I think that I could probably adapt my QG stability code (which would go very well with the nonlinear QG demo that I posted a while back). If there is an option to use a direct eigenvalue method (lapack) then that would work for 1D eigenvalue problems. Since the most unstable mode corresponds to the eigenvalue with the largest imaginary part, and I don't suspect that firedrake allows for imaginary targets, I don't imagine it would converge. That's the problem I have in FEniCS.
Colin: Thanks for the advice. For solving linear systems I know how to decompose the problem into a real and imaginary part and that would be easy enough to do. However, for an eigenvalue problem, I'm not sure how to work around this issue.
To illustrate the point, suppose that the linear system is something like
du/dt = -i*A*u
where A is a constant matrix. Formally, we can write the solution as
u = u_0 exp(-i*A*t)
and the eigenvalues that is going to grow the fastest is the one that has the largest imaginary part.
When we assume the ODE has a normal mode decomposition, exp(-i*lambda*t), we get that the eigenvalue system becomes real,
lambda*u= A*u
but the eigenvalue that corresponds to the largest growth rates is the imaginary part of the system.
We could try a solution of the form exp(sigma*t), where the most unstable mode is real (so we can use a real target) but then the system is imaginary. A bit of a catch 22 as far as I can tell.
There are two options that I think *should* work.
1) Set up the system in firedrake/fenics with real matrices, export the matrices into PETSc format, and then learn how to compute the eigenvalues directly in SLEPc. I imagine there is a steep learning curve but if I build my only SLEPc with complex numbers then that *should* work.
2) Alternatively, I could propagate the identity matrix forward for some time (how long is a bit of an issue) to find out what the propagator/fundamental matrix looks like at some later time. Then, I can compute the eigenvalues of that system, and the most unstable would be real so it is something that I can do in firedrake/fenics. One problem is that for a 2D problem, I need to integrate N^2 vectors so that would slow things down, but at least it should work.
Does anyone have any thoughts as to what might be easiest, or other clever ideas?
One goal that I have is to put together a demo to compute the linear stability of a problem and then second part that actually solves the nonlinear evolution. There are a variety of problems that I have in mind. At the moment I am working on a one-layer QG case but then I'd like to do the 3D problem that is almost identical. Would make for a very nice demo, and also help one of my PhD students who is doing this as part of his thesis right now. Any thoughts would be greatly appreciated.
Cheers, Francis