Skip to content

Boundary conditions

Each of the four domain edges carries its own boundary condition, stored as a vector of per-node entries. A single node is a mgrid::BoundaryPoint:

struct BoundaryPoint {
    ConditionType conditionType;  // dirichlet or neumann
    double        value;          // prescribed u, or prescribed normal derivative
};

ConditionType and the edge selectors leftBoundary, rightBoundary, topBoundary, bottomBoundary are plain enums in types.hpp. Two ready-made points cover the common homogeneous cases: zeroDirichletCondition and zeroNeumannCondition.

Setting conditions

mgrid::BoundaryConditions holds one Boundary vector per edge. Set a whole edge to one condition:

solution.boundaryConditions.set(mgrid::leftBoundary, mgrid::zeroNeumannCondition);

or hand it a full Boundary vector for a spatially varying profile (its length must equal the edge extent -- \(n_z\) for the left/right edges, \(n_x\) for top/bottom -- or set throws InvalidBoundaryCondition).

The bundled examples use the channel-flow pattern from Robertson & Kerr (2012), eq. 6: no-slip on the wall and floor, no-shear on the symmetry plane and the free surface.

solution.boundaryConditions.set(mgrid::leftBoundary,   mgrid::zeroNeumannCondition);   // symmetry plane  du/dx = 0
solution.boundaryConditions.set(mgrid::topBoundary,    mgrid::zeroNeumannCondition);   // free surface    du/dz = 0
solution.boundaryConditions.set(mgrid::rightBoundary,  mgrid::zeroDirichletCondition); // channel wall    u = 0
solution.boundaryConditions.set(mgrid::bottomBoundary, mgrid::zeroDirichletCondition); // channel floor   u = 0
solution.propagate_boundary_conditions();

How conditions are applied

FDArray::update_boundaries runs after every relaxation sweep and writes each edge node from its BoundaryPoint:

  • Dirichlet: the node is set to value directly.
  • Neumann: the node is set by a one-sided first-derivative stencil that reproduces the prescribed normal derivative, using the current interior values. The fourth-order form steps four cells inward:

$$ u_0 = \frac{\pm 12\,g\,h + 48u_1 - 36u_2 + 16u_3 - 3u_4}{25}, $$

where \(g\) is the prescribed derivative, \(h\) the normal spacing, and the sign follows the inward direction.

Corner ownership

The two end nodes of a horizontal (top/bottom) edge are the domain corners. update_boundaries lets the vertical (left/right) edges own them: the top/bottom sweeps skip their first and last index. A left/right condition is therefore never clobbered by an adjacent unset top/bottom edge.

Order degradation on coarse grids

The fourth-order Neumann stencil reads four points inward. On the coarsest multigrid levels the grid can be as small as \(4 \times 4\) (Settings::minimumResolution), so update_boundaries drops to a 3rd-, 2nd- or 1st-order one-sided form when fewer than four interior points are available, rather than reading out of bounds:

Interior points available Stencil order
\(\ge 4\) 4 (matches the interior operator)
3 3
2 2
1 1
0 condition ignored (a 1-wide axis cannot carry a Neumann condition)

The lower-order forms leave a small, fixed boundary-consistency error (\({\sim}10^{-3}\)) on the coarse levels that does not shrink under refinement. The fine grids in the hierarchy always fit the fourth-order form, so this does not affect the fine-grid order of accuracy -- the manufactured-solution test still measures order 2 with all-Dirichlet edges, and tests/test_boundary.cpp pins the coarse-grid Neumann path against an out-of-bounds regression.

Propagation across levels

Stack does not alias one BoundaryConditions object across grid levels. propagate_boundary_conditions takes an explicit copy onto the finest level and then builds each coarser level as a stride-2 down-sample of the next finer one -- coarse node \(k\) copies fine node \(2k\) -- which is equivalent to taking stride \(2^d\) straight from the finest edge.

This copy happens once, at construction and whenever you call propagate_boundary_conditions() again. Mutating solution.boundaryConditions after construction has no effect on the coarse levels until you re-propagate -- hence the explicit call in every example constructor after the four set calls.