Skip to main content

​
Overview

Cartesian initial value problems (IVPs) involve time-evolving PDEs on rectangular domains. Dedalus uses Fourier bases for periodic directions and Chebyshev bases for bounded directions with boundary conditions. These examples demonstrate:
  • Setting up 1D and 2D Cartesian problems
  • Implementing boundary conditions with tau terms
  • Using CFL-adaptive timestepping
  • Saving and visualizing simulation data

​
Rayleigh-Bénard Convection

File: examples/ivp_2d_rayleigh_benard/rayleigh_benard.py

​
Description

This example simulates 2D horizontally-periodic Rayleigh-Bénard convection—the classic problem of thermal convection between heated horizontal plates. When heated from below, the fluid becomes unstable and spontaneously forms convection rolls.

​
Physical Setup

  • Domain: Periodic in x, bounded in z with height Lz = 1 and width Lx = 4
  • Boundary conditions: No-slip walls with fixed temperatures (hot bottom, cold top)
  • Non-dimensionalization: Height and freefall time
  • Parameters: Ra = 2×10⁶, Pr = 1

​
Governing Equations

The Boussinesq equations for incompressible convection: 
# Continuity
problem.add_equation("trace(grad_u) + tau_p = 0")

# Buoyancy equation
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")

# Momentum equation  
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")
problem.add_equation("integ(p) = 0")  # Pressure gauge

​
Complete Code

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

# Parameters
Lx, Lz = 4, 1
Nx, Nz = 256, 64
Rayleigh = 2e6
Prandtl = 1
dealias = 3/2
stop_sim_time = 50
timestepper = d3.RK222
max_timestep = 0.125
dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=dtype)
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_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)  # First-order reduction
grad_b = d3.grad(b) + ez*lift(tau_b1)  # First-order reduction

# 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")  # Pressure gauge

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

# Initial conditions
b.fill_random('g', seed=42, distribution='normal', scale=1e-3)
b['g'] *= z * (Lz - z)  # Damp noise at walls
b['g'] += Lz - z  # Add linear background

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50)
snapshots.add_task(b, name='buoyancy')
snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity')

# CFL
CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.5, threshold=0.05,
             max_change=1.5, min_change=0.5, max_dt=max_timestep)
CFL.add_velocity(u)

# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=10)
flow.add_property(np.sqrt(u@u)/nu, name='Re')

# Main loop
try:
    logger.info('Starting main loop')
    while solver.proceed:
        timestep = CFL.compute_timestep()
        solver.step(timestep)
        if (solver.iteration-1) % 10 == 0:
            max_Re = flow.max('Re')
            logger.info('Iteration=%i, Time=%e, dt=%e, max(Re)=%f' %(solver.iteration, solver.sim_time, timestep, max_Re))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()

​
Running the Example

# Run simulation (parallel)
mpiexec -n 4 python3 rayleigh_benard.py

# Plot snapshots
mpiexec -n 4 python3 plot_snapshots.py snapshots/*.h5

​
Expected Results

  • Runtime: ~5 cpu-minutes
  • Output: Buoyancy and vorticity snapshots in snapshots/
  • Physics: Convection rolls develop from random noise, reaching a turbulent steady state

​
KdV-Burgers Equation

File: examples/ivp_1d_kdv_burgers/kdv_burgers.py

​
Description

Solves the 1D Korteweg-de Vries / Burgers equation, demonstrating the interplay between nonlinear steepening (Burgers), linear dispersion (KdV), and dissipation. This simple example runs in seconds and creates a space-time plot.

​
Equation

∂u/∂t + u∂u/∂x = a∂²u/∂x² + b∂³u/∂x³
where a is the dissipation coefficient and b is the dispersion coefficient.

​
Complete Code

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  # Dissipation
b = 2e-4  # Dispersion
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
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',1].copy()]
t_list = [solver.sim_time]
while solver.proceed:
    solver.step(timestep)
    if solver.iteration % 100 == 0:
        logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
    if solver.iteration % 25 == 0:
        u_list.append(u['g',1].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', rasterized=True, 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.tight_layout()
plt.savefig('kdv_burgers.pdf')
plt.savefig('kdv_burgers.png', dpi=200)

​
Running the Example

python3 kdv_burgers.py

​
Expected Results

  • Runtime: A few seconds (serial only)
  • Output: Space-time plot kdv_burgers.png
  • Physics: Initial soliton-like pulse evolves with dispersion and dissipation

​
Shear Flow

File: examples/ivp_2d_shear_flow/shear_flow.py

​
Description

Simulates 2D periodic incompressible shear flow with a passive tracer for visualization. Demonstrates the Kelvin-Helmholtz instability where initially smooth shear layers roll up into vortices.

​
Physical Setup

  • Domain: Fully periodic in both x and z
  • Initial condition: Two shear layers with small vertical perturbations
  • Parameters: Re = 5×10⁴, Sc = 1
  • Non-dimensionalization: Shear-layer spacing and velocity jump

​
Key Features

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

# No tau terms needed for fully periodic problem
problem = d3.IVP([u, s, p, tau_p], namespace=locals())
problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0")  # Pressure gauge

​
Initial Conditions

# Background shear
u['g'][0] = 1/2 + 1/2 * (np.tanh((z-0.5)/0.1) - np.tanh((z+0.5)/0.1))
# Match tracer to shear
s['g'] = u['g'][0]
# Add small vertical velocity perturbations
u['g'][1] += 0.1 * np.sin(2*np.pi*x/Lx) * np.exp(-(z-0.5)**2/0.01)
u['g'][1] += 0.1 * np.sin(2*np.pi*x/Lx) * np.exp(-(z+0.5)**2/0.01)

​
Running the Example

mpiexec -n 4 python3 shear_flow.py
mpiexec -n 4 python3 plot_snapshots.py snapshots/*.h5

​
Expected Results

  • Runtime: ~10 cpu-minutes
  • Physics: Kelvin-Helmholtz billows form and evolve into turbulent structures
  • Visualization: Tracer field beautifully shows vortex roll-up

Build docs developers (and LLMs) love

Get started for freeTalk to us