Skip to main content

​
Overview

The solvers module provides solver classes that build and solve the matrix systems from problem definitions. Different solver types correspond to different problem types.

​
Solver Classes

​
InitialValueSolver

InitialValueSolver(problem, timestepper, **kw)
Solver for initial value problems.

​
Parameters

  • problem (IVP): Initial value problem
  • timestepper: Timestepper class (e.g., RK222, RK443, SBDF2)
  • ncc_cutoff (float, optional): NCC expansion cutoff (default: 1e-6)
  • max_ncc_terms (int, optional): Maximum NCC terms (default: None)
  • entry_cutoff (float, optional): Matrix entry cutoff (default: 1e-12)
  • matsolver (str or class, optional): Matrix solver to use

​
Attributes

  • state (list of Fields): Problem variables
  • sim_time (float): Current simulation time
  • iteration (int): Current iteration number
  • stop_sim_time (float): Stopping simulation time
  • stop_wall_time (float): Stopping wall time
  • stop_iteration (int): Stopping iteration

​
Methods

step(dt) - Advance one timestep 
solver.step(0.01)  # Step with dt=0.01
proceed (property) - Check if solver should continue 
while solver.proceed:  
    solver.step(dt)
print_subproblem_ranks() - Print matrix ranks for diagnostics

​
Example: Heat Equation

import dedalus.public as d3  
import numpy as np

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

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

ν = 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")

# Initial condition
x = dist.local_grid(xbasis)
u['g'] = np.sin(np.pi * x)  

# Create solver
solver = problem.build_solver(d3.RK222)
solver.stop_sim_time = 1.0

# Timestepping loop
dt = 0.001
while solver.proceed:  
    solver.step(dt)
    if solver.iteration % 100 == 0:
        print(f't = {solver.sim_time:.3f}')

​
LinearBoundaryValueSolver

LinearBoundaryValueSolver(problem, **kw)
Solver for linear boundary value problems.

​
Parameters

  • problem (LBVP): Linear boundary value problem
  • ncc_cutoff (float, optional): NCC expansion cutoff
  • matsolver (str or class, optional): Matrix solver
  • Other parameters same as InitialValueSolver

​
Attributes

  • state (list of Fields): Problem variables (solution)
  • iteration (int): Solve iteration counter

​
Methods

solve(subproblems=None, rebuild_matrices=False) - Solve the BVP 
solver.solve()  # Solve for all subproblems
print_subproblem_ranks() - Print matrix ranks

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

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: ∇²u = f with u=0 on boundaries
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()

# Solution is in u['g']  
print(f"Max |u| = {np.max(np.abs(u['g']))}")

​
NonlinearBoundaryValueSolver

NonlinearBoundaryValueSolver(problem, **kw)
Solver for nonlinear BVPs using Newton-Kantorovich iteration.

​
Parameters

Same as LinearBoundaryValueSolver

​
Attributes

  • state (list of Fields): Problem variables
  • perturbations (list of Fields): Update directions
  • iteration (int): Newton iteration counter

​
Methods

newton_iteration(damping=1) - Perform one Newton iteration 
solver.newton_iteration(damping=0.5)  # With damping
print_subproblem_ranks() - Print matrix ranks

​
Example: Nonlinear Eigenvalue Problem

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

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

# Nonlinear BVP: u'' + ε*u³ = 0
ε = 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")

# Initial guess  
x = dist.local_grid(xbasis)
u['g'] = x

# Newton iteration
solver = problem.build_solver()
pert_norm = np.inf
while pert_norm > 1e-10:
    solver.newton_iteration()  
    pert_norm = sum(p.allreduce_data_norm('c', 2) for p in solver.perturbations)
    print(f'Perturbation norm: {pert_norm:.3e}')

​
EigenvalueSolver

EigenvalueSolver(problem, **kw)
Solver for linear eigenvalue problems.

​
Parameters

Same as LinearBoundaryValueSolver

​
Attributes

  • state (list of Fields): Eigenvector fields
  • eigenvalues (ndarray): Computed eigenvalues
  • eigenvectors (ndarray): Computed eigenvectors
  • eigenvalue_subproblem (Subproblem): Last solved subproblem

​
Methods

solve_dense(subproblem, rebuild_matrices=False, left=False) - Dense eigenvalue solve 
subproblem = solver.subproblems_by_group[(0,)]
solver.solve_dense(subproblem)
solve_sparse(subproblem, N, target, rebuild_matrices=False, left=False) - Sparse eigenvalue solve 
# Find 10 eigenvalues near target
solver.solve_sparse(subproblem, N=10, target=0+1j)
set_state(index, subsystem=0) - Set state to an eigenmode 
solver.set_state(0)  # Set to first eigenmode
print_subproblem_ranks(target=0) - Print matrix ranks

​
Example: Eigenvalue Problem

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

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

# EVP: u'' = λ*u with u(0)=u(1)=0
problem = d3.EVP([u, τ1, τ2], eigenvalue=λ, namespace=locals())  
problem.add_equation("dx(dx(u)) + lift(τ1,-1) + lift(τ2,+1) = λ*u")
problem.add_equation("u(x='left') = 0")
problem.add_equation("u(x='right') = 0")

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

# Extract eigenvalues
λ_vals = solver.eigenvalues
print(f"First 5 eigenvalues: {λ_vals[:5]}")

# Set state to first eigenmode  
solver.set_state(0)
x = dist.local_grid(xbasis)
import matplotlib.pyplot as plt
plt.plot(x, u['g'])
plt.title('First eigenmode')

​
Solver Options

​
Matrix Solvers

Available matrix solvers (specified via matsolver parameter):
  • ‘SuperluNaturalSpsolve’: Default sparse direct solver
  • ‘SuperluColamdSpsolve’: COLAMD ordering
  • ‘UmfpackSpsolve’: UMFPACK solver
  • ‘CsrMatrixSolver’: Generic CSR solver

​
Matrix Construction Options

  • bc_top (bool): Place BCs at top of matrix
  • tau_left (bool): Place tau columns at left
  • interleave_components (bool): Interleave tensor components
  • store_expanded_matrices (bool): Store preconditioned matrices

​
See Also

Build docs developers (and LLMs) love

Get started for freeTalk to us