Cartesian Eigenvalue Problems

Linear stability analysis and spectral problems in Cartesian geometries

Overview

Eigenvalue problems (EVPs) in Dedalus are used to compute linear stability characteristics and eigenmodes of fluid systems. These examples demonstrate:

• Setting up eigenvalue problems with the d3.EVP class
• Using sparse and dense eigensolvers
• Computing growth rates and critical parameters
• Comparing numerical results with analytical solutions

---

Rayleigh-Bénard Stability Analysis

File: examples/evp_1d_rayleigh_benard/rayleigh_benard_evp.py

Description

Computes the maximum linear growth rates for Rayleigh-Bénard convection over a range of horizontal wavenumbers. This classic stability problem determines the onset of convection and the preferred wavelength of convection rolls.

Problem Setup

The eigenvalue problem for linear perturbations around the conductive state:

# Time derivative becomes -iω
dt = lambda A: -1j*omega*A

# Eigenvalue problem
problem = d3.EVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], 
                  namespace=locals(), eigenvalue=omega)
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) - ez@u = 0")
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = 0")
problem.add_equation("b(z=0) = 0")
problem.add_equation("u(z=0) = 0")
problem.add_equation("b(z=Lz) = 0")
problem.add_equation("u(z=Lz) = 0")
problem.add_equation("integ(p) = 0")

Key Technique: Parametrized Wavenumber

To scan over horizontal wavenumbers, we use a Fourier basis with 2 modes and set the domain size to match the desired fundamental wavenumber:

def max_growth_rate(Rayleigh, Prandtl, kx, Nz, NEV=10, target=0):
    # Build Fourier basis for x with prescribed kx as the fundamental mode
    Nx = 2
    Lx = 2 * np.pi / kx
    
    xbasis = d3.ComplexFourier(coords['x'], size=Nx, bounds=(0, Lx))
    zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz))
    
    # ... build and solve EVP ...
    
    solver.solve_sparse(solver.subproblems[1], NEV, target=target)
    return np.max(solver.eigenvalues.imag)

Complete Example

rayleigh_benard_evp.py

import numpy as np
from mpi4py import MPI
import dedalus.public as d3
import logging
import matplotlib.pyplot as plt
logger = logging.getLogger(__name__)

def max_growth_rate(Rayleigh, Prandtl, kx, Nz, NEV=10, target=0):
    """Compute maximum linear growth rate."""
    Lz = 1
    Nx = 2
    Lx = 2 * np.pi / kx

    # Bases
    coords = d3.CartesianCoordinates('x', 'z')
    dist = d3.Distributor(coords, dtype=np.complex128, comm=MPI.COMM_SELF)
    xbasis = d3.ComplexFourier(coords['x'], size=Nx, bounds=(0, Lx))
    zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz))

    # Fields
    omega = dist.Field(name='omega')
    p = dist.Field(name='p', bases=(xbasis,zbasis))
    b = dist.Field(name='b', bases=(xbasis,zbasis))
    u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
    tau_p = dist.Field(name='tau_p')
    tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
    tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
    tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
    tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)

    # Substitutions
    kappa = (Rayleigh * Prandtl)**(-1/2)
    nu = (Rayleigh / Prandtl)**(-1/2)
    x, z = dist.local_grids(xbasis, zbasis)
    ex, ez = coords.unit_vector_fields(dist)
    lift_basis = zbasis.derivative_basis(1)
    lift = lambda A: d3.Lift(A, lift_basis, -1)
    grad_u = d3.grad(u) + ez*lift(tau_u1)
    grad_b = d3.grad(b) + ez*lift(tau_b1)
    dt = lambda A: -1j*omega*A

    # Problem
    problem = d3.EVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], 
                      namespace=locals(), eigenvalue=omega)
    problem.add_equation("trace(grad_u) + tau_p = 0")
    problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) - ez@u = 0")
    problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = 0")
    problem.add_equation("b(z=0) = 0")
    problem.add_equation("u(z=0) = 0")
    problem.add_equation("b(z=Lz) = 0")
    problem.add_equation("u(z=Lz) = 0")
    problem.add_equation("integ(p) = 0")

    # Solver
    solver = problem.build_solver(entry_cutoff=0)
    solver.solve_sparse(solver.subproblems[1], NEV, target=target)
    return np.max(solver.eigenvalues.imag)

if __name__ == "__main__":
    import time
    comm = MPI.COMM_WORLD

    # Parameters
    Nz = 64
    Rayleigh = 1710
    Prandtl = 1
    kx_global = np.linspace(3.0, 3.25, 50)
    NEV = 10

    # Compute growth rate over local wavenumbers
    kx_local = kx_global[comm.rank::comm.size]
    t1 = time.time()
    growth_local = np.array([max_growth_rate(Rayleigh, Prandtl, kx, Nz, NEV=NEV) 
                             for kx in kx_local])
    t2 = time.time()
    logger.info('Elapsed solve time: %f' %(t2-t1))

    # Reduce growth rates to root process
    growth_global = np.zeros_like(kx_global)
    growth_global[comm.rank::comm.size] = growth_local
    if comm.rank == 0:
        comm.Reduce(MPI.IN_PLACE, growth_global, op=MPI.SUM, root=0)
    else:
        comm.Reduce(growth_global, growth_global, op=MPI.SUM, root=0)

    # Plot growth rates from root process
    if comm.rank == 0:
        plt.figure(figsize=(6,4))
        plt.plot(kx_global, growth_global, '.')
        plt.xlabel(r'$k_x$')
        plt.ylabel(r'$\mathrm{Im}(\omega)$')
        plt.title(r'Rayleigh-Benard growth rates ($\mathrm{Ra} = %.2f, \; \mathrm{Pr} = %.2f$)' 
                  %(Rayleigh, Prandtl))
        plt.tight_layout()
        plt.savefig('growth_rates.pdf')

Running the Example

mpiexec -n 4 python3 rayleigh_benard_evp.py

Expected Results

• Runtime: A few seconds
• Output: Plot of growth rates vs wavenumber
• Physics: Near Ra = 1708, the critical Rayleigh number, growth rates transition from negative (stable) to positive (unstable)

---

Waves on a String

File: examples/evp_1d_waves_on_a_string/waves_on_a_string.py

Description

Computes eigenmodes of waves on a clamped string. This simple problem has analytical solutions, making it perfect for validating the eigenvalue solver.

Mathematical Problem

s*u + ∂²u/∂x² = 0
u(x=0) = 0
u(x=L) = 0

Analytical eigenvalues: s_n = -(nπ/L)² for n = 1, 2, 3, …

Complete Code

waves_on_a_string.py

import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# Parameters
Lx = 1
Nx = 128
dtype = np.complex128

# Bases
xcoord = d3.Coordinate('x')
dist = d3.Distributor(xcoord, dtype=dtype)
xbasis = d3.Legendre(xcoord, size=Nx, bounds=(0, Lx))

# Fields
u = dist.Field(name='u', bases=xbasis)
tau_1 = dist.Field(name='tau_1')
tau_2 = dist.Field(name='tau_2')
s = dist.Field(name='s')

# Substitutions
dx = lambda A: d3.Differentiate(A, xcoord)
lift_basis = xbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
ux = dx(u) + lift(tau_1)  # First-order reduction
uxx = dx(ux) + lift(tau_2)

# Problem
problem = d3.EVP([u, tau_1, tau_2], eigenvalue=s, namespace=locals())
problem.add_equation("s*u + uxx = 0")
problem.add_equation("u(x=0) = 0")
problem.add_equation("u(x=Lx) = 0")

# Solve
solver = problem.build_solver()
solver.solve_dense(solver.subproblems[0])
evals = np.sort(solver.eigenvalues)
n = 1 + np.arange(evals.size)
true_evals = (n * np.pi / Lx)**2
relative_error = np.abs(evals - true_evals) / true_evals

# Plot eigenvalue errors
plt.figure(figsize=(6, 4))
plt.semilogy(n, relative_error, '.')
plt.xlabel("eigenvalue number")
plt.ylabel("relative eigenvalue error")
plt.tight_layout()
plt.savefig("eigenvalue_error.pdf")
plt.savefig("eigenvalue_error.png", dpi=200)

# Plot first few eigenvectors
plt.figure(figsize=(6, 4))
x = dist.local_grid(xbasis)
for n, idx in enumerate(np.argsort(solver.eigenvalues)[:5], start=1):
    solver.set_state(idx, solver.subsystems[0])
    ug = (u['g'] / u['g'][1]).real
    plt.plot(x, ug/np.max(np.abs(ug)), label=f"n={n}")
plt.xlim(0, 1)
plt.legend(loc="lower right")
plt.ylabel(r"mode structure")
plt.xlabel(r"x")
plt.tight_layout()
plt.savefig("eigenvectors.pdf")
plt.savefig("eigenvectors.png", dpi=200)

Running the Example

python3 waves_on_a_string.py

Expected Results

• Runtime: A few seconds (serial only)
• Output:
  • eigenvalue_error.png: Shows spectral accuracy (machine precision for low modes)
  • eigenvectors.png: First five eigenmodes (sine waves)
• Validation: Eigenvalues match analytical solution to machine precision

---

Tips for Eigenvalue Problems

Sparse vs Dense Solvers

# Sparse solver: Good for finding a few eigenvalues near a target
solver.solve_sparse(subproblem, NEV=10, target=0)

# Dense solver: Computes all eigenvalues (only feasible for small problems)
solver.solve_dense(subproblem)

Choosing a Target

• target=0: Find least stable/most unstable modes
• target=1j*ω: Find modes near frequency ω
• Adjust NEV to capture enough modes

Extracting Results

# Get eigenvalues
eigenvalues = solver.eigenvalues

# Get most unstable mode
idx = np.argmax(eigenvalues.imag)
solver.set_state(idx, subsystem)

# Now fields contain the eigenmode
u_eigenmode = u['g'].copy()