Skip to main content
The generalized tau method is a system for imposing boundary conditions when solving partial differential equations using polynomial spectral methods. It consists of explicitly adding tau terms to the PDE which introduce degrees of freedom that allow the problem to be solved exactly over polynomials.

​
Overview

Most PDEs we wish to solve do not have exact polynomial solutions. Instead, we seek polynomial solutions that approximate the true solution. The tau method makes these modifications explicit in the problem specification by adding tau terms to the equations.

​
Basic Example

Consider the simple equation: ∂xu(x)−u(x)=0,x∈[0,1]u(0)=1\begin{gathered} \partial_x u(x) - u(x) = 0, \quad x \in [0, 1] \\ u(0) = 1 \end{gathered} The exact solution is u(x)=exu(x) = e^x, but we seek an approximate polynomial solution. The generalized tau method modifies the PDE to be: ∂xu(x)−u(x)+τP(x)=0\partial_x u(x) - u(x) + \tau P(x) = 0 where τ\tau is an undetermined constant and P(x)P(x) is a specified polynomial.
If P(x)P(x) is a polynomial of degree NN, then the modified equation has an exact polynomial solution uN(x)u_N(x) that is also of degree NN.

​
First-Order Form Systems

When solving a system of nonsingular PDEs, the number of tau terms and boundary conditions will generally match the total number of derivatives in the system. This is most easily counted by converting the system to first-order form.

​
Example: 2D Incompressible Hydrodynamics

Consider linearized 2D incompressible hydrodynamics with velocity u⃗=(u,v)\vec{u} = (u, v), pressure pp, and forcing f⃗=(f,g)\vec{f} = (f, g). The domain is periodic in xx and bounded in y∈[−1,1]y \in [-1, 1] with no-slip conditions. The equations in first-order form for yy-derivatives: uy−∂yu+τ1(x)P(y)=0vy−∂yv+τ2(x)P(y)=0∂xu+vy=0∂tu−ν(∂x2u+∂yuy)+∂xp+τ3(x)P(y)=f∂tv−ν(∂x2v+∂yvy)+∂yp+τ4(x)P(y)=g\begin{gathered} u_y - \partial_y u + \tau_1(x) P(y) = 0 \\ v_y - \partial_y v + \tau_2(x) P(y) = 0 \\ \partial_x u + v_y = 0 \\ \partial_t u - \nu (\partial_x^2 u + \partial_y u_y) + \partial_x p + \tau_3(x) P(y) = f \\ \partial_t v - \nu (\partial_x^2 v + \partial_y v_y) + \partial_y p + \tau_4(x) P(y) = g \end{gathered}
We have four yy-derivatives and four boundary conditions, so we need four tau terms.

​
Implementation in Dedalus v3

In Dedalus v3, equations are entered in vectorial form using substitutions: 
import dedalus.public as d3
import numpy as np

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=np.float64)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.ChebyshevT(coords['y'], size=Ny, bounds=(-1, 1))

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

# Substitutions
ex, ey = coords.unit_vector_fields(dist)
lift_basis = ybasis.derivative_basis(1)  # Chebyshev U basis
lift = lambda A: d3.Lift(A, lift_basis, -1)  # Multiply by U_{N-1}(y)
grad_u = d3.grad(u) - ey*lift(tau_u1)  # First-order reduction

# Problem
problem = d3.IVP([p, u, tau_u1, tau_u2, tau_p], namespace=locals())
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + lift(tau_u2) = f")
problem.add_equation("u(y=-1) = 0")
problem.add_equation("u(y=+1) = 0")
problem.add_equation("integ(p) = 0")  # Pressure gauge

​
Higher-Order Systems

Dedalus v3 supports equations with arbitrary differential order. We can eliminate the first-order-reduction variables to recover the original second-order equations, but containing the same tau corrections: ∂xu+∂yv−τ2(x)P(y)=0∂tu−ν(∂x2u+∂y2u)+∂xp+τ3(x)P(y)−ντ1(x)∂yP(y)=f∂tv−ν(∂x2v+∂y2v)+∂yp+τ4(x)P(y)−ντ2(x)∂yP(y)=g\begin{gathered} \partial_x u + \partial_y v - \tau_2(x) P(y) = 0 \\ \partial_t u - \nu (\partial_x^2 u + \partial_y^2 u) + \partial_x p + \tau_3(x) P(y) - \nu \tau_1(x) \partial_y P(y) = f \\ \partial_t v - \nu (\partial_x^2 v + \partial_y^2 v) + \partial_y p + \tau_4(x) P(y) - \nu \tau_2(x) \partial_y P(y) = g \end{gathered} This system has the same solution as the first-order system, but is more efficient to solve.

​
The Lift Operator

The Lift operator handles the specification of and multiplication by the tau polynomial P(y)P(y). It multiplies its argument by the specified mode/element of a selected basis. 
# For Cartesian domains with Chebyshev bases
lift_basis = ybasis.derivative_basis(1)  # Chebyshev U basis
lift = lambda A: d3.Lift(A, lift_basis, -1)  # Use highest mode U_{N-1}(y)
The choice of tau polynomial affects the conditioning and accuracy of the solution. For Chebyshev bases, using the highest mode in the Chebyshev-U basis (via the derivative basis) is recommended.

​
Disks and Balls

In disk and ball geometries, the radial dimension has only a single (outer) boundary. Second-order elliptic and parabolic equations generally only need one boundary condition and one tau term: 
# Fields for disk problem
p = dist.Field(name='p', bases=disk_basis)
u = dist.VectorField(coords, name='u', bases=disk_basis)
tau_u = dist.VectorField(coords, name='tau_u', bases=phi_basis)  # Only one tau term
tau_p = dist.Field(name='tau_p')

# Substitutions
lift = lambda A: d3.Lift(A, disk_basis, -1)

# Problem
problem = d3.IVP([p, u, tau_u, tau_p], namespace=locals())
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = f")
problem.add_equation("u(r=1) = 0")
problem.add_equation("integ(p) = 0")
For disk and ball bases, the Lift operator automatically handles the mode-dependent tau polynomials under the hood.

​
Summary

Key points for tau formulations in Dedalus v3:
  1. Create tau fields that are supported on the boundary for your problem formulation
  2. Match tau fields to boundary conditions - you need the same number and type of tau fields as boundary conditions
  3. For Cartesian/annuli/shells: Use a first-order-style implementation with first-order substitutions that include tau terms
  4. For disks/balls: Only a single tau term is needed for second-order elliptic/parabolic problems
  5. Use the Lift operator to apply tau polynomials correctly for each geometry

​
See Also

  • Gauge Conditions - Handling pressure gauges and constraint degeneracies
  • Example scripts in the Dedalus repository demonstrate tau modifications in various domains

Build docs developers (and LLMs) love

Get started for freeTalk to us