Tutorial: Boundary Conditions

This tutorial covers the different types of boundary conditions available in Tarang.jl and how to apply them effectively using the tau method with explicit tau fields.

Overview

Boundary conditions (BCs) are essential for well-posed PDE problems. Tarang.jl supports:

• Dirichlet: Specify field value at boundary
• Neumann: Specify derivative at boundary
• Robin: Linear combination of value and derivative
• Periodic: Automatic for Fourier bases

The Tau Method Approach

Tarang.jl uses an explicit tau-method approach for handling boundary conditions. Users must explicitly create tau fields and add them to equations using the lift() operator.

Why Explicit Tau Fields?

1. Clarity: The mathematical structure is visible in your code
2. Flexibility: Full control over tau placement
3. Debugging: Easy to inspect tau field values
4. Consistency: Matches the mathematical formulation

Required Steps

For any problem with non-periodic boundary conditions:

1. Create tau fields - One scalar tau DOF per scalar constraint row (vector BCs contribute one row per component)
2. Add tau fields to the problem - Include them in the field list
3. Add lift() terms to equations - Place tau contributions at specific modes
4. Specify boundary conditions - Add the algebraic constraint rows; tau columns and BC rows couple through the assembled matrix, not by a name-based link

Complete Example: Poisson Equation

Let's solve the Poisson equation $\nabla^2 u = -2$ with Dirichlet BCs u = 0 on both z walls. This is a verified, runnable LBVP; the analytic solution is u = z(Lz - z) and the solver reproduces it to machine precision.

The tau method needs one tau variable per boundary condition. Each tau is lifted into the bulk equation with lift(tau, derivative_basis(z_basis, 2), -k) and registered as a parameter via add_parameters!. The lift order index -1 targets the last coefficient row and -2 the second-to-last; a 2nd-order problem therefore uses two taus at orders -1 and -2.

using Tarang

# Create coordinates, distributor, and bases
coords = CartesianCoordinates("x", "z")
dist   = Distributor(coords; dtype=Float64, device=CPU())
x_basis = RealFourier(coords["x"]; size=4,  bounds=(0.0, 2π))   # periodic (separable)
z_basis = ChebyshevT(coords["z"];  size=16, bounds=(0.0, 1.0))  # bounded  (coupled)
dom = Domain(dist, (x_basis, z_basis))

# Field to solve for
u = ScalarField(dom, "u")

# Step 1: Create tau fields (one per BC)
# These live on the x-basis only (boundary is a line in 2D)
tau1 = ScalarField(dist, "tau1", (x_basis,), Float64)  # For BC at z=0
tau2 = ScalarField(dist, "tau2", (x_basis,), Float64)  # For BC at z=1

# Step 2: Add ALL fields to problem (including tau fields)
problem = LBVP([u, tau1, tau2])

# Step 3: Register the lifted tau terms as parameters
#         (use the 2nd-derivative basis for a 2nd-order problem)
lb2 = derivative_basis(z_basis, 2)
add_parameters!(problem; l1=lift(tau1, lb2, -1), l2=lift(tau2, lb2, -2))

# Step 4: Add equation referencing the lifted tau parameters
add_equation!(problem, "Δ(u) + l1 + l2 = -2")

# Step 5: Add boundary conditions
add_bc!(problem, "u(z=0)   = 0")
add_bc!(problem, "u(z=1.0) = 0")

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

Dirichlet Boundary Conditions

Fix the value of a field at the boundary.

Basic Setup

# Create tau fields for each Dirichlet BC
tau_T1 = ScalarField(dist, "tau_T1", (x_basis,), Float64)
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,), Float64)

# Add to problem
problem = LBVP([T, tau_T1, tau_T2])

# Register the lifted tau terms (2nd-order problem -> 2nd-derivative basis)
lb2 = derivative_basis(z_basis, 2)
add_parameters!(problem; l1=lift(tau_T1, lb2, -1), l2=lift(tau_T2, lb2, -2))

# Add equation referencing the lifted tau parameters
add_equation!(problem, "Δ(T) + l1 + l2 = source")

# Boundary conditions
add_bc!(problem, "T(z=0) = 1")   # T(z=0) = 1
add_bc!(problem, "T(z=1) = 0")   # T(z=1) = 0

No-Slip Velocity (Vector Fields)

For viscous flows at solid walls, use vector fields for compact notation:

# Vector field for velocity
u = VectorField(dist, coords, "u", (x_basis, z_basis))
p = ScalarField(dist, "p", (x_basis, z_basis))

# First-order tau fields: one in the gradient substitution, one in momentum
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,))
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,))
tau_p = ScalarField(dist, "tau_p", ())

# First-order viscous operator
ex, ez = unit_vector_fields(coords, dist)
lift_basis = derivative_basis(z_basis, 1)
τ_lift(A) = lift(A, lift_basis, -1)
grad_u = grad(u) + ez * τ_lift(tau_u1)

# Pass vector fields directly to problem
problem = IVP([u, p, tau_u1, tau_u2, tau_p])

# Add parameters and substitutions
add_parameters!(problem, nu=nu, grad_u=grad_u, τ_lift=τ_lift)

# Momentum equation (single vector equation)
add_equation!(problem, "∂t(u) - nu*div(grad_u) + ∇(p) + τ_lift(tau_u2) = -u⋅∇(u)")

# Continuity with tau_p (removes degeneracy)
add_equation!(problem, "trace(grad_u) + tau_p = 0")

# No-slip boundary conditions (vector notation)
add_bc!(problem, "u(z=0) = 0")   # No-slip bottom (all components)
add_bc!(problem, "u(z=1) = 0")   # No-slip top (all components)
add_bc!(problem, "integ(p) = 0")  # Fix the global pressure constant

Neumann Boundary Conditions

Specify the derivative (flux) at the boundary.

Basic Setup

# Tau fields for Neumann BCs
tau_T1 = ScalarField(dist, "tau_T1", (x_basis,), Float64)
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,), Float64)

problem = LBVP([T, tau_T1, tau_T2])

# Register the lifted tau terms (2nd-order problem -> 2nd-derivative basis)
lb2 = derivative_basis(z_basis, 2)
add_parameters!(problem; l1=lift(tau_T1, lb2, -1), l2=lift(tau_T2, lb2, -2))

add_equation!(problem, "Δ(T) + l1 + l2 = source")

# Neumann: specify derivative at boundary
add_bc!(problem, "∂z(T)(z=0) = 1")   # ∂T/∂z(z=0) = 1
add_bc!(problem, "∂z(T)(z=1) = 0")   # ∂T/∂z(z=1) = 0

Stress-Free Conditions

For free surfaces or slip boundaries (∂u/∂z = 0):

# Mixed: no-slip at bottom, stress-free at top
add_bc!(problem, "u_x(z=0) = 0")      # No-slip at bottom
add_bc!(problem, "∂z(u_x)(z=1) = 0")  # Stress-free at top (∂u/∂z = 0)

Robin Boundary Conditions

Linear combination: $\alpha u + \beta \frac{\partial u}{\partial n} = \gamma$

tau_T1 = ScalarField(dist, "tau_T1", (x_basis,), Float64)
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,), Float64)

problem = LBVP([T, tau_T1, tau_T2])

# Register the lifted tau terms and scalar parameters
lb2 = derivative_basis(z_basis, 2)
add_parameters!(problem; l1=lift(tau_T1, lb2, -1), l2=lift(tau_T2, lb2, -2),
                         h=10.0, k=1.0, T_amb=25.0)

add_equation!(problem, "Δ(T) + l1 + l2 = source")

# Convective heat transfer at top: h*T + k*dT/dn = h*T_ambient
add_bc!(problem, "h*T(z=1) + k*∂z(T)(z=1) = h*T_amb")

# Dirichlet at bottom
add_bc!(problem, "T(z=0) = 100")

Periodic Boundary Conditions

Periodic boundaries are automatically handled by Fourier bases - no tau fields needed!

# Fourier basis implies periodicity
x_basis = RealFourier(coords["x"], size=128, bounds=(0.0, 2π))

# No boundary conditions or tau fields needed for x-direction
# The Fourier representation automatically enforces u(x=0) = u(x=2π)

The lift() Operator

The lift(tau, derivative_basis(basis, order), -k) operator places a tau correction at a specific coefficient row of the bounded (Chebyshev/Jacobi) basis. You build the lifted terms once and register them as parameters, then reference those parameter names in the equation string:

lb2 = derivative_basis(z_basis, 2)   # 2nd-order problem -> 2nd-derivative basis
add_parameters!(problem; l1=lift(tau_u1, lb2, -1), l2=lift(tau_u2, lb2, -2))
add_equation!(problem, "Δ(u) + l1 + l2 = f")

The order index -1 lifts into the last coefficient row, -2 the second-to-last, and so on. The number of taus (and lift terms) must match the operator order in the bounded direction.

Mode Selection Guidelines

Example: Fourth-Order Problem

For a biharmonic equation (∇⁴u = f) with 4 boundary conditions:

# Four tau fields needed
tau_u1 = ScalarField(dist, "tau_u1", (x_basis,), Float64)
tau_u2 = ScalarField(dist, "tau_u2", (x_basis,), Float64)
tau_u3 = ScalarField(dist, "tau_u3", (x_basis,), Float64)
tau_u4 = ScalarField(dist, "tau_u4", (x_basis,), Float64)

problem = LBVP([u, tau_u1, tau_u2, tau_u3, tau_u4])

# Register the lifted tau terms (4th-order problem -> 4th-derivative basis)
lb4 = derivative_basis(z_basis, 4)
add_parameters!(problem; l1=lift(tau_u1, lb4, -1), l2=lift(tau_u2, lb4, -2),
                         l3=lift(tau_u3, lb4, -3), l4=lift(tau_u4, lb4, -4))

# Biharmonic equation with all four lifted tau terms
add_equation!(problem, "Δ(Δ(u)) + l1 + l2 + l3 + l4 = f")

# Clamped beam: u = 0 and du/dz = 0 at both ends
add_bc!(problem, "u(z=0) = 0")
add_bc!(problem, "∂z(u)(z=0) = 0")
add_bc!(problem, "u(z=1) = 0")
add_bc!(problem, "∂z(u)(z=1) = 0")

Complete Examples

Channel Flow (Poiseuille Flow)

using Tarang

coords = CartesianCoordinates("x", "z")
x_basis = RealFourier(coords["x"]; size=64, bounds=(0.0, 2π))
z_basis = ChebyshevT(coords["z"]; size=64, bounds=(0.0, 1.0))
dist = Distributor(coords; dtype=Float64, device=CPU())
domain = Domain(dist, (x_basis, z_basis))

# Vector velocity field
u = VectorField(domain, "u")
p = ScalarField(domain, "p")

# First-order tau fields
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,), Float64)
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,), Float64)

# Tau for pressure (removes degeneracy)
tau_p = ScalarField(dist, "tau_p", (), Float64)

# Parameters
nu = 0.01
dpdx = -1.0

# First-order viscous operator
ex, ez = unit_vector_fields(coords, dist)
lift_basis = derivative_basis(z_basis, 1)
τ_lift(A) = lift(A, lift_basis, -1)
grad_u = grad(u) + ez * τ_lift(tau_u1)

# Create problem with all fields and register substitutions
problem = IVP([u, p, tau_u1, tau_u2, tau_p])
add_parameters!(problem, nu=nu, dpdx=dpdx, ex=ex,
                         grad_u=grad_u, τ_lift=τ_lift)

# Momentum equation (vector form) - dpdx is the driving pressure gradient
add_equation!(problem, "∂t(u) - nu*div(grad_u) + ∇(p) + τ_lift(tau_u2) = -u⋅∇(u) - dpdx*ex")

# Continuity with tau_p (removes degeneracy)
add_equation!(problem, "trace(grad_u) + tau_p = 0")

# No-slip at both walls (vector notation)
add_bc!(problem, "u(z=0) = 0")
add_bc!(problem, "u(z=1) = 0")
add_bc!(problem, "integ(p) = 0")

Rayleigh-Bénard Convection

using Tarang

coords = CartesianCoordinates("x", "z")
x_basis = RealFourier(coords["x"]; size=128, bounds=(0.0, 4.0))
z_basis = ChebyshevT(coords["z"]; size=64, bounds=(0.0, 1.0))
dist = Distributor(coords; dtype=Float64, device=CPU())
domain = Domain(dist, (x_basis, z_basis))

# Vector velocity field and scalar fields
u = VectorField(domain, "u")
p = ScalarField(domain, "p")
T = ScalarField(domain, "T")

# Vector tau fields for velocity BCs
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,), Float64)
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,), Float64)

# Scalar tau fields for temperature BCs
tau_T1 = ScalarField(dist, "tau_T1", (x_basis,), Float64)
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,), Float64)

# Tau for pressure (removes degeneracy)
tau_p = ScalarField(dist, "tau_p", (), Float64)

# Parameters
Ra = 1e6   # Rayleigh number
Pr = 1.0   # Prandtl number

# First-order substitutions
ex, ez = unit_vector_fields(coords, dist)
lift_basis = derivative_basis(z_basis, 1)
τ_lift(A) = lift(A, lift_basis, -1)
grad_u = grad(u) + ez * τ_lift(tau_u1)
grad_T = grad(T) + ez * τ_lift(tau_T1)

# Create problem with all fields and register substitutions
problem = IVP([p, T, u, tau_p, tau_T1, tau_T2, tau_u1, tau_u2])
add_parameters!(problem, nu=Pr, buoy=Ra*Pr, ez=ez,
                         grad_u=grad_u, grad_T=grad_T, τ_lift=τ_lift)

# Momentum equation (vector form with buoyancy)
add_equation!(problem, "∂t(u) - nu*div(grad_u) + ∇(p) - buoy*T*ez + τ_lift(tau_u2) = -u⋅∇(u)")

# Continuity with tau_p (removes degeneracy)
add_equation!(problem, "trace(grad_u) + tau_p = 0")

# Temperature equation
add_equation!(problem, "∂t(T) - div(grad_T) + τ_lift(tau_T2) = -u⋅∇(T)")

# Boundary conditions (vector notation for velocity)
add_bc!(problem, "u(z=0) = 0")   # No-slip bottom
add_bc!(problem, "u(z=1) = 0")   # No-slip top

# Fixed temperature
add_bc!(problem, "T(z=0) = 1")   # Hot bottom
add_bc!(problem, "T(z=1) = 0")   # Cold top
add_bc!(problem, "integ(p) = 0")  # Pressure gauge

Time- and Space-Dependent Boundary Conditions

Tarang supports BCs whose value varies in time, space, or both. The BC expression string is re-evaluated at solver build time (for space-dependent) and at every time step / RK stage (for time-dependent), and the result is projected onto the appropriate Fourier mode for each per-subproblem solve.

Nothing extra is required — no add_coordinate_field! call, no set_time_variable! — the solver auto-registers coordinate arrays from the problem's bases and uses t as the default time variable.

Time-dependent BC (scalar)

add_bc!(problem, "T(z=0) = 1.0 + 0.1 * cos(2*pi*t)")

The BC value is re-evaluated at every RK stage time t + c[i]·dt, so multi-stage methods retain full formal order of accuracy even for rapidly-varying BCs. Internally this goes through update_time_dependent_bcs!(bc_manager, stage_time) followed by _apply_bc_values_to_equations!, and the resulting ConstantOperator(value) is written into equation_data[eq_idx]["F"].

Space-dependent BC (1D variation in 2D problem)

add_bc!(problem, "T(z=0) = 1.0 + 0.1 * sin(2*pi*x/4.0)")
add_bc!(problem, "T(z=Lz) = 0")

At solver build, _apply_bc_values_to_equations!(solver, 0.0) evaluates "1.0 + 0.1*sin(2*pi*x/4.0)" against the auto-registered global x coordinate array, yielding an Nx-long grid-space array. This is wrapped in an ArrayOperator and stored in equation_data[eq_idx]["F"]. At each stepper call, gather_alg_F! runs _bc_array_projection on this array — taking an unnormalized FFTW.rfft and extracting each subproblem's own Fourier-mode coefficient.

The bottom-wall temperature is thus enforced at T(x, z=0) = 1 + 0.1·sin(2πx/Lx) in grid space.

Space + time dependent BC

add_bc!(problem, "T(z=0) = 1.0 + 0.1 * sin(2*pi*x/4.0) * cos(2*pi*t)")

Combines both paths. The spatial pattern is re-projected at each stage time, so the stepper always sees the correct instantaneous BC. This is what enables oscillating boundary temperature patterns, traveling thermal waves at the wall, etc.

Two-axis spatial BCs in 3D

For a 3D problem with two periodic axes (x, y):

# 3D problem: x and y periodic, z coupled
add_bc!(problem, "T(z=0) = sin(2*pi*x/Lx) * cos(2*pi*y/Ly)")

The array is evaluated as a 2D (Nx, Ny) grid-space array and transformed via 2D FFT. Each subproblem (one per (kx, ky) mode) extracts its (kx_global, ky_global) coefficient. This works out-of-the-box — no extra configuration required.

1D BC expression in a higher-dimensional problem

If you write a 1D expression like sin(2*pi*x/Lx) in a 3D problem with two periodic axes, Tarang auto-broadcasts the 1D array along the missing axis — so the BC is treated as "uniform in y":

# 3D problem; this is constant in y by design
add_bc!(problem, "T(z=0) = sin(2*pi*x/Lx)")

The broadcast-and-FFT machinery places rfft(sin(2π x/Lx))[k_x] · Ny at the DC y mode and zero at non-DC y modes, matching the grid-space meaning exactly.

User parameters in BC expressions

Expression strings can reference user parameters registered via add_parameters!:

add_parameters!(problem, Lx=4.0, Ly=2.0, amplitude=0.1)
add_bc!(problem, "T(z=0) = 1.0 + amplitude * sin(2*pi*x/Lx)")

Or, more directly, bake the constant into the string via Julia interpolation:

Lx = 4.0
amplitude = 0.1
add_bc!(problem, "T(z=0) = 1.0 + $amplitude * sin(2*pi*x/$Lx)")

Both patterns work. The first is cleaner if you plan to sweep parameters; the second is simpler for one-off setups.

Common pitfalls

• Missing coordinate field — for a BC to reference x, the problem must have a basis with element_label == "x". This is automatic if you build bases with coords["x"], but if you have multiple coordinate systems the auto-registration may pick the wrong x. A "space-dependent BCs detected but no coordinate fields registered" warning at solver build means the auto-registration didn't find a matching basis.
• Time variable other than t — the parser and the runtime evaluator both hard-code t as the time symbol. Custom time variable names aren't fully supported; use t in your BC expressions.
• Non-constant custom function — the BC string evaluator whitelists sin, cos, tan, exp, log, sqrt, abs, sign, floor, ceil, round, rem, min, max, and hyperbolic variants. Arbitrary Julia functions in BC strings are not supported; if you need a custom function, build it into the value at simulation-script level (outside the string) and register via add_parameters!.

Validation

Checking BC Satisfaction

# After solving, verify BCs are satisfied
function check_dirichlet_bc(field, coord, location, expected_value; tol=1e-10)
    to_grid!(field)
    data = get_grid_data(field)

    # Get boundary values
    if coord == "z" && location == :left
        bc_values = data[:, 1]
    elseif coord == "z" && location == :right
        bc_values = data[:, end]
    end

    error = maximum(abs.(bc_values .- expected_value))
    @assert error < tol "BC error: $error"

    return error
end

# Usage
check_dirichlet_bc(T, "z", :left, 1.0)
check_dirichlet_bc(T, "z", :right, 0.0)

Troubleshooting

Missing or Unused Tau Fields

There is no separate automatic name-based pairing between a BC and a tau field. A missing lift column, an unused tau variable, or the wrong number of scalar tau DOFs normally appears during matrix construction or factorization as a non-square or rank-deficient system.

Solution: Create the required tau fields, include every one in the problem variables, and reference every one through a lift term:

# Wrong: No tau fields or lift terms
add_equation!(problem, "Δ(u) = f")
add_bc!(problem, "u(z=0) = 0")

# Correct: Create tau fields, register lifted tau parameters, add them to the equation
tau_u1 = ScalarField(dist, "tau_u1", (x_basis,), Float64)
tau_u2 = ScalarField(dist, "tau_u2", (x_basis,), Float64)
lb2 = derivative_basis(z_basis, 2)
add_parameters!(problem; l1=lift(tau_u1, lb2, -1), l2=lift(tau_u2, lb2, -2))
add_equation!(problem, "Δ(u) + l1 + l2 = f")
add_bc!(problem, "u(z=0) = 0")
add_bc!(problem, "u(z=1) = 0")

Over-Specified System

Problem: Too many boundary conditions cause singular matrices.

Solution: Match the number of BCs (and tau fields) to the operator order in each direction.

Under-Specified System

Problem: Not enough BCs lead to non-unique solutions.

Solution: Add appropriate boundary conditions for the problem physics.

Tau Field Dimension Mismatch

Problem: Tau field has wrong dimensionality.

Solution: Tau fields should live on the boundary, not the full domain:

# For a 2D problem with x-basis (periodic) and z-basis (non-periodic):
# Boundaries are at constant z, so tau fields live on x-basis only
tau_u = ScalarField(dist, "tau_u", (x_basis,), Float64)  # Correct: 1D on x

# NOT:
tau_u = ScalarField(dist, "tau_u", (x_basis, z_basis), Float64)  # Wrong: 2D

See Also

• Tau Method (Advanced): Mathematical details
• Problems API: Problem definition
• Bases API: Spectral basis selection
• 2D Rayleigh-Bénard Tutorial: Complete example with BCs