The Tau Method for Boundary Conditions

Tarang enforces boundary conditions using the tau method, a spectral technique that lets you solve PDEs with non-periodic boundary conditions without modifying the spectral basis. This page explains what the tau method is, how Tarang's implementation works, and how to write code that uses it correctly.

If you have used Dedalus, the approach here is almost identical: tau fields are added to the state vector, and lift() operators inject them into the equations as extra degrees of freedom to match boundary conditions.

Why Do We Need the Tau Method?

Consider the steady diffusion problem:

\[\frac{d^2 u}{dz^2} = f(z), \quad z \in [-1, 1], \quad u(-1) = u(+1) = 0.\]

In a Chebyshev spectral method, u is expanded as

\[u(z) = \sum_{n=0}^{N-1} a_n\, T_n(z),\]

and the PDE is projected onto the basis to give an N × N linear system in the coefficients $a_n$. The problem: this system has no leftover degrees of freedom to enforce the boundary conditions. All $N$ unknowns are already determined by the interior equations, but we also need two boundary constraints.

The tau method solves this by adding extra unknowns (the tau fields) that act as corrections designed to make the boundary conditions hold. The classical formulation, due to Lanczos (1938), replaces the highest-order rows of the spectral system with BC rows — but that makes the interior equation different at those rows, which is awkward for nonlinear and time-dependent problems.

Modern spectral codes (Dedalus, Tarang) use a cleaner variant: instead of replacing rows, they add a tau term to the equation itself:

\[\mathcal{L}[u] \;+\; \tau_1\,\phi_1(z) \;+\; \tau_2\,\phi_2(z) \;=\; f(z),\]

where $\phi_1, \phi_2$ are chosen basis functions that live in a "lift basis". The tau scalars $\tau_1, \tau_2$ become new unknowns in the system, and the boundary conditions become new equation rows that link them to the state. The interior of the PDE is untouched; the tau corrections only perturb the equation at a small number of high spectral modes, and for smooth solutions the required $\tau_k$ are spectrally small.

That is exactly how Tarang implements BCs.

Quick-Start Example

Here's the smallest complete verified example — a 2D Poisson problem Δu = -2 with homogeneous Dirichlet BCs on the two z walls, whose exact solution is u = z(Lz - z). The domain has one separable (Fourier) x direction and one bounded (Chebyshev) z direction; the Chebyshev direction is the coupled one that needs tau corrections:

using Tarang

coords = CartesianCoordinates("x", "z")
dist   = Distributor(coords; dtype=Float64, device=CPU())

xb = RealFourier(coords["x"]; size=4,  bounds=(0.0, 2π))    # periodic (separable)
zb = ChebyshevT(coords["z"];  size=16, bounds=(0.0, 1.0))   # bounded  (coupled)
dom = Domain(dist, (xb, zb))

# State field
u = ScalarField(dom, "u")

# Tau fields — one per boundary condition. Each lives on the COMPLEMENT of the
# state's bases: the coupled (Chebyshev) z axis is dropped, leaving (xb,).
tau1 = ScalarField(dist, "tau1", (xb,), Float64)
tau2 = ScalarField(dist, "tau2", (xb,), Float64)

# The lift basis. This equation is written as a DIRECT second-order operator
# (`Δ(u)`), so the second-derivative basis is the matching one. The first-order
# `grad_u`/`grad_T` formulation further down uses `derivative_basis(zb, 1)`.
lb2 = derivative_basis(zb, 2)

problem = LBVP([u, tau1, tau2])

# Register the lifts as named parameters, then reference them in the equation.
# lift(tau, lb2, -k) places tau at the k-th-from-last Chebyshev coefficient
# slot (-1 = last, -2 = second-to-last).
add_parameters!(problem; l1=lift(tau1, lb2, -1), l2=lift(tau2, lb2, -2))
add_equation!(problem, "Δ(u) + l1 + l2 = -2")

add_bc!(problem, "u(z=0)   = 0")
add_bc!(problem, "u(z=1.0) = 0")

solver = BoundaryValueSolver(problem)
solve!(solver)
ensure_layout!(u, :g)            # scatter writes coefficients; switch to grid

This reproduces the exact solution u = z(1 - z) to 1.4e-16, with both wall residuals at 0.0.

Notice that:

1. Tau fields (tau1, tau2) are added to the state vector, not computed from the state afterwards.
2. The lift terms are registered as named parameters via add_parameters! and then referenced by name (l1, l2) in the equation. Each lift uses the 3-arg form lift(tau, lift_basis, -k).
3. The BCs are declared with add_bc! (u(z=0) = 0) — Tarang converts each into an algebraic row that constrains the tau fields.

The solver builds one square tau subproblem per separable Fourier mode and determines u, tau1, and tau2 simultaneously in a single linear solve.

The lift Operator

lift is how a tau field enters an equation. Always use the explicit 3-argument form, with the lift basis built from derivative_basis, choosing the derivative order that matches the space the lifted term is added to:

lift(tau, derivative_basis(basis, order), -k)   # explicit lift basis

• tau is the tau field (a ScalarField or VectorField whose bases are a strict subset of the state field's bases — it's missing the coupled direction).
• derivative_basis(basis, order) is the lift basis. A direct second-order equation (Δ(u)) uses order = 2; the first-order grad_u/grad_T formulation used by the flow examples uses order = 1. Bind the result once and reuse it. See Why still pass the derivative basis? below for what this argument does — and does not — control today.
• -k is an integer mode index with wraparound semantics: -1 is the last coefficient slot of the coupled direction (lift_mode = N-1 → Julia index N), -2 is the second-to-last, 0 is the first, and so on. For a ChebyshevT state of size N, -1 places the tau value at coefficient N-1. A second-order problem uses two taus with lift orders -1 and -2.

The 2-argument short form does not work

lift(tau, n) — the basis-less short form that appears in the lift docstring — is broken in the current release. It tries to auto-detect the lift basis by scanning tau.dist.layouts, but that dictionary is keyed by (bases, dtype) pairs rather than by bases alone, so the scan hands a tuple to is_fourier_basis and dies:

MethodError: no method matching is_fourier_basis(::Tuple{RealFourier, ChebyshevT})

It fails this way regardless of the order in which you create your fields — creating the full state field first, so the distributor has already cached the complete layout, does not rescue it. Pass the lift basis explicitly instead — bind it once and reuse it:

lb1    = derivative_basis(zbasis, 1)
tau_u1 = ScalarField(dist, "tau_u1", (xbasis,), Float64)
grad_u = grad(u) + ez * lift(tau_u1, lb1, -1)

All of Tarang's shipped examples use this explicit form.

What lift(tau, basis, n) actually computes (solver view)

At the solver level, lift(tau, basis, n) resolves to a single-column sparse matrix with a 1 at row lift_mode (where lift_mode = n is wraparound-resolved: -1 → last, -2 → second-to-last, etc.). This column is a discrete delta in the coupled direction's coefficient space.

\[\text{lift}(\tau, \cdot, n) \;\longrightarrow\; \tau \cdot e_{\text{lift\_mode}}\]

where $e_{\text{lift\_mode}}$ is the unit vector with a 1 at the lift_mode-th Chebyshev-coefficient slot. When this is added to an equation's LHS, it contributes the unknown $\tau$ to exactly one coefficient row, leaving every other interior row untouched. The linear solver then chooses $\tau$ so that the BC rows hold — and because the perturbation is confined to one high-order coefficient, the interior PDE residual stays spectrally small for smooth solutions.

Implementation detail worth knowing: In Tarang's current solver, subproblem_matrix(op::Lift, sp) reads the dimension N from _subproblem_cheb_basis(sp) — the problem's own Chebyshev basis — and does not use op.basis to compute the matrix. That means lift(tau, zb, -1), lift(tau, derivative_basis(zb, 1), -1) and lift(tau, derivative_basis(zb, 2), -1) produce identical delta columns at row N-1. (Measured: all three spellings give bit-identical solutions and bit-identical tau values on the quick-start problem.) The basis argument is a semantic hint inherited from Dedalus's type system; it's retained so that future refinements — e.g. explicit basis tracking through expression trees — can hook in without breaking user code.

Why still pass the derivative basis?

Even though the explicit basis choices produce identical lift columns in the current solver, passing the intended derivative basis is worth doing:

1. Forward-compatibility: if Tarang later tracks the lift output basis in the matrix representation, your expression already identifies the intended space. Code that passed the state basis zb would silently start producing different matrices.
2. Semantic clarity: in first-order formulations like grad_u = grad(u) + ez * lift(tau_u1, lb1, -1), the lift lives alongside grad(u), which maps a ChebyshevT field into the first-derivative (ChebyshevU) space. Passing the derivative basis tells the reader what the lift is conceptually contributing, even if the linear algebra doesn't currently care.
3. Consistency with examples and tutorials: all of Tarang's shipped examples build the lift basis with derivative_basis, so sticking with it makes your code pattern-match familiar code.

Rule of thumb: bind the lift basis once at the top of your setup block and reuse it — derivative_basis(zb, 1) for a lift added to a first-order gradient substitution (which is what examples/ivp/rayleigh_benard_2d.jl does), and derivative_basis(zb, 2) for lifts added directly to a second-order equation (Δ(u), two lifts at -1 and -2). Don't invent your own lift basis.

Boundary Conditions as Algebraic Constraint Rows

When you write

add_bc!(problem, "u(z=-1) = 0")

Tarang produces a new equation row that looks like

\[\sum_{n=0}^{N-1} a_n\, T_n(-1) \;=\; 0,\]

i.e. a linear combination of u's coefficients that equals the boundary value. This row is added to the system, not substituted in. It's an algebraic equation — there is no time derivative, so it contributes nothing to the M matrix in the M·dX/dt + L·X = F formulation. In the linear-algebra picture:

• u contributes Nz rows (PDE interior) + BC rows
• tau_u1, tau_u2 contribute columns (one each, at the lift modes)
• BC rows are zero in M (no time derivative) → they're pure algebraic constraints
• The combined LHS matrix is square because (PDE rows + BC rows) = (state cols + tau cols)

This square-system condition is enforced automatically when you declare the right number of tau fields to match the number of BCs. If you miss a BC you'll get a singular / non-square system at solver-build time.

DAE-style handling in IVP steppers

For initial-value problems, BC rows have M_row = 0, which makes the full system a differential-algebraic equation (DAE) rather than a pure ODE. Tarang's subproblem stepper handles this correctly via a per-stage row override on what it classifies as "BC rows".

A note on how this classification works: sp.bc_rows is not built by scanning M for zero rows directly. Instead, build_matrices! labels rows by equation size — any equation whose per-subproblem row count is smaller than the coupled-direction basis size Nz is classified as a BC row (pressure gauge integ(p) = 0 → 1 row; wall BC T(z=0) = 1 → 1 row; vector wall BC → 2 rows; etc.). Rows from equations with eq_size ≥ Nz — the PDE interior and any Nz-sized algebraic equation like the continuity equation — are classified as "bulk".

At each IMEX-RK stage (or multistep update), the solver assembles the stage RHS from the accumulated dt·Σ(AᴱF − Aⁱ·L·X) formula for every row, then calls apply_bc_override! to replace the sp.bc_rows entries with dt·a_ii·F_BC. After the LHS solve, those rows yield exactly L_row·X = F_BC, which is the correct BC enforcement.

Why a size-based classifier is good enough: the only algebraic rows that don't get overridden are Nz-sized, F-zero equations like continuity. For them, the accumulated-RHS formula naturally produces zero on the right-hand side (since their F is zero and L*X = 0 is the target), so no override is needed. The override only matters for small algebraic rows whose F can be nonzero — i.e., BCs and gauge constraints — and those are exactly what sp.bc_rows captures.

Why the override is necessary: without it, the raw accumulated-RHS formula produces a wrong 1/γ scaling factor for inhomogeneous algebraic rows. For RK222 (γ = 1 - 1/√2, 1/γ = 2 + √2 ≈ 3.414), a BC like T(z=0) = 1 would be enforced as T(z=0) = 2+√2 instead of 1. The override is built into step_subproblem_rk! and step_subproblem_multistep!, so you don't have to do anything to enable it — but it's worth knowing the mechanism exists if you're debugging a "BC value is off by a constant factor" symptom.

First-Order Formulation (Recommended)

For time-dependent problems and anything involving the Navier–Stokes equations, the first-order formulation is the right default. Instead of writing a 2nd-order operator like Δ(u) directly, introduce an auxiliary gradient field with a tau correction. The lift basis here is the first-derivative basis, because the lift is added to grad(u), which lives in the first-derivative space:

lift_basis = derivative_basis(zbasis, 1)
τ_lift(A)  = lift(A, lift_basis, -1)

grad_u = grad(u) + ez * τ_lift(tau_u1)

Now grad_u is an "augmented gradient" that lives in the derivative (ChebyshevU) space and carries its own tau. Second derivatives become div(grad_u) rather than Δ(u). Register grad_u and τ_lift as parameters so the equation string can name them:

add_parameters!(problem, nu=nu, grad_u=grad_u, τ_lift=τ_lift)
add_equation!(problem, "∂t(u) - nu*div(grad_u) + ∇(p) + τ_lift(tau_u2) = -u⋅∇(u)")

How the two formulations compare:

The convention we use throughout Tarang's examples is:

• tau_u1 — correction added inside grad_u, one-dimensional (xbasis only)
• tau_u2 — correction added to the evolution equation directly, one-dimensional (xbasis only)

For a vector field u, both tau_u1 and tau_u2 are VectorFields too (one component per velocity component).

Worked Example: 2D Rayleigh–Bénard Convection

Here is the full first-order RBC setup, which is the examples/ivp/rayleigh_benard_2d.jl example in the repository. This is the canonical pattern to copy for any 2D channel-flow problem with non-periodic BCs.

using Tarang

# Domain and physical parameters
Lx, Lz = 4.0, 1.0
Nx, Nz = 256, 64
Rayleigh, Prandtl = 2e6, 1.0
nu   = Prandtl                 # viscous coefficient (diffusive-time scaling)
buoy = Rayleigh * Prandtl      # buoyancy forcing Ra·Pr

coords = CartesianCoordinates("x", "z")
dist   = Distributor(coords; dtype=Float64, device=CPU())
xbasis = RealFourier(coords["x"]; size=Nx, bounds=(0.0, Lx), dealias=3/2)
zbasis = ChebyshevT(coords["z"]; size=Nz, bounds=(0.0, Lz), dealias=3/2)

domain = Domain(dist, (xbasis, zbasis))

# State variables: pressure, temperature, velocity
p = ScalarField(domain, "p")
T = ScalarField(domain, "T")
u = VectorField(domain, "u")

# Tau fields. Each lives on `(xbasis,)` only — the coupled direction (z)
# is dropped because the tau correction lives at a single z mode.
#
# tau_p:  scalar, no bases — gauge for the pressure constraint
# tau_T1: gradient-substitution correction for T
# tau_T2: evolution-equation correction for T
# tau_u1, tau_u2: ditto for each velocity component (so they are VectorFields)
tau_p  = ScalarField(dist, "tau_p",  (),         Float64)
tau_T1 = ScalarField(dist, "tau_T1", (xbasis,),  Float64)
tau_T2 = ScalarField(dist, "tau_T2", (xbasis,),  Float64)
tau_u1 = VectorField(dist, coords, "tau_u1", (xbasis,), Float64)
tau_u2 = VectorField(dist, coords, "tau_u2", (xbasis,), Float64)

# First-order substitutions — FIRST-derivative lift basis
ex, ez = unit_vector_fields(coords, dist)
lift_basis = derivative_basis(zbasis, 1)
τ_lift(A) = lift(A, lift_basis, -1)

grad_u = grad(u) + ez * τ_lift(tau_u1)
grad_T = grad(T) + ez * τ_lift(tau_T1)

# Problem declaration includes ALL tau fields as state
problem = IVP([p, T, u, tau_p, tau_T1, tau_T2, tau_u1, tau_u2])

add_parameters!(problem,
    nu=nu, buoy=buoy, ez=ez,
    grad_u=grad_u, grad_T=grad_T, τ_lift=τ_lift)

# Equations
add_equation!(problem, "trace(grad_u) + tau_p = 0")
add_equation!(problem, "∂t(T) - div(grad_T) + τ_lift(tau_T2) = -u⋅∇(T)")
add_equation!(problem,
    "∂t(u) - nu*div(grad_u) + ∇(p) - buoy*T*ez + τ_lift(tau_u2) = -u⋅∇(u)")

# Boundary conditions.
# NOTE the `$Lz` interpolation: a BC string is parsed by Tarang, not by Julia,
# so it cannot see the Julia binding `Lz`. Writing "T(z=Lz) = 0" drops the BC
# and the first solve fails with a DimensionMismatch. Interpolate, or write the
# number literally.
add_bc!(problem, "T(z=0) = 1")        # hot bottom wall
add_bc!(problem, "T(z=$Lz) = 0")      # cold top wall
add_bc!(problem, "u(z=0) = 0")        # no-slip
add_bc!(problem, "u(z=$Lz) = 0")      # no-slip
add_bc!(problem, "integ(p) = 0")      # pressure gauge

solver = InitialValueSolver(problem, RK222(); dt=1e-3)

Pattern summary:

• The five add_bc! calls expand to 7 scalar constraint rows: two temperature rows, four velocity rows (two vector components at each of two walls), and one pressure-gauge row. The tau variables supply the matching seven scalar tau DOFs: one each from tau_T1 and tau_T2, four from the two-component tau_u1/tau_u2, and one from tau_p.
• Each tau field drops the coupled direction (xbasis only, not zbasis).
• tau_p is a 0-D scalar (no bases). It contributes a one-DOF candidate column during each subproblem build, but valid-mode filtering removes it at non-DC modes, leaving the actual gauge correction only at DC.
• The pressure gauge integ(p) = 0 is an algebraic constraint on the mean; it lives alongside the other BCs.

Pressure Gauge and Valid-Mode Filtering

In incompressible flow the pressure is defined only up to one global constant, so the gauge ambiguity belongs to the all-zero (DC) Fourier mode. At non-DC modes the pressure fluctuation is determined, while the domain-integral gauge equation is identically 0 = 0. Tarang handles the DC gauge and those trivial non-DC rows with two mechanisms that are sometimes confused:

1. Pressure gauge (a user-visible tau). You add tau_p to continuity and provide a gauge-fixing BC like integ(p) = 0:

add_equation!(problem, "trace(grad_u) + tau_p = 0")
add_bc!(problem, "integ(p) = 0")

tau_p is typically declared as a 0-D scalar (ScalarField(dist, "tau_p", (), Float64)). In the raw, pre-filtered matrix this gives it 1 DOF per subproblem (because a 0-D field has subproblem_field_size == 1 at every Fourier mode). After valid-mode filtering (see next section) tau_p is dropped from every non-DC subproblem — so in practice it's a single gauge-fixing scalar at the DC Fourier mode only. The non-DC subproblems don't carry it, and scatter_inputs on those modes is a no-op because the filter removed the column.

That's the picture you should hold in your head: tau_p fixes the pressure gauge at the DC mode; at every other mode, pressure is well-defined by the trace(grad_u) = 0 constraint together with the wall BCs, with no gauge ambiguity.

2. Valid-mode filtering (invisible, done by the solver). At the non-DC Fourier modes (and in some situations at the DC mode too), the integral constraint integ(p) = 0 produces a zero row in the raw LHS matrix — the integral of a non-DC Fourier mode over the domain is identically zero, so the row just says 0 = 0, which is trivially true and carries no information.

The solver detects these zero rows during matrix assembly (build_matrices!, in src/core/subsystems/subproblem_matrix_build.jl) and pairs each with a 1-DOF tau column — any variable whose per-subproblem size equals 1. That includes:

• tau_p (a 0-D gauge scalar, always vsz=1)
• tau_T1, tau_T2 (declared with (xbasis,) — at each Fourier mode, the xbasis contribution reduces to a single coefficient, so vsz=1 per subproblem)
• Components of tau_u1, tau_u2 (same story, per component)

The pairing uses a smallest-column-norm heuristic: for each zero row, the filter picks the unused 1-DOF tau column with the smallest total |L| + |M| norm — i.e., the tau that appears in the fewest or smallest-magnitude entries elsewhere. That's typically tau_p at non-DC subproblems because tau_p only contributes to continuity (one nonzero entry from +tau_p), while tau_T1/tau_T2 contribute through the full gradient substitution and so carry more weight. After pairing, both the row and the column are dropped from the filtered system, yielding a smaller square matrix that the sparse LU factorizes cleanly.

You do not have to do anything to enable this. Just declare tau_p and the integ(p) = 0 BC, and the valid-mode filter takes care of the rest. The take-away is: if you ever see a "singular pencils" or "Non-square filtered system" warning, it usually means a BC is missing, a tau field is missing, or the number of small-eq-size algebraic rows doesn't match the number of 1-DOF tau columns.

Number of Tau Terms

The number of scalar tau DOFs must match the number of scalar constraint rows (including gauge conditions). A vector-valued add_bc! call contributes one scalar row per component, and a vector tau field contributes the corresponding component-wise tau DOFs:

In a 2D problem (x-Fourier, z-coupled) each equation carries tau corrections only in the z direction — not "2 per direction". The Fourier x direction is periodic and needs no tau terms. Only the coupled z direction contributes.

Vector equations need vector tau fields. In RBC, u is a 2-component velocity, so tau_u1 and tau_u2 are VectorFields (one scalar per component per Fourier mode).

Time- and Space-Dependent BCs

Tarang supports boundary conditions whose value varies in time, space, or both. They are refreshed by the stepper:

• Time-dependent BCs (e.g. T(z=0) = sin(t)): the RK stepper re-evaluates the BC value at each stage time t + c[i]*dt, so multi-stage methods retain their full formal order of accuracy for rapidly-varying BCs.
• Space-dependent BCs (e.g. T(z=0) = sin(2*pi*x/4.0)): at solver-build time the BC expression is evaluated on the global coordinate grid, and the resulting array is projected onto the Fourier modes via an unnormalized FFTW.rfft. Each subproblem picks its own mode from the cached coefficient array.
• Space+time BCs (e.g. T(z=0) = sin(2*pi*x/4.0) * cos(2*pi*t)): combined — re-projected on every stage.

You don't need to register coordinate fields manually. The solver auto-registers global grid arrays for every Fourier/Chebyshev axis under its element label ("x", "y", "z", ...), so BC string expressions can reference those coordinates directly:

add_bc!(problem, "T(z=0) = 1 + 0.1*sin(2*pi*x/4.0)")   # Lx = 4.0, written literally
add_bc!(problem, "T(z=0) = sin(t)")                     # time-dependent

Both are enforced to machine precision — measured max|T(z=0) − target| of 2.2e-16 (space) and 2.3e-16 (time) after 5 RK222 steps of a diffusion problem.

Under MPI, every rank evaluates the BC expression on the full (global) grid and computes a local FFT — no inter-rank communication is needed because all ranks produce identical coefficient arrays.

Inspecting BC Satisfaction

After a solve, you can verify that the BCs are actually being enforced:

# After solve!()
ensure_layout!(u, :g)
u_data = get_grid_data(u)

# For Chebyshev-T, the Gauss-Lobatto grid includes the endpoints, and the grid
# runs from z_min to z_max.
# u_data[:, 1]    is the solution at z = z_min (first Cheb grid point)
# u_data[:, end]  is the solution at z = z_max (last Cheb grid point)
@assert maximum(abs.(u_data[:, 1]))   < 1e-10  "Bottom BC not satisfied"
@assert maximum(abs.(u_data[:, end])) < 1e-10  "Top BC not satisfied"

For a vector field, index a component: get_grid_data(u.components[2])[:, 1].

After a BVP solve!, the tau fields are scattered back with the corrections the solver chose, so you can read them to see how hard the tau method had to work. For smooth, well-resolved solutions they are spectrally small:

tau_data = get_coeff_data(tau2)          # ScalarField tau; for a VectorField use
@info "tau2 magnitude: $(maximum(abs.(tau_data)))"   # tau2.components[i]

In the quick-start problem the exact solution z(1 - z) is a quadratic that the Chebyshev basis represents exactly, so both taus come back as 0.0 — the tau correction has nothing to do. Swap the constant forcing for a non-polynomial one (a grid field holding sin(7z), passed in with add_parameters!) at Nz=16 and they become nonzero but tiny: max|tau1| = 2.1e-8, max|tau2| = 6.5e-8.

For time-dependent problems, check BC satisfaction inside a callback:

run!(solver;
     stop_iteration=30,
     callbacks=[on_interval(10) do s
         ensure_layout!(T, :g)
         T_data = get_grid_data(T)
         @info "T(z=0) residual: $(maximum(abs.(T_data[:, 1] .- 1.0)))"
     end])

Common Pitfalls

1. Wrong number of tau fields

If you have two wall BCs on u but only declare one tau field in the equation, the matrix system is under-determined. You'll see Warning: Matrix is not square: rows=70, cols=65 followed by Warning: Non-square filtered system: group=(3, nothing), valid_eqn=14, valid_var=13 at solver-build time (the exact counts depend on your resolution), and the solve then dies with a DimensionMismatch.

Fix: count the scalar constraint rows — including gauge conditions like integ(p) = 0, and counting one row per component for a vector BC — and declare one tau DOF per row.

2. Missing lift() term in the equation

A tau field declared in problem.variables but never referenced in any equation contributes a zero column to the LHS matrix, making the system rank-deficient.

Fix: every tau field must appear inside exactly one lift() (or τ_lift() substitution) on the implicit (left-hand) side of some equation.

3. Wrong lift basis — and what it does not cause

Writing lift(tau, zb, -1) (the state basis), or using the first-derivative basis where the second-derivative basis is meant, is a style problem, not a numerical one: as explained above, the current solver ignores op.basis and sizes the delta column from the problem's own Chebyshev basis, so these spellings build identical matrices and produce the same solution to the last bit. Nothing diverges at high Nz because of this choice — the lift basis is not a conditioning control.

Fix: still pass the basis that names the space you are lifting into — derivative_basis(zb, 1) for a first-order gradient substitution, derivative_basis(zb, 2) for a direct second-order equation. It says what you mean, matches every shipped example, and stays correct if Tarang later makes the lift basis load-bearing.

4. Tau field bases don't match

A tau field must live on the complement of the state field's bases — the state's bases minus the coupled direction. For a 2D state on (xbasis, zbasis), the tau lives on (xbasis,); for a 0-D gauge like tau_p it's ().

If you accidentally declare tau_u1 = ScalarField(dist, "tau_u1", (xbasis, zbasis)) (same bases as u), the lift call itself still constructs — the mistake is not caught there. It surfaces at solver build as Warning: Matrix is not square: rows=70, cols=180 (the tau contributes a whole Chebyshev spectrum of columns instead of one), and the solve then throws a DimensionMismatch.

Fix: drop the Chebyshev axis when declaring tau fields.

5. Non-square system at DC mode

Most non-square warnings at the DC Fourier mode come from a missing gauge BC like integ(p) = 0. The valid-mode filter can only drop a zero row if there's a paired 0-D tau column (like tau_p) to drop along with it.

Fix: make sure every PDE with a gauge ambiguity has both a tau_* field and a corresponding integ() BC.

6. BC F value not reaching the stepper

For inhomogeneous BCs like T(z=0) = 1, the stepper has to carry the value 1 through the stage RHS assembly. This happens automatically through gather_alg_F! + the apply_bc_override! path in step_subproblem_rk!, but only if the BC equation has a non-zero F expression in equation_data[eq_idx]["F"].

Symptom of the bug: max|T| decays to zero over time even though T(z=0) = 1 is declared, OR max|T| sticks at 1/γ = 2+√2 ≈ 3.414 (the classical 1/γ scaling factor from RK222's implicit coefficient).

If you see this: confirm you're on a recent version — this was a regression fixed after the subproblem-architecture rewrite. The apply_bc_override! override is what enforces L_row·X = F_BC at every stage regardless of accumulation history.

Historical Note

The tau method was introduced by Cornelius Lanczos in 1938 as an approximation technique: rather than solving a PDE exactly, he sought polynomial approximations that satisfied the PDE with a small residual (the "tau error"). It was refined for spectral methods by Steven Orszag, David Gottlieb, and others in the 1970s–80s, and adapted to modern lift-based formulations by Keaton Burns et al. in the Dedalus project in the 2010s.

The name "tau" (τ) comes from Lanczos's notation for the residual/correction terms introduced when truncating the polynomial expansion and enforcing boundary conditions.

References

Textbooks

1. Canuto, C., Hussaini, M. Y., Quarteroni, A., & Zang, T. A. (2006). Spectral Methods: Fundamentals in Single Domains. Springer. — Rigorous treatment of tau and Galerkin methods.

2. Boyd, J. P. (2001). Chebyshev and Fourier Spectral Methods (2nd ed.). Dover. — Very readable; freely available online.

3. Trefethen, L. N. (2000). Spectral Methods in MATLAB. SIAM. — Practical, code-oriented.

4. Peyret, R. (2002). Spectral Methods for Incompressible Viscous Flow. Springer. — Detailed treatment of Navier–Stokes with spectral methods.

Key papers

1. Lanczos, C. (1938). "Trigonometric interpolation of empirical and analytical functions." Journal of Mathematics and Physics, 17(1–4), 123–199. — Original tau method paper.

2. Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Brown, B. P. (2020). "Dedalus: A flexible framework for numerical simulations with spectral methods." Physical Review Research, 2, 023068. — Modern lift-based tau method as implemented in Tarang and Dedalus.

3. Orszag, S. A. (1971). "Accurate solution of the Orr-Sommerfeld stability equation." Journal of Fluid Mechanics, 50(4), 689–703. — Classic application to hydrodynamic stability.

See Also

• Boundary Conditions Tutorial: step-by-step BC examples
• Bases: spectral bases (Chebyshev, Fourier, Legendre, Jacobi)
• Solvers: using IVP / LBVP / NLBVP solvers
• API: Problems: programmatic API for adding equations and BCs
• 2D RBC Tutorial: complete Rayleigh–Bénard convection walkthrough