Quick Start

Overview

This tutorial will walk you through your first Dedalus simulation step by step. We’ll solve the 2D Poisson equation with mixed boundary conditions - a simple but complete example that demonstrates the core workflow.

This example should run in just a few seconds and requires no MPI (single process). It’s perfect for getting started!

Prerequisites

Before starting, make sure you have:

Dedalus installed (see Installation)
Basic familiarity with Python and PDEs
A text editor or Python IDE

Ensure threading is disabled for best performance: export OMP_NUM_THREADS=1

The Problem

We’ll solve the 2D Poisson equation on a rectangular domain:

\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)

With boundary conditions:

$u(y=0) = g(x)$ (Dirichlet condition at bottom)
$\frac{\partial u}{\partial y}(y=L_y) = h(x)$ (Neumann condition at top)
Periodic in $x$

Complete Example Code

Create a file called poisson.py with the following code:

poisson.py

"""
Dedalus script solving the 2D Poisson equation with mixed boundary conditions.
This demonstrates solving a 2D Cartesian linear boundary value problem.
"""

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

# Parameters
Lx, Ly = 2*np.pi, np.pi
Nx, Ny = 256, 128
dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(0, Ly))

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

# Forcing
x, y = dist.local_grids(xbasis, ybasis)
f = dist.Field(bases=(xbasis, ybasis))
g = dist.Field(bases=xbasis)
h = dist.Field(bases=xbasis)
f.fill_random('g', seed=40)
f.low_pass_filter(shape=(64, 32))
g['g'] = np.sin(8*x) * 0.025
h['g'] = 0

# Substitutions
dy = lambda A: d3.Differentiate(A, coords['y'])
lift_basis = ybasis.derivative_basis(2)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

# Problem
problem = d3.LBVP([u, tau_1, tau_2], namespace=locals())
problem.add_equation("lap(u) + lift(tau_1,-1) + lift(tau_2,-2) = f")
problem.add_equation("u(y=0) = g")
problem.add_equation("dy(u)(y=Ly) = h")

# Solver
solver = problem.build_solver()
solver.solve()

# Gather global data
x = xbasis.global_grid(dist, scale=1)
y = ybasis.global_grid(dist, scale=1)
ug = u.allgather_data('g')

# Plot
if dist.comm.rank == 0:
    plt.figure(figsize=(6, 4))
    plt.pcolormesh(x.ravel(), y.ravel(), ug.T, cmap='viridis', 
                   shading='gouraud', rasterized=True)
    plt.gca().set_aspect('equal')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title("Randomly forced Poisson equation")
    plt.colorbar(label='u')
    plt.tight_layout()
    plt.savefig('poisson.png', dpi=200)
    print("Plot saved to poisson.png")

Step-by-Step Explanation

Import Dedalus

Import the Dedalus public interface and other required packages:

import dedalus.public as d3
import numpy as np
import matplotlib.pyplot as plt

The d3 module contains all the main Dedalus classes and functions for version 3.

Define Domain Parameters

Set the domain size and resolution:

Lx, Ly = 2*np.pi, np.pi  # Domain size
Nx, Ny = 256, 128         # Number of grid points
dtype = np.float64        # Data type for fields

Higher resolution gives more accurate results but takes longer to solve. Start with these values and experiment.

Set Up Coordinate System and Bases

Create the coordinate system and spectral bases:

# Define 2D Cartesian coordinates
coords = d3.CartesianCoordinates('x', 'y')

# Create distributor for MPI parallelization
dist = d3.Distributor(coords, dtype=dtype)

# Define spectral bases
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(0, Ly))

RealFourier: For periodic directions (efficient for real-valued fields)
Chebyshev: For non-periodic directions with boundaries
Distributor: Handles MPI parallelization automatically

Create Fields

Define the solution field and tau terms for boundary conditions:

# Solution field on both bases
u = dist.Field(name='u', bases=(xbasis, ybasis))

# Tau terms for boundary conditions
tau_1 = dist.Field(name='tau_1', bases=xbasis)
tau_2 = dist.Field(name='tau_2', bases=xbasis)

For a scalar Laplacian with two boundaries, we need two tau terms. These enforce the boundary conditions using the tau method.

Set Up Forcing and Boundary Data

Create the right-hand side and boundary condition fields:

x, y = dist.local_grids(xbasis, ybasis)

# Random forcing function
f = dist.Field(bases=(xbasis, ybasis))
f.fill_random('g', seed=40)
f.low_pass_filter(shape=(64, 32))

# Boundary conditions
g = dist.Field(bases=xbasis)
g['g'] = np.sin(8*x) * 0.025

h = dist.Field(bases=xbasis)
h['g'] = 0

The 'g' notation refers to grid space (as opposed to coefficient space).

Define Operators

Set up differential operators and lifting:

# Differentiation operator
dy = lambda A: d3.Differentiate(A, coords['y'])

# Lifting for tau terms
lift_basis = ybasis.derivative_basis(2)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

Lifting raises the tau terms to the appropriate basis for the PDE.

Formulate the Problem

Create the linear boundary value problem:

problem = d3.LBVP([u, tau_1, tau_2], namespace=locals())
problem.add_equation("lap(u) + lift(tau_1,-1) + lift(tau_2,-2) = f")
problem.add_equation("u(y=0) = g")
problem.add_equation("dy(u)(y=Ly) = h")

Notice how the equations are written symbolically, just like mathematical notation! Dedalus handles the discretization automatically.

lap(u): Laplacian operator $\nabla^2 u$
lift(...): Incorporates tau terms into the equation
u(y=0): Boundary evaluation

Build and Solve

Build the solver and solve the linear system:

solver = problem.build_solver()
solver.solve()

This constructs the sparse linear system and solves it using efficient direct methods.

Extract and Visualize Results

Gather the solution and create a plot:

# Get global grids and solution
x = xbasis.global_grid(dist, scale=1)
y = ybasis.global_grid(dist, scale=1)
ug = u.allgather_data('g')

# Plot (only on rank 0 for MPI)
if dist.comm.rank == 0:
    plt.figure(figsize=(6, 4))
    plt.pcolormesh(x.ravel(), y.ravel(), ug.T, cmap='viridis')
    plt.colorbar(label='u')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title("Poisson equation solution")
    plt.savefig('poisson.png', dpi=200)

Running the Example

Run the script from the command line:

# Run on a single process
python3 poisson.py

Expected Output

You should see:

Console output: Minimal logging from Dedalus
Image file: poisson.png showing the solution field

The solution should show a smooth 2D field satisfying the boundary conditions with the random forcing pattern.

The script should complete in just a few seconds. If it takes much longer, check that threading is disabled (OMP_NUM_THREADS=1).

Key Concepts Demonstrated

This example showcases several important Dedalus features:

Symbolic Equations

Write PDEs naturally using lap(), grad(), div() operators

Multiple Bases

Combine Fourier and Chebyshev bases for different directions

Boundary Conditions

Enforce Dirichlet and Neumann conditions using tau method

MPI Ready

Same code works in serial or parallel with no changes

Next Example: Time Evolution

Now let’s try an initial value problem (IVP) that evolves in time. Here’s the 1D KdV-Burgers equation:

kdv_burgers.py

"""
Dedalus script simulating the 1D Korteweg-de Vries / Burgers equation.
"""

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

# Parameters
Lx = 10
Nx = 1024
a = 1e-4  # Diffusion coefficient
b = 2e-4  # Dispersion coefficient
dealias = 3/2
stop_sim_time = 10
timestepper = d3.SBDF2
timestep = 2e-3
dtype = np.float64

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

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

# Substitutions
dx = lambda A: d3.Differentiate(A, xcoord)

# Problem: dt(u) + u*dx(u) = a*dx(dx(u)) + b*dx(dx(dx(u)))
problem = d3.IVP([u], namespace=locals())
problem.add_equation("dt(u) - a*dx(dx(u)) - b*dx(dx(dx(u))) = - u*dx(u)")

# Initial conditions
x = dist.local_grid(xbasis)
n = 20
u['g'] = np.log(1 + np.cosh(n)**2/np.cosh(n*(x-0.2*Lx))**2) / (2*n)

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

# Main loop
u_list = [u['g'].copy()]
t_list = [solver.sim_time]

while solver.proceed:
    solver.step(timestep)
    if solver.iteration % 100 == 0:
        logger.info(f'Iteration={solver.iteration}, Time={solver.sim_time:.2e}')
    if solver.iteration % 25 == 0:
        u_list.append(u['g'].copy())
        t_list.append(solver.sim_time)

# Plot
plt.figure(figsize=(6, 4))
plt.pcolormesh(x.ravel(), np.array(t_list), np.array(u_list), 
               cmap='RdBu_r', shading='gouraud', clim=(-0.8, 0.8))
plt.xlim(0, Lx)
plt.ylim(0, stop_sim_time)
plt.xlabel('x')
plt.ylabel('t')
plt.title(f'KdV-Burgers, (a,b)=({a},{b})')
plt.colorbar(label='u')
plt.tight_layout()
plt.savefig('kdv_burgers.png', dpi=200)
print("Plot saved to kdv_burgers.png")

Run it with:

python3 kdv_burgers.py

Key Differences for Time Evolution

IVP instead of LBVP: Initial Value Problem for time-dependent equations
Timestepper: Use d3.SBDF2 (semi-implicit backward differentiation)
Main loop: Iterate with solver.step() until solver.proceed is False
Dealiasing: Set dealias=3/2 to prevent aliasing errors in nonlinear terms

Common Patterns

Problem Types

IVP
LBVP
EVP

Initial Value Problem - time evolution:

problem = d3.IVP([u, v], namespace=locals())
problem.add_equation("dt(u) = ...")
solver = problem.build_solver(d3.RK222)
while solver.proceed:
    solver.step(dt)

Linear Boundary Value Problem - elliptic equations:

problem = d3.LBVP([u, tau], namespace=locals())
problem.add_equation("lap(u) + ... = f")
problem.add_equation("u(y=0) = 0")
solver = problem.build_solver()
solver.solve()

Eigenvalue Problem - stability analysis:

problem = d3.EVP([u, tau], eigenvalue='omega', namespace=locals())
problem.add_equation("omega*u = ...")
solver = problem.build_solver()
solver.solve_dense(solver.subproblems[0])

Useful Operators

# Differential operators
d3.div(u)          # Divergence
d3.grad(u)         # Gradient
d3.lap(u)          # Laplacian
d3.curl(u)         # Curl
d3.Differentiate(u, coord)  # Directional derivative

# Vector operations
u @ v              # Dot product
u @ d3.grad(v)     # Advection u·∇v
d3.cross(u, v)     # Cross product
d3.skew(u)         # Skew-symmetric gradient

# Integration
d3.integ(u)        # Volume integral
d3.Average(u, coord)  # Average over coordinate

Tips and Best Practices

Choosing resolution

Start with moderate resolution and increase until results converge:

Monitor spectral coefficients to check if they decay to zero
Use field.allgather_data('c') to examine coefficients
Increase resolution by factors of 2 for testing

Handling nonlinear terms

For equations with nonlinear terms like

u \frac{\partial u}{\partial x}

Always use dealiasing: dealias=3/2
Write as: u*dx(u) or u@grad(u) for vectors
Dedalus automatically handles the dealiasing

Saving data

Use the built-in analysis system for large simulations:

analysis = solver.evaluator.add_file_handler('snapshots', sim_dt=0.1)
analysis.add_task(u, name='velocity')
analysis.add_task(d3.div(u), name='divergence')

Data is saved to HDF5 files that can be read with h5py or xarray.

Performance optimization

Always set OMP_NUM_THREADS=1
Use MPI for parallelism: mpiexec -n N
Profile with: python3 -m cProfile -o output.prof script.py
Monitor with: flow.add_property() for runtime diagnostics

Where to Go Next

Example Gallery

Browse complete examples including fluid dynamics, eigenvalue problems, and more

Core Concepts

Deep dive into bases, operators, boundary conditions, and tau method

API Reference

Complete documentation of all classes and functions

Community

Join the mailing list to ask questions and share your work

Troubleshooting

If you encounter issues:

Check threading: echo $OMP_NUM_THREADS should return 1
Verify installation: python3 -m dedalus test
Check MPI: Run with -n 1 first, then try parallel
Examine logs: Dedalus uses Python logging for diagnostics
Ask for help: Post to the mailing list

Remember: Most performance issues are due to threading not being disabled. Always check OMP_NUM_THREADS=1!

Get Started

Core Concepts

Advanced Topics

Examples

Overview

Prerequisites

The Problem

Complete Example Code

Step-by-Step Explanation

Running the Example

Expected Output

Key Concepts Demonstrated

Symbolic Equations

Multiple Bases

Boundary Conditions

MPI Ready

Next Example: Time Evolution

Key Differences for Time Evolution

Common Patterns

Problem Types

Useful Operators

Tips and Best Practices

Where to Go Next

Example Gallery

Core Concepts

API Reference

Community

Troubleshooting

Build docs developers (and LLMs) love

Get Started

Core Concepts

Advanced Topics

Examples

​Overview

​Prerequisites

​The Problem

​Complete Example Code

​Step-by-Step Explanation

​Running the Example

​Expected Output

​Key Concepts Demonstrated

Symbolic Equations

Multiple Bases

Boundary Conditions

MPI Ready

​Next Example: Time Evolution

​Key Differences for Time Evolution

​Common Patterns

​Problem Types

​Useful Operators

​Tips and Best Practices

​Where to Go Next

Example Gallery

Core Concepts

API Reference

Community

​Troubleshooting

Build docs developers (and LLMs) love

Overview

Prerequisites

The Problem

Complete Example Code

Step-by-Step Explanation

Running the Example

Expected Output

Key Concepts Demonstrated

Next Example: Time Evolution

Key Differences for Time Evolution

Common Patterns

Problem Types

Useful Operators

Tips and Best Practices

Where to Go Next

Troubleshooting