Flow boundary conditions

Supported boundary conditions:

1. Free slip $\frac{\partial u_i}{\partial x_i}=0$ at the boundary $\mathrm{\Gamma }$

2. No slip $u_i=0$ at the boundary $\mathrm{\Gamma }$

3. Periodic velocity or displacement across paired boundaries

Defining the boundary conditions

We have two ways of defining the boundary condition formulations: - VelocityBoundaryConditions, and - DisplacementBoundaryConditions. The first one is used for the velocity-pressure formulation, and the second one is used for the displacement-pressure formulation. The flow boundary conditions can be switched on and off by setting them as true or false at the appropriate boundaries. Valid boundary names are left and right, top and bot, and for the 3D case, front and back.

For example, if we want to have free free-slip in every single boundary in a 2D simulation, we need to instantiate VelocityBoundaryConditions or DisplacementBoundaryConditions as:

bcs = VelocityBoundaryConditions(;
    no_slip      = (left=false, right=false, top=false, bot=false),
    free_slip    = (left=true, right=true, top=true, bot=true),
)
bcs = DisplacementBoundaryConditions(;
    no_slip      = (left=false, right=false, top=false, bot=false),
    free_slip    = (left=true, right=true, top=true, bot=true),
)

The equivalent for the 3D case would be:

bcs = VelocityBoundaryConditions(;
    no_slip      = (left=false, right=false, top=false, bot=false, front=false, back=false),
    free_slip    = (left=true, right=true, top=true, bot=true, front=true, back=true),
)
bcs = DisplacementBoundaryConditions(;
    no_slip      = (left=false, right=false, top=false, bot=false, front=false, back=false),
    free_slip    = (left=true, right=true, top=true, bot=true, front=true, back=true),
)

The face convention is left/right = x-min/x-max, front/back = y-min/y-max, and bot/top = z-min/z-max in 3D. Flow ghost values are written with the staggered-grid reflection rules appropriate to each component. The kernels run on the selected ParallelStencil backend, including CUDA and AMDGPU extensions.

Periodic faces are set through the periodic keyword and must be enabled in pairs (left/right, front/back, or bot/top):

bcs = VelocityBoundaryConditions(;
    no_slip      = (left=false, right=false, top=false, bot=false),
    free_slip    = (left=false, right=false, top=true,  bot=true),
    periodic     = (left=true,  right=true,  top=false, bot=false),
)

A face can carry at most one of no_slip, free_slip, and periodic; the constructor throws otherwise. A periodic top face is also incompatible with free_surface=true. Faces where all three are false are left untouched by flow_bcs!, which is how a prescribed velocity or displacement is imposed: the caller writes those boundary and ghost values itself.

Periodic boundaries in the Stokes solvers

The two coincident faces of a periodic direction are a single physical plane, so they share one momentum unknown. That extra unknown only exists if the containers are built from the boundary conditions, which is what the bcs argument of StokesArrays is for:

bcs    = VelocityBoundaryConditions(; periodic = (left = true, right = true, top = false, bot = false))
stokes = StokesArrays(backend, ni, bcs)   # NOT StokesArrays(backend, ni)

Building stokes without them leaves the seam out of the momentum residual, and its velocity stays pinned to whatever the initial guess put there. The solvers check the two against each other and throw if they disagree, so the boundary conditions have to be constructed first.

Periodic flow boundaries are currently implemented for solve_DYREL! in 2D and 3D, and for solve_VariationalDYREL! in 2D. solve! and solve_VariationalStokes! reject them rather than silently leaving the seam unsolved.

For the variational solver the volume fractions have to be periodic too: on the Vx and vertex arrays of a RockRatio, index 1 and index end along a periodic direction are the same plane and must carry the same fraction, otherwise the reduced space itself is not periodic.

Two further restrictions are checked at solve time: the periodic direction must have uniform grid spacing (the seam face spans the wrap, so one spacing value has to describe it), and the run must be on a single MPI rank. In particular, do not pass periodx/periody/periodz to init_global_grid for these models: ImplicitGlobalGrid then reports a global grid shrunk by the halo overlap, which is the count the residual norms are normalised by. The periodicity is carried entirely by the velocity boundary conditions.

For what a periodic direction changes inside the solver — the seam momentum row, the wrapping stencils, and the preconditioner — see Periodic boundaries in the Stokes solvers: implementation notes.

For a thermal periodic model there is no momentum row involved: configure the faces on TemperatureBoundaryConditions as shown below and, for a distributed run, mirror the directions in ImplicitGlobalGrid (periodx, periody, periodz) and update the halos after applying the boundary conditions.

Prescribing the velocity/displacement boundary conditions

Normally, one would prescribe the velocity/displacement boundary conditions by setting the velocity/displacement field at the boundary through the application of a background strain rate εbg. Depending on the formulation, the velocity/displacement field is set as follows for the 2D case:

Velocity formulation

stokes.V.Vx .= PTArray(backend)([ x*εbg for x in xvi[1], _ in 1:ny+2]) # Velocity in x direction
stokes.V.Vy .= PTArray(backend)([-y*εbg for _ in 1:nx+2, y in xvi[2]]) # Velocity in y direction

For a pure-shear background field, the package provides a backend-aware, ParallelStencil implementation:

pureshear_bc!(stokes, xci, xvi, εbg)

In 2D it sets Vx = εbg*x and Vy = -εbg*y; in 3D it sets Vx = εbg*x, Vy = εbg*y, and Vz = -εbg*z. Each component uses the vertex coordinates of its own direction. The backend is inferred from stokes. Only the staggered interior ranges are initialized; ghost layers remain available for the configured boundary conditions.

For an xy simple-shear background field, use:

simpleshear_bc!(stokes, xci, xvi, γbg)

This sets Vx = γbg*y and sets the other velocity components to zero: Vy = 0 in 2D, and Vy = Vz = 0 in 3D. The coordinate y is taken from the cell-center coordinates xci[2], while xvi is retained in the API for consistency with pureshear_bc!. As with pure shear, only staggered-grid interior values are initialized; ghost layers are left untouched.

Make sure to apply the set velocity to the boundary conditions. You do this by calling the flow_bcs! function,

flow_bcs!(stokes, flow_bcs)

and then applying the velocities to the halo

update_halo!(@velocity(stokes)...)

Displacement formulation

stokes.U.Ux .= PTArray(backend)([ x*εbg*lx*dt for x in xvi[1], _ in 1:ny+2]) # Displacement in x direction
stokes.U.Uy .= PTArray(backend)([-y*εbg*ly*dt for _ in 1:nx+2, y in xvi[2]]) # Displacement in y direction
flow_bcs!(stokes, flow_bcs)

Make sure to initialize the displacement according to the extent of your domain. Here, lx and ly are the domain lengths in the x and y directions, respectively. Also for the displacement formulation it is important that the displacement is converted to velocity before updating the halo. This can be done by calling the displacement2velocity! function.

displacement2velocity!(stokes, dt) # convert displacement to velocity
update_halo!(@velocity(stokes)...)

Thermal Boundary Conditions

Thermal boundary conditions are collected in TemperatureBoundaryConditions. The same type is used in 2D and 3D.

Supported thermal boundary conditions:

1. No flux $\frac{\partial T}{\partial x_i}=0$ at the boundary $\mathrm{\Gamma }$

2. Constant temperature on the outer boundary $T=T_{\mathrm{\Gamma }}$ at the boundary $\mathrm{\Gamma }$

3. Constant heat flux in the pseudo-transient diffusion kernels $q_T=q_{\mathrm{\Gamma }}$ at the boundary $\mathrm{\Gamma }$

4. Periodic temperature $T({\mathrm{\Gamma }}_i)=T({\mathrm{\Gamma }}_j)$ across paired boundaries

5. Mask-based Dirichlet values inside the domain $T=f(x_i)$ at selected points in $\mathrm{\Omega }$

Face Names

In 2D, boundary tuples use left, right, top, and bot. In 3D, add front and back:

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (
        left = true,
        right = true,
        front = true,
        back = true,
        top = false,
        bot = false,
    ),
)

Faces omitted from no_flux, constant_flux, constant_value, or periodic are treated as inactive for that condition. The dimensionality is inferred from the longest tuple you provide, and the defaults are four-face 2D tuples: pass a complete six-face tuple to obtain a 3D boundary-condition set. A tuple with any other number of faces is rejected.

No-Flux Boundaries

Use no_flux to copy the adjacent interior temperature into the ghost layer. For example, this applies no-flux boundaries on the left and right sides of a 2D domain:

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (left = true, right = true, top = false, bot = false),
)

thermal_bcs!(thermal, thermal_bc)

Constant-Value Boundaries

Use constant_value for fixed-temperature outer boundaries. These values are applied by thermal_bcs! through the ghost-cell relation Tghost = 2 * Tboundary - Tinterior.

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (left = true, right = true, top = false, bot = false),
    constant_value = (left = false, right = false, top = 273.0, bot = 1573.0),
)

thermal_bcs!(thermal, thermal_bc)

If no_flux and constant_value are both active on the same face, thermal_bcs! applies constant_value first and no_flux second.

Periodic Boundaries

Use periodic to copy the opposite interior temperature into the ghost layer. For example, this applies periodic temperature boundaries on the left and right sides of a 2D domain:

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (left = false, right = false, top = false, bot = false),
    periodic = (left = true, right = true, top = false, bot = false),
)

thermal_bcs!(thermal, thermal_bc)

In 3D, include front and back when those faces should also be periodic:

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (
        left = false,
        right = false,
        front = false,
        back = false,
        top = false,
        bot = false,
    ),
    periodic = (
        left = true,
        right = true,
        front = true,
        back = true,
        top = false,
        bot = false,
    ),
)

Periodic faces must be paired by direction and cannot also carry no_flux, constant_flux, or constant_value.

For a single-process run, periodic in the boundary-condition object is sufficient. For a distributed run, also configure the same directions in ImplicitGlobalGrid and exchange halos after applying boundary conditions:

igg = IGG(init_global_grid(nx, ny, 1; periodx = true)...)
thermal_bcs!(thermal, thermal_bc)
update_halo!(thermal.T)

Use periodx, periody, and periodz for the left/right, front/back, and bot/top pairs, respectively. The same grid topology is used by periodic Stokes velocity halos.

Constant-Flux Boundaries

Use constant_flux to prescribe flux values in the pseudo-transient heat diffusion solver. These values are consumed by the PT compute_flux! kernels, not by thermal_bcs!.

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (
        left = true,
        right = true,
        front = true,
        back = true,
        top = false,
        bot = false,
    ),
    constant_flux = (
        left = false,
        right = false,
        front = false,
        back = false,
        top = 0.0,
        bot = 0.03,
    ),
)

Mask-Based Dirichlet Conditions

Use dirichlet for fixed values inside the domain, selected by a mask:

thermal_bc = TemperatureBoundaryConditions(;
    no_flux = (left = true, right = true, top = false, bot = false),
    dirichlet = (; constant = 273.0, mask = mask),
)