Lie-Poisson integrators for

quasi-geostrophic equations

Klas Modin

Collaborator: Milo Viviani

Quasi-geostrophic equation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

3D Navier-Stokes equation

Shallow water equation

Quasi-geostrophic equation

Quasi-geostrophic equation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

Jule Charney

3D Navier-Stokes equation

Shallow water equation

Quasi-geostrophic equation

Quasi-geostrophic equation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

3D Navier-Stokes equation

Shallow water equation

Quasi-geostrophic equation

horizontal scale \(\gg\) vertical scale

Quasi-geostrophic equation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

3D Navier-Stokes equation

Shallow water equation

Quasi-geostrophic equation

horizontal scale \(\gg\) vertical scale

Quasi-geostrophic equation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

3D Navier-Stokes equation

Shallow water equation

Quasi-geostrophic equation

Coriolis + pressure \(\gg\) inertial forces

Vorticity formulation

\dot v + \nabla_v v = -\nabla p - 2 \tilde\Omega\times v
\operatorname{div} v = 0

Apply curl to \(v\)

\dot\omega - \{\psi,\omega \} = 0
\Delta\psi = \omega - f
\dot\omega - L_v\omega = 0
\omega: S^2\to \mathbb{R} \qquad \psi: S^2\to \mathbb{R}

level-sets of \(\omega\)

Geometry of QGE

Lie-Poisson system on \(\mathfrak{X}_\mu(S^2)^* \simeq C^\infty_0(S^2) \)

\(G=\mathrm{Diff}_\mu(S^2)\)

\(T_e^*G\simeq\mathfrak g^*\)

Casimir functions:

\displaystyle\mathcal C_f(\omega) = \int_{S^2}f(\omega)\mu

Finite-dim (weak) co-adjoint orbits:

\displaystyle\omega = \sum_{k=1}^N \Gamma_k \delta_{q^k}

Statistical mechanics theories for 2D Euler

Idea by Onsager (1949):

  • Approximate \(\omega\) by PV for large \(N\)
  • Find invariant measure and presume ergodicity

Miller (1990) and Robert & Sommeria (1991): (MRS)

  • Minimize microcanonical entropy under energy and Casimir constraints

Is the MRS prediction correct?

No!

2D Euler equations are not ergodic

...but perhaps MRS is "generically" correct

Flow ergodic except at "KAM islands"

Poincaré section of finite dimensional Hamiltonian system

Ultimate problem

for geometric integration

To test MRS we need to:

  • Run long simulations
  • Preserve the Casimirs
    (energy + enstrophy alone not enough)
  • Preserve the Lie-Poisson structure

(criterion in MRS)

On \(\mathbb{T}^2\) such discretization exists (sine-bracket)

[Zeitlin 1991, McLachlan 1993]

based on quantization theory by Hoppe (1989)

Numerical simulations support MRS on \(\mathbb{T}^2\)

[Abramov & Majda 2003]

A torus is not a sphere

MRS generally assumed valid also on \(S^2\)

However, non-structure preserving simulations by Dritschel, Qi, & Marston (2015) contradict MRS on \(S^2\)

DQM simulation yield persistent unsteadiness

Our mission: construct trustworthy discretization on \(S^2\)

2D Euler to isospectral flow via Berezin-Toeplitz quantization

C^\infty(M)\ni f \mapsto T^N_f \in \mathfrak{g}_N

Exists if \(M\) compact quantizable Kähler manifold

Idea: approximate Poisson algebra with matrix algebras

\{f,g\} \to [T^N_f,T^N_g] \quad N\to\infty
\displaystyle\dot \omega = \left\{\Delta^{-1}\omega,\omega \right\}
\displaystyle\dot W = [\Delta_N^{-1}W,W]

From 2D Euler

To isospectral

\omega \mapsto W

Lie-Poisson isospectral flows

\dot W = [B(W),W]

Let \(B\colon\mathfrak{g}\to\mathfrak{g}\)

isospectral flow

Analytic function \(f\) yields first integral

C_f(W) := \mathrm{tr}(f(W))

Casimir function

Hamiltonian case

B(W) = \nabla H(W)^\dagger

Hamiltonian function

Note: Non-canonical Poisson structure (Lie-Poisson)

Explicit B-T quantization on \(S^2\)

C^\infty(\mathbb S^2)\ni \omega\mapsto W \in \mathfrak{su}(N)

[Hoppe, 1989]

  • Express \(\omega\) in spherical harmonics expansion \[ \omega = \sum_{l=1}^\infty \sum_{m=-l}^l \omega^{lm}Y_{lm}\]
  • Truncate at \(l_{\it max}=N-1\)
  • For fixed \(m\), linear map between \((\omega^{lm})_{l=1}^{N-m})\) and \(m\):th diagonal of \(W\)
  • Gives \(N\) linear maps

Complicated coefficients, expressed by Wigner 3-j symbols of very high order

~2 weeks to compute coefficients for \(N=1025\)

Discrete \(S^2\) Laplacian on \(\mathfrak{su}(N)\)

  • "Magic" formula [Hoppe & Yau, 1998]
    \[\Delta_N W =\frac{N^2-1}{2}\left([X^N,[X^N,W]]- \frac{1}{2}[X_+^N,[X_-^N,W]]- \frac{1}{2}[X_-^N,[X_+^N,W]] \right)  \]

banded matrices

\displaystyle\dot W = [\Delta_N^{-1}W,W]

Recall

What is \(\Delta_N\) and how compute \(\Delta_N^{-1}W\) ?

(Naive approach requires \(O(N^3)\) operations with large constant)

\(O(N^2)\) operations

  • \(\Delta_N\) admits sparse \(LU\)-factorization with \(O(N^2\) non-zeros

Spatial discretization obtained!

\displaystyle \dot W = [\Delta_N^{-1}W,W]

Note: corresponds to

\(N^2\) spherical harmonics

\(O(N^2)\) operations

\(O(N^3)\) operations

Isospectral flow \(\Rightarrow\) discrete Casimirs

\displaystyle C^N_f(W) = \operatorname{tr}(f(W))

In progress: spatial convergence

Conjecture [M. & Viviani, 2020]
Fixed interval \([0,T]\), constant \(C=C(T,\omega_0)\) s.t. \[ \lVert\omega(t,\cdot)-T_N^{-1}(W^N(t))\rVert_\infty \leq C/N^2 \]

Classical global existence and uniqueness in \(L^\infty\) setting

\(\Longrightarrow\) \(L^\infty\)-norm conserved (it's a Casimir)

Toeplitz-Berezin quantization theory gives

\lVert T_N^{-1}([T_N(\omega_1),T_N(\omega_2)]) - \{\omega_1,\omega_2 \} \rVert_{\infty} = O(1/N^2)

Time discretization

Aim: numerical integrator that is

  • isospectral, \(W_{k}\to W_{k+1}\) an isospectral map
    necessary to preserve Casimirs
     
  • symplectic, \(W_{k}\to W_{k+1}\) a Lie-Poisson map \(\mathfrak{su}(N)^*\to\mathfrak{su}(N)^*\)
    necessary to (nearly) preserve energy and phase space structure

What about symplectic Runge-Kutta methods (SRK)?

  • Not Lie-Poisson preserving!
  • Not isospectral!

Isospectral Symplectic

Runge-Kutta methods

[M. & Viviani 2019]

\dot W = [B(W),W]

Given \(s\)-stage Butcher tableau \((a_{ij},b_i)\) for SRK

\begin{aligned} X_i &= - (W_n + \sum_{j=1}^s a_{ij}X_j)h B(\tilde W_i) \\ Y_i &= hB(\tilde W_i)(W_k + \sum_{j=1}^sa_{ij}Y_j \\ K_{ij} &= h B(\tilde W_i)(\sum_{j'=1}^s(a_{ij'}X_{j'}+a_{jj'}K_{ij'})) \\ \tilde W_i &= W_k + \sum_{j=1}^s a_{ij}(X_j + Y_j + K_{ij} \\ W_{k+1} &= W_k + \sum_{i=1}^s [hB(\tilde W_i),\tilde W_i] \end{aligned}

Theorem: method is isospectral and Lie-Poisson preserving on any reductive Lie algebra

Numerical results for QGE

Zonal jet and vortex structures on Jupiter

Copyright: NASA, Cassini Imaging Team

Simulation of unstable Rossby-Haurwitz wave

Thank you!

M. & Viviani

A Casimir preserving scheme for long-time simulation of spherical ideal hydrodynamics
J. Fluid Mech., 2020

 

M. & Viviani
Lie-Poisson methods for isospectral flows

Found. Comp. Math., 2020