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.

For a distributed run, configure the same directions in ImplicitGlobalGrid (periodx, periody, periodz) and update the velocity 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),
)