Inviscid inertial modes in the sphere

The evolution equation of the velocity $\mathbf{u}$ is given by the momentum equation

\[\frac{\partial\mathbf{u}}{\partial t} + 2\Omega \mathbf{e}_z\times\mathbf{u} = -\frac{1}{\rho}\nabla p,\]

satisfying $\mathbf{u}\cdot\mathbf{n} = 0$ at $r=1$ (the surface of the sphere).

By projecting this equation onto poloidal and toroidal basis vectors, we eliminate the pressure $p$.

Assuming $\mathbf{u}(\mathbf{r},t) = \mathbf{u}(\mathbf{r}) \exp(\lambda t)$, the momentum equation reduces to an eigen problem

\[\lambda \mathbf{B}\mathbf{x} = \mathbf{A}\mathbf{x},\]

where the eigen vector $\mathbf{x}$ contain the spectral coefficients of the basis elements and

\[B_{ij} = \int \mathbf{u}_i \cdot \mathbf{u}_j\,\mathrm{d}V,\]

and

\[A_{ij} = \int \mathbf{u}_i \cdot \left(2\Omega\mathbf{e}_z\times\mathbf{u}_j\right)\,\mathrm{d}V.\]

For moderate polynomial degrees we can quickly solve this.

Solve using Limace.jl

First, load the packages:

using Limace
using LinearAlgebra, SparseArrays, CairoMakie, GeoMakie
using Limace.Discretization: spectospat
CairoMakie.activate!(; type="png", px_per_unit=2)

Define a polynomial truncation degree N.

N = 8
8

Create inviscid velocity basis. Here, we include all azimuthal wave numbers m = -N:N (we can also consider each m individually for the inviscid inertial modes).

u = Inviscid(N)
Basis{Inviscid, Sphere}(8, -8:8, 0:0, InviscidBC(), Sphere(1.0), Dict{Symbol, Float64}())

Assemble the matrix $\mathbf{A}$ of the Coriolis operator.

A = Limace.coriolis(u)
276×276 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 632 stored entries:
⎡⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⢄⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⡄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠦⣀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢣⡀⠀⠀⠀⠈⠳⠠⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⠄⠀⠀⠀⠀⠈⠢⣀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢤⡀⠀⠀⠀⠀⠙⢤⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠡⢄⡀⠀⠀⠀⠙⠢⡄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⠤⣀⠀⠀⠈⠳⣄⠀⎥
⎢⠀⠀⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠂⠀⠀⠈⠓⎥
⎢⢤⠀⠀⠉⠰⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⢑⠀⠀⠀⠘⠦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠑⢄⠀⠀⠀⠈⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⢦⡀⠀⠀⠀⠁⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠢⡀⠀⠀⠀⠁⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⠈⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⣄⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡀⠀⠀⠣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⎦

The inertial matrix (sometimes called mass matrix) is just the unit operator for the Inviscid basis, as it is orthonormal. Due to this orthonormality, the eigen problem reduces to a standard eigenvalue problem. We call Matrix(A) to convert the sparse matrix A to a dense one. eigen gives the dense eigenvalue spectrum. The eigenvalues λ have zero real part, i.e. no viscous damping (you can type λ by typing \lambda<tab>).

λ, x = eigen(Matrix(A));

The eigenvectors are stored in a matrix x, so that λ[i]*x[:,i] ≈ A*x[:,i].

λ[1]*x[:,1] ≈ A*x[:,1]
true

We can do the same using the high-level interface through a LimaceProblem object:

problem = LimaceProblem([u], [Limace.Inertial(u), Limace.Coriolis(u)]);
Limace.assemble!(problem);
Limace.solve!(problem; method=:dense);
λ2, x2 = problem.sol.values, problem.sol.vectors;
per = sortperm(λ,by=abs);
per2 = sortperm(λ2,by=abs);
(λ[per] ≈ λ2[per2]) && (x[:,per] ≈ x2[:,per2])
true

We can compare some of the numerically calculated eigenvalues λ to the analytical equation given by Zhang et al. (2001) for equatorially symmetric inertial modes.

function λ_analytical(m, n)
	sm = sign(m)
	m = abs(m)
	return -sm*2 / (m + 2) * (√(1 + m * (m + 2) / (n * (2n + 2m + 1))) - 1) * im
end;

all(m->any(isapprox(λ_analytical(m,1)), λ), 1:N-1)
true

Since the Inviscid basis is complete in the Cartesian polynomials, all values are exactly as the analytical value (up to machine precision). For higher values of n, the analytical equation zhang(m, n) is only an approximation and one should compare to the roots of the univariate polynomial that gives the eigenvalues of the inertial modes.

Visualization

We can plot the velocity of the modes, for example on the surface. The azimuthal and latitudinal velocity of the mode corresponding to zhang(3,1) is shown below.

r, nθ, nϕ = 1.0, 100, 200

imode = findmin(x->abs(x-λ_analytical(3,1)), λ)[2];

ur, uθ, uϕ, θ, ϕ = spectospat(x[:,imode], u, r, nθ, nϕ);

Here, we have used Limace.Discretization.spectospat to convert the spectral coefficients of the mode to a spatial representation on a grid with nθ and nϕ points in the latitudinal and azimuthal direction, respectively. The radial coordinate is fixed at r=1.0, i.e. the surface of the sphere.

let
	lons, lats = rad2deg.(ϕ).-180, rad2deg.(θ);

	fig = Figure(backgroundcolor = :transparent)
	kwargs = (; dest="+proj=moll", yticklabelsvisible=false, xticklabelsvisible=false, xticksvisible=false, yticksvisible=false)
	ax1 = GeoAxis(fig[1,1]; kwargs..., title=L"u_\theta")
	ax2 = GeoAxis(fig[1,2]; kwargs..., title=L"u_\phi")
	for ax in (ax1, ax2)
		hidedecorations!(ax)
	end

	_uθ = real.(uθ[1,:,:].+eps())
	umax = maximum(abs,_uθ)
	surface!(ax1,lons, lats,_uθ'; colorrange=(-umax,umax), colormap=:balance,shading=NoShading)

	_uϕ = real.(uϕ[1,:,:].+eps())
	umax = maximum(abs,_uϕ)
	surface!(ax2,lons, lats,_uϕ'; colorrange=(-umax,umax), colormap=:balance, shading=NoShading)
	fig
end
Example block output

This page was generated using Literate.jl.