Documentation of the main branch

These pages follow Mallard's main branch, which is ahead of release 0.8.0: some features they describe come with the next release.

Design: characteristic boundary conditions and sponge layers¶

Status: implemented (nscbc_outlet, nscbc_inlet, [[sponges]]).

Flame, jet and wake simulations need open boundaries that let acoustic waves, vortices and flames leave without reflecting. They must also hold a mean pressure or inflow state. p_out and upt impose their state exactly and reflect almost all of an acoustic wave. extrapolation holds nothing, and (see below) still reflects. The Navier-Stokes characteristic boundary conditions (NSCBC) of Poinsot & Lele (1992) let outgoing waves leave unchanged and make the incoming wave a weak relaxation toward the target.

The incoming wave of a face¶

In NSCBC, the normal derivatives at a boundary point are split into the amplitudes L_k = lambda_k dw_k/dn of the characteristic waves. Here w_k is the characteristic variable and lambda_k its speed along the outward normal n. Outgoing amplitudes come from the interior; incoming ones are prescribed, and the boundary point advances in time by these LODI relations.

Mallard has no boundary points: a boundary face's flux comes from a Riemann solver, with the reconstructed interior state W_l on one side and an exterior state W_g on the other. The outgoing acoustic wave is that of W_l:

outgoing w+ = p + Z u_n, incoming w- = p - Z u_n, with Z = rho c
c^2 = gamma p / rho, with the face's own gamma for mixtures (its frozen surrogate)

The incoming acoustic wave is the face's own boundary point. Each characteristic face keeps a pressure p_b and normal velocity u_b. Once per step, at the first Runge-Kutta stage, its w- = p_b - Z u_b advances by a forward-Euler step of the LODI relation:

dw-/dt = -K (p_b - p_t) - beta T-            (outlets)
dw-/dt = K Z (u_b - u_n,t) - beta T-         (inlets, Poinsot & Lele's L_5)
K = sigma c (1 - M^2) / L

Here p_b and u_b are first refreshed from the current outgoing wave. K dt is capped at 1. The face's w- goes into the exterior state at every stage of the step.

The face keeps p_b and u_b rather than w- itself: pressure and velocity are continuous at contacts, w- is not. A hot front leaving through an outlet changes Z by a factor of two, and a stored w- would then put an O(Z u) pressure jump on the face.

p_b and u_b start from the interior face state at the first step. They are the only state of the boundary, and restart files carry them (format version 4), so restarted runs reproduce uninterrupted ones bitwise. Each face is keyed by the global id of its cell and its local face in that cell, and the faces are written in key order, so a file written on any number of ranks continues on any other. A file without some face of the run (version 3 and earlier, or one written with that boundary not characteristic) starts every face afresh from the solution.

An inlet's target velocity is per face, and may vary in time (synthetic turbulence, synthetic_inflow.md). Its change over the step enters the incoming wave as the target's own incoming wave, -Z (u_n,t(t + dt) - u_n,t(t)) (Guézennec & Poinsot 2009), so that the face follows the target without reflecting outgoing waves.

This is the LODI relation of Poinsot & Lele itself. For a plane wave at normal incidence, dw-/dt = -K p' gives the reflection coefficient R = -1 / (1 - 2 i omega / K) (Selle, Nicoud & Poinsot 2004). Low frequencies reflect, so the mean pressure is held; high frequencies leave. The pulse test checks this.

Stateless formulations do not work¶

The first implementation was stateless. It shifted the incoming wave of a transmissive state by -(V/A) L / |lambda|, so that the boundary cell itself obeyed dw-/dt = -L. Two problems:

Steady offsets. That cell's w- is frozen except for L. Any steady source in the cell is balanced by K (p - p_t), and the outlet settles away from its target by source / K:
viscous stress in a Poiseuille flow: 20% of the flow rate
heat release of a flame at the outlet: the flame stayed anchored there, or the domain kept pressurizing
Instability. Taking the incoming wave from the face's own state instead differences the incoming wave downwind, and blows up.

With its own boundary point, the face lets the boundary cell relax toward w- at the rate |lambda| / h. Sources in the cell are carried inward as waves, as in the interior, and the face's pressure settles at p_t.

Exterior state¶

Exterior state on a subsonic face:

Wave	Source in the exterior state
Outgoing acoustic	`W_l`
Incoming acoustic	the face's `p_b`, `u_b`, combined with the local `Z`
Entropy, tangential velocity and composition at outflow	`W_l`, with an isentropic density change
Incoming convective waves (backflow at outlets)	the transmissive state `W_e` of `extrapolation` (the reconstructed state on the interior image face), not a copy of the boundary cell

Special cases:

Supersonic outflow (u_n >= c): the exterior state is W_l, and the face follows the interior.
Supersonic inflow at inlets (u_n <= -c): the exterior state is the target.

Low-Mach correction. Characteristic faces take the low-Mach correction (Thornber et al.) like interior faces. Without it, the boundary cell's two sides have different dissipation. In a steady shear flow, this held the outlet pressure away from its target.

Reconstruction next to the boundary¶

A wave that leaves through an open boundary must leave the boundary cell's reconstruction intact. Zero-gradient ghosts break this: the mirror cells of TENO-E stencils and the ghost of the MUSCL least-squares fit force the normal derivative toward zero, and the reconstructed states jump by O(h) at the boundary cell's inner face. With exact upwinding, such a jump in an outgoing wave does nothing. But the default low-Mach correction scales the velocity jump by z = max(M, 0.1), which mixes it into the incoming wave.

A normal acoustic pulse leaving through extrapolation (400 cells over the domain, the pulse 20 cells wide, default low-Mach correction) reflects by these fractions:

Reconstruction	Quadrilaterals	Triangles	Jittered mixed
MUSCL	15%	5%	0.6%
TENO5	25%	4%	14%

MUSCL boundary cells therefore leave out the ghost of a characteristic face. They use the vertex-neighbor least-squares fit without it, the fit MUSCL already uses on tetrahedra. On triangles the face-neighbor fit without the ghost is exactly determined and unstable.

TENO-E keeps its mirror cells across characteristic faces. Its one-sided stencils without them are unstable in single precision, and with the face's own incoming wave the mirrors cost little:

Reconstruction	`extrapolation`	`nscbc_outlet`, `sigma = 0`
MUSCL	15%	0.2% (1% with the ghost)
TENO5	34%	1%

(the pulse reflected on quadrilaterals, 200 cells over the domain, pulse 10 cells wide)

Everything else treats the face like extrapolation:

the image-face values and gradients of the viscous fluxes, or zero normal derivatives where no image exists (the outflow viscous conditions of Poinsot & Lele)
the ghosts of viscous gradients (without the characteristic face's ghost, as for MUSCL) and of scalar reconstruction

Transverse terms¶

T- holds the transverse terms of the incoming wave's equation. These are the tangential convection and the divergence of the tangential velocity, which normal derivatives do not see:

T- = u_t . grad p + rho c^2 div_t u_t - rho c (u_t . grad u_n)

With beta = 1, the face evolves under them in full, as in Poinsot & Lele's original condition. At low Mach number this pushes the low pressure of a leaving vortex core back toward p_t, which sends a pressure wave upstream and sets off a transverse slosh between walls.

Yoo et al. (2005), Yoo & Im (2007) and Lodato, Domingo & Vervisch (2008) cancel them in the incoming wave except for a fraction beta, equal to the Mach number. Granet et al. (2010) show this matters most at low Mach number. The default is the local beta = M; beta = 0 drops the terms.

The tangential derivatives are a least-squares fit in the plane of the face, as finite-difference NSCBC differentiates along the boundary:

Data. The fit uses the states of the cells of the neighboring characteristic faces: those with nearly the same normal sharing an edge (3D) or a node (2D), also across periodic seams, at the nearest periodic image of their centers. A one-sided fit at a seam biased the transverse terms there: with turbulence leaving through an outlet periodic across, the outlet's pressure along the seam fell and the run diverged after a few flow-through times. Edge faces have no fit (below).
Normal derivatives. A point along the normal sets them to zero.
Why not the boundary cell's own gradient. It picks up the normal variation of waves crossing the boundary. On triangles it also picks up the grid-scale odd-even velocity pattern as a transverse divergence.
Why not reconstructed face states. With TENO, they made the correction unstable.

Stability: stencil, edges and reversed flow¶

The incoming wave integrates -beta T- every step, and only the weak relaxation K pulls it back (K = 0.02 for L = 12 and sigma = 0.25). So any mode for which the transverse terms have the wrong sign grows. Where the interior acts on it only weakly, the outlet's pressure drifts by up to beta T- / K, the flow there reverses, and the run diverges. Turbulence leaving through an outlet found three such modes:

Faces sharing only a corner (3D). On a structured outlet, a face and its eight neighbors sharing a node hold a grid-scale odd-even pattern across one direction, say v(y) (-1)^k along z. In their fit, two of the three columns that set dv/dy have the opposite sign, so the fit returns -dv/dy / 3: a divergence of the wrong sign. The incoming wave then amplifies the pattern, a checkerboard in the outlet's last cell layer. The fit therefore uses the faces sharing an edge (in 2D, a node). For this pattern they give the right sign.
Edges: faces sharing a node with a wall, symmetry plane or other non-characteristic boundary. There the fit is one-sided. A pressure rise at such a face pushes fluid away from the wall only, so the tangential velocity falls with distance from the wall. The fit reads this as a negative divergence; with the wall's v = 0 the true one is positive. The incoming wave then raises the face's pressure further. These faces get no transverse terms.
Reversed faces of outlets (u_n < 0). The transverse terms let structures leave without the outlet pushing their pressure back. A reversed face lets nothing out. Its incoming convective waves come from the transmissive state, so beta T- integrates the spreading of the reversed jet itself, whose pressure then rises with the jet. These faces have beta = 0: they only relax their pressure toward p_t, like an inflow at held pressure.

Test case. A 3D box [0, 2] x [0, 1]^2:

nscbc_inlet at u = 0.2 (c = 1), and nscbc_outlet with L = 12
symmetry planes on the sides, or no-slip walls with cells clustered toward them
mu = 2e-4, MUSCL without a limiter
decaying eddies A sin(2 pi x) cos(2 pi y) cos(2 pi z) (plus smaller modes), with A up to six times the mean flow, crossing the outlet with reversed flow

The table gives the largest Mach number and the lowest pressure at the end (p_t = 0.714):

Case	Before	Now	`beta = 0`
`48 x 24 x 24`, symmetry, `A = 0.4`, t = 25	diverges at t = 11.6	0.31, 0.728	0.30, 0.703 (t = 20)
`48 x 24 x 24`, walls, `A = 0.4`, t = 25	diverges at t = 18.5	0.35, 0.711	0.38, 0.704
`48 x 24 x 24`, symmetry, `A = 0.8`, t = 25	diverges at t = 6.0	0.43, 0.804	0.41, 0.707
`48 x 24 x 24`, symmetry, `A = 1.2`, t = 25	diverges at t = 0.7	0.50, 0.878	0.58, 0.668
`96 x 48 x 48`, walls, `A = 0.4`, t = 30	diverges at t = 10.5	0.51, 0.710	0.53, 0.696
`96 x 48 x 48`, symmetry, `A = 0.8`, t = 30	diverges at t = 1.7	0.61, 0.810	0.45, 0.679

Each of the three parts is needed:

Without the reversed-face rule (and still with the corner neighbors), A = 1.2 diverged at t = 13.2.
With only the reversed-face rule, both A = 0.4 cases stayed bounded, but at a Mach number of 0.55 to 0.58 instead of 0.3 to 0.35.
With the corner neighbors, the finer walled box grew the checkerboard: its lowest pressure fell to 0.607 by t = 29.

The mean pressure still drifts, as it does with beta = M on any outlet with a small K (the vortex test below removes it): beta T- averaged over the outlet need not vanish, and only K holds the mean. In the boxes with symmetry sides it rose by 2 to 23% of p_t by t = 25. A larger sigma or a shorter L holds it closer.

The spatially developing channel of the synthetic-inflow work (Re_tau = 180, outlet L = 25 h) showed the same signatures before it diverged:

In the last cell layer, the pressure fluctuation next to both walls was two to three times that of the layer inside. In mid-channel they were equal.
Where it diverged, near the top wall, the last layer held a checkerboard.

Alternatives that did not hold:

Constant beta = 0.2 instead of the local Mach number: it still diverged, at t = 16.4 instead of 11.6.
Backflow stabilization, as for incompressible outlets (Esmaily Moghadam et al. 2011): a traction rho min(u_n, 0) u / 2 on reversed faces, which removes the kinetic energy they let in. It only delayed the divergence, to t = 21.1 instead of 18.5.
Zero tangential velocity entering reversed faces: no effect.
Mirroring the boundary cell across walls into the fit, instead of dropping edge faces: bounded, but the mean pressure of the symmetry boxes rose by 6 to 53% of p_t.

Inlet (`nscbc_inlet`)¶

The target is u_t, T_t and, for mixtures, the composition.

Acoustic wave: the normal velocity relaxes through the face's incoming wave, as above.
Entropy and tangential velocity: these waves are convective (lambda = u_n), so the relaxation dphi/dt = -K_phi (phi - phi_t) is a blend of the transmissive state toward the target:

phi_g = phi_e + min(1, h K_phi / |u_n|) (phi_t - phi_e),   K_phi = sigma_phi c / L

for phi = T or u_t, with h twice the distance from the boundary cell's centroid to the face.

Without sigma_T or sigma_t, T and u_t are imposed exactly. This does not reflect acoustic waves: in linear theory, the reflected acoustic wave depends only on the incoming acoustic condition. - Ghost density: p_g / (R_t T_g), with R_t the gas constant of the target composition. - Composition: imposed exactly, as for upt. It is a convective wave and does not reflect acoustics. A relaxed composition would need the ghost thermodynamics [gamma, e0] of a blended composition at every evaluation. The flow block uses the target's [gamma, e0] on inflow faces.

Multicomponent and reacting flows¶

Species are convective waves. They leave through outlets from the interior, and enter through inlets at the target composition (mass-flux upwinding). Acoustics uses each face's frozen gamma.

Reaction terms. Yoo & Im add the reaction source terms to the LODI system, so that heat release next to a boundary is not mistaken for an incoming wave. Mallard Strang-splits chemistry from the flow:

A reaction half-step changes the boundary cell at constant volume.
The cell relaxes toward the face's incoming wave and sends the pressure rise inward as a wave, as any interior cell does.
The face's own state sees the heat release only through the outgoing wave.

So no separate reaction term is needed. In the flame example (an H2/air flame blown out through the outlet), the pressure stays within 0.7% of the target while the flame leaves, and within 0.3% across the domain. The largest deviation is a uniform dip while the exit velocity falls from 68 to 10 m/s, relaxed on the outlet's time scale.

Sponge layers (`[[sponges]]`)¶

A sponge adds S = -sigma_s(x) (U - U_ref(x)) to the conservative variables, and the matching term to each species partial density (Bodony 2006; Mani 2012).

strength (1/s) is an expression in x, y, z, zero outside the layer. A ramp such as A ((x - x0)/d)^2 keeps its own reflections small.
The reference state is u, p, T (plus X or Y), each a number or an expression.
Strength and reference are evaluated once at cell centroids.
Several sponges add up.
The source is explicit, so the time step is capped at 1 / max sigma_s when a sponge is present.

A sponge works with any boundary condition. The usual setup is a sponge in front of an nscbc_outlet or farfield.

Validation (`test/nscbc_test.cpp`, `test/boundary3d_test.cpp`, MPI tests)¶

Acoustic pulse at normal incidence, default low-Mach correction:
sigma = 0 reflects under 0.25% with MUSCL (quadrilaterals, triangles and jittered mixed meshes) and 1% with TENO5 (quadrilaterals), against 0.6 to 34% for extrapolation.
sigma = 20 reflects as dw-/dt = -K p' predicts, within 3%.
p_out reflects 97 to 99%.
Vortex (M = 0.25, swirl 0.05) leaving through an outlet between walls, against a domain three times longer. RMS pressure error at t = 3:
default: 0.05 times the vortex's pressure deficit, on quadrilaterals, triangles and jittered meshes
beta = 1 (Poinsot & Lele's original condition): 1.0 times (1.2 to 1.3 before the edges dropped their transverse terms)
Poiseuille flow driven by a body force through an outlet: the profile stays within 1.2% and the pressure within 3.5% of the force balance, on quadrilaterals, triangles and jittered meshes.
Gas mixtures:
an acoustic pulse in air leaves (under 2% reflected, against over 90% for p_out)
a hot H2/N2 front leaves without disturbing the pressure by more than 5% of its dynamic pressure
the flame example above
3D: a plane pulse leaves a box of hexahedra or tetrahedra (under 3%).
Eddies reversing the flow at the outlet (the box of "Stability", 24 x 12 x 12 cells, symmetry sides, A = 0.8): bounded to t = 12. Before these changes it diverged at t = 9.8.
MPI: runs with characteristic boundaries and sponges are bitwise equal to serial ones (MUSCL and TENO, 2D and 3D).
Single precision: the same tests pass with round-off noise of about 1e-4 of the pressure.
Unchanged behavior: results with the existing conditions are bitwise unchanged; the examples were compared with the previous build.