.. highlight:: rst .. _sec:tutorialBFSRecycling: 3D backward facing step with recycling-plane inflow =================================================== .. _sec:TUTO_BFSREC::Intro: Introduction ------------ This tutorial sets up a non-reactive 3D turbulent backward-facing step (BFS) flow using the `recycling-plane inflow` capability described in :doc:`LMeXControls`. 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 .. _sec:TUTO_BFSREC::PrepStep: .. include:: Tutorials_SettingUp.rst Case setup ---------- Geometry, grid and boundary conditions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The computational domain is a 0.08 × 0.02 × 0.02 :math:`m^3` box with the streamwise direction (:math:`x`) aligned with the bulk flow and a periodic spanwise direction (:math:`z`). The step occupies the lower-left quarter of the :math:`x`-:math:`y` plane and is extruded across the full periodic :math:`z` range. The physical setup here is identical that used in the BFSFlame tutorial. Boundary conditions are ``Inflow`` on :math:`x_{lo}`, ``Outflow`` on :math:`x_{hi}`, no-slip walls on :math:`y`, and ``Interior`` (periodic) on :math:`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 :math:`\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 :math:`H = 0.01` m, equal to the post-step channel half- height, giving an expansion ratio of 2:1. The bulk inlet velocity is :math:`U_{\text{bulk}} = 10` m/s, so the flow-through time is :math:`\tau = L_x / U_{\text{bulk}} = 8` ms and the step Reynolds number is :math:`\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 :math:`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 :math:`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 :doc:`LMeXControls` 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 :math:`x` and located at :math:`x_{\text{plane}} = 0.065` m, which is :math:`5.5\,H` downstream of the step end and :math:`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 :math:`\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 :math:`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: #. The dt-shrink ramp during the first ~50 steps: ``dt`` starts at :math:`10^{-2}` of the convective CFL estimate and grows by the ``dt_change_max`` factor each step until it reaches the CFL limit (:math:`\sim\!8\times 10^{-6}` s). #. 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. #. 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 :math:`W\tau` ending at candidate time :math:`t^*`: if the trailing window :math:`[t^* - W\tau,\, t^*]` and leading window :math:`[t^*,\, t^* + W\tau]` agree within a tolerance, declare onset at :math:`t^*`. With :math:`W=4` and a 10\% tolerance, this case typically identifies onset at :math:`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 :math:`x_{\text{plane}} = 0.065` configuration on 25,000 timesteps (:math:`\sim 23\,\tau`) is: * stationary onset at :math:`4\,\tau`, robustness 16 / 31 (52\%), * :math:`\langle \mathrm{KE} \rangle = 1.04 \times 10^{-3}` J with CoV 37\%, * :math:`\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 :math:`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.** :math:`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.** :math:`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 :math:`8.5\,\tau` and a post-onset window-pair pass rate of only 5\%. * **Plane far downstream (e.g.** :math:`x = 0.05,\, 0.065` **).** The plane samples post-transition broadband turbulence with weak large-scale coherence. The recycling loop delay :math:`(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 :math:`4\,\tau`, with ``runP25`` (:math:`x = 0.065`, :math:`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 :math:`x_{\text{plane}} = 0.065` (only 5 mm from the outflow at :math:`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.