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?

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

Required Steps

For any problem with non-periodic boundary conditions:

Create tau fields - One per boundary condition
Add tau fields to the problem - Include them in the field list
Add lift() terms to equations - Place tau contributions at specific modes
Specify boundary conditions - Link each BC to its tau field

Complete Example: Poisson Equation

Let's solve the Poisson equation $\nabla^2 u = f$ with Dirichlet BCs:

using Tarang

# Create coordinates and bases
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))

# Create distributor and fields
dist = Distributor(coords)
u = ScalarField(dist, "u", (x_basis, z_basis))
f = ScalarField(dist, "f", (x_basis, z_basis))  # Source term

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

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

# Step 3: Add source term parameter
add_parameters!(problem, f=f)

# Step 4: Add equation with lift() operators (string format)
add_equation!(problem, "Δ(u) + lift(tau_u1, -1) + lift(tau_u2, -2) = f")

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

# Solve
solver = BoundaryValueSolver(problem)
solve!(solver)

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,))
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,))

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

# Add equation with lift terms (string format)
add_equation!(problem, "Δ(T) + lift(tau_T1, -1) + lift(tau_T2, -2) = 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))

# Vector tau fields for BCs at each wall
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,))  # Wall at z=0
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,))  # Wall at z=1
tau_p = ScalarField(dist, "tau_p", ())

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

# Add parameters
add_parameters!(problem, nu=nu)

# Momentum equation (single vector equation)
add_equation!(problem, "∂t(u) - nu*Δ(u) + ∇(p) + lift(tau_u2, -2) = -u⋅∇(u)")

# Continuity with tau_p (removes degeneracy)
add_equation!(problem, "div(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)

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,))
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,))

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

add_equation!(problem, "Δ(T) + lift(tau_T1, -1) + lift(tau_T2, -2) = 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,))
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,))

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

# Add parameters
add_parameters!(problem, h=10.0, k=1.0, T_amb=25.0)

add_equation!(problem, "Δ(T) + lift(tau_T1, -1) + lift(tau_T2, -2) = 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π)

Mixing Periodic and Non-Periodic

When using Fourier (periodic) and Chebyshev (non-periodic) bases together, only create tau fields and boundary conditions for the non-periodic directions.

The lift() Operator

The lift() operator places tau corrections at specific spectral modes. In the string equation format:

"Δ(u) + lift(tau_u1, -1) + lift(tau_u2, -2) = f"

The tau field name in the lift() operator should match the tau field name you created.

Mode Selection Guidelines

Operator Order	Number of BCs	Number of lift() terms
1st (∂/∂z)	1	1
2nd (∂²/∂z²)	2	2
4th (∇⁴)	4	4

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,))
tau_u2 = ScalarField(dist, "tau_u2", (x_basis,))
tau_u3 = ScalarField(dist, "tau_u3", (x_basis,))
tau_u4 = ScalarField(dist, "tau_u4", (x_basis,))

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

# Biharmonic equation with all four lift terms
add_equation!(problem, "Δ(Δ(u)) + lift(tau_u1, -1) + lift(tau_u2, -2) + lift(tau_u3, -3) + lift(tau_u4, -4) = 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)

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

# Vector tau fields for velocity BCs at each wall
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,))  # Wall at z=0
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,))  # Wall at z=1

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

# Parameters
nu = 0.01
dpdx = -1.0

# Create problem with all fields
problem = IVP([u, p, tau_u1, tau_u2, tau_p])

# Add parameters
add_parameters!(problem, nu=nu, dpdx=dpdx)

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

# Continuity with tau_p (removes degeneracy)
add_equation!(problem, "div(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")

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)

# Vector velocity field and scalar fields
u = VectorField(dist, coords, "u", (x_basis, z_basis))
p = ScalarField(dist, "p", (x_basis, z_basis))
T = ScalarField(dist, "T", (x_basis, z_basis))

# Vector tau fields for velocity BCs
tau_u1 = VectorField(dist, coords, "tau_u1", (x_basis,))  # Wall at z=0
tau_u2 = VectorField(dist, coords, "tau_u2", (x_basis,))  # Wall at z=1

# Scalar tau fields for temperature BCs
tau_T1 = ScalarField(dist, "tau_T1", (x_basis,))  # BC at z=0
tau_T2 = ScalarField(dist, "tau_T2", (x_basis,))  # BC at z=1

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

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

# Create problem with all fields
problem = IVP([u, p, T, tau_u1, tau_u2, tau_T1, tau_T2, tau_p])

# Add parameter substitutions
add_parameters!(problem, Ra=Ra, Pr=Pr)

# Momentum equation (vector form with buoyancy)
# ez is the unit vector in z-direction
add_equation!(problem, "∂t(u) - Pr*Δ(u) + ∇(p) - Ra*Pr*T*ez + lift(tau_u2, -2) = -u⋅∇(u)")

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

# Temperature equation
add_equation!(problem, "∂t(T) - Δ(T) + lift(tau_T2, -2) = -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

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 Tau Fields

Error: ArgumentError: Missing tau field specifications for boundary conditions

Solution: Create tau fields and include lift() terms in your equations:

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

# Correct: Create tau fields and add lift terms
tau_u1 = ScalarField(dist, "tau_u1", (x_basis,))
tau_u2 = ScalarField(dist, "tau_u2", (x_basis,))
add_equation!(problem, "Δ(u) + lift(tau_u1, -1) + lift(tau_u2, -2) = f")
add_bc!(problem, "u(z=0) = 0")
add_bc!(problem, "u(z=1) = 0")

# 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,))  # Correct: 1D on x

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

Tutorial: Boundary Conditions

Overview

The Tau Method Approach

Why Explicit Tau Fields?

Required Steps

Complete Example: Poisson Equation

Dirichlet Boundary Conditions

Basic Setup

No-Slip Velocity (Vector Fields)

Neumann Boundary Conditions

Basic Setup

Stress-Free Conditions

Robin Boundary Conditions

Periodic Boundary Conditions

The lift() Operator

Mode Selection Guidelines

Example: Fourth-Order Problem

Complete Examples

Channel Flow (Poiseuille Flow)

Rayleigh-Bénard Convection

Time- and Space-Dependent Boundary Conditions

Time-dependent BC (scalar)

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

Space + time dependent BC

Two-axis spatial BCs in 3D

1D BC expression in a higher-dimensional problem

User parameters in BC expressions

Common pitfalls

Validation

Checking BC Satisfaction

Troubleshooting

Missing Tau Fields

Over-Specified System

Under-Specified System

Tau Field Dimension Mismatch

See Also