3D backward facing step with recycling-plane inflow

Introduction

This tutorial sets up a non-reactive 3D turbulent backward-facing step (BFS) flow using the recycling-plane inflow capability described in PeleLMeX controls. The mean inflow is prescribed in the usual way through the Inflow boundary condition; the new machinery samples the velocity field on a cross-section in the interior of the domain and re-injects the fluctuating part onto the inflow ghost cells every timestep, providing the inlet with self-consistent turbulent forcing without resorting to a synthetic-turbulence library or a precomputed plotfile.

The case demonstrates:

  • setting up a 3D BFS in PeleLMeX with embedded boundaries and a spanwise periodic direction,

  • enabling the recycling-plane inflow,

  • seeding the initial condition with a small-amplitude perturbation to trigger the shear-layer instability,

  • monitoring domain-integrated kinetic energy and enstrophy via the do_temporals machinery to detect statistically stationary turbulence,

  • the practical effect of the source-plane position on convergence speed and statistical quality.

Note that the results here are colored strongly by the fact that the low spatial resolution is unable to produce any significant turbulence that would otherwise be generated by the boundary along the top of the step. The resolved turbulence source in this case is rather produced by the domain expansion past the step, which generates a large-scale coherent vortex that sheds KH fluctuations. If wall boundary-layer turbulence is instead of interest, dramatically higher resolution in space and time is required along the boundary and the sample plane should be located on the step.

The case files referenced below are located in:

Exec/RegTests/EB_BFS_Recycling

Setting-up your environment

Getting a functioning environment in which to compile and run PeleLMeX is the first step of this tutorial. Follow the steps listed below to get to this point:

  1. The first step is to get PeleLMeX and its dependencies. To do so, use a recursive git clone:

    git clone --recursive --shallow-submodules --single-branch https://github.com/AMReX-Combustion/PeleLMeX.git

    The --shallow-submodules and --single-branch flags are recommended for most users as they substantially reduce the size of the download by skipping extraneous parts of the git history. Developers may wish to omit these flags in order download the complete git history of PeleLMeX and its submodules, though standard git commands may also be used after a shallow clone to obtain the skipped portions if needed.

  2. Move into the Exec folder containing your tutorial. To do so:

    cd PeleLMeX/Exec/RegTests/<CaseName>

    where <CaseName> is the name of your tutorial, e.g. HotBubble, FlameSheet, EB_BackwardStepFlame, EB_FlowPastCylinder, or TripleFlame.

You’re good to go!

Note

The makefile system is set up such that default paths are automatically set to the submodules obtained with the recursive git clone, however advanced users can set their own dependencies in the GNUmakefile for each case by updating the top-most lines as follows:

PELE_HOME     = <path_to_PeleLMeX>
AMREX_HOME        = <path_to_MyAMReX>
AMREX_HYDRO_HOME  = <path_to_MyAMReXHydro>
PELE_PHYSICS_HOME = <path_to_MyPelePhysics>
SUNDIALS_HOME     = <path_to_MySUNDIALS>

or directly through shell environment variables (using bash for instance):

export PELE_HOME=<path_to_PeleLMeX>
export AMREX_HOME=<path_to_MyAMReX>
export AMREX_HYDRO_HOME=<path_to_MyAMReXHydro>
export PELE_PHYSICS_HOME=<path_to_MyPelePhysics>
export SUNDIALS_HOME=<path_to_MySUNDIALS>

Note that using the first option will overwrite any environment variables you might have previously defined when using this GNUmakefile.

Case setup

Geometry, grid and boundary conditions

The computational domain is a 0.08 × 0.02 × 0.02 \(m^3\) box with the streamwise direction (\(x\)) aligned with the bulk flow and a periodic spanwise direction (\(z\)). The step occupies the lower-left quarter of the \(x\)-\(y\) plane and is extruded across the full periodic \(z\) range. The physical setup here is identical that used in the BFSFlame tutorial. Boundary conditions are Inflow on \(x_{lo}\), Outflow on \(x_{hi}\), no-slip walls on \(y\), and Interior (periodic) on \(z\).

#---------------------- DOMAIN DEFINITION ------------------------
geometry.is_periodic = 0 0 1
geometry.coord_sys   = 0
geometry.prob_lo     = -0.01  -0.01   0.00
geometry.prob_hi     =  0.07   0.01   0.02

#---------------------- BC FLAGS ---------------------------------
peleLM.lo_bc = Inflow   NoSlipWallAdiab  Interior
peleLM.hi_bc = Outflow  NoSlipWallAdiab  Interior

The base grid is 256 × 64 × 64 with isotropic cell size \(\Delta x = 3.125 \times 10^{-4}\) m. AMR is disabled in this tutorial; max_grid_size = 32 is used to obtain enough boxes to scale onto 32 MPI ranks.

#---------------------- AMR CONTROL ------------------------------
amr.n_cell          = 256 64 64
amr.max_level       = 0
amr.blocking_factor = 16
amr.max_grid_size   = 32

The step is described as an axis-aligned box, extruded across the spanwise direction:

#---------------------- EB SETUP ---------------------------------
eb2.geom_type = box
eb2.box_lo =    -0.01  -0.01  0.00
eb2.box_hi =     0.01   0.00  0.02
eb2.box_has_fluid_inside = 0
eb2.small_volfrac = 1.0e-4

The step height is \(H = 0.01\) m, equal to the post-step channel half- height, giving an expansion ratio of 2:1. The bulk inlet velocity is \(U_{\text{bulk}} = 10\) m/s, so the flow-through time is \(\tau = L_x / U_{\text{bulk}} = 8\) ms and the step Reynolds number is \(\mathrm{Re}_H = U_{\text{bulk}} H / \nu \approx 6700\).

Problem specifications and initial-condition seeding

The user parameters are gathered in MyProbParm:

struct MyProbParm : public ProbParmDefault
{
  amrex::Real T_mean              = 298.0;
  amrex::Real meanFlowMag         = 10.0;
  amrex::Real turb_seed_amplitude = 0.05;
};

  • T_mean : inlet and initial gas temperature

  • meanFlowMag : inlet \(x\)-velocity magnitude

  • turb_seed_amplitude : amplitude of the initial-condition velocity perturbation, expressed as a fraction of meanFlowMag

Recycling on its own does not produce turbulence: if the velocity at the source plane has no fluctuation around its running mean, the injected contribution is zero and the inlet stays steady. To break this trivial fixed point the initial condition adds a small-amplitude, deterministic multi-mode perturbation to all three velocity components, tapered to zero at the wall-normal boundaries:

const amrex::Real ywall = 1.0 - std::pow(2.0 * y / Ly, 2.0);
const amrex::Real wy = std::max(amrex::Real(0.0), ywall);
state(i, j, k, VELX) = U + A * wy * pu * norm;
state(i, j, k, VELY) = 0.0 + A * wy * pv * norm;
state(i, j, k, VELZ) = 0.0 + A * wy * pw * norm;

where A = turb_seed_amplitude * U and pu, pv, pw are sums of three sinusoidal triads with incommensurate wavenumbers and phase offsets (see pelelmex_prob.H for the full expression). The shear-layer instability behind the step amplifies the perturbation into developed turbulence within a few flow-through times; once that turbulence reaches the source plane, the recycling loop closes.

The standard Inflow boundary condition supplies the prescribed mean profile through bcnormal (uniform \(U_{\text{bulk}}\) here); the recycling layer adds the zero-mean fluctuation on top of it.

Recycling-plane inflow controls

The recycling feature is enabled with five inputs (see PeleLMeX controls for the full reference):

#---------------------- Recycling-plane inflow --------------------
peleLM.use_inlet_from_plane      = 1
peleLM.inlet_plane_dir           = 0
peleLM.inlet_plane_position      = 0.065
peleLM.inlet_plane_warmup_steps  = 500
# peleLM.inlet_plane_avg_window  = 1.0e-4   # uncomment for EMA

The source plane is normal to \(x\) and located at \(x_{\text{plane}} = 0.065\) m, which is \(5.5\,H\) downstream of the step end and \(0.5\,H\) from the outflow boundary. The choice of plane position has a strong effect on convergence and statistical quality; this is discussed below.

The 500-sample warmup gives the running mean roughly one flow-through to settle before any fluctuation is injected. With the cumulative mean (no avg_window set) the running average is \(\langle u \rangle_n = \frac{1}{n}\sum_i u_i\); specifying inlet_plane_avg_window switches to an exponential moving average with relaxation timescale equal to the window. The cumulative form is appropriate for one-shot stationary- statistics studies; the EMA form is preferred for cases where the bulk flow is slowly evolving.

Note

The fluctuation is injected on every velocity component on the inflow ghost layer. The implementation requires that all components share the same ext_dir boundary status on the inflow face; mixed BCs on different velocity components are not supported and will abort the simulation.

Numerical parameters

#---------------------- PeleLM CONTROL ---------------------------
peleLM.v               = 1
peleLM.sdc_iterMax     = 2
peleLM.num_init_iter   = 1
peleLM.chem_integrator = "ReactorNull"

amr.cfl              = 0.4
amr.dt_shrink        = 0.01

The case is non-reactive (ReactorNull skips chemistry source evaluation entirely) but the full low-Mach algorithm is otherwise active, so the recycling machinery exercises the same code path it would in a reacting simulation.

num_init_iter = 1 is used because successive initial pressure iterations on this 3D EB geometry can drive the MAC projection divergent on the first real step; a single initial iteration is sufficient to establish a consistent pressure field.

Building the executable

The case uses the non-reactive air chemistry model with Fuego EOS and Simple transport:

DIM             = 3
USE_EB          = TRUE
USE_MPI         = TRUE
Chemistry_Model = air
Eos_Model       = Fuego
Transport_Model = Simple

Build third-party libraries first, then the executable:

make -j4 TPL
make -j4

This produces PeleLMeX3d.gnu.MPI.ex.

Running the simulation

The simulation is non-trivial in cost: each timestep takes around 1.6 s on 32 MPI ranks (~0.05 s/rank/step at this resolution), and reaching fully converged stationary statistics requires of order \(20\,\tau \approx 25\,000\) timesteps, so plan for roughly 8–11 hours of wall time on 32 ranks for a single case.

Launch the run, redirecting stdout to a log file:

mpiexec -n 32 ./PeleLMeX3d.gnu.MPI.ex inputs.3d-recycle peleLM.inlet_plane_position=0.065 > run.log 2>&1

Three things in particular are worth watching during the run:

  1. The dt-shrink ramp during the first ~50 steps: dt starts at \(10^{-2}\) of the convective CFL estimate and grows by the dt_change_max factor each step until it reaches the CFL limit (\(\sim\!8\times 10^{-6}\) s).

  2. The first 500 steps where the warmup gate is active: the snapshot path is exercised every step (and the running mean accumulates) but no fluctuation is yet injected.

  3. The transition to active recycling once the warmup completes: this is visible in the temporals/tempState log as a change of slope in the integrated kinetic energy.

Diagnostics and stationarity detection

The case enables PeleLMeX’s temporal-output machinery to produce a time history of the domain-integrated kinetic energy and enstrophy:

peleLM.do_temporals      = 1
peleLM.temporal_int      = 10
peleLM.do_extremas       = 1
peleLM.do_mass_balance   = 1

This writes one CSV row every 10 steps to temporals/tempState with columns iter, time, dt, kinEnergy, enstrophy, pressure, fuelConsumption, heatRelease (the last two are zero for this non-reactive case).

A simple stationarity check is to compare the mean of kinEnergy over two adjacent windows of width \(W\tau\) ending at candidate time \(t^*\): if the trailing window \([t^* - W\tau,\, t^*]\) and leading window \([t^*,\, t^* + W\tau]\) agree within a tolerance, declare onset at \(t^*\). With \(W=4\) and a 10% tolerance, this case typically identifies onset at \(t^* \approx 4\tau\) and shows post-onset robustness (fraction of subsequent window-pairs that pass the same test) above 50%. The script analyze_stationarity.py shipped with the case reports both numbers, along with the running mean and standard deviation of kinEnergy and enstrophy over the stationary regime.

A representative final-run summary for the \(x_{\text{plane}} = 0.065\) configuration on 25,000 timesteps (\(\sim 23\,\tau\)) is:

  • stationary onset at \(4\,\tau\), robustness 16 / 31 (52%),

  • \(\langle \mathrm{KE} \rangle = 1.04 \times 10^{-3}\) J with CoV 37%,

  • \(\langle \mathrm{enstrophy} \rangle = 4.0 \times 10^{2}\) with CoV 18%.

The kinetic-energy CoV reflects coherent shear-layer shedding in the post-step recovery zone; the enstrophy CoV is much tighter because the small-scale vorticity averages more uniformly across the domain.

Choice of source-plane position

The source-plane location is the single most influential parameter for this feature. The tutorial input shipped in-tree uses \(x_{\text{plane}} = 0.065\). The trends summarized below come from four additional out-of-tree runs in which only peleLM.inlet_plane_position was changed, taking values from {0.05, 0.035, 0.005, 0.0095} to sweep the plane from near the step to farther downstream. The main trends observed in those separate runs are:

  • Plane upstream of the step edge (e.g. \(x = 0.005,\, 0.0095\) ). The flow at the source plane is laminar boundary-layer development of the inflow itself: the running mean equals the instantaneous field to within numerical noise, no fluctuation is injected, and the recycling effectively disables itself. Domain- integrated KE has a coefficient of variation of order 1%, but the flow simply runs on whatever turbulence the IC perturbation seeded through the natural BFS shedding — recycling is not contributing.

  • Plane near the step end (e.g. \(x = 0.035,\, 2.5\,H\) downstream). The plane samples the strongly coherent shear-layer roll-up. Recycling that signal closes a tight feedback loop with the shear-layer mode and produces a large amplitude transient that overshoots and oscillates, with stationary onset pushed out to \(8.5\,\tau\) and a post-onset window-pair pass rate of only 5%.

  • Plane far downstream (e.g. \(x = 0.05,\, 0.065\) ). The plane samples post-transition broadband turbulence with weak large-scale coherence. The recycling loop delay \((x_{\text{plane}} - x_{\text{inlet}})/U_{\text{bulk}}\) is now comparable to a flow-through time, so the injected perturbation is decorrelated from the resulting shear-layer dynamics. Stationary onset is reached around \(4\,\tau\), with runP25 (\(x = 0.065\), \(5.5\,H\) downstream) giving the highest post-onset robustness of the five cases.

The general guidance, then, is to place the source plane in a region of developed broadband turbulence rather than near a region of strong coherent shedding, and to ensure the loop delay is on the order of an integral timescale of the flow. The geometry-specific bounds (sufficient distance from the geometric feature, but not so close to the outflow that the BC interferes with the recycled signal) are the practical constraints; in this case runP25 lands on a usable balance.

Note

Even at \(x_{\text{plane}} = 0.065\) (only 5 mm from the outflow at \(x = 0.07\)) the outflow boundary condition does not detectably contaminate the recycled signal; the gentle outflow extrapolation in fact mildly de-correlates the streamwise structure, which in this setup is helpful. With more reactive or compressibility-sensitive physics the user should re-evaluate that conclusion.

Restarting and continuing the simulation

The running mean of the source plane is written into the checkpoint header (the line beginning with RecyclingPlane:) and the per-level mean fields are written under recycle_mean in each level’s checkpoint directory. On restart, the storage is rebuilt against the current run’s grids and the previously-accumulated mean is ParallelCopy’d back into place, so a checkpoint/restart of an already-stationary simulation continues without reseeding the running mean.

Restart is triggered the usual way:

amr.restart = chk20000
amr.max_step = 30000

Older checkpoints written before the recycling feature was added do not contain the RecyclingPlane block; the read path silently falls back to “no saved mean” and the running mean reseeds from the next snapshot. This makes the feature backwards-compatible with pre-existing checkpoints.