Skip to main content

​
Overview

The problems module provides classes for defining different types of PDE problems: initial value problems (IVPs), linear boundary value problems (LBVPs), nonlinear boundary value problems (NLBVPs), and eigenvalue problems (EVPs).

​
Problem Types

​
InitialValueProblem (IVP)

InitialValueProblem(variables, time='t', namespace=None)
Initial value problems of the form: M·dt(X) + L·X = F(X, t)

​
Parameters

  • variables (list of Fields): Problem variables to solve for
  • time (str or Field, optional): Time variable name or field (default: ‘t’)
  • namespace (dict, optional): Namespace for parsing equations (recommended: locals())

​
Constraints

  • LHS must be linear in variables
  • LHS must be time-independent
  • LHS must be first-order in time derivatives
  • RHS must not contain time derivatives

​
Methods

add_equation(equation, condition=“True”) - Add equation to problem 
problem.add_equation("dt(u) - ν*lap(u) = -u@grad(u)")
build_solver(timestepper) - Create solver 
solver = problem.build_solver(d3.RK222)
build_EVP(eigenvalue, backgrounds, perturbations) - Convert to EVP

​
Example: Heat Equation

import dedalus.public as d3
import numpy as np

# Setup  
coord = d3.Coordinate('x')
dist = d3.Distributor(coord)
xbasis = d3.ChebyshevT(coord, size=64, bounds=(0, 1))

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

# Problem
ν = 0.01
problem = d3.IVP([u, τ1, τ2], namespace=locals())
problem.add_equation("dt(u) - ν*dx(dx(u)) + lift(τ1,-1) + lift(τ2,+1) = 0")
problem.add_equation("u(x='left') = 0")  
problem.add_equation("u(x='right') = 0")

# Solve
solver = problem.build_solver(d3.RK222)
solver.stop_sim_time = 10.0
while solver.proceed:
    solver.step(dt)

​
Example: Navier-Stokes

import dedalus.public as d3
import numpy as np

# Setup 2D periodic box
coords = d3.CartesianCoordinates('x', 'y')  
dist = d3.Distributor(coords)
xbasis = d3.RealFourier(coords['x'], size=256, bounds=(0, 2*np.pi))
ybasis = d3.RealFourier(coords['y'], size=256, bounds=(0, 2*np.pi))

# Fields
u = dist.VectorField(coords, bases=(xbasis, ybasis), name='u')  
p = dist.Field(bases=(xbasis, ybasis), name='p')

# Problem
ν = 1e-4
problem = d3.IVP([u, p], namespace=locals())
problem.add_equation("dt(u) + grad(p) - ν*lap(u) = -u@grad(u)")  
problem.add_equation("div(u) = 0")

# Timestepping
solver = problem.build_solver(d3.RK443)

​
LinearBoundaryValueProblem (LBVP)

LinearBoundaryValueProblem(variables, namespace=None)
Linear boundary value problems: L·X = F

​
Parameters

  • variables (list of Fields): Problem variables
  • namespace (dict, optional): Namespace for parsing

​
Constraints

  • LHS must be linear in variables
  • RHS must be independent of variables

​
Methods

add_equation(equation, condition) - Add equation build_solver() - Create solver

​
Example: Poisson Equation

import dedalus.public as d3
import numpy as np

# Setup
coord = d3.Coordinate('x')
dist = d3.Distributor(coord)
xbasis = d3.ChebyshevT(coord, size=128, bounds=(0, 1))  

# Fields
u = dist.Field(bases=xbasis, name='u')
τ1 = dist.Field(name='τ1')
τ2 = dist.Field(name='τ2')
f = dist.Field(bases=xbasis, name='f')

# Source term  
x = dist.local_grid(xbasis)
f['g'] = -np.sin(2*np.pi*x)

# Problem
problem = d3.LBVP([u, τ1, τ2], namespace=locals())
problem.add_equation("dx(dx(u)) + lift(τ1,-1) + lift(τ2,+1) = f")  
problem.add_equation("u(x='left') = 0")
problem.add_equation("u(x='right') = 0")

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

​
NonlinearBoundaryValueProblem (NLBVP)

NonlinearBoundaryValueProblem(variables, namespace=None)
Nonlinear boundary value problems: G(X) - H(X) = 0 Solved using Newton-Kantorovich iteration.

​
Methods

add_equation(equation, condition) - Add equation build_solver() - Create solver

​
Example

import dedalus.public as d3

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

# Nonlinear problem
ε = 0.1
problem = d3.NLBVP([u, τ1, τ2], namespace=locals())
problem.add_equation("dx(dx(u)) + ε*u**3 + lift(τ1,-1) + lift(τ2,+1) = 0")  
problem.add_equation("u(x='left') = 0")
problem.add_equation("u(x='right') = 1")

# Solve with Newton iteration
solver = problem.build_solver()
pert_norm = np.inf
while pert_norm > 1e-10:  
    solver.newton_iteration()
    pert_norm = sum(pert.allreduce_data_norm('c', 2) for pert in solver.perturbations)

​
EigenvalueProblem (EVP)

EigenvalueProblem(variables, eigenvalue, namespace=None)
Linear eigenvalue problems: λ·M·X + L·X = 0

​
Parameters

  • variables (list of Fields): Eigenvector fields
  • eigenvalue (Field): Eigenvalue field
  • namespace (dict, optional): Namespace

​
Constraints

  • LHS must be linear in variables
  • LHS must be affine in eigenvalue
  • RHS must be zero

​
Methods

add_equation(equation, condition) - Add equation build_solver() - Create EVP solver

​
Example: Orr-Sommerfeld

import dedalus.public as d3
import numpy as np

# Setup
coord = d3.Coordinate('z')  
dist = d3.Distributor(coord)
zbasis = d3.ChebyshevT(coord, size=128, bounds=(-1, 1))

# Fields
v = dist.Field(bases=zbasis, name='v')
η = dist.Field(bases=zbasis, name='η')  
λ = dist.Field(name='λ')

# Base flow
U = dist.Field(bases=zbasis)
z = dist.local_grid(zbasis)
U['g'] = 1 - z**2  # Poiseuille flow

# Parameters
Re = 5000
k = 1.0  # Wavenumber

# Orr-Sommerfeld equation
problem = d3.EVP([v, η], eigenvalue=λ, namespace=locals())  
problem.add_equation("λ*(dz(dz(v)) - k**2*v) + 1j*k*U*(dz(dz(v)) - k**2*v) - 1j*k*dz(dz(U))*v - (1/Re)*(dz(dz(dz(dz(v)))) - 2*k**2*dz(dz(v)) + k**4*v) + lift(η,-1) = 0")
problem.add_equation("v(z=-1) = 0")
problem.add_equation("v(z=+1) = 0")  

# Solve EVP
solver = problem.build_solver()
subproblem = solver.subproblems_by_group[(0,)]
solver.solve_dense(subproblem)

# Extract eigenvalues  
λ_vals = solver.eigenvalues

​
Aliases

IVP = InitialValueProblem
LBVP = LinearBoundaryValueProblem  
NLBVP = NonlinearBoundaryValueProblem
EVP = EigenvalueProblem

​
Equation Parsing

Equations can be specified as:
  1. Strings (parsed using namespace): 
    problem.add_equation("dt(u) - ν*lap(u) = 0")
  2. Expression tuples: 
    problem.add_equation((d3.dt(u) - ν*d3.lap(u), 0))

​
See Also

Build docs developers (and LLMs) love

Get started for freeTalk to us