r"""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 :meth:`Solver.run` and the configuration parameters in :class:`Solver` that
control each stage.
Overview of one time step
--------------------------
:meth:`Solver.run` performs the following operations each step:
1. **Cache flush and boundary conditions**:
:meth:`~ember.grid.Grid.update_cached_conserved` refreshes the pressure and
temperature caches left stale by the previous step's in-place march, then
:meth:`~ember.grid.Grid.update_bconds` updates boundary patch targets
(mass-flow throttle, radial equilibrium, mixing-plane exchange) and
:meth:`~ember.grid.Grid.apply_bconds` imposes those boundary conditions by
modifying :attr:`~ember.block.Block.conserved_nd`.
2. **NaN check**: :meth:`~ember.grid.Grid.check_nan` aborts the run early if
the conserved state has gone non-finite. Sets a flag on the returned
:class:`~ember.convergence_history.ConvergenceHistory` and leaves the
invalid field in place for inspection.
3. **Source terms**: :meth:`~ember.grid.Grid.update_sources` rebuilds the
viscous forces in :attr:`~ember.block.Block.F_body_nd`. With the
Runge--Kutta integrator, :attr:`~Solver.n_stage` ``> 0``, this runs once
every full step with the source terms held constant for all stages. When
using the scree integrator, :attr:`~Solver.n_stage` ``== 0``, the source
terms are recomputed every few steps to save cost.
4. **Update time step**: :meth:`~ember.grid.Grid.update_timestep` computes
the time step and stores it pre-divided by cell volume in
:attr:`~ember.block.Block.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
:attr:`~Solver.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**: :meth:`~ember.grid.Grid.update_residual` calculates the
unintegrated net-flow residual, with optional implicit residual smoothing
, :attr:`~Solver.sf_resid`.
6. **Convergence logging**: every :attr:`~Solver.n_step_log` steps,
:meth:`~ember.convergence_history.ConvergenceHistory.record_convergence`
and :meth:`~ember.convergence_history.ConvergenceHistory.format_message`
record and print a :class:`~ember.convergence_history.ConvergenceHistory`
row using the current residual.
7. **March**: advance the solution with the selected integrator -- Denton's
scree march, :func:`scree_step`, or Jameson multi-stage Runge--Kutta,
:func:`rk_step` -- optionally accelerated by multigrid.
8. **Smoothing**: :meth:`~ember.grid.Grid.smooth` applies a
constant-coefficient blended second- and fourth-order filter to the
post-march :attr:`ember.block.Block.conserved_nd` field, to provide artificial dissipation and
suppress odd-even decoupling.
9. **Pseudotime averaging**: over the final :attr:`~Solver.n_step_avg`
steps, :meth:`~ember.grid.Grid.accumulate_avg` accumulates the conserved
state into a running average, which
:meth:`~ember.grid.Grid.finalise_average` uses to replace the
instantaneous state once the run completes. With
:attr:`~Solver.n_step_avg` of 0 or 1 the whole stage is skipped and the
final state stands as the solution.
.. _march-schemes:
Time integrators
-----------------
:attr:`~Solver.n_stage` selects one of two integrators, applied to every
block each step:
:func:`scree_step` -- :attr:`~Solver.n_stage` ``== 0``
Implements Denton's basic scree scheme (two steps forward, one step back)
.. math::
\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, :math:`(d\mathcal{U}/dt)_{n-1}`, is kept in
:attr:`~ember.block.Block.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.
:func:`rk_step` -- :attr:`~Solver.n_stage` ``>= 1``
-- Jameson's classic Runge--Kutta scheme. Every full step
snapshots the starting conserved state
:math:`\mathcal{U}_n`, saving it to :attr:`~ember.block.Block.store`.
Then, :attr:`~Solver.n_stage` substages call
:func:`advance_rk_stage_mg`, each stage :math:`k` marching relative to the starting snapshot but with the residual freshly evaluated on the most recent stage:
.. math::
\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)}
for :math:`m =` :attr:`~Solver.n_stage` stages with coefficients
:math:`\alpha_k = 1 / (m - k + 1)`, so the final stage takes the full step
:math:`\alpha_m = 1`.
Both integrators share the :ref:`multigrid` acceleration and the constant-coefficient
:ref:`smoothing` step; they differ only in the march formula and in how many
residual evaluations occur per step. At the end of the step, :attr:`~ember.block.Block.conserved_nd` contains the advanced solution.
.. _cfl:
CFL number
----------
The integrators of the previous section scale each cell's residual by
:attr:`~Solver.cfl` and the local volumetric timestep
:math:`\Delta t_\mathrm{vol}` set by :meth:`~ember.grid.Grid.update_timestep`,
so the :math:`\Delta t` appearing in the march formulae is
:math:`\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):
.. math::
\Delta t_\mathrm{vol} = \frac{1}{\max(\lambda_\mathrm{conv},\,\lambda_\mathrm{diff})}
For each direction :math:`d \in \{i, j, k\}`, with :math:`\mathbf{S}_d` the mean
of the two opposing face-area vectors, :math:`\mathbf{V}_\mathrm{rel}` the
relative-frame velocity, :math:`a` the speed of sound, :math:`\mu_t` the
turbulent viscosity, :math:`\rho` the density and :math:`\mathcal{V}` the cell
volume, the directional convective radius and the two combined radii are
.. math::
\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}}
Taking the max over directions (rather than Blazek's sum) makes
:attr:`~Solver.cfl` the true 1D Courant limit, so it stays
aspect-ratio-independent for both limits and a single value scales them
consistently; :attr:`~Solver.fac_visc` tightens only the diffusion radius
so the viscous march tolerates the same CFL as the inviscid one. Setting
:math:`\mu_t = 0` drops :math:`\lambda_\mathrm{diff}` and recovers the bare
convective limit. The 4-stage Runge--Kutta march is stable up to
:math:`\mathrm{cfl} \approx 2\sqrt{2}`; the scree scheme is stable up to :math:`\mathrm{cfl} \approx 0.6` on a uniform mesh.
:attr:`~Solver.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 (:attr:`~Solver.sf_resid`)
damps high-frequency residual content and so tolerates a substantially higher
CFL for a given scheme.
.. _smoothing:
Smoothing
---------
A constant-coefficient blend of second- and fourth-difference operators, :meth:`ember.grid.Grid.smooth`, is applied to :attr:`ember.block.Block.conserved_nd` after
each step to suppress odd--even decoupling and high-frequency content
introduced by the march and multigrid corrections.
:attr:`~Solver.sf4` and :attr:`~Solver.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.
:attr:`~Solver.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 :attr:`~Solver.sf4` and :attr:`~Solver.sf2`
keep their meaning across the switch.
.. _multigrid:
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:
- :attr:`~Solver.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 :meth:`Solver.run` raises before
marching.
- :attr:`~Solver.fac_mgrid` -- scaling on the coarse correction; 0 also
disables multigrid.
- :attr:`~Solver.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)``).
:attr:`~Solver.sf_resid` additionally drives implicit residual smoothing
(Jameson IRS) on the fine-grid residual via
:meth:`~ember.grid.Grid.update_residual`, independent of whether multigrid is
enabled.
**Full multigrid startup.** :meth:`Solver.run_fmg` runs the same solver
coarse-to-fine as a startup schedule, rather than within a single step. It
builds :attr:`~Solver.n_levels` progressively-halved grids
(:meth:`~ember.grid.Grid.resample`), solves the coarsest first, then
prolongs each converged solution onto the next finer grid as its initial
guess (:meth:`~ember.grid.Grid.interp_from_grid`) and calls :meth:`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 :attr:`~Solver.n_levels`, identical to calling :meth:`Solver.run`
directly on the finest grid. With ``n_levels <= 0`` it reduces to a single call
to :meth:`Solver.run`.
**Soft start.** :meth:`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:
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: :meth:`ember.grid.Grid.update_sources` zeroes it and adds the viscous
terms, then :meth:`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-coupling:
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:
- :class:`~ember.patch.InletPatch` and :class:`~ember.patch.OutletPatch` take
one under-relaxed step of the characteristic condition per timestep, scaled
by :attr:`~ember.patch.NonReflectingPatch.sigma`; see
:attr:`~ember.solver.Solver.rf_inlet` and :attr:`~ember.solver.Solver.rf_outlet`.
- :class:`~ember.patch.MixingPatch` takes the same under-relaxed
characteristic step as the inlet and outlet, scaled by its own
:attr:`~ember.patch.NonReflectingPatch.sigma`; see
:attr:`~ember.solver.Solver.rf_mix`.
- :class:`~ember.mixing_communicator.MixingCommunicator` relaxes the
mixing-plane target exchanged between adjacent blocks with the patches'
``rf_exchange``, separately from either side's own step.
:attr:`~ember.solver.Solver.mix_reflective` replaces that whole exchange,
and both sides' characteristic steps, with a direct mixed-out state.
- :class:`~ember.patch.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=...)``.
:meth:`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); :meth:`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-and-averaging:
Logging, averaging, and convergence history
---------------------------------------------
Convergence diagnostics are recorded into a
:class:`~ember.convergence_history.ConvergenceHistory` every
``Solver.n_step_log`` steps: mean residual, mass flow / stagnation
enthalpy / entropy at row interfaces, and outlet throttle state
(:meth:`ember.grid.Grid.get_convergence`).
Pseudotime averaging of the conserved variables accumulates over the final
``Solver.n_step_avg`` steps of the run
(:meth:`ember.grid.Grid.accumulate_avg`). On completion, the time-averaged
state replaces the instantaneous state
(:meth:`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.
"""
import logging
from abc import ABC, abstractmethod
from dataclasses import dataclass, replace
import ember
from ember import util
from ember.convergence_history import ConvergenceHistory
from ember.grid import DivergenceError
logger = logging.getLogger(__name__)
def _log_rss(what, *args):
"""Debug-log process RSS and its high-water mark at a labelled point.
``what`` is a printf-style label, formatted with ``args`` only when debug
logging is on. The high-water mark comes from the kernel (``VmHWM``), so it
catches transient allocations that came and went between two probes --
which is the whole point of scattering these through the march.
"""
if not logger.isEnabledFor(logging.DEBUG):
return
rss, peak = util.rss_bytes()
logger.debug(
"rss: %s -- %.1f MB resident, %.1f MB peak",
what % args if args else what,
rss / 1024**2,
peak / 1024**2,
)
[docs]
class BaseSolver(ABC):
"""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 :class:`Solver` and external-solver
plugins (e.g. ``ember.plugins.ts.TS3Solver``) implement this contract so
they are drop-in interchangeable.
"""
[docs]
@abstractmethod
def run(self, grid):
"""Solve on ``grid`` in place; return a ``ConvergenceHistory``."""
raise NotImplementedError
[docs]
@abstractmethod
def soft(self):
"""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``:
.. code-block:: python
solver.soft().run(grid) # robust startup, in place
hist = solver.run(grid) # production settings from there
"""
raise NotImplementedError
[docs]
@dataclass(frozen=True)
class Solver(BaseSolver):
"""Configuration for the explicit time-marching solver.
Also the entry point: build one with the parameters below and call
:meth:`run` (or :meth:`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 :meth:`run_fmg` already derives its
per-level configurations with :func:`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 :attr:`n_step` --- see :meth:`__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 :attr:`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 :attr:`sf2_T` says
otherwise."""
sf2_T: float | None = None
"""Weight on the temperature limb of the adaptive sensor.
None takes :attr:`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 :attr:`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
:attr:`sf2` and :attr:`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
:func:`ember.block_util.wall_yplus` for consistent y+."""
sf_resid: float = 1.0
"""Implicit residual smoothing factor. Applied to the fine residual by
:meth:`~ember.grid.Grid.update_residual` (``sf``) and, on both integrators,
to the coarse block-restricted residual of the multigrid correction
(:func:`advance_rk_stage_mg`'s ``sf_irs`` for RK, :func:`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 (:func:`scree_step` and :func:`rk_step`)."""
fac_mgrid: float = 0.2
"""Scaling factor on multigrid corrections. Honored by both integrators
(:func:`scree_step` and :func:`rk_step`)."""
expon_mgrid: float = 1.414
"""Base of the per-level multigrid decay, ``coef_l ~ expon_mgrid**-(l-1)``.
Honored by both integrators (:func:`scree_step` and :func:`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
(:attr:`~ember.patch.NonReflectingPatch.sigma`) on every
:class:`~ember.patch.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 :attr:`rf_inlet`, for every :class:`~ember.patch.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 :attr:`rf_inlet`, for every :class:`~ember.patch.MixingPatch`. This
is each side's own characteristic relaxation; :attr:`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
:class:`~ember.patch.MixingPatch`. Read from the patches by
:class:`~ember.mixing_communicator.MixingCommunicator` at each exchange.
As :attr:`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
:class:`~ember.patch.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 :attr:`rf_inlet`, the default is imposed and None leaves each
plane's own value alone."""
mix_endwall_rf: float | None = 0.1
"""Fraction of :attr:`rf_exchange` applied at the endwall itself, rising
linearly to all of it :attr:`mix_endwall_span` in from the wall. Relative
rather than absolute, unlike SU2's, so it keeps its meaning when
:attr:`rf_exchange` is retuned. Must lie in ``(0, 1]``; ignored while
:attr:`mix_endwall_span` is zero."""
mix_reflective: bool | None = False
r"""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 :cite:t:`Saxer1993`: the cross-plane
mismatch is split by direction of propagation, relaxed onto a target with
:attr:`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. :attr:`rf_mix` and
:attr:`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 :attr:`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."""
def __post_init__(self):
"""Reject averaging windows the march cannot honour, and unknown wall laws.
:attr:`n_step_avg` counts the steps at the end of the march that
:meth:`~ember.grid.Grid.accumulate_avg` sums into the pseudotime
average, dividing each sample by it. A negative window has no meaning.
One longer than the march is worse than meaningless: the loop can only
add ``n_step`` samples but the accumulator still divides by
``n_step_avg``, so the solution comes back scaled by
``n_step / n_step_avg`` --- silently, because a scaled field is every
bit as finite as the real one and neither the divergence check nor the
convergence history has any way to notice. Caught at construction
rather than after a march has been paid for.
An unknown :attr:`wall_law` is caught here for the same reason.
"""
ember.block_util._wall_law_code(self.wall_law)
if self.n_step_avg < 0:
raise ValueError(f"n_step_avg must be >= 0, got {self.n_step_avg}.")
if self.n_step_avg > self.n_step:
raise ValueError(
f"n_step_avg={self.n_step_avg} exceeds n_step={self.n_step}: the "
"march cannot accumulate more samples than it takes steps, so the "
"averaged solution would come back scaled by n_step/n_step_avg."
)
[docs]
def run(self, grid):
"""Drive ``grid`` through ``n_step`` steps in place; return the history.
The public :class:`BaseSolver` entry point for the stage-by-stage march.
"""
return _run(grid, self)
[docs]
def run_fmg(self, grid):
"""Full-multigrid startup on ``grid`` in place.
Returns a list of per-level :class:`ConvergenceHistory`, coarsest first.
Not part of the :class:`BaseSolver` contract (plugins have no FMG
analogue).
"""
return _run_fmg(grid, self)
[docs]
def soft(self):
"""Return a copy detuned for a robust start; see :meth:`BaseSolver.soft`.
Returns
-------
Solver
A new robust configuration.
"""
return replace(
self,
cfl=1.5,
sf2=0.02,
sf4=0.01,
dampin=10.0,
fac_mgrid=0.0,
n_step_avg=1,
adaptive_smoothing=False,
# A cold start's mixing planes need the faster exchange: at 0.02
# a turbine stage's soft pass diverged at the hub by step 34 that
# survives at 0.05, while the production march is steadier at the
# lower value (turbigen sweep9 0094 and 0017).
rf_exchange=0.05,
)
[docs]
def 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."""
# Preconditions: dt_vol_nd populated and cached P/T consistent with
# conserved_nd on entry. The caller invalidates caches and applies boundary
# conditions between steps, and smooths once on the post-step state
# (Grid.smooth, shared with the RK path); no smoothing happens here.
#
# fac_mgrid == 0 scales every coarse correction to identically zero, so
# collapse to the plain no-MG dispatch (which also makes sf_irs inert)
# rather than paying restrict/prolong per level for a guaranteed-zero push.
n_levels_eff = n_levels if fac_mgrid > 0.0 else 0
for block in grid:
ni, nj, nk = block.shape
cell_shape = (ni - 1, nj - 1, nk - 1, 5)
# block.store is sized to the nodal shape (shared with the RK path) but
# the Denton residual history is cell-shaped, so carve a leading
# cell-shaped view of it. It is a separate, persistent buffer from the
# arena and so is carved on its own.
store_cell = util.carve_view(block.store, cell_shape)
# residual_nd is read here (a cache hit from the loop's get_convergence, or
# a fresh evaluation) BEFORE the kernel runs. Evaluating residual_nd
# borrows block.scratch for both the IRS work vector and the rolling
# flow buffers, and the kernels below reuse the same arena, so it must
# be fully materialised first -- passing it as an argument guarantees
# that ordering.
if n_levels_eff > 0:
# Denton block-sum multigrid over the scheme-agnostic engine
# (mg_restrict_levels), like advance_rk_stage_mg's coarse path at
# alpha=1 but for the scree fine term: every level uses the
# Denton-lagged q = 2*residual - store, and the correction is
# injected into the fine cell increment, fused with its cell->node
# scatter, so a rolling two-plane buffer serves. sf_irs > 0 picks
# the coarse-IRS kernel; the fine term arrives already smoothed by
# the caller. ONE carve: rbuf and coarse scratch reach one call.
rbuf, *mg_bufs = util.carve_view(
block.scratch,
(ni - 1, nj - 1, 5, 2),
*mg_coarse_shapes(ni, nj, nk, n_levels_eff),
)
kernel = (
ember.fortran.scree_mg_irs
if sf_irs > 0.0
else ember.fortran.scree_mg_noirs
)
kernel(
cons=block.conserved_nd,
residual=block.residual_nd,
store=store_cell,
dt_vol=block.dt_vol_nd,
vol=block.vol_nd,
cfl=cfl,
fmgrid=fac_mgrid,
expon_mgrid=expon_mgrid,
sf_irs=sf_irs,
n_levels=n_levels_eff,
rbuf=rbuf,
rfac=block.damp_rfac,
dampin=dampin,
**dict(zip(MG_COARSE_NAMES, mg_bufs)),
)
else:
# Multigrid off: fine term only, no coarse scratch. Forms
# q = 2*residual - store in store, builds the increment cfl*dt_vol*q,
# rolls the history (store <- residual) and frozen-scatters the
# increment straight onto cons -- bypassing the setters so the P/T
# cache stays frozen. Untouched by fac_mgrid/n_levels. This path
# still materialises the full-volume increment, as rk_plain does.
tmp = util.carve_view(block.scratch, cell_shape)
ember.fortran.scree_plain(
cons=block.conserved_nd,
residual=block.residual_nd,
store=store_cell,
dt_vol=block.dt_vol_nd,
cfl=cfl,
tmp=tmp,
rfac=block.damp_rfac,
dampin=dampin,
)
def _mg_coarse_scratch_sizes(ni, nj, nk, n_levels, np=5):
"""Element counts for the hier2 kernels' flat packed scratch.
``n_corr`` sizes ``corr_all``, which holds every coarse level's scaled
correction back-to-back (coarsest level first, where the collapse starts) --
the sum of per-level element counts. Each slot is compact and exactly the
size the coarse-residual smoother wants, so the restriction gathers straight
into it and scales it in place; there is no separate coarse-residual buffer.
``triw`` is reused per level rather than packed, so it only needs the
largest (level-1) IRS Thomas-coefficient size. Both are carved once per call
from ``block.scratch``, never reallocated -- see
:func:`advance_rk_stage_mg` and :func:`scree_step`.
"""
n_corr = 0
for lvl in range(1, n_levels + 1):
b = 2**lvl
n_corr += ((ni - 1) // b) * ((nj - 1) // b) * ((nk - 1) // b) * np
nc1i, nc1j, nc1k = (ni - 1) // 2, (nj - 1) // 2, (nk - 1) // 2
n_tri = 2 * (nc1i + nc1j + nc1k) if n_levels > 0 else 0
return n_corr, n_tri
# The multigrid kernels' scratch arguments, in the order mg_coarse_shapes
# returns their shapes. Kept together so the two cannot drift apart.
MG_COARSE_NAMES = ("dtblk", "rawbuf", "sdt", "sv", "corr_all", "triw")
[docs]
def 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.
"""
nc1i, nc1j, nc1k = (ni - 1) // 2, (nj - 1) // 2, (nk - 1) // 2
n_corr, n_tri = _mg_coarse_scratch_sizes(ni, nj, nk, n_levels)
return (
(nc1i, nc1j, nc1k),
(nc1i, nc1j, nc1k, 5),
(nc1i, nc1j, nc1k),
(nc1i, nc1j, nc1k),
(n_corr,),
(n_tri,),
)
[docs]
def advance_rk_stage_mg(
grid, alpha, cfl, fac_mgrid, n_levels, expon_mgrid=2.0, sf_irs=0.0, dampin=0.0
):
r"""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.
"""
# fac_mgrid == 0 makes the coarse loop a no-op; collapse to no-MG dispatch.
n_levels_eff = n_levels if fac_mgrid > 0.0 else 0
for block in grid:
ni, nj, nk = block.shape
if n_levels_eff > 0:
# Multigrid-on RK wrappers over the scheme-agnostic engine
# (mg_restrict_levels). sf_irs > 0 selects the coarse-IRS kernel,
# otherwise the plain _noirs kernel, which enters no smoothing code
# at all; the two share the engine and differ only in the smoother
# passed. mg_collapse_levels leaves the total in corr_all's finest
# slot and the fine scatter injects it alongside the fine term, so
# a rolling two-plane buffer serves in place of a full-volume
# increment. One carve, as in scree_step: both reach one call.
rbuf, *mg_bufs = util.carve_view(
block.scratch,
(ni - 1, nj - 1, 5, 2),
*mg_coarse_shapes(ni, nj, nk, max(n_levels_eff, 0)),
)
kernel = (
ember.fortran.rk_mg_irs if sf_irs > 0.0 else ember.fortran.rk_mg_noirs
)
kernel(
cons=block.conserved_nd,
snapshot=block.store,
residual=block.residual_nd,
dt_vol=block.dt_vol_nd,
vol=block.vol_nd,
alpha=alpha,
cfl=cfl,
fmgrid=fac_mgrid,
expon_mgrid=expon_mgrid,
sf_irs=sf_irs,
n_levels=n_levels_eff,
rbuf=rbuf,
rfac=block.damp_rfac,
dampin=dampin,
**dict(zip(MG_COARSE_NAMES, mg_bufs)),
)
else:
# Multigrid off: plain Jameson RK fine-term stage, no coarse scratch.
tmp = util.carve_view(block.scratch, (ni - 1, nj - 1, nk - 1, 5))
ember.fortran.rk_plain(
cons=block.conserved_nd,
snapshot=block.store,
residual=block.residual_nd,
dt_vol=block.dt_vol_nd,
alpha=alpha,
cfl=cfl,
tmp=tmp,
rfac=block.damp_rfac,
dampin=dampin,
)
[docs]
@util.profile
def rk_step(grid, conf):
"""Advance every block one Jameson multi-stage RK step in place."""
# Snapshot the bconds-consistent step-top state U_n into block.store (the
# caller's residual build did not touch conserved_nd); every stage marches
# off this frozen snapshot. Smoothing is not applied here -- Grid.smooth runs
# once on the post-step state, shared with the scree path.
for block in grid:
block.store[...] = block.conserved_nd
for i_stage in range(conf.n_stage):
# Stage coefficient alpha_k = 1/(n_stage - k); the final stage takes
# the full step. advance_rk_stage_mg marches off the snapshot and folds
# in the Denton block-sum multigrid correction (an empty coarse loop at
# conf.n_levels == 0, i.e. a plain RK stage). A nonzero conf.sf_resid
# smooths both the fine residual and the coarse block-restricted one.
# Stage 0 reuses the step-top P/T flush and residual; later stages
# march off a changed conserved_nd, so P/T and bconds are refreshed
# below before the next advance.
alpha = 1.0 / (conf.n_stage - i_stage)
advance_rk_stage_mg(
grid,
alpha,
conf.cfl,
conf.fac_mgrid,
conf.n_levels,
expon_mgrid=conf.expon_mgrid,
sf_irs=conf.sf_resid,
dampin=conf.dampin,
)
grid.update_cached_conserved()
grid.apply_bconds()
# The residual just feeds the next stage's advance. The final stage has
# no next stage, and nothing reads residual_nd between here and the next
# step's top-of-loop rebuild (smooth/accumulate_avg/get_convergence do
# not, and the pre-march update_residual recomputes it first), so skip
# the redundant final rebuild.
if i_stage < conf.n_stage - 1:
grid.update_residual(sf=conf.sf_resid)
def _apply_bcond_relaxation(grid, conf):
"""Push the configured boundary settings onto the grid's patches.
The relaxation factors live on the patch rather than on the solver, so that
they survive a restart and so that two planes of a multi-row grid can be
damped differently; :attr:`~ember.solver.Solver.mix_reflective` is here for
a different reason -- the patch and the communicator read it from places
that never see a Solver. Either way the solver has the last word: a
configured value is imposed on every patch it names, so a run follows from
its own configuration rather than from whatever a setup script left behind.
Pass None to opt out of that for one setting and keep what the patches
carry.
Called per level from :func:`_run`, so the full-multigrid chain configures
each of its separately resampled grids.
"""
for sigma, patches in (
(conf.rf_inlet, grid.patches.inlet),
(conf.rf_outlet, grid.patches.outlet),
(conf.rf_mix, grid.patches.mixing),
):
if sigma is not None:
for patch in patches:
patch.sigma = sigma
# The communicator reads this from the patches at every exchange, so a
# cached communicator picks it up without being rebuilt.
if conf.rf_exchange is not None:
for patch in grid.patches.mixing:
patch.rf_exchange = conf.rf_exchange
# Read by the communicator and by each side's own correction, and so held
# on the patch for the same reason as rf_exchange.
if conf.mix_endwall_span is not None and not 0.0 <= conf.mix_endwall_span < 0.5:
raise ValueError(
f"mix_endwall_span must lie in [0, 0.5), got {conf.mix_endwall_span}: "
"it is measured in from each endwall, so a half would meet in the middle."
)
if conf.mix_endwall_rf is not None and not 0.0 < conf.mix_endwall_rf <= 1.0:
raise ValueError(
f"mix_endwall_rf must lie in (0, 1], got {conf.mix_endwall_rf}: it is "
"the fraction of rf_exchange left at the wall."
)
for name, value in (
("endwall_span", conf.mix_endwall_span),
("endwall_rf", conf.mix_endwall_rf),
):
if value is not None:
for patch in grid.patches.mixing:
setattr(patch, name, value)
# Read by the patch itself (set_block_avg, step, apply, _calc_reference)
# as well as by the communicator, none of which see a Solver -- hence the
# stamp here rather than an argument. Imposed on both sides of every plane
# from one value, which is what check_match's agreement requirement needs.
if conf.mix_reflective is not None:
for patch in grid.patches.mixing:
patch._reflective = conf.mix_reflective
def _validate_mg(grid, n_levels):
"""Raise if any block cannot be evenly divided into the coarsest MG blocks.
Denton block-sum multigrid with ``n_levels`` coarse levels groups cells into
blocks of linear size up to ``2**n_levels``. We require exact division (no
runt edge blocks), so every cell dimension ``ni-1, nj-1, nk-1`` must be a
multiple of the coarsest block size. Since the sizes are powers of two,
divisibility by the coarsest implies it for every finer level, so one check
per dimension suffices. No-op when ``n_levels <= 0``.
Also caps ``n_levels`` at :data:`ember.block.MAX_MG_LEVELS`, because the
shared scratch arena is sized for that depth and a deeper hierarchy would
carve past its end. Checked here rather than at allocation so a bad
configuration is refused before marching, alongside the divisibility rule
it belongs with.
"""
if n_levels <= 0:
return
if n_levels > ember.block.MAX_MG_LEVELS:
raise ValueError(
f"multigrid n_levels={n_levels} exceeds the maximum "
f"{ember.block.MAX_MG_LEVELS}. Block.scratch is one arena sized for "
"that many levels (ember.block._scratch_len), so a deeper hierarchy "
"would carve past the end of it. Raise MAX_MG_LEVELS if you need "
"one, and the arena grows to match."
)
b_coarse = 2**n_levels
for i_block, block in enumerate(grid):
ni, nj, nk = block.shape
for name, n_cell in (("i", ni - 1), ("j", nj - 1), ("k", nk - 1)):
if n_cell % b_coarse != 0:
raise ValueError(
f"Block {i_block} has {n_cell} cells in {name} but "
f"multigrid n_levels={n_levels} needs a multiple of the "
f"coarsest block size {b_coarse}."
)
def _validate_throttle(grid):
"""Raise if more than one outlet patch is throttled to a mass flow.
Each :class:`~ember.patch.OutletPatch` runs its own controller on the mass
flow through its own face, so two throttles on one exit would each see most
of the same error and each apply the whole correction for it, over-throttling
by roughly the number of patches. Splitting a target between them is not the
fix either: the split is only known once the answer is.
An exit spread over several blocks therefore has to prescribe pressure on
all but one of its patches, or be closed by a single patch. Checked at the
start of a run rather than in :meth:`~ember.patch.OutletPatch.set_throttle`,
which sees one patch and cannot know what the rest of the grid carries.
"""
throttled = [p for p in grid.patches.outlet if p.mdot_target is not None]
if len(throttled) > 1:
labels = ", ".join(repr(p.label) for p in throttled)
raise ValueError(
f"{len(throttled)} outlet patches are throttled ({labels}); at most "
"one may be. Prescribe P on the others with set_P, or clear their "
"throttle with set_throttle(None)."
)
@util.profile
def _run(grid, conf):
"""Drive a grid through ``n_step`` explicit time-marching steps.
See the module-level `Overview of one time step`_ for the stage-by-stage
ordering within each step.
``dt_vol`` is relaxed in place by :meth:`~ember.grid.Grid.update_timestep`
(the kernel blends ``rf*new + (1-rf)*old``); the initial call before the
step loop uses ``rf=1.0`` to seed the uninitialised buffer, ``rf=0.2``
thereafter.
Returns
-------
ConvergenceHistory
The recorded history, already trimmed to the steps it logged, so every
row holds data and no ``isfinite`` masking is needed. If the march blew
up, the step loop breaks early,
:attr:`~ember.convergence_history.ConvergenceHistory.diverged` is True, and ``grid`` keeps the
invalid field for inspection (the pseudotime average is not finalised,
since it would overwrite ``conserved_nd`` with a buffer that
``accumulate_avg`` may never have populated).
"""
# Fail fast if the grid cannot be evenly blocked for multigrid.
_validate_mg(grid, conf.n_levels)
_validate_throttle(grid)
_apply_bcond_relaxation(grid, conf)
# Initialise timesteps
grid.update_timestep(rf=1.0, fac_visc=conf.fac_visc)
# One record per log step: i_step % n_step_log == 0 fires ceil(n_step /
# n_step_log) times over range(n_step); floor division would under-allocate.
n_log = -(-conf.n_step // conf.n_step_log)
hist = ConvergenceHistory.from_grid(n_log, grid)
_log_rss(
"march start, n_levels=%d, shape(s) %s",
conf.n_levels,
[block.shape for block in grid],
)
for i_step in range(conf.n_step):
#
# We overwrote conserved_nd in place last step
# So flush the cache to recalculate P and T
grid.update_cached_conserved()
# Throttle/radial equilibrium targets, and the mixing-plane exchange
grid.update_bconds(cfl=conf.cfl)
grid.apply_bconds()
try:
grid.check_nan()
except DivergenceError as err:
logger.error("Solver diverged at step %d: %s", i_step, err)
hist.diverged = True
break
# Refresh source terms every step for the multi-stage RK march, which
# re-evaluates the residual each substep and needs current sources; the
# scree march reuses one residual per step, so lag the expensive viscous
# pass to every fifth step there. Recompute the timestep every step
# regardless so dt_vol tracks the flow through fast transients: a lagged
# dt_vol is sized for an already-stale state and overshoots the
# stability limit during a cold start, and the timestep refresh is cheap.
n_step_source = 5 if conf.n_stage == 0 else 1
add_sources = i_step % n_step_source == 0
if add_sources:
grid.update_sources(conf.inviscid, conf.wall_law)
_log_rss("step %d after update_sources", i_step)
# The timestep pass also finishes the body force on a refresh step (the
# polar source and the SFD force, onto the viscous part just built) and
# advances the SFD low-pass filter on the dt_vol it has just computed.
# The filter advances every step regardless of the source cadence: its
# dt is a per-step increment, so advancing it on the scree march's
# every-fifth refresh would stretch the time constant fivefold. The
# force reads the filter before this step's update, the usual explicit
# SFD coupling. At the default zero gain the filter is never touched.
grid.update_timestep(
rf=0.2,
fac_visc=conf.fac_visc,
add_sources=add_sources,
gain_filt=conf.gain_filt,
cfl=conf.cfl,
delta_filt=conf.delta_filt,
)
# Prepare the residual
grid.update_residual(sf=conf.sf_resid)
_log_rss("step %d after update_residual", i_step)
# Convergence logging of the pre-march state
if i_step % conf.n_step_log == 0:
hist.record_convergence(i_step, grid.get_convergence())
logger.info(
"%s",
hist.format_message(n_step=conf.n_step),
)
_log_rss("step %d after get_convergence", i_step)
# Take a step with the selected integrator. Both reuse the first
# residual evaluated above, RK then recalculates each substep
if conf.n_stage == 0:
scree_step(
grid,
conf.cfl,
fac_mgrid=conf.fac_mgrid,
expon_mgrid=conf.expon_mgrid,
n_levels=conf.n_levels,
sf_irs=conf.sf_resid,
dampin=conf.dampin,
)
else:
rk_step(grid, conf)
_log_rss("step %d after integrator", i_step)
# Smooth the post-step conserved solution
grid.smooth(
conf.sf4 * conf.cfl,
conf.sf2 * conf.cfl,
adaptive=conf.adaptive_smoothing,
sf2_T=None if conf.sf2_T is None else conf.sf2_T * conf.cfl,
)
_log_rss("step %d after smooth", i_step)
# Pseudotime average over the last n_step_avg steps. A window of 0 or
# 1 needs no accumulator: the mean of one sample is that sample, which
# conserved_nd already holds. Skipping it keeps Block.conserved_avg_nd
# -- a full nodal five-component buffer, 19.3 MB at a 1M-cell block --
# from ever being allocated.
if conf.n_step_avg > 1 and i_step >= (conf.n_step - conf.n_step_avg):
grid.accumulate_avg(conf.n_step_avg)
# Copy the final average back into the primary storage. Skipped on
# divergence, where the loop broke before accumulate_avg ran and this would
# overwrite the invalid conserved_nd with a zeroed buffer, and when no
# averaging was asked for, where it would write zeros over a good solution.
# The cache invalidation is NOT optional on those paths: the integrators
# and smooth write conserved_nd through the frozen-pressure path, which
# does not bump the conserved versions, so the last cached P/T can be
# marked current while the field underneath has moved on.
if not hist.diverged:
if conf.n_step_avg > 1:
grid.finalise_average()
else:
grid.update_cached_conserved()
_log_rss("march end, n_levels=%d", conf.n_levels)
# A completed march logs on every one of the ceil(n_step / n_step_log) rows
# that from_grid allocated, so this only bites when the loop broke early:
# the caller never sees the unwritten NaN tail a divergence leaves behind.
return hist.trim()
def _run_fmg(grid, conf):
"""Full-multigrid startup: solve coarse-to-fine, prolonging each guess.
``conf.n_levels`` is the single grid-hierarchy depth. With ``n_levels == 0``
this is exactly :func:`run(grid, conf) <run>`. Otherwise it builds
``n_levels`` grids successively halved by :meth:`~ember.grid.Grid.resample`
at factor ``0.5``, solves the coarsest, and
:meth:`~ember.grid.Grid.interp_from_grid`'s the solution up onto each finer grid
as its initial guess. Sequencing level ``i`` (``0`` = coarsest) is marched
with in-step Denton block-sum multigrid depth ``i`` -- the same grid
hierarchy ``n_levels`` already names -- so the coarsest runs plain and the
finest runs at full ``n_levels``, identical to :meth:`Solver.run` on the finest.
The single validation ``_validate_mg(grid, n_levels)`` is sufficient for the
whole chain: level ``i`` holds ``N / 2**(n_levels - i)`` cells per dimension
and its depth-``i`` march needs those divisible by ``2**i``, which reduces to
``N`` (the finest cell count) divisible by ``2**n_levels`` at every level.
That divisibility also keeps every ``resample(0.5)`` node-coincident with the
finer grid (an exact subset of its nodes), so the coarse coordinates are the
true geometry, not interpolated.
Every level runs the same ``conf`` apart from ``n_levels`` (fixed
``n_step`` per level). ``grid`` is the finest level and is mutated in place
to carry the final solution, matching :meth:`Solver.run`.
Parameters
----------
grid : Grid
Finest grid, already carrying its initial guess. Mutated in place.
conf : Solver
Solver configuration; ``conf.n_levels`` sets the hierarchy depth.
Returns
-------
list of ConvergenceHistory
Per-level histories, coarsest first, finest last.
"""
_validate_mg(grid, conf.n_levels)
if conf.n_levels <= 0:
return [_run(grid, conf)]
# Build finest -> coarsest, then reverse. resample carries the already-set
# fine guess down, so the coarsest starts from the coarsened cold start.
_log_rss("fmg: before building the coarse chain")
chain = [grid]
for _ in range(conf.n_levels):
chain.append(chain[-1].resample(0.5))
_log_rss(
"fmg: resampled level with shape(s) %s",
[block.shape for block in chain[-1]],
)
chain.reverse() # chain[-1] is the original `grid`
histories = []
for i, level_grid in enumerate(chain): # i == in-step MG depth for this mesh
if i > 0:
level_grid.interp_from_grid(chain[i - 1]) # prolong previous solution
_log_rss("fmg: interpolated onto level %d", i)
logger.info(
"FMG level %d/%d, shape(s) %s",
i,
conf.n_levels,
[block.shape for block in level_grid],
)
histories.append(_run(level_grid, replace(conf, n_levels=i)))
_log_rss("fmg: finished level %d", i)
return histories