Skip to main content
When entering a system of PDEs in Dedalus, the left-hand side (LHS) of the equations is parsed into a sparse linear system. For the solver to succeed, this LHS matrix must be square and nonsingular - it must have a unique solution for any possible values on the right-hand side (RHS).

​
The Problem

The linear system must constrain all degrees of freedom of the variables, including gauge freedoms. Failure to do so results in singular matrices that cannot be inverted.

​
Example: Incompressible Flow in a Periodic Domain

Consider incompressible hydrodynamics in a fully periodic domain discretized with Fourier bases: 
import dedalus.public as d3
import numpy as np

coords = d3.CartesianCoordinates('x', 'y', 'z')
dist = d3.Distributor(coords, dtype=np.float64)
bases = (d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx)),
         d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly)),
         d3.RealFourier(coords['z'], size=Nz, bounds=(0, Lz)))

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

# Problem - INCORRECT: This will produce singular matrices!
problem = d3.IVP([p, u], namespace=locals())
problem.add_equation("div(u) = 0")
problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
This formulation produces square matrices but they are not nonsingular. The pressure gauge is undetermined - any constant can be added to the pressure without affecting the equations.

​
Understanding the Degeneracy

Two related issues exist:
  1. Underdetermined variable: Any constant can be added to the pressure
  2. Degenerate constraint: For the mean k⃗=0\vec{k} = 0 Fourier mode, the divergence equation becomes "0 = 0"
While this seems consistent, the system must be solvable for any RHS, not just the one we entered. If we entered "div(u) = 1", the mean mode would give "0 = 1" - clearly impossible.

​
The Solution: Gauge Conditions with Tau Variables

We expand the system by adding a spatially-constant variable (a tau variable) to the divergence equation to absorb the degeneracy, and impose the pressure gauge as a separate equation: 
# Fields
p = dist.Field(name='p', bases=bases)
u = dist.VectorField(coords, name='u', bases=bases)
tau_p = dist.Field(name='tau_p')  # Constant scalar field

# Problem - CORRECT
problem = d3.IVP([p, u, tau_p], namespace=locals())
problem.add_equation("div(u) + tau_p = 0")  # Tau absorbs degeneracy
problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
problem.add_equation("integ(p) = 0")  # Fix pressure gauge
We’ve added one degree of freedom (tau_p) and one constraint (integ(p) = 0), so the system is still square. Now it’s also nonsingular since the mean pressure is fixed.

​
How It Works

The tau variable mechanism works elegantly:
  • The mean pressure is fixed by the gauge condition
  • The degeneracy is lifted - the tau variable absorbs/acquires the mean value of any possible RHS
  • The mean divergence of velocity is always zero, as required by the periodic discretization

​
Gauge Conditions in Bounded Domains

Similar modifications work for other types of gauges and geometries. For incompressible hydrodynamics in a bounded domain, we still need this modification so the integral of the divergence equation is compatible with specified inflow boundary conditions. 
coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.float64)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz))

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

# Substitutions for tau method
ez = coords.unit_vector_fields(dist)[1]
lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
grad_u = d3.grad(u) + ez*lift(tau_u1)

# Problem
problem = d3.IVP([p, u, tau_p, tau_u1, tau_u2], namespace=locals())
problem.add_equation("trace(grad_u) + tau_p = 0")  # Tau for compatibility
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + lift(tau_u2) = - u@grad(u)")
problem.add_equation("u(z=0) = 0")
problem.add_equation("u(z=Lz) = 0")
problem.add_equation("integ(p) = 0")  # Pressure gauge

​
Interpretation with Inflow

From the modified equation, we can see:
  • If the prescribed net inflow is nonzero, the tau variable acquires a corresponding nonzero value
  • The velocity then has a spatially uniform convergence equal to this tau value
  • For properly specified boundary conditions with no net inflow, the tau variable is zero and velocity is divergence-free

​
General Pattern

The gauge condition pattern applies to various situations:
1
Identify the gauge freedom
2
Determine which variables have underdetermined components (e.g., mean pressure, total angular momentum)
3
Find the degenerate constraint
4
Identify which equation becomes "0 = 0" or is otherwise degenerate
5
Add a tau variable
6
Introduce a constant scalar field to absorb the degeneracy in the affected equation
7
Add the gauge equation
8
Specify an additional constraint that fixes the gauge (e.g., integ(p) = 0, integ(L) = L0)

​
Common Gauge Conditions

tau_p = dist.Field(name='tau_p')
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0")
Fixes the mean pressure to zero.
tau_T = dist.Field(name='tau_T')
problem.add_equation("dt(T) - kappa*lap(T) + tau_T = source")
problem.add_equation("integ(T) = T_mean")
Maintains a specified mean temperature when there’s no natural reference.
tau_L = dist.Field(name='tau_L')
problem.add_equation("dt(L) - div(stress) + tau_L = 0")
problem.add_equation("integ(L) = L0")
Fixes total angular momentum in a rotating reference frame.

​
Key Takeaways

Always Check Matrix Properties

The LHS matrix must be square and nonsingular - gauge freedoms break singularity even if dimensions match.

Use Tau Variables

Constant tau fields elegantly absorb constraint degeneracies while maintaining system structure.

Fix All Gauges

Every gauge freedom requires both a tau variable to lift degeneracy and an equation to fix the gauge.

Check Compatibility

Tau variables ensure integral constraints are compatible with boundary conditions and forcing.

​
See Also

  • Tau Method - Understanding tau variables for boundary conditions
  • Example scripts in the Dedalus repository demonstrate gauge conditions in various domains

Build docs developers (and LLMs) love

Get started for freeTalk to us