A variational principle based ideal MHD stability solver and optimizer

[rahulgaur104]

Finite-n stability is important for tokamaks and non-poloidally omnigenous stellarators. Currently, there are no open source codes to solve this problem, much less perform optimization with it.

After a preliminary attempt discretizing and working with the ideal MHD force operator in #1791 , I realized that it would be easier to discretize the energy principle instead. Therefore, this PR creates a variational principle based eigenvalue solver and optimizer.

After discretization, the generalized eigenvalue problem becomes

$$A x = \lambda B x $$

where B is symmetric positive definite and A is symmetric. With some efficient linear algebra we cholesky decompose B and convert the problem above into a standard eigenvalue problem.

$$\tilde{A} x = \lambda x$$

Addresses #1790

TODOs for the current PR:

• Test different basis for representing radial derivatives. Fourier is fine for toroidal and poloidal derivatives. Finite-difference, or Legendre can be used for radial.
• Decide the API to use for creating differentiation matrices and ensure that the matrices are not recalculated at each iteration. They only have to be built once and the optimizer should not track their gradients. Tests needed to ensure that.
  Since the fundamental matrices themselves are pretty small (< 100 x 100 ~ few MBs), they can be passed via transforms.
• Ensure the mass matrix is symmetric positive definite and the stiffness matrix without the instability term is close to symmetric positive definite (no or very small negative eigenvalues). This can then be corrected using Gershgorin circle theorem by adding a small positive value to the diagonal before adding the instability drive and eigendecomposing.
• Effcient decomposition of the mass matrix. Mass matrix is N x N, so vanilla Cholesky will be O(N^3) which is prohibitive for large N ~ 100 k. However, the structure of the mass matrix allows us to rearrange it so that the cost of doing Cholesky becomes O(N).
• Enforce incompressibility and make sure it works without polluting the eigenspectrum. Tests needed($\mathbf{\gamma \rightarrow \infty}$ should be equivalent to imposing incompressibility). Observation: Eigenfunctions are smooth when incompressibility is not enforced so perhaps we should optimize without imposing incompressibility (and compare with the incompressible modes post-optimization)
• New plotting routine needed to plot the eigenfunction and the physical quantities associated with it. Probably needs to be merged in plot_section.
• Compare and benchmark with another finite-n stability solver, possibly TERPSICHORE.
• Compare with NIMSTELL Landreman-Buller-Drevlak QH interchange example from Patil et al. Optional.
• Test the eigenvalue gradient calculation and test the resolution limits that can be fit on a single A100 GPU.
• Try eigendecomposing with scipy.eigh on a CPU but gradient calculation on a GPU using the Hellman-Feynman formula.
• ITo study the mode familyN (mod NFP) = 0, we can solve the stability problem in a single-field period, reducing the matrix size $$1/NFP^2$$. Add a separate condition to study this mode family. Future work.

Future PRs

• Write a custom, matrix-free eigenvalue solver to allow high-resolution stability calculations. This will save a lot of memory but will take longer. Future PR.
• External mode stability calculation. Future work.

Checkpoints (only for me):

commit 7166c2d
commit 28e41c8

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB