Skip to main content

​
Overview

Dedalus supports four main problem types for solving PDEs: Initial Value Problems (IVP), Eigenvalue Problems (EVP), Linear Boundary Value Problems (LBVP), and Nonlinear Boundary Value Problems (NLBVP). Each problem type is defined in dedalus/core/problems.py.

​
Initial Value Problems (IVP)

Time-dependent problems that evolve from initial conditions.

​
Problem Structure

IVPs support equations of the form: 
M·∂X/∂t + L·X = F(X, t)
Where:
  • M, L: Linear operators (time-independent)
  • X: Problem variables
  • F: Possibly nonlinear RHS (no explicit time derivatives)
  • t: Time variable

​
Creating an IVP

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), dealias=3/2)
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1), dealias=3/2)

# Define 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
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)

# Create IVP with time variable
problem = d3.IVP(
    [p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2],
    time='t',
    namespace=locals()
)

​
Parameters

  • variables: List of Field objects to solve for
  • time: Time variable name (default: ‘t’) or Field object
  • namespace: Dictionary for parsing string equations (recommended: locals())

​
Adding Equations

# Parameters
kappa = 1e-3
nu = 1e-3
Lz = 1

# Grids and substitutions
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)

# Add equations
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*ez + lift(tau_u2) = -u@grad(u)")

# Boundary conditions
problem.add_equation("b(z=0) = Lz")
problem.add_equation("u(z=0) = 0")
problem.add_equation("b(z=Lz) = 0")
problem.add_equation("u(z=Lz) = 0")

# Pressure gauge
problem.add_equation("integ(p) = 0")

​
Equation Constraints

IVP LHS Requirements:
  • Must be linear in problem variables
  • Must be time-independent
  • Must be first-order in time derivatives
IVP RHS Requirements:
  • Can be nonlinear in problem variables
  • Can depend on time
  • Cannot contain explicit time derivatives

​
Eigenvalue Problems (EVP)

Linear eigenvalue problems for stability analysis.

​
Problem Structure

EVPs support equations of the form: 
λ·M·X + L·X = 0
Where:
  • λ: Eigenvalue
  • M, L: Linear operators
  • X: Eigenvector (problem variables)

​
Creating an EVP

import numpy as np
import dedalus.public as d3

coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.complex128)  # Usually complex

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

# Define perturbation fields
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
b = dist.Field(name='b', bases=(xbasis, zbasis))
p = dist.Field(name='p', bases=(xbasis, zbasis))

# Tau terms
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)
tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
tau_p = dist.Field(name='tau_p')

# Create eigenvalue field
eigenvalue = dist.Field(name='λ')

# Create EVP
problem = d3.EVP(
    [u, b, p, tau_u1, tau_u2, tau_b1, tau_b2, tau_p],
    eigenvalue,
    namespace=locals()
)

​
Adding Equations

# Parameters
Ra = 1e6
Pr = 1
kappa = (Ra * Pr)**(-1/2)
nu = (Ra / Pr)**(-1/2)

# Substitutions
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)

# EVP equations (λ appears explicitly)
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("λ*b - kappa*div(grad_b) + lift(tau_b2) = 0")
problem.add_equation("λ*u - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = 0")

# Boundary conditions
problem.add_equation("b(z=0) = 0")
problem.add_equation("u(z=0) = 0")
problem.add_equation("b(z=1) = 0")
problem.add_equation("u(z=1) = 0")
problem.add_equation("integ(p) = 0")

​
Equation Constraints

EVP Requirements:
  • LHS must be linear in problem variables
  • LHS must be affine in the eigenvalue (can have constant term)
  • RHS must be zero

​
Converting IVP to EVP

For time-independent IVPs, you can automatically generate an EVP: 
# Create IVP
ivp = d3.IVP([u, b, p, ...], namespace=locals())
ivp.add_equation("dt(u) + grad(p) = 0")
# ... add more equations ...

# Convert to EVP for stability analysis
evp = ivp.build_EVP(eigenvalue=λ, 
                     backgrounds=[U, B, P],  # Background state
                     perturbations=[u, b, p])  # Perturbation fields
This converts dt(X) → λ*X and linearizes the RHS around the background state.

​
Linear Boundary Value Problems (LBVP)

Steady-state linear problems.

​
Problem Structure

LBVPs support equations of the form: 
L·X = F
Where:
  • L: Linear operator
  • X: Problem variables
  • F: Inhomogeneous term (independent of X)

​
Creating an LBVP

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))

# Define fields
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
p = dist.Field(name='p', bases=(xbasis, zbasis))

# Tau terms
tau_p = dist.Field(name='tau_p')
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)

# Create LBVP
problem = d3.LBVP(
    [u, p, tau_p, tau_u1, tau_u2],
    namespace=locals()
)

​
Adding Equations

# Parameters
nu = 1e-3
U0 = 1.0  # Forcing

# Substitutions  
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)

# Stokes flow equations
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("-nu*div(grad_u) + grad(p) + lift(tau_u2) = 0")

# Boundary conditions
problem.add_equation("u(z=0) = 0")
problem.add_equation("u(z=1) = U0*ex")  # Driven lid
problem.add_equation("integ(p) = 0")

​
Equation Constraints

LBVP Requirements:
  • LHS must be strictly linear in problem variables
  • RHS must be independent of problem variables
  • Cannot have affine terms on LHS (use RHS instead)

​
Nonlinear Boundary Value Problems (NLBVP)

Steady-state nonlinear problems solved via Newton iteration.

​
Problem Structure

NLBVPs support equations of the form: 
G(X) = H(X)
Solved via Newton-Kantorovich method: 
F(X) = G(X) - H(X) = 0
dF(X[n])·dX = -F(X[n])
X[n+1] = X[n] + dX

​
Creating an NLBVP

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))

# Define fields
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
T = dist.Field(name='T', bases=(xbasis, zbasis))
p = dist.Field(name='p', bases=(xbasis, zbasis))

# Tau terms
tau_p = dist.Field(name='tau_p')
tau_T1 = dist.Field(name='tau_T1', bases=xbasis)
tau_T2 = dist.Field(name='tau_T2', bases=xbasis)
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)

# Create NLBVP
problem = d3.NLBVP(
    [u, T, p, tau_p, tau_T1, tau_T2, tau_u1, tau_u2],
    namespace=locals()
)

​
Adding Equations

# Parameters
Ra = 1e6
Pr = 1
kappa = (Ra * Pr)**(-1/2)
nu = (Ra / Pr)**(-1/2)

# Substitutions
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_T = d3.grad(T) + ez*lift(tau_T1)

# Nonlinear steady-state equations
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("-kappa*div(grad_T) + lift(tau_T2) = -u@grad(T)")  # Nonlinear
problem.add_equation("-nu*div(grad_u) + grad(p) - T*ez + lift(tau_u2) = -u@grad(u)")  # Nonlinear

# Boundary conditions
problem.add_equation("T(z=0) = 1")
problem.add_equation("u(z=0) = 0")
problem.add_equation("T(z=1) = 0")
problem.add_equation("u(z=1) = 0")
problem.add_equation("integ(p) = 0")

​
Equation Constraints

NLBVP Requirements: No linearity restrictions - equations can be arbitrarily nonlinear. The Frechet differential is computed symbolically.

​
Boundary Conditions

All problem types support boundary conditions via interpolation.

​
Dirichlet Conditions

# Fixed value at boundary
problem.add_equation("u(z=0) = 0")  # No-slip at bottom
problem.add_equation("T(z=1) = 273")  # Fixed temperature at top

​
Neumann Conditions

# Fixed derivative at boundary
problem.add_equation("dz(T)(z=0) = 0")  # Insulating bottom

​
Robin Conditions

# Mixed condition
problem.add_equation("dz(T)(z=0) + alpha*T(z=0) = 0")  # Newton cooling

​
Stress Conditions

# Stress-free boundary
problem.add_equation("dz(u)(z=1) + ez*tau_u1 = 0")

​
Complete Examples

​
IVP: 2D Rayleigh-Bénard Convection

import numpy as np
import dedalus.public as d3

# Parameters
Lx, Lz = 4, 1
Nx, Nz = 256, 64
Rayleigh = 2e6
Prandtl = 1

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

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

# 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_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)

# Problem
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*ez + lift(tau_u2) = -u@grad(u)")

problem.add_equation("b(z=0) = Lz")
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")

# Initial conditions
b.fill_random('g', seed=42, distribution='normal', scale=1e-3)
b['g'] *= z * (Lz - z)
b['g'] += Lz - z

# Solver
solver = problem.build_solver(d3.RK222)
solver.stop_sim_time = 50

​
EVP: Stability Analysis

import numpy as np
import dedalus.public as d3

# Setup for EVP
coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.complex128)

kx = 3.117  # Horizontal wavenumber
xbasis = d3.ComplexFourier(coords['x'], size=1, bounds=(0, 2*np.pi/kx))
zbasis = d3.ChebyshevT(coords['z'], size=64, bounds=(0, 1))

# Perturbation fields
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
b = dist.Field(name='b', bases=(xbasis, zbasis))
p = dist.Field(name='p', bases=(xbasis, zbasis))

tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)
tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
tau_p = dist.Field(name='tau_p')

# Eigenvalue
eigenvalue = dist.Field(name='λ')

# Problem
problem = d3.EVP([u, b, p, tau_u1, tau_u2, tau_b1, tau_b2, tau_p],
                  eigenvalue, namespace=locals())

# EVP equations...
# (same as IVP but with dt(...) → λ*...)

# Solver
solver = problem.build_solver()
solver.solve_dense(solver.subproblems[0])
print(f"Eigenvalues: {solver.eigenvalues}")

​
Best Practices

Problem Selection:
  • Use IVP for time evolution problems
  • Use EVP for linear stability analysis
  • Use LBVP for linear steady-state problems
  • Use NLBVP for nonlinear steady-state solutions
Namespace Management: Always pass namespace=locals() when creating problems to ensure string equations can access all variables and operators in scope.
First-Order Reduction: For problems with Chebyshev/Jacobi bases and boundary conditions, use first-order reduction by defining auxiliary gradient variables with lifted tau terms.

​
See Also

Build docs developers (and LLMs) love

Get started for freeTalk to us