Skip to main content

​
Overview

Dedalus supports native polar coordinate problems on disks and annuli using specialized bases:
  • DiskBasis: For problems on the full disk (r ∈ [0, R])
  • AnnulusBasis: For problems on an annular region (r ∈ [R_i, R_o])
These geometries are ideal for rotating fluid systems, centrifugal effects, and pipe flows.

​
Disk Libration

File: examples/ivp_disk_libration/libration.py

​
Description

Simulates the librational instability in a disk—a phenomenon that occurs when a cylinder of fluid undergoes oscillatory (librational) forcing. The simulation solves the incompressible Navier-Stokes equations linearized around a background librating flow.

​
Physical Setup

  • Geometry: Unit disk (radius = 1)
  • Background flow: Oscillating azimuthal velocity matching libration
  • Boundary conditions: No-slip at r = 1
  • Non-dimensionalization: Disk radius and librational frequency
  • Parameters: Ekman = 1/(2×20²), Ro = 40

​
Key Features

Disk basis setup: 
coords = d3.PolarCoordinates('phi', 'r')
dist = d3.Distributor(coords, dtype=dtype)
disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dealias=dealias, dtype=dtype)
edge = disk.edge
Background librating flow: 
# Construct oscillating background using Bessel functions
from scipy.special import jv
u0_real = dist.VectorField(coords, bases=disk)
u0_imag = dist.VectorField(coords, bases=disk)
u0_real['g'][0] = Ro * np.real(jv(1, (1-1j)*r/np.sqrt(2*Ekman)) / 
                                jv(1, (1-1j)/np.sqrt(2*Ekman)))
u0_imag['g'][0] = Ro * np.imag(jv(1, (1-1j)*r/np.sqrt(2*Ekman)) / 
                                jv(1, (1-1j)/np.sqrt(2*Ekman)))
t = dist.Field()
u0 = np.cos(t) * u0_real - np.sin(t) * u0_imag
Linearized equations: 
problem = d3.IVP([p, u, tau_u, tau_p], time=t, namespace=locals())
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - u@grad(u0) - u0@grad(u)")
problem.add_equation("u(r=1) = 0")
problem.add_equation("integ(p) = 0")

​
Complete Code

import numpy as np
import dedalus.public as d3
from scipy.special import jv
import logging
logger = logging.getLogger(__name__)

# Parameters
Nphi, Nr = 32, 128
Ekman = 1 / 2 / 20**2
Ro = 40
dealias = 3/2
stop_sim_time = 50
timestepper = d3.SBDF2
timestep = 1e-3
dtype = np.float64

# Bases
coords = d3.PolarCoordinates('phi', 'r')
dist = d3.Distributor(coords, dtype=dtype)
disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dealias=dealias, dtype=dtype)
edge = disk.edge

# Fields
u = dist.VectorField(coords, name='u', bases=disk)
p = dist.Field(name='p', bases=disk)
tau_u = dist.VectorField(coords, name='tau_u', bases=edge)
tau_p = dist.Field(name='tau_p')

# Substitutions
phi, r = dist.local_grids(disk)
nu = Ekman
lift = lambda A: d3.Lift(A, disk, -1)

# Background librating flow
u0_real = dist.VectorField(coords, bases=disk)
u0_imag = dist.VectorField(coords, bases=disk)
u0_real['g'][0] = Ro * np.real(jv(1, (1-1j)*r/np.sqrt(2*Ekman)) / jv(1, (1-1j)/np.sqrt(2*Ekman)))
u0_imag['g'][0] = Ro * np.imag(jv(1, (1-1j)*r/np.sqrt(2*Ekman)) / jv(1, (1-1j)/np.sqrt(2*Ekman)))
t = dist.Field()
u0 = np.cos(t) * u0_real - np.sin(t) * u0_imag

# Problem
problem = d3.IVP([p, u, tau_u, tau_p], time=t, namespace=locals())
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - u@grad(u0) - u0@grad(u)")
problem.add_equation("u(r=1) = 0")
problem.add_equation("integ(p) = 0")

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

# Initial conditions
u.fill_random('g', seed=42, distribution='standard_normal')
u.low_pass_filter(scales=0.25)

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.1, max_writes=20)
snapshots.add_task(-d3.div(d3.skew(u)), scales=(4, 1), name='vorticity')
scalars = solver.evaluator.add_file_handler('scalars', sim_dt=0.01)
scalars.add_task(d3.integ(0.5*u@u), name='KE')

# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=100)
flow.add_property(u@u, name='u2')

# Main loop
try:
    logger.info('Starting main loop')
    while solver.proceed:
        solver.step(timestep)
        if (solver.iteration-1) % 100 == 0:
            max_u = np.sqrt(flow.max('u2'))
            logger.info("Iteration=%i, Time=%e, dt=%e, max(u)=%e" 
                        %(solver.iteration, solver.sim_time, timestep, max_u))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()

​
Running the Example

mpiexec -n 4 python3 libration.py
mpiexec -n 4 python3 plot_disk.py snapshots/*.h5
python3 plot_scalars.py scalars/*.h5

​
Expected Results

  • Runtime: ~20 cpu-minutes
  • Physics: Librational forcing excites instability, leading to exponential growth of vorticity
  • Output: Vorticity snapshots and kinetic energy time series

​
Annulus Centrifugal Convection

File: examples/ivp_annulus_centrifugal_convection/centrifugal_convection.py

​
Description

Simulates 2D centrifugal convection in an annulus—convection driven by radial buoyancy forces in a rotating system. This geometry is relevant to Taylor-Couette flows and astrophysical disks.

​
Physical Setup

  • Geometry: Annulus with inner radius R_i and outer radius R_o
  • Radii ratio: η = 3 (R_o = 3×R_i)
  • Boundary conditions: No-slip walls, fixed temperatures
  • Parameters: Ra = 10⁶, Pr = 1
  • Buoyancy: Radial gravity proportional to radius

​
Key Features

Annulus basis setup: 
# Derived parameters
eta = 3
Ri = 2 / (1 + eta)
Ro = 2 * eta / (1 + eta)

# Bases
coords = d3.PolarCoordinates('phi', 'r')
dist = d3.Distributor(coords, dtype=dtype)
annulus = d3.AnnulusBasis(coords, shape=(Nphi, Nr), radii=(Ri, Ro), 
                          dealias=dealias, dtype=dtype)
edge = annulus.outer_edge
Radial gravity: 
# Radial vector field
rvec = dist.VectorField(coords, bases=annulus.radial_basis)
rvec['g'][1] = r

# Gravity proportional to radius
g = rvec * 2 * (eta - 1) / (eta + 1)
First-order formulation with lifting: 
lift_basis = annulus.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
grad_u = d3.grad(u) + rvec*lift(tau_u1)
grad_b = d3.grad(b) + rvec*lift(tau_b1)

problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + b*g + lift(tau_u2) = - u@grad(u)")
problem.add_equation("b(r=Ri) = 0")
problem.add_equation("u(r=Ri) = 0")
problem.add_equation("b(r=Ro) = 1")
problem.add_equation("u(r=Ro) = 0")
problem.add_equation("integ(p) = 0")

​
Running the Example

mpiexec -n 4 python3 centrifugal_convection.py
mpiexec -n 4 python3 plot_polar.py snapshots/*.h5
python3 plot_scalars.py scalars/*.h5

​
Expected Results

  • Runtime: ~10 cpu-minutes
  • Physics: Convection develops in the radial direction, forming plumes that rise/fall
  • Output: Vorticity and buoyancy fields in polar coordinates

​
Pipe Flow Eigenvalue Problem

File: examples/evp_disk_pipe_flow/pipe_flow.py

​
Description

Solves the linear stability eigenvalue problem for pipe flow (Poiseuille flow in a cylinder). This classic problem demonstrates how to handle axial wavenumbers in the disk geometry.

​
Physical Setup

  • Geometry: Unit disk (pipe cross-section)
  • Background flow: Parabolic axial velocity w₀ = 1 - r²
  • Parameters: Re = 10⁴, k_z = 1 (axial wavenumber), m = 5 (azimuthal mode)
  • Boundary conditions: No-slip at r = 1

​
Key Technique: Axial Derivatives

# Axial derivative as multiplication by ik_z
dz = lambda A: 1j*kz*A

# Eigenvalue problem
problem = d3.EVP([u, w, p, tau_u, tau_w, tau_p], eigenvalue=s, namespace=locals())
problem.add_equation("div(u) + dz(w) = 0")
problem.add_equation("dt(u) + w0*dz(u) + grad(p) - (1/Re)*(lap(u)+dz(dz(u))) + lift(tau_u) = 0")
problem.add_equation("dt(w) + w0*dz(w) + u@grad(w0) + dz(p) - (1/Re)*(lap(w)+dz(dz(w))) + lift(tau_w) = 0")
problem.add_equation("u(r=1) = 0")
problem.add_equation("w(r=1) = 0")

​
Running the Example

python3 pipe_flow.py

​
Expected Results

  • Runtime: A few seconds
  • Output:
    • Slowest decaying eigenvalue printed to console
    • Plot of eigenfunction structure pipe_eigenfunctions.png
  • Physics: All modes are stable (negative real part of eigenvalue)
  • Validation: Results compare with reference values from Vasil et al. (2016)

​
Polar Coordinate Tips

​
Accessing Edges

disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dtype=dtype)
edge = disk.edge  # For boundary conditions

annulus = d3.AnnulusBasis(coords, shape=(Nphi, Nr), radii=(Ri, Ro), dtype=dtype)
inner_edge = annulus.inner_edge
outer_edge = annulus.outer_edge

​
Radial Vector Fields

Often needed for lifting in first-order formulations: 
rvec = dist.VectorField(coords, bases=basis.radial_basis)
rvec['g'][1] = r  # Set radial component

​
Coordinate Conversion

For plotting in Cartesian coordinates: 
x, y = coords.cartesian(phi, r)
plt.pcolormesh(x, y, field['g'])

Build docs developers (and LLMs) love

Get started for freeTalk to us