Composite elements: one element, many atoms

A composite element presents as one element but expands at build time into a small graph of atomic elements joined by internal edges.

Source notebook: examples/flow/composite_elements.ipynb. This page is executed from it at build time.

A composite element presents as a single element but expands, at build time, into a small graph of atomic elements joined by internal edges. The solver, Jacobian and perturbation layers never see a composite: at solve time the expanded graph is indistinguishable from a hand-built one (no new kernels). The user’s declared node/edge ids are preserved (internals append at the tail), so existing edge-indexed analysis keeps working unchanged.

This behaviour actually causes the coefficient matrix to be less optimal in terms of its band structure, flagged to be addressed in a future update. The present sparse solver from scipy already does an internal re-ordering, hence the actual band structure does not really matter there.

Two families of composite elements are shown here:

import warnings
import numpy as np
import plotly.graph_objects as go

import nefes
from nefes.elements import catalog as cat
from nefes.elements import auto_refine
from nefes.plotting import use_nefes_theme, COLORWAY

use_nefes_theme()

GAMMA, R = 1.4, 287.0
GAS = nefes.perfect_gas(R, GAMMA)
P0, T0 = 101325.0, 300.0

1: the orifice round-trips a hand-placed iac + sac

The orifice is an isentropic contraction to the throat followed by a sudden expansion (Borda-Carnot). cat.orifice(throat_area) is exactly that assembly behind one element. We solve it both ways and confirm the mean flow on the shared inlet/outlet edges and the compact scattering matrix are identical.

The throat area represents the effective area of the orifice, e.g. A_{\text{eff}} = C_d A_{\text{geo}}, where C_d is the associated discarge coefficient. This element can, in principle, be used to represent arbitrary orifice-like flow devices without being subjected to the low-Mach restrictions of the generalized loss element.

A1, AT, A2 = 3.0e-3, 1.0e-3, 3.0e-3
PT_IN = 140000.0

#  (a) hand-built: inlet - iac - sac - outlet
hand = nefes.Network(
    GAS,
    nodes=[
        cat.total_pressure_inlet(PT_IN, T0),
        cat.isentropic_area_change(),
        cat.sudden_area_change(),
        cat.pressure_outlet(P0, T0),
    ],
    edges=[(0, 1, A1), (1, 2, AT), (2, 3, A2)],
)
sh = hand.solve()
assert sh.converged, (sh.residual_norm, sh.print_residuals())

#  (b) composite: inlet - orifice - outlet
orifice_net = nefes.Network(
    GAS,
    nodes=[
        cat.total_pressure_inlet(PT_IN, T0),
        cat.orifice(AT, name="orifice"),
        cat.pressure_outlet(P0, T0),
    ],
    edges=[(0, 1, A1), (1, 2, A2)],
)
sc = orifice_net.solve()
assert sc.converged, (sc.residual_norm, sc.print_residuals())

probe = np.array([120.0])
S_hand = sh.perturbation_response(probe).acoustic_scattering_matrix(0, 2)[0]
S_comp = sc.perturbation_response(probe).acoustic_scattering_matrix(0, 1)[0]
throat = sc.composite("orifice").throat  # throat edge id, without knowing the expanded layout
print(f"inlet Mach   hand {sh.edge(0)['M']:.5f}   composite {sc.edge(0)['M']:.5f}")
print(f"throat Mach  hand {sh.edge(1)['M']:.5f}   composite {sc.edge(throat)['M']:.5f}")
print(f"outlet Mach  hand {sh.edge(2)['M']:.5f}   composite {sc.edge(1)['M']:.5f}")
print(f"scattering matrix max |difference|: {np.max(np.abs(S_hand - S_comp)):.2e}")
inlet Mach   hand 0.19703   composite 0.19703
throat Mach  hand 0.95158   composite 0.95158
outlet Mach  hand 0.26416   composite 0.26416
scattering matrix max |difference|: 4.44e-15

Result. Identical to solver tolerance: the composite is a pure build-time convenience. solution.composite('orifice').throat reads the throat state directly, without knowing the expanded edge layout.

2: the Helmholtz resonator as one element

cat.helmholtz_resonator(volume, neck_length, neck_area) wraps the tee + neck duct + cavity side-branch behind one element (the TODO’s first user of composites). It resonates at the lumped Helmholtz frequency, shorting the main line into a transmission-loss peak.

V, AN, LN, AM, LM = 1.0e-3, 5.0e-4, 0.02, 3.0e-3, 0.05
freqs = np.linspace(50.0, 1100.0, 1100)

net = nefes.Network(
    GAS,
    nodes=[
        cat.total_pressure_inlet(P0, T0, name="inlet"),
        cat.duct(LM, name="upstream"),
        cat.helmholtz_resonator(V, LN, AN, name="resonator"),
        cat.duct(LM, name="downstream"),
        cat.pressure_outlet(P0, T0, name="outlet"),
    ],
    edges=[(0, 1, AM), (1, 2, AM), (2, 3, AM), (3, 4, AM)],
)
sol = net.solve()
assert sol.converged, (sol.residual_norm, sol.print_residuals())

e_in = net.edge_between("inlet", "upstream")
e_out = net.edge_between("downstream", "outlet")

resp = sol.perturbation_response(freqs)
with warnings.catch_warnings():  # inlet/outlet straddle the internal tee (expected)
    warnings.simplefilter("ignore")
    tau = resp.acoustic_scattering_matrix(e_in, e_out)[:, 1, 0]
tl = -20.0 * np.log10(np.abs(tau))
c = sol.field("c")[e_in]
f0 = c * np.sqrt(AN / (V * LN)) / (2.0 * np.pi)

fig = go.Figure()
fig.add_trace(go.Scatter(x=freqs, y=tl, name="transmission loss", line=dict(color=COLORWAY[0], width=2)))
fig.add_vline(x=f0, line=dict(color=COLORWAY[4], width=1, dash="dot"), annotation_text="f0", annotation_position="top")
fig.update_layout(
    title="helmholtz_resonator: one element, a Helmholtz peak",
    xaxis_title="Frequency, Hz",
    yaxis_title="Transmission Loss, dB",
    width=820,
    height=440,
)
fig.show()
print(f"analytic f0 = {f0:.1f} Hz   Nefes peak = {freqs[int(np.argmax(tl))]:.1f} Hz")
analytic f0 = 276.3 Hz   Nefes peak = 275.5 Hz

Result. The composite lands the resonance on the analytic Helmholtz frequency, exactly as the hand-built tee + neck + cavity would.

3: fanno_pipe: a discretized pipe resolves the Mach rise

A long pipe is Fanno flow: wall friction accelerates the subsonic flow toward the exit, so a single lumped friction coefficient misses the gradient. cat.fanno_pipe(length, diameter, friction, N) chains N pipe atoms; as N grows the chain converges to the true Fanno solution.

L, D, F, MDOT = 8.0, 0.05, 0.03, 0.55
area = np.pi * D**2 / 4.0


def fanno_net(N):
    return nefes.Network(
        GAS,
        nodes=[cat.mass_flow_inlet(MDOT, T0), cat.fanno_pipe(L, D, F, N, name="pipe"), cat.pressure_outlet(P0, T0)],
        edges=[(0, 1, area), (1, 2, area)],
    )


def solved_fanno(N):
    "Solve the fanno_pipe network at N segments and gate on convergence before returning."
    s = fanno_net(N).solve()
    assert s.converged, (s.residual_norm, s.print_residuals())
    return s


fig = go.Figure()
for N, col in zip((2, 8, 32), (COLORWAY[0], COLORWAY[2], COLORWAY[1])):
    s = solved_fanno(N)
    cv = s.composite("pipe")
    edges = [0] + list(cv.internal_edges) + [1]  # inlet, interior (flow order), outlet
    x = np.linspace(0.0, L, len(edges))
    M = s.field("M")
    fig.add_trace(go.Scatter(x=x, y=[M[e] for e in edges], mode="lines+markers", name=f"N = {N}", line_color=col))
fig.update_layout(
    title="fanno_pipe: the Mach profile resolves as N grows",
    xaxis_title="Axial position, m",
    yaxis_title="Mach number",
    width=820,
    height=440,
)
fig.show()

#  probe the resolving quantities: the inlet Mach and the integrated friction loss
ref = auto_refine(
    build=solved_fanno,
    n_start=16,
    probe=lambda s: {"M_in": float(s.field("M")[0]), "pt_drop": float(s.field("p_t")[0] - s.field("p_t")[1])},
)
print(f"auto-refinement to N={ref.n_final}: quantities change at most {ref.worst:.2e} (converged: {ref.converged})")
auto-refinement to N=32: quantities change at most 6.32e-05 (converged: True)

Result. Each segment carries a locally-uniform mean state, so the chain captures the continuous Mach rise: and a doubling of N barely moves the resolving quantities (the inlet Mach, the integrated loss), which is the built-in convergence check (auto_refine).

4: tapered_duct: a con-di nozzle chokes at its true throat

cat.tapered_duct(table) discretizes a continuously area-varying passage into compact area-change + duct segments. The standard input is a table of (x, A) pairs: the axial positions x [m] set the station spacing (which may be non-uniform, e.g. clustered near a throat) and the total length is inferred from them; a callable A(x) is also accepted. Each segment carries a real duct, so the taper propagates waves through its interior. A converging-diverging profile chokes at its narrowest station: the throat: which the composite’s throat edge identifies.

A_in, A_th, A_out, Ln = 3.0e-3, 1.5e-3, 3.0e-3, 0.4
N = 16
half = N // 2
areas = list(np.linspace(A_in, A_th, half + 1)) + list(np.linspace(A_th, A_out, N - half + 1))[1:]
table = list(zip(np.linspace(0.0, Ln, len(areas)), areas))  # (x, A) pairs; length inferred from x

net = nefes.Network(
    GAS,
    nodes=[cat.total_pressure_inlet(108000.0, T0), cat.tapered_duct(table, name="nozzle"), cat.pressure_outlet(P0, T0)],
    edges=[(0, 1, A_in), (1, 2, A_out)],
)
sol = net.solve()
assert sol.converged, (sol.residual_norm, sol.print_residuals())
cv = sol.composite("nozzle")

#  order the composite's edges along the axis (contract then expand)
edges = [0] + list(cv.internal_edges) + [1]
xs = np.linspace(0.0, Ln, len(edges))
M = [sol.edge(e)["M"] for e in edges]
A = [sol.edge(e)["area"] for e in edges]

fig = go.Figure()
fig.add_trace(
    go.Scatter(x=xs, y=M, mode="lines+markers", name="Mach", line=dict(color=COLORWAY[0], width=2), yaxis="y1")
)
fig.add_trace(
    go.Scatter(
        x=xs, y=np.array(A) * 1e3, name="area [10⁻³ m²]", line=dict(color="#8896a5", width=2, dash="dash"), yaxis="y2"
    )
)
fig.update_layout(
    title="tapered_duct: the con-di throat is the fastest, narrowest edge",
    xaxis_title="Axial position, m",
    yaxis=dict(title="Mach number"),
    yaxis2=dict(title="area, 10⁻³ m²", overlaying="y", side="right", showgrid=False),
    width=820,
    height=440,
)
fig.show()
print(
    f"throat edge {cv.throat}:  area {sol.edge(cv.throat)['area']:.2e} (min), "
    f"Mach {sol.edge(cv.throat)['M']:.4f} (max, subsonic)"
)
throat edge 16:  area 1.50e-03 (min), Mach 0.9065 (max, subsonic)

Result. The resolved Mach peaks at the geometric throat: the min-area station: with the isentropic-area-change complementarity engaging on exactly that segment. Pushing the inlet pressure further drives the throat to M = 1, the v1 (subsonic) choke limit.

5: sudden_contraction: the resolved vena contracta

A sudden contraction necks to a vena contracta of area cc * downstream_area, where the static pressure bottoms out, then re-expands with loss. cat.sudden_contraction(cc=...) resolves that state (the downstream area is read off the attached outflow edge) (isentropic to the vena contracta + Borda re-expansion), so the minimum static pressure and the compressible loss are exact: unlike the lumped O(M²) cc-loss head.

A_BIG, A_SMALL, CC = 4.0e-3, 1.0e-3, 0.62


def contraction(pt_in):
    net = nefes.Network(GAS)
    i = net.add(cat.total_pressure_inlet(pt_in, T0))
    sc = net.add(cat.sudden_contraction(cc=CC, name="contr"))
    o = net.add(cat.pressure_outlet(P0, T0))
    a = net.connect(i, sc, A_BIG)
    b = net.connect(sc, o, A_SMALL)
    sol = net.solve()
    assert sol.converged, (sol.residual_norm, sol.print_residuals())
    return sol, a, b


def sac_loss(pt_in):  # the incompressible O(M^2) reference
    net = nefes.Network(GAS)
    i = net.add(cat.total_pressure_inlet(pt_in, T0))
    s = net.add(cat.sudden_area_change(cc=CC))
    o = net.add(cat.pressure_outlet(P0, T0))
    a = net.connect(i, s, A_BIG)
    b = net.connect(s, o, A_SMALL)
    sol = net.solve()
    assert sol.converged, (sol.residual_norm, sol.print_residuals())
    pt = sol.field("p_t")
    return (pt[a] - pt[b]) / pt[a], float(sol.edge(b)["M"])


#  the static-pressure dip at the vena contracta
sol, a, b = contraction(130000.0)
th = sol.composites[0].throat
fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x=["inlet", "vena contracta", "outlet"],
        y=[sol.edge(a)["p"] / 1e5, sol.edge(th)["p"] / 1e5, sol.edge(b)["p"] / 1e5],
        mode="lines+markers",
        line=dict(color=COLORWAY[0], width=2),
    )
)
fig.update_layout(
    title="sudden_contraction: static pressure bottoms at the vena contracta",
    yaxis_title="Static pressure, bar",
    width=820,
    height=400,
)
fig.show()

#  compressible accuracy: the resolved loss diverges from the O(M^2) model as Mach rises
pts = np.linspace(102000.0, 130000.0, 12)
m_out, loss_sac, loss_res = [], [], []
for p in pts:
    ls, m = sac_loss(p)
    loss_sac.append(ls)
    m_out.append(m)
    s, aa, bb = contraction(p)
    pt = s.field("p_t")
    loss_res.append((pt[aa] - pt[bb]) / pt[aa])
fig2 = go.Figure()
fig2.add_trace(
    go.Scatter(
        x=m_out,
        y=loss_sac,
        name="sudden_area_change  (incompressible, O(M²))",
        line=dict(color="#8896a5", width=2, dash="dash"),
    )
)
fig2.add_trace(
    go.Scatter(
        x=m_out, y=loss_res, name="sudden_contraction  (resolved vena contracta)", line=dict(color=COLORWAY[0], width=2)
    )
)
fig2.update_layout(
    title="Compressible loss correction grows with Mach",
    xaxis_title="Downstream Mach",
    yaxis_title="Total-pressure loss fraction",
    width=820,
    height=440,
)
fig2.show()

Result. The vena contracta carries the lowest static pressure in the element: the cavitation-relevant minimum a lumped loss cannot report. The two loss models agree at low Mach and separate as the flow compresses, the resolved composite being exact within the v1 subsonic range.

Takeaway

composite class expands to reads back
orifice, lossy_nozzle 1 macro iac (+ iac) + sac throat state
helmholtz_resonator 1 macro tee + neck duct + cavity resonance
sudden_contraction 1 macro isentropic + Borda vena-contracta state
fanno_pipe 2 discretization N × pipe atom Mach profile, converges in N
tapered_duct 2 discretization N × (iac + duct) from an (x, A) table choke at the true throat

Every recipe rides the same element-independent expander; the solved network is byte-identical to a hand build, and solution.composite(name) projects the result back to the user-facing element.

Export for Nemo

Composite elements serialize to the UI format as the single node the user specified (an orifice saves as one Orifice node carrying its throatArea, never its expanded internals) and the UI’s Composite elements palette offers the same six elements, so a saved case round-trips both ways. Save with sol.to_yaml(path) (mean-flow results embedded) or network.to_yaml(path) (topology only), then open the file in Nemo; loading a UI-saved case back is nefes.load_case. The cell below exports the composite orifice line from Section 1 directly.

import os, tempfile

#  composites serialize as themselves: this case carries a single Orifice node
_out = os.path.join(tempfile.mkdtemp(), "composite_elements.yaml")
sc.to_yaml(_out)  # embeds the mean-flow results; use orifice_net.to_yaml(_out) for topology only
print("saved case:", _out)
saved case: /tmp/tmpkphngazg/composite_elements.yaml
Back to top