Solver

Configure and run the Navier–Stokes solver.

ember integrates the compressible Navier–Stokes equations to a steady state using an explicit pseudo-time march, accelerated by multigrid and residual smoothing. The finite-volume discretisation is second-order accurate in space on structured multi-block grids; cell-centred residuals are distributed equally back to nodal conserved variables. This page describes the solver loop driven by Solver.run() and the configuration parameters in Solver that control each stage.

Overview of one time step

Solver.run() performs the following operations each step:

1. Cache flush and boundary conditions: update_cached_conserved() refreshes the pressure and temperature caches left stale by the previous step’s in-place march, then update_bconds() updates boundary patch targets (mass-flow throttle, radial equilibrium, mixing-plane exchange) and apply_bconds() imposes those boundary conditions by modifying conserved_nd.

2. NaN check: check_nan() aborts the run early if the conserved state has gone non-finite. Sets a flag on the returned ConvergenceHistory and leaves the invalid field in place for inspection.

3. Source terms: update_sources() rebuilds the viscous forces in F_body_nd. With the Runge–Kutta integrator, n_stage > 0, this runs once every full step with the source terms held constant for all stages. When using the scree integrator, n_stage == 0, the source terms are recomputed every few steps to save cost.

4. Update time step: update_timestep() computes the time step and stores it pre-divided by cell volume in dt_vol_nd. On the steps where stage 3 ran, the same pass adds the cell sources to the body force: the polar coordinates source term needed to balance the radial momentum equation and, when gain_filt is nonzero, the selective-frequency-damping force. When the gain is nonzero it also advances the SFD low-pass filter one step on the time step just computed, every step regardless of the source cadence; the force reads the filter before that update.

5. Residual: update_residual() calculates the unintegrated net-flow residual, with optional implicit residual smoothing , sf_resid.

6. Convergence logging: every n_step_log steps, record_convergence() and format_message() record and print a ConvergenceHistory row using the current residual.

7. March: advance the solution with the selected integrator – Denton’s scree march, scree_step(), or Jameson multi-stage Runge–Kutta, rk_step() – optionally accelerated by multigrid.

8. Smoothing: smooth() applies a constant-coefficient blended second- and fourth-order filter to the post-march ember.block.Block.conserved_nd field, to provide artificial dissipation and suppress odd-even decoupling.

9. Pseudotime averaging: over the final n_step_avg

   steps, accumulate_avg() accumulates the conserved state into a running average, which finalise_average() uses to replace the instantaneous state once the run completes. With n_step_avg of 0 or 1 the whole stage is skipped and the final state stands as the solution.

Time integrators

n_stage selects one of two integrators, applied to every block each step:

scree_step() – n_stage == 0 Implements Denton’s basic scree scheme (two steps forward, one step back)

\[\mathcal{U}_{n+1} = \mathcal{U}_n + \left[ 2 \left(\frac{d\mathcal{U}}{dt}\right)_n - \left(\frac{d\mathcal{U}}{dt}\right)_{n-1} \right] \Delta t\]

The residual from the previous step, \((d\mathcal{U}/dt)_{n-1}\), is kept in store between steps. This scheme is only first-order accurate in pseudotime, unlike Adams–Bashforth, but experience shows it requires less artificial dissipation and is more robust.

rk_step() – n_stage >= 1 – Jameson’s classic Runge–Kutta scheme. Every full step snapshots the starting conserved state \(\mathcal{U}_n\), saving it to store. Then, n_stage substages call advance_rk_stage_mg(), each stage \(k\) marching relative to the starting snapshot but with the residual freshly evaluated on the most recent stage:

\[\begin{split}\mathcal{U}^{(0)} &= \mathcal{U}_n \\ \mathcal{U}^{(k)} &= \mathcal{U}_n + \alpha_k\,\Delta t \left(\frac{d\mathcal{U}}{dt}\right)^{(k-1)}, \qquad k = 1, \dots, m \\ \mathcal{U}_{n+1} &= \mathcal{U}^{(m)}\end{split}\]

for \(m =\) n_stage stages with coefficients \(\alpha_k = 1 / (m - k + 1)\), so the final stage takes the full step \(\alpha_m = 1\).

Both integrators share the Multigrid acceleration and the constant-coefficient Smoothing step; they differ only in the march formula and in how many residual evaluations occur per step. At the end of the step, conserved_nd contains the advanced solution.

CFL number

The integrators of the previous section scale each cell’s residual by cfl and the local volumetric timestep \(\Delta t_\mathrm{vol}\) set by update_timestep(), so the \(\Delta t\) appearing in the march formulae is \(\mathrm{cfl}\,\Delta t_\mathrm{vol}\) per cell. The timestep is the reciprocal of the larger of a convective and a turbulent-diffusion spectral radius (a max-of-directional-radii variant of the JST/Blazek definition):

\[\Delta t_\mathrm{vol} = \frac{1}{\max(\lambda_\mathrm{conv},\,\lambda_\mathrm{diff})}\]

For each direction \(d \in \{i, j, k\}\), with \(\mathbf{S}_d\) the mean of the two opposing face-area vectors, \(\mathbf{V}_\mathrm{rel}\) the relative-frame velocity, \(a\) the speed of sound, \(\mu_t\) the turbulent viscosity, \(\rho\) the density and \(\mathcal{V}\) the cell volume, the directional convective radius and the two combined radii are

\[\begin{split}\Lambda_d &= \left| \mathbf{V}_\mathrm{rel} \cdot \mathbf{S}_d \right| + a \left\| \mathbf{S}_d \right\| \\ \lambda_\mathrm{conv} &= \max_d \Lambda_d \\ \lambda_\mathrm{diff} &= f_\mathrm{visc}\, \frac{\mu_t}{\rho}\, \frac{\max_d \left\| \mathbf{S}_d \right\|^2}{\mathcal{V}}\end{split}\]

Taking the max over directions (rather than Blazek’s sum) makes cfl the true 1D Courant limit, so it stays aspect-ratio-independent for both limits and a single value scales them consistently; fac_visc tightens only the diffusion radius so the viscous march tolerates the same CFL as the inviscid one. Setting \(\mu_t = 0\) drops \(\lambda_\mathrm{diff}\) and recovers the bare convective limit. The 4-stage Runge–Kutta march is stable up to \(\mathrm{cfl} \approx 2\sqrt{2}\); the scree scheme is stable up to \(\mathrm{cfl} \approx 0.6\) on a uniform mesh.

cfl is a single constant applied uniformly to every cell for the entire run – there is no per-cell adaptive CFL field and no tolerance-driven backoff inside the solver. A larger CFL converges faster but risks divergence; implicit residual smoothing (sf_resid) damps high-frequency residual content and so tolerates a substantially higher CFL for a given scheme.

Smoothing

A constant-coefficient blend of second- and fourth-difference operators, ember.grid.Grid.smooth(), is applied to ember.block.Block.conserved_nd after each step to suppress odd–even decoupling and high-frequency content introduced by the march and multigrid corrections. sf4 and sf2 are coefficients on the fourth- and second-difference terms, each scaled by the run’s cfl to make the effective dissipation independent of the time step.

adaptive_smoothing swaps the constant second-order factor for a JST normalised-curvature sensor on pressure and temperature, clipping the fourth-order factor against it as max(sf4 - sf2n, 0). Fourth-order dissipation then switches off inside a shock, where it would otherwise oscillate, and second-order dissipation appears only where the sensor fires. The stencil family and the biased boundary closures are shared with the constant-coefficient kernel, so sf4 and sf2 keep their meaning across the switch.

Multigrid

In-step Denton block-sum multigrid. Coarse-grid corrections are computed in place within a single march call over one grid – coarse block-sum corrections are folded directly into the fine-grid increment before it is scattered onto the conserved state – rather than as a classical restrict/prolong V-cycle across separate coarse grids. Both integrators honor two knobs:

• n_levels – number of coarse levels; 0 disables multigrid. Each block’s cell counts must be an exact multiple of the coarsest block size 2**n_levels in every direction, or Solver.run() raises before marching.

• fac_mgrid – scaling on the coarse correction; 0 also disables multigrid.

• expon_mgrid – base of the per-level geometric decay applied on top of fac_mgrid, so successively coarser levels are damped further (coef_l ~ expon_mgrid**-(l-1)).

sf_resid additionally drives implicit residual smoothing (Jameson IRS) on the fine-grid residual via update_residual(), independent of whether multigrid is enabled.

Full multigrid startup. Solver.run_fmg() runs the same solver coarse-to-fine as a startup schedule, rather than within a single step. It builds n_levels progressively-halved grids (resample()), solves the coarsest first, then prolongs each converged solution onto the next finer grid as its initial guess (interp_from_grid()) and calls Solver.run() again with the in-step multigrid depth set to that level’s index – so the coarsest level runs with no in-step multigrid and the finest runs at the full requested n_levels, identical to calling Solver.run() directly on the finest grid. With n_levels <= 0 it reduces to a single call to Solver.run().

Soft start. BaseSolver.soft() returns a detuned copy of a solver – a lower CFL, heavier second-difference smoothing, the increment limiter on, no coarse correction – whose converged answer is an initial guess the production settings can survive. The detuned values are one fixed known-robust configuration, not a factor on what came in, and the march length and the model choices are left alone: it is the same problem on the same grid. It is a configuration transform, not a march, so the schedule is written at the call site (solver.soft().run(grid) then solver.run(grid), the grid carrying the guess between them in place). It is compulsory on the interface because only each solver knows what detuning means for its own settings, and orthogonal to the grid hierarchy above: solver.soft().run_fmg(grid) softens every level of a full-multigrid startup.

Body forces and viscous model

The cell-centred body-force buffer (block.F_body_nd) accumulates all source terms before they are added to the residual. It is rebuilt in two passes: ember.grid.Grid.update_sources() zeroes it and adds the viscous terms, then ember.grid.Grid.update_timestep() with add_sources=True adds the cell sources from the cell averages its timestep walk already forms:

• Viscous shear stresses and heat flux, computed unless Solver.inviscid is set. The viscous pass is phased across the whole grid (tau/q on every block, then a periodic-seam halo exchange, then face-flux accumulation) so block-to-block periodic interfaces stay consistent.

• A polar (axisymmetric) source term to balance the cylindrical coordinate metric.

• An optional selective-frequency-damping (SFD) force when Solver.gain_filt is nonzero.

The mixing-length turbulent viscosity uses a fixed turbulent Prandtl number of 1.0 and is evaluated from the absolute-frame vorticity magnitude. No-slip walls take their shear stress from a wall function, a skin-friction law in the Reynolds number of the wall-adjacent cell, chosen by Solver.wall_law. Solver.fac_visc multiplies the turbulent-diffusion timestep radius independently of this, tightening the viscous stability limit to recover the inviscid stable CFL where needed.

Boundary patches and inter-block coupling

Inlet, outlet, and mixing-plane patches each relax their own state towards a target every step, with their own relaxation factor rather than a single solver-wide setting:

• InletPatch and OutletPatch take one under-relaxed step of the characteristic condition per timestep, scaled by sigma; see rf_inlet and rf_outlet.

• MixingPatch takes the same under-relaxed characteristic step as the inlet and outlet, scaled by its own sigma; see rf_mix.

• MixingCommunicator relaxes the mixing-plane target exchanged between adjacent blocks with the patches’ rf_exchange, separately from either side’s own step. mix_reflective replaces that whole exchange, and both sides’ characteristic steps, with a direct mixed-out state.

• OutletPatch relaxes its spanwise radial-equilibrium profile separately, via set_adjustment(rf=...), and damps its mass-flow throttle separately again, via the dimensionless gains of set_throttle(mdot_target, Kp=..., Ki=...).

ember.grid.Grid.update_bconds() advances the slowly-varying boundary targets once per step (mixing-plane exchange, characteristic mean state, outlet throttle and spanwise target); ember.grid.Grid.apply_bconds() then imposes the full set of physical boundary conditions and closes periodic seams every time it is called, including between Runge–Kutta substages.

Logging, averaging, and convergence history

Convergence diagnostics are recorded into a ConvergenceHistory every Solver.n_step_log steps: mean residual, mass flow / stagnation enthalpy / entropy at row interfaces, and outlet throttle state (ember.grid.Grid.get_convergence()).

Pseudotime averaging of the conserved variables accumulates over the final Solver.n_step_avg steps of the run (ember.grid.Grid.accumulate_avg()). On completion, the time-averaged state replaces the instantaneous state (ember.grid.Grid.finalise_average()) – skipped if the run diverged, so the invalid field is preserved for inspection rather than overwritten by a partially-accumulated average, and skipped when Solver.n_step_avg is 0 or 1, where the mean of one sample is the sample and no accumulator is allocated.

class solver.BaseSolver

Common interface for in-place flow solvers.

A solver is constructed from its configuration and run with solver.run(grid), which advances grid in place and returns the convergence history. Both the built-in Solver and external-solver plugins (e.g. ember.plugins.ts.TS3Solver) implement this contract so they are drop-in interchangeable.

abstractmethod run(grid)

Solve on grid in place; return a ConvergenceHistory.

abstractmethod soft()

Return a copy of this solver configured for a soft start.

A soft configuration trades convergence rate for robustness: the same solver, detuned enough to survive an initial guess too far from the solution for the production settings to tolerate. It solves the same problem on the same grid, so its answer is a legitimate – if expensive to reach and slack – initial guess for a run of self:

solver.soft().run(grid)   # robust startup, in place
hist = solver.run(grid)   # production settings from there

class solver.Solver(n_step, n_step_log=10, n_step_avg=1, cfl=5.0, sf4=0.008, sf2=0.002, sf2_T=None, adaptive_smoothing=False, inviscid=False, fac_visc=1.0, wall_law='fit', sf_resid=1.0, gain_filt=0.0, delta_filt=1.0, n_stage=0, n_levels=3, fac_mgrid=0.2, expon_mgrid=1.414, dampin=0.0, rf_inlet=0.05, rf_outlet=0.05, rf_mix=0.01, rf_exchange=0.01, mix_endwall_span=0.0, mix_endwall_rf=0.1, mix_reflective=False)

Configuration for the explicit time-marching solver.

Also the entry point: build one with the parameters below and call run() (or run_fmg()) to march a grid in place.

Frozen, because a solver is a set of parameters rather than a thing with state: nothing here is written after construction, the march keeps its working state on the grid, and run_fmg() already derives its per-level configurations with dataclasses.replace(). Settings that would otherwise be adjusted in place are made by building another one.

n_step: int

n_step_log: int = 10

Number of steps between convergence log messages.

n_step_avg: int = 1

Number of steps at the end of the march to average the solution over.

0 and 1 both mean no averaging: a one-sample mean is the sample, so the march skips the accumulator entirely and leaves the final state as the solution. Must not exceed n_step — see __post_init__().

cfl: float = 5.0

Constant CFL number for the march

sf4: float = 0.008

Fourth-order smoothing factor.

sf2: float = 0.002

Second-order smoothing factor.

Constant everywhere unless adaptive_smoothing is set, in which case it is the ceiling a JST shock sensor scales: the weight on the pressure sensor, and on the temperature sensor too unless sf2_T says otherwise.

sf2_T: float | None = None

Weight on the temperature limb of the adaptive sensor.

None takes sf2, so both limbs share one weight, as they always have; zero leaves the sensor on pressure alone, as classical JST has it. Read only when adaptive_smoothing is set.

adaptive_smoothing: bool = False

Drive the second-order smoothing with a JST shock sensor.

Replaces the constant second-order factor with a normalised-curvature sensor on pressure and temperature, and clips the fourth-order factor against it, so fourth-order dissipation switches off inside a shock and second-order dissipation appears only where it is needed. Same stencil family and same boundary closures as the constant-coefficient kernel, so sf2 and sf4 keep their meaning.

inviscid: bool = False

Skip viscous terms in the sources evaluation.

fac_visc: float = 1.0

Multiplier on the turbulent-diffusion timestep radius; >1 tightens the viscous limit to recover the inviscid stable CFL.

wall_law: str = 'fit'

Skin-friction law at no-slip walls, "fit" or "reichardt".

Both give the wall shear stress from the Reynolds number of the wall-adjacent cell, on the velocity at the first off-wall node. "fit" is a curve fit in ln(Re), switching to the laminar cf = 2/Re at y+ = 11.3; it follows the log law to a few percent above y+ = 30 but reads up to 25% low in the buffer layer just above the switch. "reichardt" inverts Reichardt’s law, one expression from the viscous sublayer through the buffer layer to the log layer, at the cost of two Newton steps per wall face. Pass the same law to ember.block_util.wall_yplus() for consistent y+.

sf_resid: float = 1.0

Implicit residual smoothing factor. Applied to the fine residual by update_residual() (sf) and, on both integrators, to the coarse block-restricted residual of the multigrid correction (advance_rk_stage_mg()’s sf_irs for RK, scree_step()’s sf_irs for scree). The coarse smoothing needs n_levels > 0 to have any effect.

gain_filt: float = 0.0

Selective frequency damping gain.

delta_filt: float = 1.0

Selective frequency damping filter width (higher is smoother).

n_stage: int = 0

Number of time integration stages per step. 0 for scree, >=1 for RK.

n_levels: int = 3

Number of coarse multigrid levels; 0 disables multigrid. Honored by both integrators (scree_step() and rk_step()).

fac_mgrid: float = 0.2

Scaling factor on multigrid corrections. Honored by both integrators (scree_step() and rk_step()).

expon_mgrid: float = 1.414

Base of the per-level multigrid decay, coef_l ~ expon_mgrid**-(l-1). Honored by both integrators (scree_step() and rk_step()).

dampin: float = 0.0

Negative-feedback change limiter (multall’s DAMP); 0 disables it.

Applied to the ASSEMBLED increment – fine term plus injected coarse multigrid correction – immediately before it reaches the nodes. Each cell’s increment is soft-clipped to dU / (1 + |dU| / (dampin * mean|dU|)), per conserved variable, so cells far above their block mean are pulled back towards dampin times it.

rf_inlet: float | None = 0.05

Characteristic under-relaxation (sigma) on every InletPatch. Imposed on every such patch at the start of the run, so the default overrides a value the patches carried in; pass None to leave whatever they already hold.

rf_outlet: float | None = 0.05

As rf_inlet, for every OutletPatch. This is the characteristic relaxation only; the spanwise radial-equilibrium profile has its own, set via set_adjustment(rf=...).

rf_mix: float | None = 0.01

As rf_inlet, for every MixingPatch. This is each side’s own characteristic relaxation; rf_exchange is the separate factor on the cross-plane exchange between them.

rf_exchange: float | None = 0.01

Relaxation of the cross-plane mismatch on every MixingPatch. Read from the patches by MixingCommunicator at each exchange. As rf_inlet, the default is imposed and None leaves each plane’s own value alone.

mix_endwall_span: float | None = 0.0

Fraction of span, measured in from hub and casing, over which every MixingPatch relaxes the exchange more gently and ramps out the harmonic part of its own correction, down to the pitchwise mean alone at the wall.

The endwall treatment of SU2’s GILES_EXTRA_RELAXFACTOR. The exchange is linearised about one state per span station, and an endwall station is the first where the two sides stop being near one: a corner separation can leave them streaming away from the plane in opposite directions, with a mean that cancels to nothing between them. Zero is off; must be below one half. As rf_inlet, the default is imposed and None leaves each plane’s own value alone.

mix_endwall_rf: float | None = 0.1

Fraction of rf_exchange applied at the endwall itself, rising linearly to all of it mix_endwall_span in from the wall. Relative rather than absolute, unlike SU2’s, so it keeps its meaning when rf_exchange is retuned. Must lie in (0, 1]; ignored while mix_endwall_span is zero.

mix_reflective: bool | None = False

Run every mixing plane as a reflective one, imposing the mixed-out state directly instead of the characteristic exchange.

The default plane is the whole of Saxer and Giles [2]: the cross-plane mismatch is split by direction of propagation, relaxed onto a target with rf_exchange, and the target drives only the mean mode of a boundary condition that stays non-reflecting to the harmonics. This replaces all of that with the simplest thing that couples two rows – at every span station the conserved variables of both faces are set to the average of the two sides’ circumferential means, the whole face reset to it on every application.

What is given up is accuracy at the plane. Every pitchwise harmonic is annihilated at the boundary node rather than absorbed, so the plane reflects; the state imposed is the area average of the conserved variables, which preserves each row’s mass flow exactly but not its momentum or energy flux, so it is not the flux-conserving mixed-out state a loss audit would want; and with no under-relaxation anywhere the two rows are yanked onto their common mean every stage. rf_mix and rf_exchange address machinery this switches off, and so do nothing while it is set.

What is not given up is conservation. The inviscid face flow is built from the four boundary nodes of each face quad and the face area vector alone, so two faces carrying the same pitch-uniform state pass identical mass, meridional momentum, angular momentum and energy per unit annulus, whatever their blade counts, pitchwise resolutions and rotational speeds.

As rf_inlet, this is imposed on every mixing patch of every level at the start of the run, so the default overrides whatever a grid was pickled with, and both sides of every plane necessarily agree; pass None to leave each plane as it is. A grid with no mixing plane is unaffected.

run(grid)

Drive grid through n_step steps in place; return the history.

The public BaseSolver entry point for the stage-by-stage march.

run_fmg(grid)

Full-multigrid startup on grid in place.

Returns a list of per-level ConvergenceHistory, coarsest first. Not part of the BaseSolver contract (plugins have no FMG analogue).

soft()

Return a copy detuned for a robust start; see BaseSolver.soft().

Returns:

A new robust configuration.

Return type:

Solver

solver.scree_step(grid, cfl, fac_mgrid=0.0, expon_mgrid=2.0, n_levels=0, sf_irs=0.0, dampin=0.0)

Advance every block one Denton scree step in place.

solver.mg_coarse_shapes(ni, nj, nk, n_levels)

Shapes of the multigrid kernels’ seven scratch buffers, in MG_COARSE_NAMES order.

Separated from the carve so a CALLER can fold these into the single util.carve_view that also carves its own buffers. That matters because the multigrid scratch and the caller’s increment buffer are live at the same time and come from the same arena (Block.scratch): carving them together is what makes them provably disjoint rather than disjoint by convention. It is also what the arena sizing is computed from, so the sizing and the carve cannot disagree.

solver.advance_rk_stage_mg(grid, alpha, cfl, fac_mgrid, n_levels, expon_mgrid=2.0, sf_irs=0.0, dampin=0.0)

One Jameson RK stage, optionally with Denton block-sum multigrid.

Each stage marches every block off its step-start snapshot (block.store, seeded by the caller) using the previous stage’s residual, assembling in one pass a cell-centred increment – fine RK term plus injected coarse – thus:

dU_cell = alpha*cfl*dt_vol*residual                          (fine)
        + sum_l  inject_l( coef_l * dt_coarse_l * restrict_l(residual) )
cons    = snapshot + cell_to_node(dU_cell)

n_levels counts coarse levels only; at the default 0 the coarse loop is empty and this is a plain Jameson RK step. Level l has block size b = 2**l and coef_l = alpha*cfl*fac_mgrid/b**2 * expon_mgrid**-(l-1), so coarser levels damp further. dt_coarse_l is the volume-weighted HARMONIC mean of dt_vol, the block wanting 1/<Lambda> per cell.

Restriction is hierarchical, only level 1 reading the fine grid; prolongation is injection, the transpose of the block-sum restriction on any mesh, collapsed in place inside corr_all so the correction rides the fine term’s scatter, differing from it only at block faces.

The per-block body is one fused Fortran kernel (rk_plain with no coarse levels, else rk_mg_irs/rk_mg_noirs over the shared engine mg_restrict_levels), taking all its scratch from block.scratch in ONE carve_view, which is what makes the pieces disjoint. sf_irs (default 0, off) smooths the coarse residual as the caller already smoothed the fine.

The scatter reads block.store and writes conserved_nd directly (frozen pressure, bypassing the P/T cache). Unmasked: grid.apply_bconds re-imposes the targets between stages. Assumes dt_vol_nd and residual_nd are populated.

solver.rk_step(grid, conf)

Advance every block one Jameson multi-stage RK step in place.