Periodic, Symmetric, and Robin Boundary Conditions: Advanced FiPy Patterns

Key Takeaways

• Periodic boundaries in FiPy are most cleanly implemented with PeriodicGrid objects or Gmsh periodic mesh commands rather than manual face coupling.
• Symmetric boundaries are zero-flux Neumann conditions applied to symmetry planes. They are not a separate FiPy category, but a physical use of standard boundary condition patterns.
• Robin boundaries combine value and flux constraints. They can be implemented with ImplicitRobinSource, ImplicitSourceTerm, or a FaceVariable divergence approach.
• The most common Robin boundary pitfalls are sign errors from incorrect flux direction conventions and variable aliasing when multiple Python names reference the same FaceVariable object.
• For time-dependent periodic boundaries, use conditional expressions such as t % 10 rather than Python functions. Python functions do not re-evaluate automatically as the simulation time variable changes.

If you have worked through basic FiPy tutorials, you already know that boundary conditions are where theory meets the edges of the mesh. Dirichlet, Neumann, Robin, and periodic boundaries may look simple in theory, but implementation details matter in real simulations.

This article covers advanced FiPy patterns for periodic domains, symmetry planes, and Robin boundaries. These are common in unit-cell models, axisymmetric problems, convective heat transfer, and mass exchange simulations.

The examples below show how to implement each boundary type correctly, how to choose between approaches, and how to avoid silent failures that can waste hours of debugging.

Periodic Boundary Conditions: Coupling the Edges

Periodic boundaries glue opposite sides of the domain together. A solution leaving one side reappears on the opposite side. This pattern is common in unit-cell simulations, representative volume elements, materials science, and simplified transport models where the real domain repeats.

The Simple Approach: PeriodicGrid Objects

FiPy provides built-in periodic grid objects that handle face coupling automatically. For 1D domains, the setup is simple:

from fipy import PeriodicGrid1D, CellVariable, DiffusionTerm, TransientTerm

nx = 50
dx = 1.0
mesh = PeriodicGrid1D(nx=nx, dx=dx)

C = CellVariable(name="concentration", mesh=mesh, value=0.0)

# No explicit boundary condition needed.
# FiPy handles periodic coupling through the mesh.
eq = TransientTerm() == DiffusionTerm(coeff=1.0)

# Initialize with a non-uniform perturbation.
x = mesh.cellCenters[0]
C.setValue(1.0, where=(x > 0.7) & (x < 0.8))

eq.sweep(var=C)

For 2D domains, FiPy provides PeriodicGrid2D and oriented variants:

from fipy import PeriodicGrid2D, CellVariable, DiffusionTerm

nx, ny = 50, 50
dx = dy = 1.0

# Periodic on all sides
mesh = PeriodicGrid2D(nx=nx, ny=ny, dx=dx, dy=dy)

# Oriented variants can be used when needed:
# PeriodicGrid2DLeftRight: periodic only in x
# PeriodicGrid2DTopBottom: periodic only in y

The oriented variants are useful when you need periodicity in one direction but different boundary conditions in the other. For example, a transport problem may need periodic inflow and outflow in one direction and no-flux walls in the other.

Complex Geometries: Gmsh Periodic Meshes

When the domain is not a simple line or rectangle, use Gmsh to define periodic relationships at the mesh level. The important commands are Periodic Line and Periodic Surface.

// In a Gmsh .geo file
Periodic Line {1} = {2} + 10;
Periodic Surface {3} = {4};

Then import the mesh in Python:

from fipy import GmshImporter2D

mesh = GmshImporter2D("mesh.geo")

The advantage is that periodic coupling is defined in the mesh, not manually in the solver. This is more robust for arbitrary geometries because the solver works with a mesh that already contains the periodic relationships.

Manual Face Coupling: When You Need Control

For custom cases, you may want to manually map faces and apply constraints. This gives full control but requires careful face indexing.

from fipy import CellVariable, Grid1D

nx = 50
dx = 1.0
mesh = Grid1D(nx=nx, dx=dx)

C = CellVariable(mesh=mesh, value=0.0)

left_faces = mesh.facesLeft
right_faces = mesh.facesRight

# Manual periodic-style coupling
C.constrain(C.faceValue[mesh.facesLeft.value], where=mesh.facesRight)

Manual coupling is more fragile than using PeriodicGrid. Prefer built-in periodic meshes unless you need behavior they cannot provide.

Time-Dependent Periodic Boundaries

Time-dependent periodic boundaries can be tricky. A common mistake is to define a Python function and expect FiPy to re-evaluate it automatically:

# Wrong: this evaluates once and then stays fixed
def sw_func(t):
    if t < 5:
        return 0.0
    elif t < 10:
        return 1.0
    else:
        return sw_func(t - 10)

valueRight = sw_func(t)

The function result is evaluated once. As t changes, the boundary value does not automatically update.

Use conditional expressions and FiPy-compatible variables instead:

# Correct pattern: expression can re-evaluate as t changes
t_mod = t % 10
valueRight = ((t_mod >= 5) & (t_mod < 10)) * 1.0

C.constrain(valueRight, where=mesh.facesRight)

The key point is that FiPy variables behave like symbolic expressions. Build time-dependent conditions with FiPy-compatible expressions, logical masks, and array-style operations rather than Python if statements.

When Periodic Boundaries Fail Silently

A common silent failure is using a uniform initial condition on a periodic domain:

C = CellVariable(mesh=mesh, value=1.0)

eq = TransientTerm() == DiffusionTerm(coeff=1.0)
eq.solve(var=C, dt=0.1)

The result never changes because the field has no gradient. Periodic boundaries do not create gradients by themselves.

Always initialize periodic domains with a non-uniform condition, such as a perturbation, step function, or Gaussian blob.

Symmetric Boundary Conditions: Zero Flux on Symmetry Planes

Symmetric boundaries exploit symmetry so you can simulate only part of the domain. This is common in axisymmetric heat transfer, crystallographic models, and flows with symmetry planes.

In FiPy, symmetric boundaries are implemented as zero-flux Neumann conditions on the symmetry plane.

from fipy import Grid2D, CellVariable, DiffusionTerm, TransientTerm

nx, ny = 50, 25
dx = dy = 0.1

mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)

T = CellVariable(name="temperature", mesh=mesh, value=100.0)

# Symmetry plane at the top: zero flux
T.faceGrad.constrain(0.0, where=mesh.facesTop)

# Bottom boundary: fixed temperature
T.constrain(100.0, where=mesh.facesBottom)

# Left and right: insulated
T.faceGrad.constrain(0.0, where=mesh.facesLeft)
T.faceGrad.constrain(0.0, where=mesh.facesRight)

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

The physical meaning is that heat or mass does not cross the symmetry plane because the omitted half of the domain mirrors the simulated half.

Axisymmetric Cylindrical Coordinates

Axisymmetric problems require extra care because the radial coordinate changes the diffusion operator. In cylindrical coordinates, the radial diffusion term includes a geometric factor:

(1/r) ∂(r ∂C/∂r) / ∂r

FiPy can represent this type of behavior, but the setup requires careful coefficient construction and mesh design.

from fipy import CellVariable, FaceVariable, DiffusionTerm, TransientTerm

r = mesh.cellCenters[0]

# Example structure for a radial coefficient
axisymmetric_coeff = FaceVariable(mesh=mesh, value=1.0)

eq = TransientTerm() == DiffusionTerm(coeff=axisymmetric_coeff)

The axis at r = 0 is singular because of the 1/r factor. Avoid placing cell centers exactly at r = 0. Start the mesh at a small positive radius, use cell-centered radial coordinates carefully, or use a mesh design that avoids the singularity.

Robin Boundary Conditions: Mixed Value-Flux Constraints

Robin boundary conditions specify a linear combination of the variable and its normal derivative. They combine a value constraint with a flux constraint.

A common example is Newton’s law of cooling:

k ∂T/∂n + h(T – T∞) = 0

Here, heat flux at the boundary is proportional to the difference between the surface temperature and the ambient temperature.

The ImplicitRobinSource Approach

ImplicitRobinSource is often the cleanest approach for convective heat transfer and mass exchange.

import fipy as fp

nx = 50
dx = 1.0
mesh = fp.Grid1D(nx=nx, dx=dx)

T = fp.CellVariable(name="temperature", mesh=mesh, value=20.0)

# Physical parameters
k = 1.0
h = 0.5
T_inf = 10.0
T_left = 100.0

# Equation setup
eq = fp.TransientTerm() == fp.DiffusionTerm(coeff=k)

# Left boundary: fixed temperature
bc_left = fp.FixedValue(mesh.facesLeft, T_left)

# Right boundary: convective heat transfer
bc_right = fp.ImplicitRobinSource(
    mesh.facesRight,
    coeff=h,
    value=h * T_inf
)

dt = 0.1
for step in range(100):
    eq.solve(var=T, boundaryConditions=[bc_left, bc_right], dt=dt)

The ImplicitRobinSource approach adds an implicit source or sink contribution to the solver matrix. This can be more stable than explicit boundary handling, especially for heat exchange or mass transfer problems.

The FaceVariable Divergence Approach

For spatially varying coefficients or coupled multi-physics boundaries, a FaceVariable divergence approach gives more control.

from fipy import CellVariable, FaceVariable, Grid1D
from fipy import DiffusionTerm, ImplicitSourceTerm, TransientTerm

nx = 100
dx = 1.0
mesh = Grid1D(nx=nx, dx=dx)

C = CellVariable(mesh=mesh, value=0.0)

D_coeff = 2.0
h_coeff = 1.0
C_ambient = 0.5

# Step 1: define face coefficient
D_face = FaceVariable(mesh=mesh, value=D_coeff)

# Step 2: zero out coefficient at the boundary as a precaution
D_face[..., mesh.facesLeft.value] = 0.0

# Step 3: formulate Robin-like boundary contribution
robin_flux = (h_coeff * (C_ambient - C)).divergence

eq = (
    TransientTerm()
    == DiffusionTerm(coeff=D_face)
    + robin_flux
    - ImplicitSourceTerm(coeff=0.0)
)

This approach requires attention to three details:

1. Zero out the diffusion coefficient at the target boundary face.
2. Formulate the Robin flux as a divergence source.
3. Add the contribution to the equation consistently.

Zeroing the boundary coefficient is often treated as a safe precaution because it prevents the default no-flux behavior from conflicting with the custom Robin contribution.

Common Robin Boundary Pitfalls

Sign Errors in the Flux Direction

The sign convention can flip depending on whether the flux is treated as an inward source or outward sink.

# Wrong sign
RobinCoeff = maskSurfaces * D * n / (-dPf.dot(a) + b)

# Corrected sign pattern
RobinCoeff = maskSurfaces * D * n / (dPf.dot(a) + b)

Verify limiting cases:

• If h → 0, the boundary should behave like zero-flux Neumann.
• If h → ∞, the boundary should behave like fixed-value Dirichlet.

Variable Aliasing

Variable aliasing happens when multiple Python names point to the same FaceVariable object.

# Wrong: all names reference the same object
D = FaceVariable(mesh=mesh, value=0.9)
Gamma = D
b = D

Gamma.setValue(0.0, where=maskSurfaces)

Changing Gamma also changes D and b. This can silently break the boundary condition.

Use a non-mutating expression instead:

# Better: create a new masked expression
Gamma = D * ~maskSurfaces

This creates a new expression that is D away from the boundary and zero on the boundary, without mutating the original object.

Convergence Issues

Robin boundaries can make steady-state solves unstable or difficult to converge. A practical workaround is to introduce time dependence even for a problem that should reach steady state. Transient evolution can help the solver approach the correct solution more gradually.

Spatially Varying and Patch-Based Boundaries

Boundary conditions can be applied only to selected regions by using spatial masks. This is useful when different parts of the same boundary have different physical behavior.

from fipy import Grid2D, CellVariable, DiffusionTerm

nx, ny = 80, 80
dx = dy = 0.1
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)

T = CellVariable(mesh=mesh, value=0.0)

X, Y = mesh.faceCenters

# Left half of the boundary
mask_left = X < 0.5 * nx * dx

# Fixed value on one patch
T.constrain(1.0, where=mesh.facesLeft & mask_left)

# Zero flux on the other patch
T.faceGrad.constrain(0.0, where=mesh.facesLeft & ~mask_left)

This patch boundary approach is important for problems with mixed physical conditions, such as a wall that is partly heated and partly insulated.

What We Recommend: Choosing the Right Approach

Periodic Boundaries

1. Use PeriodicGrid for simple structured 1D, 2D, or 3D domains.
2. Use Gmsh periodic meshes for complex geometries because mesh-level coupling is more robust than manual face constraints.
3. Initialize with non-uniform conditions because uniform values produce no evolution in pure diffusion problems.
4. Use conditional expressions such as t % N for time-dependent periodic boundaries, not Python functions.

Symmetric Boundaries

1. Apply zero-flux Neumann conditions with faceGrad.constrain(0.0) on the symmetry plane.
2. For axisymmetric problems, handle the r = 0 singularity carefully through mesh design.
3. Verify symmetry by checking that the solution gradient is zero at the symmetry plane.

Robin Boundaries

1. Use ImplicitRobinSource for convective heat transfer and mass exchange.
2. Use the FaceVariable divergence approach for spatially varying coefficients or coupled multi-physics boundaries.
3. Zero out coefficients at the boundary face as a precaution when using custom flux formulations.
4. Watch for sign errors and verify limiting cases.
5. Watch for variable aliasing. Prefer masked expressions such as D * ~mask instead of mutating shared variables with setValue().

Comparison: When to Use Each Approach

Troubleshooting: Common Boundary Condition Failures

The Solution Does Not Change

If a periodic domain produces the same value everywhere, check the initial condition. A uniform initial condition has no gradient, so pure diffusion has nothing to evolve.

The Boundary Value Is Wrong

For Robin boundaries, check:

• The sign of the flux direction.
• Whether you are constraining faces rather than cells.
• Whether the diffusion coefficient was zeroed at the boundary in custom formulations.

The Solver Diverges

For Robin boundaries, try:

• Reducing dt and increasing resolution.
• Adding transient evolution even for problems intended to reach steady state.
• Using PETSc solvers instead of simple sweeps for complex boundary setups.

The Boundary Does Not Update

For time-dependent periodic boundaries, the likely cause is a Python function used in place of a FiPy expression. Use Variable objects and logical masks instead of if and elif statements.

Conclusion

Periodic, symmetric, and Robin boundary conditions are where FiPy’s theory meets the edges of the mesh. Periodic grids handle face coupling automatically. Symmetry planes are zero-flux Neumann conditions. Robin boundaries combine value and flux constraints through ImplicitRobinSource or a FaceVariable divergence approach.

The most important habit is to verify limiting cases. For Robin boundaries, check that h → 0 gives no-flux behavior and h → ∞ gives fixed-value behavior. For periodic domains, verify that the initial condition is not uniform. For symmetric boundaries, confirm that the gradient at the symmetry plane is zero.

These checks catch most boundary condition errors before they become long debugging sessions.

Related Guides

• Boundary Conditions: Theory and Implementation in FiPy — The fundamentals of Dirichlet, Neumann, and mixed boundaries.
• Working Through Your First FiPy Example — Mesh setup, variables, diffusion equations, and boundary conditions.
• Understanding FiPy’s Core Architecture — Cell versus face variables, boundary condition constraints, and matrix assembly.
• FiPy GitHub Discussions — Community examples and troubleshooting.

Next Steps

1. Start with PeriodicGrid1D for a simple 1D periodic domain.
2. Add time dependence with conditional expressions such as t % N.
3. For Robin boundaries, begin with ImplicitRobinSource and verify limiting cases.
4. Move to Gmsh meshes for complex geometries once structured cases work.

For help implementing boundary conditions in specific research models, visit our homepage for consulting resources and FiPy development support.