Skip to main content

Overview

Dedalus uses spectral basis functions to represent field data in the spatial domain. The framework supports various 1D basis types (Fourier, Chebyshev, Jacobi) and multidimensional compound bases (Disk, Sphere, Shell, Ball).

One-Dimensional Bases

All 1D basis classes are defined in dedalus/core/basis.py.

Fourier Bases

For periodic domains.

RealFourier

Complex Fourier modes with real-valued data in grid space.
import numpy as np
import dedalus.public as d3

coords = d3.CartesianCoordinates('x')
xbasis = d3.RealFourier(coords['x'], size=256, bounds=(0, 2*np.pi), dealias=3/2)
Parameters:
  • coord: Coordinate object
  • size: Number of modes
  • bounds: Domain bounds (default: [0, 2π])
  • dealias: Dealiasing factor for grid size (default: 1)
  • library: Transform library (‘fftw’ or ‘matrix’, default: ‘fftw’)

ComplexFourier

Complex Fourier modes with complex-valued data.
xbasis = d3.ComplexFourier(coords['x'], size=256, bounds=(0, 2*np.pi), dealias=3/2)
Mode ordering:
  • Odd size N: [0, 1, …, k_max, -k_max, …, -1]
  • Even size N: [0, 1, …, k_max, 1-k_max, …, -1]
  • k_max = N // 2

Fourier Factory

Automatically select based on dtype.
# Returns RealFourier for real dtype
xbasis = d3.Fourier(coords['x'], size=256, bounds=(0, 2*np.pi), 
                     dtype=np.float64, dealias=3/2)

# Returns ComplexFourier for complex dtype
xbasis = d3.Fourier(coords['x'], size=256, bounds=(0, 2*np.pi), 
                     dtype=np.complex128, dealias=3/2)

Jacobi Polynomials

General Jacobi polynomial basis P_n^(a,b)(x) on [-1, 1].
zbasis = d3.Jacobi(coords['z'], size=64, bounds=(0, 1), a=0, b=0, 
                    a0=-1/2, b0=-1/2, dealias=3/2)
Parameters:
  • coord: Coordinate object
  • size: Number of polynomials
  • bounds: Physical bounds (mapped from [-1, 1])
  • a, b: Jacobi parameters for the basis
  • a0, b0: Grid parameters (default: same as a, b)
  • dealias: Dealiasing factor
  • library: Transform library (‘matrix’ default, ‘fftw_dct’ for a0=b0=-1/2)
Properties:
  • Native bounds: [-1, 1]
  • Affine change of variables maps to problem bounds
  • Mass = sqrt(∫ P_n^(a,b)² (1-x)^a (1+x)^b dx)

Special Cases of Jacobi

Legendre

Legendre polynomials (a=0, b=0).
zbasis = d3.Legendre(coords['z'], size=64, bounds=(0, 1), dealias=3/2)

Chebyshev

Chebyshev polynomials of the first kind T_n(x) (α=0).
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1), dealias=3/2)
# Equivalent to:
zbasis = d3.Chebyshev(coords['z'], size=64, bounds=(0, 1), dealias=3/2)
Properties:
  • Ultraspherical with α=0
  • Jacobi with a=b=-1/2
  • Efficient DCT transforms available

ChebyshevU

Chebyshev polynomials of the second kind U_n(x) (α=1).
zbasis = d3.ChebyshevU(coords['z'], size=64, bounds=(0, 1), dealias=3/2)

Ultraspherical

Ultraspherical (Gegenbauer) polynomials C_n^(α)(x).
zbasis = d3.Ultraspherical(coords['z'], size=64, bounds=(0, 1), 
                            alpha=1, alpha0=-1/2, dealias=3/2)
Parameters:
  • alpha: Ultraspherical parameter for basis
  • alpha0: Grid parameter (default: same as alpha)
  • Converts to Jacobi with a = b = alpha - 1/2

Multidimensional Bases

Compound bases for problems with special geometries.

DiskBasis

For problems on a disk in polar coordinates.
import dedalus.public as d3

coords = d3.PolarCoordinates('phi', 'r')
dist = d3.Distributor(coords, dtype=np.float64)

basis = d3.DiskBasis(coords, shape=(256, 64), radius=1, dealias=(3/2, 3/2), 
                      dtype=np.float64)
Parameters:
  • coordsys: PolarCoordinates object
  • shape: (n_phi, n_r) - number of modes
  • radius: Disk radius
  • dealias: Dealiasing factors (phi, r)
  • dtype: Data type
  • k: Radial mode offset (default: 0)
  • alpha: Jacobi α parameter (default: 0)
Structure:
  • Azimuthal: Fourier modes in phi
  • Radial: Jacobi polynomials in r
  • Automatically handles regularity at r=0

AnnulusBasis

For problems on an annulus.
basis = d3.AnnulusBasis(coords, shape=(256, 64), radii=(0.5, 1), 
                         dealias=(3/2, 3/2), dtype=np.float64)
Parameters:
  • coordsys: PolarCoordinates object
  • shape: (n_phi, n_r)
  • radii: (inner_radius, outer_radius)
  • Other parameters same as DiskBasis

SphereBasis

For problems on a spherical surface (S2).
coords = d3.S2Coordinates('phi', 'theta')
dist = d3.Distributor(coords, dtype=np.float64)

basis = d3.SphereBasis(coords, shape=(256, 128), radius=1, dealias=(3/2, 3/2), 
                        dtype=np.float64)
Parameters:
  • coordsys: S2Coordinates object
  • shape: (n_phi, n_theta)
  • radius: Sphere radius
  • dealias: Dealiasing factors
  • dtype: Data type
  • k: Spin order offset (default: 0)
  • alpha: Regularization parameter (default: 0)
Structure:
  • Uses spin-weighted spherical harmonics
  • Azimuthal: Fourier in phi
  • Meridional: Jacobi in theta with special regularity

BallBasis

For problems in a solid sphere.
coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=np.float64)

basis = d3.BallBasis(coords, shape=(256, 128, 64), radius=1, 
                      dealias=(3/2, 3/2, 3/2), dtype=np.float64)
Parameters:
  • coordsys: SphericalCoordinates object
  • shape: (n_phi, n_theta, n_r)
  • radius: Ball radius
  • dealias: Dealiasing factors
  • dtype: Data type
  • k: Spin-radial offset (default: 0)
  • alpha: Jacobi α parameter (default: 0)
Structure:
  • Combines S2 structure with radial Jacobi
  • Handles regularity at r=0 and on pole axis

ShellBasis

For problems in a spherical shell.
basis = d3.ShellBasis(coords, shape=(256, 128, 64), radii=(0.5, 1), 
                       dealias=(3/2, 3/2, 3/2), dtype=np.float64)
Parameters:
  • coordsys: SphericalCoordinates object
  • shape: (n_phi, n_theta, n_r)
  • radii: (inner_radius, outer_radius)
  • Other parameters same as BallBasis

Basis Properties

Common Properties

All bases share these properties:
# Dimensionality
print(basis.dim)  # Number of dimensions

# Shape
print(basis.shape)  # Mode counts per dimension

# Bounds
print(basis.bounds)  # Physical domain bounds

# Volume
print(basis.volume)  # Domain measure

# Coordinate system
print(basis.coordsys)  # Associated coordinate system

# Grid shape at scale
print(basis.grid_shape(scales=(3/2, 3/2)))  # Grid points per dimension

Grids and Modes

import numpy as np
import dedalus.public as d3

coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.float64)

xbasis = d3.RealFourier(coords['x'], size=128, bounds=(0, 4))
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1))

# Global grids (all grid points)
x_global = xbasis.global_grid(dist, scale=1)
z_global = zbasis.global_grid(dist, scale=1)

# Local grids (on this MPI process)
x_local = xbasis.local_grid(dist, scale=1)
z_local = zbasis.local_grid(dist, scale=1)

# Or get both at once from distributor
x, z = dist.local_grids(xbasis, zbasis, scales=1)

# Native grids (before coordinate mapping)
x_native = xbasis._native_grid(scale=1)

# Mode numbers
modes = dist.local_modes(xbasis)

Dealiasing

Dealiasing prevents aliasing errors in nonlinear terms.
# Create basis with dealiasing
xbasis = d3.RealFourier(coords['x'], size=128, bounds=(0, 4), dealias=3/2)

# Grid has 3/2 times as many points
print(xbasis.grid_shape(scales=(1,)))    # (128,)
print(xbasis.grid_shape(scales=(3/2,)))  # (192,)

# Multidimensional
basis = d3.DiskBasis(coords, shape=(128, 64), radius=1, 
                     dealias=(3/2, 3/2), dtype=np.float64)
Guidelines:
  • Use dealias=3/2 for quadratic nonlinearities
  • Use dealias=2 for cubic nonlinearities
  • Higher dealias for higher-order nonlinearities

Derivatives and Conversions

Derivative Bases

# Get basis for derivative
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1), a0=-1/2, b0=-1/2)
zbasis_deriv = zbasis.derivative_basis(order=1)

print(zbasis)         # ChebyshevT (a=b=-1/2)
print(zbasis_deriv)   # Jacobi with a=b=1/2

Change of Variables

Bases use affine coordinate transformations:
# Basis defined on [0, 1]
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1))

# Native grid on [-1, 1]
native_grid = zbasis._native_grid(scale=1)

# Transformed to [0, 1]
problem_grid = zbasis.COV.problem_coord(native_grid)

# Stretch factor affects derivatives
stretch = zbasis.COV.stretch  # (1-0) / (1-(-1)) = 0.5

Complete Example

import numpy as np
import dedalus.public as d3
from mpi4py import MPI

# Parameters
Lx, Lz = 4, 1
Nx, Nz = 256, 64
Rayleigh = 2e6
Prandtl = 1
dealias = 3/2
dtype = np.float64

# Coordinate system
coords = d3.CartesianCoordinates('x', 'z')

# Distributor
comm = MPI.COMM_WORLD
dist = d3.Distributor(coords, comm=comm, dtype=dtype)

# Bases
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias)
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias)

# Fields
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 terms for boundary conditions (no basis in z)
tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
tau_b2 = dist.Field(name='tau_b2', bases=xbasis)

# Substitutions
x, z = dist.local_grids(xbasis, zbasis)
kappa = (Rayleigh * Prandtl)**(-1/2)

# Derivative basis for first-order reduction
lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)

# First-order formulation
grad_b = d3.grad(b) + d3.CartesianCoordinates('z').unit_vector_fields(dist)[1] * lift(tau_b1)

print(f"x-basis: {xbasis}")
print(f"z-basis: {zbasis}")
print(f"z-derivative basis: {lift_basis}")
print(f"Dealiasing factor: {dealias}")
print(f"Grid shape (scale=1): {xbasis.grid_shape((1,))[0]} x {zbasis.grid_shape((1,))[0]}")
print(f"Grid shape (scale={dealias}): {xbasis.grid_shape((dealias,))[0]} x {zbasis.grid_shape((dealias,))[0]}")

Best Practices

Basis Selection:
  • Use Fourier for periodic directions
  • Use Chebyshev for non-periodic boundaries
  • Use Legendre for weak boundary conditions
  • Use compound bases (Disk, Sphere) for complex geometries
Boundary Conditions: Chebyshev and Jacobi bases require explicit boundary conditions via tau terms. The number of tau terms must match the order of the highest derivative in the problem.
Resolution: Spectral convergence requires sufficient resolution for the smoothest scales in the solution. For turbulent problems, resolution should capture the dissipation scale with some safety margin.

See Also

Build docs developers (and LLMs) love