r"""Functions for integration and averaging over 2D blocks.
This module implements averaging methods for reducing spatially-varying 2D flow fields
to representative scalar quantities, essential for comparing CFD results with mean-line design
points and experimental measurements. The module provides area averaging, mass-flux averaging,
and mixed-out averaging methods following turbomachinery conventions. As discussed by
:cite:t:`Cumpsty2005`, dimensional reduction from 2D to 0D inherently loses information, requiring
careful selection of which flow properties to conserve. :cite:t:`Burdett2022` demonstrate
that mixed-out averages minimize sensitivity to streamwise location variations. All averaging
functions operate on 2D Block objects and support both absolute and relative reference frames
for turbomachinery with rotating components, enabling consistent performance metric extraction
across blade rows.
"""
import numpy as np
from ember import util
from ember import block_util
from ember import perturbation
from ember import fluxes as ember_fluxes
from ember import set_iterative
def _get_axes(axes, triangulated):
"""Helper to validate axes argument."""
if triangulated:
if axes is not None:
raise ValueError("For triangulated grids, axes must be None.")
axes = (0,)
else:
axes = (0, 1) if axes is None else axes
if len(axes) > 2 or any(ax not in (0, 1) for ax in axes) or len(axes) == 0:
raise ValueError(
"For structured grids, axes must be a tuple of two integers from (0, 1), got {axes}."
)
return axes
def _get_dA(block):
"""Return face area array, dispatching to dA_tri or dA_quad based on block.triangulated."""
return block.dA_tri if block.triangulated else block.dA_quad
def _get_dA_nd(block):
"""Nondimensional counterpart of :func:`_get_dA`."""
return block.dA_tri_nd if block.triangulated else block.dA_quad_nd
def _dot_conserved(flux, dA, axes):
"""Specialized dot product for conserved variable fluxes."""
return np.sum(
np.einsum(
"...ij,...i->...j",
flux,
dA,
),
axis=axes,
)
def _integrate_scalar(scalar_face, dA_face, axes):
"""Helper function to integrate scalar over faces.
Agnostic to quadrilateral vs triangulated faces."""
assert scalar_face.ndim == dA_face.ndim - 1
assert dA_face.shape[-1] == 3
assert scalar_face.shape == dA_face.shape[:-1]
return np.sum(scalar_face * util.vecnorm(dA_face), axis=axes)
def _integrate_vector(vector_face, dA_face, axes):
"""Helper function to integrate vector flux over faces.
Agnostic to quadrilateral vs triangulated faces."""
return np.sum(util.dot(vector_face, dA_face), axis=axes)
def _node_to_face_2d(nodal_data):
r"""Average nodal values to face centres for 2D data.
For a :math:`(n_i, n_j)` array of nodal values, the face-centred value is
.. math::
\bar{q}_{i,j} = \tfrac{1}{4}\bigl(
q_{i,j} + q_{i+1,j} + q_{i+1,j+1} + q_{i,j+1}\bigr)
Parameters
----------
nodal_data : Array, shape (ni, nj, ...)
Values at grid nodes.
Returns
-------
Array, shape (ni-1, nj-1, ...)
Values averaged to face centres.
"""
return 0.25 * (
nodal_data[:-1, :-1, ...]
+ nodal_data[1:, :-1, ...]
+ nodal_data[1:, 1:, ...]
+ nodal_data[:-1, 1:, ...]
)
def _node_to_face(nodal_data, triangulated):
"""Distribute nodal data to faces."""
if triangulated:
# Average over triangle vertices
return np.mean(nodal_data, axis=1)
else:
# Use 2D node_to_face distribution
return _node_to_face_2d(nodal_data)
[docs]
def flow_mass(block, axes=None):
r"""Integrate mass flow through faces of a 2D block.
Calculates the mass flow over the block faces,
.. math::
\dot{m} = \int \rho \mathbf{V}\cdot\mathrm{d}\mathbf{A} \,.
Parameters
----------
block : Block, shape (ni, nj) or (ntri, 3)
2D structured or triangulated block.
axes : tuple of int, default (0, 1)
For structured grids, axes over which to sum the mass flow;
for triangulated grids, should be None to sum over all faces.
Returns
-------
mass_flow: float
Mass flow rate through the block.
"""
mass_flux_face = ember_fluxes.get_mass_flux(block)
axes = _get_axes(axes, block.triangulated)
return _integrate_vector(mass_flux_face, _get_dA(block), axes)
[docs]
def flow_conserved(block, axes=None):
r"""Integrate conserved flows through faces of a 2D block.
Calculates the conserved flows over the block faces,
.. math::
\int \mathcal{F}\cdot\mathrm{d}\mathbf{A} \,,
where the conserved flux tensor :math:`\mathcal{F}` carries the fluxes of
mass, axial momentum, radial momentum, angular momentum and stagnation
enthalpy (energy),
.. math::
\mathcal{F} = \rho \mathbf{V}
\begin{bmatrix} 1 \\ V_x \\ V_r \\ r V_\theta \\ h_0 \end{bmatrix}
+ p
\begin{bmatrix}
\mathbf{0} \\ \mathbf{e}_x \\ \mathbf{e}_r
\\ r\,\mathbf{e}_\theta \\ \Omega r\,\mathbf{e}_\theta
\end{bmatrix} \,.
Parameters
----------
block : Block, shape (ni, nj) or (ntri, 3)
2D structured or triangulated block.
axes : tuple of int, default (0, 1)
For structured grids, axes over which to sum the flows;
for triangulated grids, should be None to sum over all faces.
Returns
-------
flow_conserved : Array shape (5,)
Integrated conserved flows
"""
return (
flow_conserved_nd(block, axes) * ember_fluxes.flux_ref(block) * block.L_ref**2
)
[docs]
def flow_conserved_nd(block, axes=None):
r"""Nondimensional counterpart of :func:`flow_conserved`, shape (5,).
Scaled by ``fluxes.flux_ref(block) * block.L_ref**2``, this is
:func:`flow_conserved`. Working directly in this form lets :func:`mix_out`
run its Newton iteration without any reference-scale bookkeeping, since
the Jacobian from :mod:`ember.perturbation` is nondimensional too.
"""
axes = _get_axes(axes, block.triangulated)
return _dot_conserved(ember_fluxes.get_flux_nd(block), _get_dA_nd(block), axes)
[docs]
def mass_average(scalar_node, block, axes=None):
r"""Take mass-weighted average of a 2D nodal scalar field.
Calculates the mass-weighted average of the scalar field :math:`\phi`,
.. math::
\bar{\phi} = \frac{\int \phi\, \rho \mathbf{V}\cdot\mathrm{d}\mathbf{A}}
{\int \rho \mathbf{V}\cdot\mathrm{d}\mathbf{A}} \,.
Parameters
----------
scalar_node : Array, shape (ni, nj) or (ntri, 3)
Scalar field values at grid nodes
block : Block, shape (ni, nj) or (ntri, 3)
2D structured or triangulated block
axes : tuple of int, default (0, 1)
For structured grids, axes over which to average;
for triangulated grids, should be None to average over all faces
Returns
-------
avg_scalar : float
Mass-weighted average value
Raises
------
ValueError
If the net mass flux through the block is zero
"""
axes = _get_axes(axes, block.triangulated)
mass_flux_vector = ember_fluxes.get_mass_flux(block) # Shape (..., 3)
scalar_face = _node_to_face(
scalar_node, block.triangulated
) # Shape (...,) - scalar at faces
# Calculate scalar mass flow through each face: phi * (rho V . dA)
mass_flow_scalar = util.dot(
mass_flux_vector, _get_dA(block)
) # Shape (...,) - scalar mass flow
scalar_mass_flow = (
scalar_face * mass_flow_scalar
) # Shape (...,) - scalar-weighted mass flow
numerator = np.sum(scalar_mass_flow, axis=axes)
denominator = flow_mass(block, axes)
# Check for zero net mass flux. `np.any` because a partial reduction --
# `axes=(1,)` to average over the pitch and keep the span, say -- leaves an
# array here rather than a scalar, and a bare `<` on one of those raises
# about an ambiguous truth value instead of reporting the flux.
if np.any(np.abs(denominator) < 1e-14):
raise ValueError(
"Net mass flux through the block is zero. "
"Mass averaging requires non-zero net mass flux. "
"Consider using area_average() instead."
)
return numerator / denominator
[docs]
def mass_band(block, lo, hi, axis=0):
r"""Extract the part of a structured 2D cut carrying a band of mass flow.
Integrates the mass flow strip by strip along `axis`, summing across the
other axis, to give a cumulative mass fraction :math:`\psi` running from
zero at the first node to one at the last. The band is the region
:math:`\mathrm{lo} \le \psi \le \mathrm{hi}`, so ``lo=0.45, hi=0.55`` on a
hub-to-casing axis is the stream surface carrying the central tenth of the
flow.
The two ends of the band sit at fractional node positions inside the
strips where :math:`\psi` first reaches `lo` and `hi`. The data are
interpolated linearly to those two positions and the original nodes
between them are copied through, so the band keeps the resolution of the
cut it came from. Each end is placed by bisection on the mass flow of the
partial strip it bounds, integrated just as :func:`flow_mass` integrates
the band, so the band carries :math:`\mathrm{hi} - \mathrm{lo}` of the flow
to round-off --- not merely to within the variation of flux and face area
across a strip, which a linear interpolation of :math:`\psi` would leave.
The first crossing is taken rather than an inverse interpolation of
:math:`\psi`, which need not be monotonic: a strip of backflow in a corner
elsewhere on the cut does not stop the band being found.
Parameters
----------
block : Block, shape (ni, nj)
Structured 2D block.
lo, hi : float
Mass fractions bounding the band, :math:`0 \le \mathrm{lo} <
\mathrm{hi} \le 1`.
axis : int, default 0
Axis along which the mass fraction runs. The band spans the whole of
the other axis.
Returns
-------
band : Block, shape (nband, nj) or (ni, nband)
New block holding the band, with the metadata of `block`.
Raises
------
ValueError
If `block` is triangulated or not 2D, if `axis` is not 0 or 1, if the
fractions are out of order or out of range, or if the net mass flux
through the block is zero.
"""
if block.triangulated or block.ndim != 2:
raise ValueError("mass_band requires a structured 2D block.")
if axis not in (0, 1):
raise ValueError(f"axis must be 0 or 1, got {axis}.")
if not 0.0 <= lo < hi <= 1.0:
raise ValueError(f"Need 0 <= lo < hi <= 1, got lo={lo}, hi={hi}.")
# Mass through each strip between adjacent nodes along `axis`
dm = np.asarray(flow_mass(block, axes=(1 - axis,)), dtype=float)
total = dm.sum()
if np.abs(total) < 1e-14:
raise ValueError(
"Net mass flux through the block is zero, so it has no mass fraction."
)
# Dividing by a signed total makes the fraction run zero to one whichever
# way the flow crosses the cut. The last node is pinned to one, so that
# round-off cannot leave hi=1 short of it.
psi = np.concatenate([[0.0], np.cumsum(dm)]) / total
psi[-1] = 1.0
# Rows along the band axis, so both ends are found and built the same way
# whichever axis the fraction runs along.
data = np.moveaxis(block._data, axis, 0)
def fraction(p0, p1):
"""Mass fraction between fractional positions p0 < p1 in one strip."""
return flow_mass(_interpolate_rows(block, data, [p0, p1], axis)) / total
# The low end: the partial strip from it up to the next node must carry
# what that node's fraction exceeds `lo` by.
k = int(np.argmax(psi >= lo))
if psi[k] == lo:
pos_lo = float(k)
else:
need = psi[k] - lo
pos_lo = _bisect(lambda p: fraction(p, k) < need, k - 1, k)
# The high end: the partial strip up to it must carry what `hi` exceeds the
# start of that strip by. The strip starts at the low end if both ends fall
# in the same one.
k = int(np.argmax(psi >= hi))
if psi[k] == hi:
pos_hi = float(k)
else:
start = max(k - 1, pos_lo)
need = hi - max(psi[k - 1], lo)
pos_hi = _bisect(lambda p: fraction(start, p) >= need, start, k)
inner = np.arange(np.floor(pos_lo) + 1, np.ceil(pos_hi))
pos = np.concatenate([[pos_lo], inner, [pos_hi]])
return _interpolate_rows(block, data, pos, axis)
def _bisect(past, a, b, n=60):
"""Return the point in [a, b] where the predicate `past` turns true.
`past` is false at `a` and true at `b`. Sixty halvings of a unit interval
reach the spacing of doubles, so the result is as exact as the predicate.
"""
for _ in range(n):
p = 0.5 * (a + b)
if past(p):
b = p
else:
a = p
return 0.5 * (a + b)
def _interpolate_rows(block, data, pos, axis):
"""Return a block of `data` rows linearly interpolated to positions `pos`.
`data` has its interpolation axis first; the result has it back at `axis`.
A position on the last node is taken from the last strip rather than past
it.
"""
pos = np.asarray(pos, dtype=float)
j0 = np.minimum(np.floor(pos).astype(int), data.shape[0] - 2)
w = (pos - j0)[:, None, None]
rows = np.moveaxis((1.0 - w) * data[j0] + w * data[j0 + 1], 0, axis)
out = block.empty(shape=rows.shape[:-1])
out._data = rows
return out
[docs]
def area_average(scalar_node, block, axes=None):
r"""Take area-weighted average of a 2D nodal scalar field.
Calculates the area-weighted average of the scalar field :math:`\phi`,
.. math::
\bar{\phi} = \frac{\int \phi\, \mathrm{d}A}{\int \mathrm{d}A} \,.
Parameters
----------
scalar_node : Array, shape (ni, nj) or (ntri, 3)
Scalar field values at grid nodes
block : Block, shape (ni, nj) or (ntri, 3)
2D structured or triangulated block
axes : tuple of int, default (0, 1)
For structured grids, axes over which to average;
for triangulated grids, should be None to average over all faces
Returns
-------
avg_scalar : float
Area-weighted average value
"""
axes = _get_axes(axes, block.triangulated)
scalar_face = _node_to_face(scalar_node, block.triangulated)
numerator = _integrate_scalar(scalar_face, _get_dA(block), axes)
denominator = _integrate_scalar(np.ones_like(scalar_face), _get_dA(block), axes)
return numerator / denominator
[docs]
def total_area(block):
r"""Compute total vector area of a 2D block.
Calculates the total vector area as the integral of the face area vectors,
.. math::
\mathbf{A} = \int \mathrm{d}\mathbf{A} \,.
Parameters
----------
block : Block, shape (ni, nj) or (ntri, 3)
2D structured or triangulated block
Returns
-------
A : Array, shape (3,)
Total area of the cut [m^2] in polar coordinates (Ax, Ar, At)
"""
return total_area_nd(block) * block.L_ref**2
[docs]
def total_area_nd(block):
r"""Nondimensional counterpart of :func:`total_area`, shape (3,).
Total vector area divided by :math:`L_\mathrm{ref}^2`.
"""
assert block.ndim == 2
axes = _get_axes(None, block.triangulated)
return np.sum(_get_dA_nd(block), axis=axes)
[docs]
def mix_out(block, AR=1.0):
r"""Mix out a 2D cut to uniformity, optionally through a contracted area.
The mixed-out state is the uniform flow that, passed through the total
area :math:`\mathbf{A} = \int \mathrm{d}\mathbf{A}`, carries the same
conserved flows as the non-uniform cut. Its conserved variables
:math:`\mathcal{U}` are found by solving
.. math::
\mathcal{F}(\mathcal{U})\cdot\mathbf{A}
= \int \mathcal{F}\cdot\mathrm{d}\mathbf{A} \,,
for the five conserved flows (mass, axial and radial momentum, angular
momentum and energy), where :math:`\mathcal{F}` is the flux tensor of
:func:`flow_conserved`. The five equations are solved iteratively by
Newton steps on :math:`\mathcal{U}`. Mixing to uniformity generates
entropy, so the result has higher entropy than the original state.
The optional area ratio ``AR`` then contracts (``AR<1``) or expands the
uniform state isentropically from :math:`\mathbf{A}` to ``AR`` times
:math:`\mathbf{A}`, conserving mass, stagnation enthalpy, entropy and
angular momentum :math:`r V_\theta` (at fixed radius) while holding the
pitch angle :math:`\beta`. This second step is reversible, so the mixing
loss is independent of ``AR`` and ``AR=1`` recovers the plain mix-out. The
contraction stays on the mixed-out sub/supersonic branch and raises
:class:`RuntimeError` if it would choke.
Parameters
----------
block : Block, shape (ni, nj) or (ntri, 3)
2D block, can be structured or triangulated.
AR : float, default 1.0
Area ratio for the isentropic contraction applied after mixing out.
``AR=1`` retains the constant-area mix-out; ``AR<1`` contracts the
uniform state to area ``AR * A``.
Returns
-------
mix : Block, shape ()
New scalar block with mixed-out uniform state.
"""
# Calculate total area and conserved quantities. The whole Newton loop
# below runs nondimensionally, matching the ND Jacobian from
# ember.perturbation, so no reference scales appear in it at all.
A = total_area_nd(block)
flow = flow_conserved_nd(block)
# Ensure that mass flow is positive
if flow[0] <= 0.0:
A *= -1.0
flow *= -1.0
# Do not allow significant projected area in theta direction
assert np.all(np.abs(A[2]) < 1e-6 * np.linalg.norm(A[:2])), (
f"Block has significant projected area in theta direction: A={A}"
)
A = A[:2] # Drop theta component
A_ref = np.linalg.norm(A)
# Allocate a scalar block for the mixed-out state
mix = block.empty()
# Mixed out coordinates
rmix = np.sqrt(np.mean(np.array([block.r.min(), block.r.max()]) ** 2))
xmix = util.extent(block.x).mean()
tmix = util.extent(block.t).mean()
mix.set_x(xmix)
mix.set_r(rmix)
mix.set_t(tmix)
# Initial guess for conserved variables simple mean, written direct on the
# ND storage. This skips set_conserved's density/finite validation, which
# a mean of valid input data cannot trip, and the loop's own mix.rho check
# still guards against Newton-step divergence below.
mix.conserved_nd[...] = block.conserved_nd.mean(axis=(0, 1))
mix.update_cached_conserved()
# Get absolute tolerance for flows, on the same ND scale as the flows
rho_ref = mix.rho_nd
V_ref = mix.V_nd
rhoV_ref = rho_ref * V_ref
rhoVsq_ref = rho_ref * V_ref**2
de_ref = rho_ref * V_ref**3
atol = (
np.array(
[
rhoV_ref,
rhoVsq_ref,
rhoVsq_ref,
rhoVsq_ref * block.r_mid_nd,
de_ref,
]
)
* A_ref
* 1e-4
)
# Iteratively adjust conserved variables to match total flow
rf = 0.1
max_iter = 10000
n_stall = 50
err_flow = np.inf
best_err = np.inf
best_iter = 0
for niter in range(max_iter):
# Calculate current fluxes and flows (xr system)
flux_mix = ember_fluxes.get_flux_nd(mix)[:2, :]
flow_mix = _dot_conserved(flux_mix, A, axes=())
err_flow = flow - flow_mix
err_flux = err_flow / A_ref
# Stop once the residual stops improving, not the moment it first
# drops below atol: breaking at the crossing makes the answer depend on
# which iterate lands inside the tolerance ball first, and so on input
# ordering. Iterating to the float32 fixed point removes that. atol
# survives below as the scale that makes the five residuals comparable
# and as the post-loop divergence check.
err_scaled = np.max(np.abs(err_flow) / atol)
if err_scaled < best_err * (1.0 - 1e-6):
best_err = err_scaled
best_iter = niter
elif niter - best_iter > n_stall:
break
# Resolve to interface-aligned velocities. Beta moves every iteration
# as mix's state does, so the matrix can't be hoisted out of the loop
# -- but building both directions from it here, once, means the
# to/from pair below shares one sine and cosine rather than each
# re-deriving its own.
Beta = mix.Beta
rot_to, rot_from = util.rotation_matrices(np.radians(Beta))
block_util.resolve_to_interface(mix, rot_to)
# Calculate Jacobian of conserved/flux transformation (nondimensional)
f2c = perturbation.flux_to_conserved(mix)
# Apply the ND Jacobian to the (already ND) flux error to get the ND
# conserved correction directly
dcons_nd = util.matvec(f2c, err_flux)
# Apply relaxation to avoid overshoot, direct on the ND storage
mix.conserved_nd[...] += rf * dcons_nd
mix.update_cached_conserved()
if mix.rho < 0.0:
print(
f" NEGATIVE DENSITY at iter {niter}: rho={mix.rho:.6g}, dcons_nd={dcons_nd}"
)
raise Exception("Negative density")
# Resolve back to physical velocities
block_util.resolve_from_interface(mix, rot_from)
if (np.abs(err_flow) >= atol).any():
print(f" FAILED after {niter + 1} iters: err_flow={err_flow}, atol={atol}")
print(f" final conserved: {mix.conserved}")
raise RuntimeError(
f"Failed to converge mixing after {niter + 1} iterations, err_flow={err_flow}, atol={atol}"
)
# Optionally contract the uniform mixed-out state isentropically from area
# A to AR*A. This is a reversible second step that leaves the mixing loss
# (computed above at the true area) unchanged.
if AR != 1.0:
# Only subsonic meridional flow is supported by the contraction solver.
if mix.Vm >= mix.a:
raise NotImplementedError(
"Isentropic contraction (AR != 1) is only supported for "
f"subsonic meridional flow (got Vm={mix.Vm:.6g} >= a={mix.a:.6g})."
)
# Meridional mass flux through the contracted area AR*A. Mass is
# conserved, so rho*Vm scales as 1/AR.
rhoVm_target = mix.rho * mix.Vm / AR
set_iterative.set_ho_s_rhoVm_Vt_Beta(
mix,
mix.ho,
mix.s,
rhoVm_target,
mix.Vt,
Beta=mix.Beta,
)
return mix