Boundary Conditions: Theory and Implementation in FiPy

Boundary conditions are the “rules at the edges” of a simulation domain. If your PDE describes what happens inside the mesh, boundary conditions describe what the world outside the mesh is doing to it (or not doing to it). In practice, they often determine whether your solution is physically meaningful, stable, and reproducible.

If you’ve ever run a diffusion or heat equation and got a curve that drifts oddly, flattens too fast, or refuses to change at all, boundary conditions are one of the first places to check. In FiPy, boundary conditions are usually implemented with constraints placed on boundary faces. The good news: once you understand the few common patterns, you can implement most real-world cases with small, readable code changes.

What Boundary Conditions Do in a PDE Model

A PDE alone typically has infinitely many solutions. Boundary conditions reduce those possibilities to the one that matches your physical setup. They also influence conservation (mass/energy), the long-time steady state, and how quickly the transient solution evolves.

Two important ideas to keep in mind:

• Boundary conditions are not optional. Even “do nothing” at the boundary is a choice (often interpreted as zero-flux).
• The same PDE can behave completely differently under different boundary assumptions (e.g., fixed temperature vs insulated walls).

The Most Common Types of Boundary Conditions

Dirichlet boundary condition (fixed value)

Dirichlet means the variable is pinned to a prescribed value on the boundary. In heat transfer, this is a wall held at a fixed temperature. In diffusion, it can represent a boundary in contact
with a large reservoir at fixed concentration. Conceptually: “the field equals this value at the boundary.”

Neumann boundary condition (fixed flux or fixed gradient)

Neumann specifies the flux (or equivalently the derivative/gradient normal to the boundary). The most common special case is zero-flux (insulated boundary, no mass/heat crosses the boundary).

Conceptually: “the derivative/flux equals this value at the boundary.”

Robin boundary condition (mixed: value + flux)

Robin mixes value and flux in one relationship. A classic example is convective heat transfer at a wall: heat flux is proportional to the difference between the wall temperature and an ambient temperature.

Conceptually: “flux depends on how far the boundary value is from a reference.”

Periodic boundary condition

Periodic means the left and right boundaries are “glued together”: the solution repeats. This is common in idealized unit-cell simulations, repeating materials, and simplified transport problems.

Mixed and piecewise boundary conditions

Many real problems use different boundary types on different sides (e.g., left boundary fixed value, right boundary zero-flux). In FiPy, this is normal and usually implemented by applying different constraints on different face sets.

A Practical Mapping: Type → Meaning → FiPy Implementation

The table below summarizes the most common boundary choices and how they typically look in FiPy. The exact details can vary by equation and term, but this gives you a reliable starting point.

Baseline Example: 1D Diffusion with Different Boundary Conditions

We’ll use the standard diffusion PDE in 1D:
∂φ/∂t = D ∂²φ/∂x².
In FiPy this is often written as:

TransientTerm() == DiffusionTerm(coeff=D)

The boundary behavior is controlled by what we constrain on the boundary faces.

Setup: mesh, variable, equation

# Save as: bc_fipy_diffusion_1d.py

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm
import matplotlib.pyplot as plt

nx = 200
dx = 1.0 / nx
mesh = Grid1D(nx=nx, dx=dx)

phi = CellVariable(name="phi", mesh=mesh, value=0.0)

x = mesh.cellCenters[0]
phi.setValue(1.0, where=(x < 0.3))  # a block/step-like initial condition

D = 1.0
eq = TransientTerm() == DiffusionTerm(coeff=D)

Case A: Dirichlet boundaries (fixed values at both ends)

This forces the field to remain pinned at the boundaries. For diffusion, this often represents boundaries connected
to reservoirs.

# Dirichlet: left boundary fixed to 1.0, right boundary fixed to 0.0
phi.constrain(1.0, where=mesh.facesLeft)
phi.constrain(0.0, where=mesh.facesRight)

What to expect: over time the solution tends toward a steady gradient-like profile consistent with the fixed ends.

Case B: Neumann zero-flux boundaries (insulated ends)

A common “default” physical boundary for diffusion is no-flux: nothing crosses the boundary. In FiPy, a typical way to express a fixed gradient at faces is constraining var.faceGrad. Zero-flux is a zero gradient normal to the boundary.

# Neumann: zero gradient (no-flux) at both ends
phi.faceGrad.constrain(0.0, where=mesh.facesLeft)
phi.faceGrad.constrain(0.0, where=mesh.facesRight)

What to expect: the total “mass” (integral of φ) tends to be conserved for pure diffusion, and the profile eventually approaches a uniform value (the average), because diffusion smooths until nothing drives gradients.

Case C: Mixed boundary (Dirichlet on left, no-flux on right)

This is a realistic scenario: one side in contact with a fixed reservoir, the other insulated.

# Mixed:
phi.constrain(1.0, where=mesh.facesLeft)          # fixed value on the left
phi.faceGrad.constrain(0.0, where=mesh.facesRight) # no-flux on the right

What to expect: the domain is “fed” by the left boundary while the right boundary blocks outflow, often causing a different steady state than the symmetric cases.

Time loop and plotting (works for all cases)

dt = 1e-4
steps = 3000

snapshots = []
times = []
capture_at = {0, 50, 200, 800, 2000, 3000}

for step in range(steps + 1):
    if step in capture_at:
        snapshots.append(phi.value.copy())
        times.append(step * dt)
    eq.solve(var=phi, dt=dt)

plt.figure()
for y, t in zip(snapshots, times):
    plt.plot(x.value, y, label=f"t={t:.4f}")

plt.xlabel("x")
plt.ylabel("phi")
plt.title("Diffusion with Boundary Conditions (FiPy)")
plt.legend()
plt.show()

How to Think About Boundary Conditions in FiPy

A useful mental model is: FiPy solves for cell-centered values, but boundary conditions are applied on faces. When you constrain phi on boundary faces, you are enforcing the boundary value. When you constrain phi.faceGrad, you are controlling the gradient (and therefore flux for diffusion).

If your PDE includes diffusion, convection, reaction, sources, or multiple coupled variables, boundary conditions can interact with those terms. In many cases, the “right” implementation is the one that makes your conservation checks and expected limiting behavior come out correctly.

Robin Boundaries in Practice: The Idea and a FiPy Strategy

Robin boundaries often appear in heat and mass transfer when the boundary exchanges with an environment. A typical form is:
-D ∂φ/∂n = h (φ - φ_env)
where h is a transfer coefficient and φ_env is an ambient/reference value.

In many numerical implementations, Robin conditions can be handled by adding a boundary “sink/source” contribution that depends on the boundary value. In FiPy, a common strategy is to represent this effect through an additional term that acts on cells adjacent to the boundary (or through specialized boundary condition helpers, when used).

The key practical point: Robin is not “just a Dirichlet” and not “just a Neumann.” It is a relationship between the two. When you implement it, validate by checking limiting cases:

• If h → 0, the boundary should behave like no-flux (Neumann zero).
• If h → ∞, the boundary should behave like fixed value at φ = φ_env (Dirichlet).

Boundary Conditions in 2D: What Changes, What Stays the Same

In 2D, you still constrain faces. The difference is that you have more boundary segments: left, right, top, bottom (and possibly more if the geometry is complex).

from fipy import Grid2D, CellVariable, TransientTerm, DiffusionTerm

nx, ny = 80, 80
dx = dy = 1.0 / nx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)

phi = CellVariable(mesh=mesh, value=0.0)
D = 1.0
eq = TransientTerm() == DiffusionTerm(coeff=D)

# Example: left boundary fixed to 1, other boundaries insulated
phi.constrain(1.0, where=mesh.facesLeft)
phi.faceGrad.constrain(0.0, where=mesh.facesRight)
phi.faceGrad.constrain(0.0, where=mesh.facesTop)
phi.faceGrad.constrain(0.0, where=mesh.facesBottom)

The same validation habits apply: you should know what “steady state” looks like for your chosen boundaries, and you should run a simple scenario first before layering on complexity.

A Quick Validation Checklist

Boundary conditions are easy to write and surprisingly easy to get subtly wrong. These checks catch many issues early:

• Does the solution approach a reasonable steady state consistent with the boundaries? (e.g., fixed ends produce a stable gradient-like profile)
• For pure diffusion with insulated boundaries, is the average value conserved over time (approximately)?
• If you swap left/right boundary types, does the solution transform in a way you’d expect?
• If you reduce dt and increase nx, does the solution converge (changes become smaller)?

Troubleshooting: Common Boundary-Condition Pitfalls

The solution “sticks” and barely evolves

• dt may be extremely small relative to your spatial resolution and diffusion coefficient.
• Your initial condition might already satisfy the boundary constraints and be near steady state.
• You might be constraining more than you intended (e.g., pinning multiple sides to constants).

The boundary values are not what you expect

• Confirm you are constraining boundary faces (e.g., mesh.facesLeft), not cells.
• Check ordering: apply constraints after creating the variable and before solving.
• If you reuse a variable across experiments, reset constraints (or recreate the variable) to avoid “leftover” constraints.

Unstable or noisy behavior near boundaries

• Reduce dt and check whether the instability disappears.
• Increase resolution (nx, ny) to reduce discretization artifacts.
• Make sure the boundary condition matches the physics of the term (flux conditions should match the modeled flux).

How to Choose the Right Boundary Condition

A fast way to choose is to answer one sentence for each boundary:

• Is the boundary held at a known value (Dirichlet)?
• Is the boundary insulated / no transport crosses (Neumann zero-flux)?
• Is the boundary exchanging with an environment at a rate (Robin)?
• Is the domain repeating (Periodic)?

If you can’t answer confidently, start with the simplest physically defensible option and add complexity only after you can explain the differences in results. Boundary conditions are a modeling decision, not just a syntax detail.

Conclusion

Boundary conditions are where theory meets the edges of your mesh. In FiPy, they’re typically expressed by placing constraints on boundary faces: fixed values via var.constrain(), and fixed gradients/fluxes viavar.faceGrad.constrain() (with more advanced behaviors built from the same ideas).

If you take one habit forward, make it this: always predict the qualitative behavior before you run the model. If the simulation violates that prediction, inspect boundary conditions first, then time step and discretization. Once your boundary logic is solid, FiPy becomes much easier to scale from toy examples to research-grade models.