Overview
Fields are the fundamental data containers in Dedalus, representing scalar, vector, and tensor quantities on domains. Operators act on fields to build equation expressions.
Field Types
All field classes are defined in dedalus/core/field.py.
Creating Fields
Field (Generic)
General field with arbitrary tensor signature.
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))
# Scalar field (no tensorsig)
f = dist.Field(name='f', bases=(xbasis, zbasis))
# Field without spatial bases (constant)
const = dist.Field(name='const', dtype=np.float64)
ScalarField
Scalar field with no tensor indices.
# Equivalent to Field with no tensorsig
s = dist.ScalarField(name='s', bases=(xbasis, zbasis))
VectorField
Vector field with one tensor index.
# Vector in Cartesian coordinates
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
print(u.tensorsig) # (coords,)
print(u['g'].shape) # (2, Nx, Nz) for 2D coordinates
TensorField
Tensor field with two or more indices.
# Rank-2 tensor (e.g., stress tensor)
T = dist.TensorField((coords, coords), name='T', bases=(xbasis, zbasis))
print(T.tensorsig) # (coords, coords)
print(T['g'].shape) # (2, 2, Nx, Nz) for 2D coordinates
Field Properties
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))
u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
# Basic properties
print(u.name) # 'u'
print(u.domain) # Domain object
print(u.dist) # Distributor object
print(u.dtype) # np.float64
print(u.tensorsig) # (coords,)
# Data access
print(u.data) # Current layout data
print(u.layout) # Current layout object
print(u['g']) # Grid space data
print(u['c']) # Coefficient space data
# Shape information
print(u.data.shape) # Shape in current layout
print(u.global_shape()) # Global shape
print(u.local_shape()) # Local shape on this process
Accessing Field Data
Layout System
Fields can be in grid space (‘g’) or coefficient space (‘c’).
# Access grid space data
f['g'] = np.sin(2*np.pi*x/Lx) * np.cos(np.pi*z/Lz)
# Access coefficient space data
f['c'][0, 0] = 1.0
# Check current layout
print(f.layout) # Layout object
# Change layout explicitly
f.change_layout('c') # Transform to coefficient space
f.change_layout('g') # Transform to grid space
Indexing
# Scalar field
f['g'][i, j] = value # Grid space at point (i, j)
f['c'][kx, kz] = value # Coefficient at mode (kx, kz)
# Vector field
u['g'][0, i, j] = value # x-component at (i, j)
u['g'][1, i, j] = value # z-component at (i, j)
# Tensor field
T['g'][0, 1, i, j] = value # T_xz component at (i, j)
Slicing
# Get x-component of vector field
ux = u['g'][0]
# Get slice at z=constant
f_slice = f['g'][:, k]
# Set entire field
f['g'][:] = 0 # Zero out
f['g'] = other_field['g'] # Copy
Setting Field Values
From Arrays
import numpy as np
# Set directly
f['g'] = np.sin(2*np.pi*x/Lx)
# Set components
u['g'][0] = np.cos(2*np.pi*x/Lx)
u['g'][1] = 0
Using fill_random
# Random noise
f.fill_random('g', seed=42, distribution='normal', scale=1e-3)
# Random with specific distribution
f.fill_random('g', seed=42, distribution='uniform', low=-1, high=1)
Preset Scales
Control grid resolution for evaluation.
# Set scales for dealiasing
f.preset_scales(3/2)
# Access data at specified scale
print(f['g'].shape) # Shape with 3/2 dealiasing
Arithmetic Operations
Fields support standard arithmetic operations.
# Scalar operations
f2 = 2 * f + 1
f3 = f * g # Nonlinear term
f4 = f / g
f5 = f ** 2
# Vector operations
w = u + v # Vector addition
w = 2 * u # Scalar multiplication
# Dot product (returns scalar)
kinetic = u @ u # u·u
# Vector-vector multiplication (outer product, returns tensor)
T = u * v # Returns tensor with u_i * v_j
Differential Operators
All operators are defined in dedalus/core/operators.py.
Gradient
Compute the gradient of a scalar field.
import dedalus.public as d3
# Gradient of scalar field (returns vector)
grad_f = d3.grad(f)
# Or use string parsing
problem.add_equation("grad(f) = 0")
- Input: Scalar field
- Output: Vector field
- Components: ∇f = (∂f/∂x, ∂f/∂y, ∂f/∂z)
Divergence
Compute the divergence of a vector field.
# Divergence of vector field (returns scalar)
div_u = d3.div(u)
# In equations
problem.add_equation("div(u) = 0") # Incompressibility
- Input: Vector field
- Output: Scalar field
- Value: ∇·u = ∂u_x/∂x + ∂u_y/∂y + ∂u_z/∂z
Curl
Compute the curl of a vector field.
# Curl of vector field (returns vector in 3D, scalar in 2D)
vorticity = d3.curl(u)
# In equations
problem.add_equation("w - curl(u) = 0")
- Input: Vector field
- Output: Vector field (3D) or Scalar field (2D)
- 2D: ∇×u = ∂u_y/∂x - ∂u_x/∂y
- 3D: (∇×u)_i = ε_ijk ∂u_k/∂x_j
Laplacian
Compute the Laplacian of a field.
# Laplacian (works on scalars, vectors, tensors)
lap_f = d3.lap(f)
lap_u = d3.lap(u)
# In equations
problem.add_equation("dt(f) - kappa*lap(f) = 0") # Diffusion
- ∇²f = ∂²f/∂x² + ∂²f/∂y² + ∂²f/∂z²
- Acts component-wise on vectors and tensors
Trace
Compute the trace of a tensor.
# Trace of rank-2 tensor (returns scalar)
trace_T = d3.trace(T)
# Trace of gradient (equivalent to divergence)
problem.add_equation("trace(grad(u)) = 0") # div(u) = 0
Skew-Symmetric Part
Compute the skew-symmetric part of a tensor.
# Skew-symmetric part
skew_grad_u = d3.skew(d3.grad(u))
# Vorticity from skew(grad(u))
omega = -d3.div(d3.skew(d3.grad(u))) # Equivalent to curl in 2D
Transpose
Transpose tensor components.
# Transpose tensor
T_transpose = d3.TransposeComponents(T)
# Or use property
T_transpose = T.T
Advanced Operators
Interpolation
Interpolate field to a specific location.
# Interpolate to z=0
f_bottom = d3.Interpolate(f, z=0)
# Or use parentheses syntax
f_bottom = f(z=0)
# In equations (boundary conditions)
problem.add_equation("f(z=0) = 1")
Integration
Integrate field over a coordinate.
# Integrate over z
f_integrated = d3.Integrate(f, coords['z'])
# In equations (pressure gauge)
problem.add_equation("integ(p) = 0")
Averaging
Compute average over a coordinate.
# Average over x
f_avg = d3.Average(f, coords['x'])
# In string parsing
problem.add_equation("u_mean - Average(u, x) = 0")
Lifting
Lift boundary data to the domain interior.
# Lift tau term for first-order reduction
lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
# Use in gradient
grad_u = d3.grad(u) + ez * lift(tau_u1)
Component Access
Unit Vectors
# Get coordinate unit vectors
ex, ey, ez = coords.unit_vector_fields(dist)
# Use in equations
problem.add_equation("dt(u) + grad(p) - b*ez = -u@grad(u)")
Vector Components
# Extract specific component
ux = u @ ex # x-component
uz = u @ ez # z-component
# Or use indexing
ux_data = u['g'][0] # x-component data
Radial/Angular Components
For curvilinear coordinates.
# Radial component in spherical coordinates
u_r = d3.RadialComponent(u)
# Angular component
u_angular = d3.AngularComponent(u)
# Azimuthal component
u_phi = d3.AzimuthalComponent(u)
Complete Example
import numpy as np
import dedalus.public as d3
# Setup
Lx, Lz = 4, 1
Nx, Nz = 256, 64
Rayleigh = 2e6
Prandtl = 1
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 terms (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)
# Parameters
kappa = (Rayleigh * Prandtl)**(-1/2)
nu = (Rayleigh / Prandtl)**(-1/2)
# Grids and unit vectors
x, z = dist.local_grids(xbasis, zbasis)
ex, ez = coords.unit_vector_fields(dist)
# First-order reduction
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())
# Equations using operators
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")
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)
print(f"Velocity field: {u}")
print(f"Velocity shape: {u['g'].shape}")
print(f"Gradient of u: {grad_u}")
print(f"Divergence: {d3.div(u)}")
print(f"Laplacian of b: {d3.lap(b)}")
Best Practices
Operator Composition:
Compose operators to build complex expressions:
div(grad(f)) for Laplaciantrace(grad(u)) for divergence-div(skew(grad(u))) for vorticity
Grid vs Coefficient Space:
Nonlinear operations must be performed in grid space. Dedalus automatically transforms fields as needed, but explicit change_layout('g') can avoid repeated transforms.
First-Order Reduction:
For problems with boundary conditions, use first-order reduction by defining auxiliary gradient variables and lifting tau terms. This ensures proper boundary treatment.
See Also