Skip to content

qlroms.intrusive_qlroms.pinball — fluidic pinball (FEniCSx)

Incompressible Navier–Stokes past three cylinders, solved with a dolfinx fractional-step scheme on Taylor–Hood P2/P1. The intrusive member uses mass-orthonormal velocity POD modes plus supremizers, UFL assembly of the projected momentum and pressure-Poisson equations, and pressure elimination through the PPE, so it has the same closed quadratic form as the KS members.

dolfinx is needed only to generate snapshots, run the FOM, or build a model on a cache miss: pinball, pinball.rom and pinball.config import without it, and cached models reload without it. Everything that touches dolfinx lives in aux_fenics and solver.

Fluidic pinball: velocity magnitude on the P2 nodes, three cylinders and the wake probes

\(|\mathbf{u}|\) at the first test snapshot of the \(Re = 90\) campaign (interpolated from the P2 velocity nodes of the Taylor–Hood mesh); white disks are the three cylinders, crosses the near-wake probes FOM observes (obs_idx).

At a glance

Name One-liner
config.PODGalerkinODEOperators Dolfinx-free container of one chart's reduced momentum + PPE operators; closed_quadratic_operators() eliminates the pressure. .Psi (the chart's pressure POD modes, centered on .pressure_center) plus .recover_pressure(a) reconstruct the physical pressure field from a reduced velocity coordinate; None on a model cached before this existed. .correct_divergence(a) is the minimal-norm projection of a onto the chart's own reduced-divergence constraint \(\mathbf{D}\bm{a}+\bm{d}=0\).
fom.FOM The solver as a dynamodels discrete map (state = stacked \([\bm{u}; p]\) dofs, obs_idx sensors, serial, \(m = 1\)).
rom.ROM The case's GalerkinROM: mass-orthonormal velocity chart (incl. supremizers) + pressure-eliminated \((\bm{b},\mathbf{A},\mathsf{B})\).
rom.build_local_model(Xv, Xp, *, K, r_velocity, r_pressure, reynolds, dt, save_dir, mesh_dir, ...) Cluster → weighted POD + supremizers → UFL Galerkin + PPE elimination → qlGalerkin; cached, dolfinx only on a miss. The snapshots are clustered and the states assigned to charts with the velocity mass matrix (cluster_metric="mass", the default; "euclidean" reproduces the earlier dof-distance partition), and for \(K>1\) the exact atlas is built in that metric. pressure_centroid (default "poisson_from_velocity_centroid", a Poisson solve consistent with the PPE the fluctuation obeys; "empirical" is the plain cluster mean) sets the point the pressure is expanded around -- the pressure POD basis and the supremizer enrichment are centered on that same choice. The cache name carries both the clustering metric (\(K>1\) only) and the pressure centroid (always).
rom.assemble_galerkin_operators(...) The UFL assembly of one chart's operators.
rom.free_run_divergence_free(model, x0, n_steps) Alias of qlroms.free_run: every pinball qlGalerkin step is predict (one ETDRK4 step) / correct (PODGalerkinODEOperators.correct_divergence), in step and forecast alike -- kills the reduced-divergence drift a plain free_run accumulates over a long closed loop, but is a kinematic projection, not a stabiliser: it does not by itself fix instability from an unstable eigenvalue of \(\mathbf{A}_{\text{eff}}\). Needs model.pressure_ops.
solver.PinballFOM, solver.PinballFOMParameters The fractional-step solver: run() for snapshot campaigns (MPI-capable), step() / get_state() / set_state() for the FOM adapter. wall_bc selects the channel side walls: "no_slip" (default, physical rigid walls) or "far_field" (walls reuse the inlet's \(\bm{u}=(1,0)\) profile, how the fom_re130-style production campaigns were generated); obstacles are always no-slip. Saved in snapshot_metadata.npz and in FOM's cache filename.
aux_fenics Mesh labels, loader and gmsh generator, taylor_hood_spaces, assemble_mass_matrix, dirichlet_bcs, supremizer and pressure-centroid solves.
build_pinball.load_snapshots Snapshot loading — see qlroms.intrusive_qlroms.
build_pinball.load_mass_matrix(kind, mesh_dir, n_dof) The FEM mass matrix (velocity or pressure) read without dolfinx from the file build_local_model or run_pinball_fom.py --export-mass writes next to the mesh; checks the dof count against the snapshots.

qlroms.intrusive_qlroms.pinball.config

pinball/config.py -- PODGalerkinODEOperators, the dolfinx-free container of one chart's reduced momentum + pressure-Poisson operators assembled by pinball.rom (UFL); cached models carry a list of them in model.pressure_ops, so this class must stay importable without dolfinx and keep its module path.

PODGalerkinODEOperators dataclass

Reduced continuous-time operators for one affine POD chart (pinball/rom.py):

a_dot = f + A a + B(a, a) + P b A_pr b = p0 - lap(p_c) + P1 a + P2(a, a)

Source code in qlroms/intrusive_qlroms/pinball/config.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
@dataclass
class PODGalerkinODEOperators:
    """Reduced continuous-time operators for one affine POD chart (pinball/rom.py):

        a_dot = f + A a + B(a, a) + P b
        A_pr b = p0 - lap(p_c) + P1 a + P2(a, a)
    """

    f: np.ndarray
    A: np.ndarray
    B: np.ndarray
    P: np.ndarray
    A_pr: np.ndarray
    pressure_rhs_const: np.ndarray
    pressure_rhs_linear: np.ndarray
    pressure_rhs_quadratic: np.ndarray
    pressure_center_laplacian: np.ndarray
    D: np.ndarray
    d: np.ndarray
    pressure_center: np.ndarray
    pressure_centroid_mode: str
    # Pressure POD modes (Np, r_pressure), Mw_p-orthonormal, centered on pressure_center
    # (see pinball.rom.assemble_galerkin_operators): p ~= pressure_center + Psi @ b. None
    # for a model cached before this field existed -- rebuild it to recover pressure fields.
    Psi: np.ndarray | None = None

    @property
    def r_velocity(self) -> int:
        return int(self.f.shape[0])

    @property
    def r_pressure(self) -> int:
        return int(self.A_pr.shape[0])

    def pressure_rhs(self, a: np.ndarray) -> np.ndarray:
        a = np.asarray(a, dtype=float).reshape(-1)
        return (
            self.pressure_rhs_const
            - self.pressure_center_laplacian
            + self.pressure_rhs_linear @ a
            + np.einsum("ljk,j,k->l", self.pressure_rhs_quadratic, a, a, optimize=True)
        )

    def solve_pressure(self, a: np.ndarray) -> np.ndarray:
        if self.A_pr.size == 0:
            return np.zeros((0,), dtype=float)
        return _solve_square_or_lstsq(self.A_pr, self.pressure_rhs(a))

    def rhs(self, a: np.ndarray, b: np.ndarray | None = None) -> np.ndarray:
        a = np.asarray(a, dtype=float).reshape(-1)
        if b is None:
            b = self.solve_pressure(a)
        return (
            self.f
            + self.A @ a
            + np.einsum("ijk,j,k->i", self.B, a, a, optimize=True)
            + self.P @ np.asarray(b, dtype=float).reshape(-1)
        )

    def divergence_residual(self, a: np.ndarray) -> np.ndarray:
        return self.D @ np.asarray(a, dtype=float).reshape(-1) + self.d

    @cached_property
    def _divergence_projection(self) -> tuple[np.ndarray, np.ndarray] | None:
        """(P, q) folding the minimal-norm projection onto D a + d = 0 into ONE affine map,
        correct_divergence(a) = P @ a + q -- D and d never change, so the (small, r_pressure
        x r_pressure) linear solve behind the projection is worth doing only once, offline,
        like every other per-chart operator this class hands the online step (e.g.
        QuadROM.etdrk4_rom), rather than repeating it every call in a forecast loop.
        None when the chart carries no pressure modes."""
        if self.D.size == 0:
            return None
        # K = D^T (D D^T)^{-1}, via K^T = (D D^T)^{-1} D (solve, not explicit inverse)
        K = _solve_square_or_lstsq(self.D @ self.D.T, self.D).T
        P = np.eye(self.D.shape[1]) - K @ self.D
        return P, -(K @ self.d)

    def correct_divergence(self, a: np.ndarray) -> np.ndarray:
        """Minimal-norm correction of a onto this chart's reduced-divergence constraint
        D a + d = 0 (divergence_residual(a_corrected) ~= 0): the corrector half of a
        predict (one ETDRK4 step of the closed quadratic ODE) / correct (this) forecast --
        see pinball.rom.free_run_divergence_free. Because Phi is Mw-orthonormal, the plain
        Euclidean least-squares correction on a IS already the minimal correction in the
        model's own physical (Mw) norm: ||Phi c||_Mw = ||c||_2 for any coefficient vector c.
        A no-op when the chart carries no pressure modes (r_pressure == 0)."""
        a = np.asarray(a, dtype=float).reshape(-1)
        if self._divergence_projection is None:
            return a
        P, q = self._divergence_projection
        return P @ a + q

    def recover_pressure(self, a: np.ndarray, b: np.ndarray | None = None) -> np.ndarray:
        """Physical pressure field p ~= pressure_center + Psi @ b (b solved from the reduced
        velocity coordinate a via the PPE, solve_pressure, if not given directly)."""
        if self.Psi is None:
            raise RuntimeError(
                "no pressure POD modes stored on this model (Psi is None) -- pressure "
                "recovery needs a model built after this field was added; rebuild it.")
        if b is None:
            b = self.solve_pressure(a)
        return self.pressure_center + self.Psi @ np.asarray(b, dtype=float).reshape(-1)

    def closed_quadratic_operators(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Pressure-eliminated (f_eff, A_eff, B_eff): substitute the PPE solution
        b(a) = A_pr^{-1}(p0 - lap(p_c) + P1 a + P2(a,a)) into the momentum equation."""
        if self.A_pr.size == 0 or self.P.size == 0:
            return self.f.copy(), self.A.copy(), self.B.copy()

        p0 = _solve_square_or_lstsq(self.A_pr, self.pressure_rhs_const - self.pressure_center_laplacian)
        p1 = _solve_square_or_lstsq(self.A_pr, self.pressure_rhs_linear)
        p2 = _solve_square_or_lstsq(
            self.A_pr, self.pressure_rhs_quadratic.reshape(self.r_pressure, -1)
        ).reshape(self.r_pressure, self.r_velocity, self.r_velocity)

        f_eff = self.f + self.P @ p0
        A_eff = self.A + self.P @ p1
        B_eff = self.B + np.einsum("il,ljk->ijk", self.P, p2, optimize=True)
        return f_eff, A_eff, B_eff

correct_divergence(a)

Minimal-norm correction of a onto this chart's reduced-divergence constraint D a + d = 0 (divergence_residual(a_corrected) ~= 0): the corrector half of a predict (one ETDRK4 step of the closed quadratic ODE) / correct (this) forecast -- see pinball.rom.free_run_divergence_free. Because Phi is Mw-orthonormal, the plain Euclidean least-squares correction on a IS already the minimal correction in the model's own physical (Mw) norm: ||Phi c||_Mw = ||c||_2 for any coefficient vector c. A no-op when the chart carries no pressure modes (r_pressure == 0).

Source code in qlroms/intrusive_qlroms/pinball/config.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
def correct_divergence(self, a: np.ndarray) -> np.ndarray:
    """Minimal-norm correction of a onto this chart's reduced-divergence constraint
    D a + d = 0 (divergence_residual(a_corrected) ~= 0): the corrector half of a
    predict (one ETDRK4 step of the closed quadratic ODE) / correct (this) forecast --
    see pinball.rom.free_run_divergence_free. Because Phi is Mw-orthonormal, the plain
    Euclidean least-squares correction on a IS already the minimal correction in the
    model's own physical (Mw) norm: ||Phi c||_Mw = ||c||_2 for any coefficient vector c.
    A no-op when the chart carries no pressure modes (r_pressure == 0)."""
    a = np.asarray(a, dtype=float).reshape(-1)
    if self._divergence_projection is None:
        return a
    P, q = self._divergence_projection
    return P @ a + q

recover_pressure(a, b=None)

Physical pressure field p ~= pressure_center + Psi @ b (b solved from the reduced velocity coordinate a via the PPE, solve_pressure, if not given directly).

Source code in qlroms/intrusive_qlroms/pinball/config.py 
116
117
118
119
120
121
122
123
124
125
def recover_pressure(self, a: np.ndarray, b: np.ndarray | None = None) -> np.ndarray:
    """Physical pressure field p ~= pressure_center + Psi @ b (b solved from the reduced
    velocity coordinate a via the PPE, solve_pressure, if not given directly)."""
    if self.Psi is None:
        raise RuntimeError(
            "no pressure POD modes stored on this model (Psi is None) -- pressure "
            "recovery needs a model built after this field was added; rebuild it.")
    if b is None:
        b = self.solve_pressure(a)
    return self.pressure_center + self.Psi @ np.asarray(b, dtype=float).reshape(-1)

closed_quadratic_operators()

Pressure-eliminated (f_eff, A_eff, B_eff): substitute the PPE solution b(a) = A_pr^{-1}(p0 - lap(p_c) + P1 a + P2(a,a)) into the momentum equation.

Source code in qlroms/intrusive_qlroms/pinball/config.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
def closed_quadratic_operators(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Pressure-eliminated (f_eff, A_eff, B_eff): substitute the PPE solution
    b(a) = A_pr^{-1}(p0 - lap(p_c) + P1 a + P2(a,a)) into the momentum equation."""
    if self.A_pr.size == 0 or self.P.size == 0:
        return self.f.copy(), self.A.copy(), self.B.copy()

    p0 = _solve_square_or_lstsq(self.A_pr, self.pressure_rhs_const - self.pressure_center_laplacian)
    p1 = _solve_square_or_lstsq(self.A_pr, self.pressure_rhs_linear)
    p2 = _solve_square_or_lstsq(
        self.A_pr, self.pressure_rhs_quadratic.reshape(self.r_pressure, -1)
    ).reshape(self.r_pressure, self.r_velocity, self.r_velocity)

    f_eff = self.f + self.P @ p0
    A_eff = self.A + self.P @ p1
    B_eff = self.B + np.einsum("il,ljk->ijk", self.P, p2, optimize=True)
    return f_eff, A_eff, B_eff

qlroms.intrusive_qlroms.pinball.fom

Pinball full-order model as a dynamodels.Model (discrete map).

FOM wraps the dolfinx fractional-step solver (pinball.solver.PinballFOM) behind the dynamodels Model protocol, observing the velocity-dof indices obs_idx (psi[obs_idx]).

Serial and single-member by design: the state is the stacked [velocity; pressure] dof vector and m = 1 (FOM ensembles are compute campaigns, not a Model feature). dolfinx is imported only inside _build_solver, so the module stays importable without a FEniCSx install (tests stub _build_solver).

FOM

Bases: Model

Fractional-step pinball Navier-Stokes solver as a dynamodels discrete map.

Parameters:

Name Type Description Default
psi0 ndarray

Initial stacked [velocity; pressure] state, (Nu + Np,); zeros = from rest.

 None
dt float

Output step; defaults to dt_solver (coarser -> the integrator interpolates).

 None
reynolds float

Reynolds number -- a declared dynamodels parameter (enters alpha0 and the dataset filename). SET it, don't assimilate it; synced to the solver's viscosity Constant at every time_step.

 60.0
mesh_dir path - like

Pregenerated mesh directory (see pinball.aux_fenics); None uses the loader default.

 None
obs_idx array_like of int

Sensor indices into the stacked state (velocity dofs; a qlROM observes the same sensors through its H = Phi[obs_idx] rows).

 (0,)
dt_solver float

Time step of the fractional-step scheme itself.

 0.01
wall_bc (no_slip, far_field)

Channel side-wall condition: "no_slip" (default) are physical rigid walls; "far_field" reuses the inlet's u=(1,0) profile there instead, the large-domain approximation the fom_re130-style production campaigns were actually generated with. Obstacles are always no-slip regardless. See pinball.solver.PinballFOMParameters.wall_bc.

 "no_slip"
Source code in qlroms/intrusive_qlroms/pinball/fom.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
class FOM(Model):
    """Fractional-step pinball Navier-Stokes solver as a dynamodels discrete map.

    Parameters
    ----------
    psi0 : np.ndarray, optional
        Initial stacked [velocity; pressure] state, ``(Nu + Np,)``; zeros = from rest.
    dt : float, optional
        Output step; defaults to `dt_solver` (coarser -> the integrator interpolates).
    reynolds : float
        Reynolds number -- a declared dynamodels parameter (enters alpha0 and the
        dataset filename). SET it, don't assimilate it; synced to the solver's
        viscosity Constant at every `time_step`.
    mesh_dir : path-like, optional
        Pregenerated mesh directory (see pinball.aux_fenics); None uses the loader default.
    obs_idx : array_like of int
        Sensor indices into the stacked state (velocity dofs; a qlROM observes the
        same sensors through its H = Phi[obs_idx] rows).
    dt_solver : float
        Time step of the fractional-step scheme itself.
    wall_bc : {"no_slip", "far_field"}
        Channel side-wall condition: "no_slip" (default) are physical rigid walls;
        "far_field" reuses the inlet's u=(1,0) profile there instead, the large-domain
        approximation the fom_re130-style production campaigns were actually generated
        with. Obstacles are always no-slip regardless. See
        pinball.solver.PinballFOMParameters.wall_bc.
    """

    t_transient = 0.0
    t_CR = 0.0            # set at construction (defaults to 100 output steps)

    params = ['reynolds']
    # Structural parameters: declaring them lets ntsa.respawn / Model.reset_model
    # rebuild a FOM; dt_solver/wall_bc need the scalar CLASS default below to land in
    # the filename encoding (Model.filename only encodes scalar class-level defaults).
    fixed_params = ['mesh_dir', 'obs_idx', 'dt_solver', 'wall_bc']

    reynolds = 60.0
    dt_solver = 0.01
    wall_bc = "no_slip"

    def __init__(self, psi0=None, dt=None, reynolds=60.0, mesh_dir=None,
                 obs_idx=(0,), dt_solver=0.01, wall_bc="no_slip", **model_kwargs):
        if wall_bc not in WALL_BC_MODES:
            raise ValueError(f"wall_bc must be one of {WALL_BC_MODES}, got {wall_bc!r}.")
        self.reynolds = float(reynolds)
        self.mesh_dir = mesh_dir
        self.dt_solver = float(dt_solver)
        self.wall_bc = wall_bc
        self.obs_idx = np.asarray(obs_idx, dtype=int).reshape(-1)
        self.Nq = len(self.obs_idx)

        self.solver = self._build_solver()
        if getattr(self.solver, "comm", None) is not None and self.solver.comm.size > 1:
            raise RuntimeError("FOM (the dynamodels adapter) is serial; run MPI-parallel "
                               "snapshot generation through PinballFOM.run() instead.")
        self._Nu, self._Np = self.solver.num_dofs

        N = self._Nu + self._Np
        if psi0 is None:
            psi0 = np.zeros(N)     # from-rest start, as in PinballFOM.run()
        psi0 = np.asarray(psi0, dtype=np.float64)
        if psi0.shape[0] != N:
            raise ValueError(f"psi0 must have {N} rows ({self._Nu} velocity + "
                             f"{self._Np} pressure dofs), got {psi0.shape[0]}.")
        if self.obs_idx.size and (self.obs_idx.min() < 0 or self.obs_idx.max() >= N):
            raise ValueError(f"obs_idx out of range for the {N}-dof state.")

        # Re-seed the solver only when the tracked state changed under us (time_step).
        self._synced = False
        self._last_state = None

        dt_eff = float(dt) if dt is not None else self.dt_solver
        model_kwargs.setdefault('t_CR', 100 * dt_eff)
        super().__init__(psi0=psi0, dt=dt_eff,
                         integrator_class=DiscreteIntegrator, **model_kwargs)

    # ---- solver construction (the only dolfinx-touching path) ----
    def _build_solver(self):
        """Build the dolfinx fractional-step solver; requires a FEniCSx install."""
        from .aux_fenics import load_pinball_mesh
        from .solver import PinballFOM, PinballFOMParameters
        mesh, facet_tags, labels = load_pinball_mesh(self.mesh_dir)
        parameters = PinballFOMParameters(dt=self.dt_solver, num_steps=0,
                                          reynolds=self.reynolds, wall_bc=self.wall_bc)
        return PinballFOM(mesh, facet_tags, labels, parameters, snapshot_dir=None)

    # ---- discrete-map timing ----
    @property
    def _precision_t_step(self):
        """int: decimal precision of the solver's step (same rule Model applies to dt)."""
        return int(np.ceil(-np.log10(self.dt_solver) + 2))

    @property
    def dt_step(self):
        """float: solver step rounded as `Model` rounds `dt`, so dt == dt_solver
        compares equal inside `DiscreteIntegrator` (coarser dt -> interpolation)."""
        return float(np.round(self.dt_solver, self._precision_t_step))

    # ---- discrete map ----
    def time_step(self, Nt=1, **kwargs):
        """Advance the solver `Nt` fractional steps from `current_state`.

        Returns (psi, t): trajectory ``(Nt + 1, Nphi, 1)`` with ``psi[0]`` the
        current state, and the ``(Nt + 1,)`` time grid.
        """
        s = self.current_state
        if s.shape[1] != 1:
            raise NotImplementedError("the pinball FOM is single-member (m = 1).")
        state = np.asarray(s[:, 0], dtype=np.float64)

        # Re-seed only when the tracked state was changed externally (fresh model,
        # DA analysis, interpolated grid): a re-seed collapses the AB2 convection
        # history onto the given field (warm-restart approximation), while
        # back-to-back windows keep the solver's exact multistep history.
        if not (self._synced and self._last_state is not None
                and np.array_equal(self._last_state, state)):
            self.solver.set_state(state[:self._Nu], state[self._Nu:])
        self.solver.set_reynolds(self.reynolds)
        # keep the solver clock on the model clock (inlet profiles may be time-dependent)
        self.solver.time = float(self.current_time)

        psi = np.empty((Nt + 1, self.Nphi, 1))
        psi[0, :, 0] = state
        for n in range(Nt):
            self.solver.step()
            psi[n + 1, :, 0] = self.solver.get_state()
        self._last_state = psi[-1, :, 0].copy()
        self._synced = True

        prec = max(self.precision_t, self._precision_t_step)
        t = np.round(self.current_time + np.arange(Nt + 1) * self.dt_step, prec)
        return psi, t

    # ---- observables ----
    def get_observables(self, Nt=1, **kwargs):
        """Sensor values: the `obs_idx` rows of the trailing state(s)."""
        if Nt == 1:
            return self.hist[-1, self.obs_idx, :]
        return self.hist[-Nt:, self.obs_idx, :]

    @property
    def obs_labels(self):
        r"""list of str: LaTeX labels for the sensor observables, $u(x_i)$."""
        return [f'$u(x_{{{i}}})$' for i in self.obs_idx]

    # ---- copying ----
    def __deepcopy__(self, memo):
        """Deep copy everything EXCEPT the solver (PETSc objects do not deepcopy):
        the copy rebuilds it and re-seeds from current_state on the first time_step."""
        cls = self.__class__
        new = cls.__new__(cls)
        memo[id(self)] = new
        for key, value in self.__dict__.items():
            if key == 'solver':
                continue
            setattr(new, key, deepcopy(value, memo))
        new.solver = new._build_solver()
        new._synced = False
        new._last_state = None
        return new

    # ---- ensembles ----
    def init_ensemble(self, m=1, est_alpha=None, std_alpha=0.001, **kwargs):
        if m > 1:
            raise NotImplementedError(
                "FOM ensembles are full CFD campaigns; run ensembles on a fitted "
                "qlROM and keep the FOM as the truth model.")
        if est_alpha or isinstance(std_alpha, dict):
            raise NotImplementedError(
                "parameter estimation on the FOM is not supported (time_step reads the "
                "reynolds attribute, not an augmented state row): set reynolds instead.")
        return super().init_ensemble(m=m, est_alpha=[], std_alpha=std_alpha, **kwargs)

dt_step property

float: solver step rounded as Model rounds dt, so dt == dt_solver compares equal inside DiscreteIntegrator (coarser dt -> interpolation).

obs_labels property

list of str: LaTeX labels for the sensor observables, \(u(x_i)\).

time_step(Nt=1, **kwargs)

Advance the solver Nt fractional steps from current_state.

Returns (psi, t): trajectory (Nt + 1, Nphi, 1) with psi[0] the current state, and the (Nt + 1,) time grid.

Source code in qlroms/intrusive_qlroms/pinball/fom.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def time_step(self, Nt=1, **kwargs):
    """Advance the solver `Nt` fractional steps from `current_state`.

    Returns (psi, t): trajectory ``(Nt + 1, Nphi, 1)`` with ``psi[0]`` the
    current state, and the ``(Nt + 1,)`` time grid.
    """
    s = self.current_state
    if s.shape[1] != 1:
        raise NotImplementedError("the pinball FOM is single-member (m = 1).")
    state = np.asarray(s[:, 0], dtype=np.float64)

    # Re-seed only when the tracked state was changed externally (fresh model,
    # DA analysis, interpolated grid): a re-seed collapses the AB2 convection
    # history onto the given field (warm-restart approximation), while
    # back-to-back windows keep the solver's exact multistep history.
    if not (self._synced and self._last_state is not None
            and np.array_equal(self._last_state, state)):
        self.solver.set_state(state[:self._Nu], state[self._Nu:])
    self.solver.set_reynolds(self.reynolds)
    # keep the solver clock on the model clock (inlet profiles may be time-dependent)
    self.solver.time = float(self.current_time)

    psi = np.empty((Nt + 1, self.Nphi, 1))
    psi[0, :, 0] = state
    for n in range(Nt):
        self.solver.step()
        psi[n + 1, :, 0] = self.solver.get_state()
    self._last_state = psi[-1, :, 0].copy()
    self._synced = True

    prec = max(self.precision_t, self._precision_t_step)
    t = np.round(self.current_time + np.arange(Nt + 1) * self.dt_step, prec)
    return psi, t

get_observables(Nt=1, **kwargs)

Sensor values: the obs_idx rows of the trailing state(s).

Source code in qlroms/intrusive_qlroms/pinball/fom.py 
156
157
158
159
160
def get_observables(self, Nt=1, **kwargs):
    """Sensor values: the `obs_idx` rows of the trailing state(s)."""
    if Nt == 1:
        return self.hist[-1, self.obs_idx, :]
    return self.hist[-Nt:, self.obs_idx, :]

__deepcopy__(memo)

Deep copy everything EXCEPT the solver (PETSc objects do not deepcopy): the copy rebuilds it and re-seeds from current_state on the first time_step.

Source code in qlroms/intrusive_qlroms/pinball/fom.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def __deepcopy__(self, memo):
    """Deep copy everything EXCEPT the solver (PETSc objects do not deepcopy):
    the copy rebuilds it and re-seeds from current_state on the first time_step."""
    cls = self.__class__
    new = cls.__new__(cls)
    memo[id(self)] = new
    for key, value in self.__dict__.items():
        if key == 'solver':
            continue
        setattr(new, key, deepcopy(value, memo))
    new.solver = new._build_solver()
    new._synced = False
    new._last_state = None
    return new

qlroms.intrusive_qlroms.pinball.rom

pinball/rom.py -- the intrusive Galerkin equations for the pinball (Navier-Stokes) case, mirroring ks1d/rom.py: ROM (the case's single-cluster member), qlROM, the governing-equation projection (assemble_galerkin_operators) and build_local_model, which caches itself and needs dolfinx only on a cache miss (all dolfinx/ufl imports are local to the two build functions, so cached models reload and this module imports without FEniCSx). Everything generic (charts, clustering, POD, compilation, DA) comes from qlroms.

The reduced model, per affine chart (velocity modes Phi mass-orthonormal, incl. supremizers; velocity centered on the cluster mean; pressure expanded around the selected pressure centroid -- the pressure POD basis and the supremizers are built around that SAME centroid, not the empirical mean, whenever the two differ):

a_dot   = f + A a + B(a, a) + P b          (projected momentum equation)
A_pr b  = p0 - lap(p_c) + P1 a + P2(a, a)  (pressure-Poisson equation, PPE)

with the UFL assembly below providing f = -(u_c.grad u_c, phi_i) + (p_c, div phi_i) - nu (grad u_c, grad phi_i) A = -(u_c.grad phi_j + phi_j.grad u_c, phi_i) - nu (grad phi_j, grad phi_i) B = -(phi_j.grad phi_k, phi_i) P = D_r^T, D_r[l, i] = (div phi_i, psi_l) A_pr = (grad psi_k, grad psi_j), PPE right side from -(conv, grad psi_l).

Eliminating pressure through the PPE closes the velocity dynamics into the affine- quadratic form a_dot = b + A a + B(a, a) every qlroms family shares (PODGalerkinODEOperators.closed_quadratic_operators) -- a compilation of those members IS a qlGalerkin (build_local_model below; K=1 is the global case).

ROM

Bases: GalerkinROM

Pinball single-cluster member: a GalerkinROM whose chart is the cluster's mass-orthonormal velocity POD basis plus supremizers (Mw = the assembled FEM mass matrix) and whose (b, A, B) are the pressure-eliminated projected Navier-Stokes operators (PODGalerkinODEOperators.closed_quadratic_operators).

Source code in qlroms/intrusive_qlroms/pinball/rom.py 
60
61
62
63
64
class ROM(GalerkinROM):
    """Pinball single-cluster member: a GalerkinROM whose chart is the cluster's
    mass-orthonormal velocity POD basis plus supremizers (Mw = the assembled FEM mass
    matrix) and whose (b, A, B) are the pressure-eliminated projected Navier-Stokes
    operators (PODGalerkinODEOperators.closed_quadratic_operators)."""

assemble_galerkin_operators(Umodes, Pmodes, u_mean, pressure_center, nu, pressure_centroid_mode='poisson_from_velocity_centroid')

UFL-assemble one chart's reduced momentum + PPE operators (module docstring).

Umodes must already be mass-orthonormal (incl. supremizers); u_mean is the chart centroid; the pressure expansion is around pressure_center (all dolfinx Functions).

Source code in qlroms/intrusive_qlroms/pinball/rom.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def assemble_galerkin_operators(
    Umodes,
    Pmodes,
    u_mean,
    pressure_center,
    nu: float,
    pressure_centroid_mode: str = "poisson_from_velocity_centroid",
) -> PODGalerkinODEOperators:
    """UFL-assemble one chart's reduced momentum + PPE operators (module docstring).

    Umodes must already be mass-orthonormal (incl. supremizers); u_mean is the chart
    centroid; the pressure expansion is around `pressure_center` (all dolfinx Functions).
    """
    from dolfinx.fem import assemble_scalar, form
    from ufl import div, dot, dx, grad, inner, nabla_grad

    def _scalar(expr):
        return assemble_scalar(form(expr))

    ru, rp = len(Umodes), len(Pmodes)

    # -- momentum equation: viscous stiffness, convective tensor, mean-flow terms --
    K_r = np.zeros((ru, ru))
    C = np.zeros((ru, ru, ru))
    A_lin = np.zeros((ru, ru))
    f = np.zeros(ru)
    conv_mean = dot(u_mean, nabla_grad(u_mean))
    for i, phi_i in enumerate(Umodes):
        print(f"    [rom] Momentum test mode {i + 1}/{ru}", flush=True)
        f[i] = (
            -_scalar(inner(conv_mean, phi_i) * dx)
            + _scalar(pressure_center * div(phi_i) * dx)
            - nu * _scalar(inner(grad(u_mean), grad(phi_i)) * dx)
        )
        for j, phi_j in enumerate(Umodes):
            K_r[i, j] = _scalar(inner(grad(phi_j), grad(phi_i)) * dx)
            A_lin[i, j] = -(
                _scalar(inner(dot(u_mean, nabla_grad(phi_j)), phi_i) * dx)
                + _scalar(inner(dot(phi_j, nabla_grad(u_mean)), phi_i) * dx)
            )
            for k, phi_k in enumerate(Umodes):
                C[i, j, k] = _scalar(inner(dot(phi_j, nabla_grad(phi_k)), phi_i) * dx)

    # -- pressure: Poisson stiffness, divergence coupling, PPE right side --
    A_pr = np.zeros((rp, rp))
    D_r = np.zeros((rp, ru))
    D_mean = np.zeros(rp)
    p0 = np.zeros(rp)
    p1 = np.zeros((rp, ru))
    p2 = np.zeros((rp, ru, ru))
    center_laplacian = np.zeros(rp)
    for l, psi_l in enumerate(Pmodes):
        print(f"    [rom] Pressure test mode {l + 1}/{rp}", flush=True)
        grad_psi = grad(psi_l)
        for m, psi_m in enumerate(Pmodes):
            A_pr[l, m] = _scalar(inner(grad(psi_m), grad_psi) * dx)
        D_mean[l] = _scalar(div(u_mean) * psi_l * dx)
        p0[l] = _scalar(-inner(conv_mean, grad_psi) * dx)
        center_laplacian[l] = _scalar(inner(grad(pressure_center), grad_psi) * dx)
        for j, phi_j in enumerate(Umodes):
            D_r[l, j] = _scalar(div(phi_j) * psi_l * dx)
            conv_linear = dot(u_mean, nabla_grad(phi_j)) + dot(phi_j, nabla_grad(u_mean))
            p1[l, j] = _scalar(-inner(conv_linear, grad_psi) * dx)
            for k, phi_k in enumerate(Umodes):
                p2[l, j, k] = _scalar(-inner(dot(phi_j, nabla_grad(phi_k)), grad_psi) * dx)

    Np = pressure_center.x.array.shape[0]
    Psi = (np.stack([psi.x.array.copy() for psi in Pmodes], axis=1) if Pmodes
          else np.zeros((Np, 0)))

    return PODGalerkinODEOperators(
        f=f,
        A=A_lin - nu * K_r,
        B=-C,
        P=D_r.T.copy(),
        A_pr=A_pr,
        pressure_rhs_const=p0,
        pressure_rhs_linear=p1,
        pressure_rhs_quadratic=p2,
        pressure_center_laplacian=center_laplacian,
        D=D_r,
        d=D_mean,
        pressure_center=pressure_center.x.array.copy(),
        pressure_centroid_mode=pressure_centroid_mode,
        Psi=Psi,
    )

intrusive_model_path(K, r_velocity, r_pressure, reynolds, *, Ntrain, dt, save_dir=None, random_state=1, cluster_space='pod_lossless', cluster_metric='mass', pressure_centroid='poisson_from_velocity_centroid', tag=None)

Cache file that build_local_model reads and writes for these settings.

Every argument that changes the fitted chart is part of the name, so two builds that differ in any of them never collide on the same file: Ntrain (the training window length -- a cache keyed only on (K, r_velocity, r_pressure, reynolds) would otherwise silently hand back a model fit on a different window/dt/seed than the one requested) and dt (baked into the cached model's own ETDRK4 stepping) always; random_state and cluster_space only for K > 1, where they actually affect anything (K = 1 has no partition to seed). The pressure-centroid tag (_pcpoisson / _pcemp) is always present -- it changes the pressure POD basis and supremizers of every chart, K = 1 included (see build_local_model's pressure_centroid). tag is the snapshot campaign's tag (fom_re140__, qlroms.intrusive_qlroms.build_pinball.campaign_name), appended as __<tag>: a repeated campaign at the same Reynolds number with the same settings would otherwise silently reload the model fit on the untagged campaign's snapshots. None (an untagged campaign) keeps the historical name, so existing caches still load.

Source code in qlroms/intrusive_qlroms/pinball/rom.py 
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def intrusive_model_path(K: int, r_velocity: int, r_pressure: int, reynolds: float, *,
                         Ntrain: int, dt: float, save_dir=None,
                         random_state: int = 1, cluster_space: str = "pod_lossless",
                         cluster_metric: ClusterMetric = "mass",
                         pressure_centroid: PressureCentroidMode = "poisson_from_velocity_centroid",
                         tag: str | None = None,
                         ) -> Path:
    """Cache file that `build_local_model` reads and writes for these settings.

    Every argument that changes the fitted chart is part of the name, so two builds that
    differ in any of them never collide on the same file: `Ntrain` (the training window
    length -- a cache keyed only on (K, r_velocity, r_pressure, reynolds) would otherwise
    silently hand back a model fit on a different window/dt/seed than the one requested)
    and `dt` (baked into the cached model's own ETDRK4 stepping) always; `random_state` and
    `cluster_space` only for K > 1, where they actually affect anything (K = 1 has no
    partition to seed). The pressure-centroid tag (`_pcpoisson` / `_pcemp`) is always
    present -- it changes the pressure POD basis and supremizers of every chart, K = 1
    included (see build_local_model's `pressure_centroid`). `tag` is the snapshot campaign's
    tag (fom_re140__<tag>, qlroms.intrusive_qlroms.build_pinball.campaign_name), appended as
    `__<tag>`: a repeated campaign at the same Reynolds number with the same settings would
    otherwise silently reload the model fit on the untagged campaign's snapshots. None (an
    untagged campaign) keeps the historical name, so existing caches still load."""
    save_dir = Path(save_dir) if save_dir is not None else PINBALL_MODELS
    cluster_tag = ("" if K == 1 else
                  f"_{cluster_space}_seed{random_state}_{'mass' if cluster_metric == 'mass' else 'euclid'}")
    pc_tag = "_pcpoisson" if pressure_centroid == "poisson_from_velocity_centroid" else "_pcemp"
    return (save_dir /
           f"intrusive_model_K{K}_rv{r_velocity}_rp{r_pressure}_Re{reynolds:g}"
           f"_Ntr{Ntrain}_dt{dt:g}{cluster_tag}{pc_tag}{f'__{tag}' if tag else ''}.pth")

build_local_model(Xv, Xp, *, K, r_velocity, r_pressure, reynolds, dt, save_dir=None, mesh_dir=None, pressure_centroid='poisson_from_velocity_centroid', random_state=1, cluster_space='pod_lossless', cluster_metric='mass', tag=None)

Offline intrusive build (cached): cluster velocity snapshots (qlroms.fit_clusters), per cluster do mass-weighted POD (velocity + pressure), supremizer enrichment, mass-re-orthonormalization, and UFL Galerkin operator assembly; compile the pressure-eliminated members into a qlGalerkin (K=1 is the global case).

Parameters:

Name Type Description Default
Xv, Xp

aligned (Nu, Nt) velocity and (Np, Nt) pressure snapshot matrices -- chronological, the caller's training window.

 required
K, r_velocity, r_pressure

clusters and POD ranks. Each chart's final velocity dimension is r_velocity + (r_pressure + 1) supremizers.

 required
reynolds, dt

physics + snapshot spacing (dt sizes the ETDRK4 stepping).

 required
save_dir

cache directory (default qlroms.utils.paths.PINBALL_MODELS); the cache file (see intrusive_model_path) is checked FIRST, so a cached model reloads without dolfinx.

 None
mesh_dir

pinball mesh directory (default PINBALL_MESHES), loaded only on a cache miss.

 None
pressure_centroid PressureCentroidMode

"poisson_from_velocity_centroid" (default: a Poisson solve driven by the velocity centroid, consistent with the same PPE the fluctuation obeys) or "empirical" (the cluster's plain pressure mean). The pressure POD basis and the supremizer enrichment are built around whichever centroid is chosen, not always the empirical mean, so this changes the fitted chart materially -- part of the cache filename (intrusive_model_path).

 'poisson_from_velocity_centroid'
random_state, cluster_space

forwarded to the clustering.

 required
cluster_metric ClusterMetric

"mass" (default) clusters the snapshots and assigns states to charts with the FEM velocity mass matrix, the same inner product as the POD and the transition maps; for K > 1 the exact atlas is built too. "euclidean" reproduces the earlier plain dof-distance partition and affiliation (for comparisons only); transitions stay exact in the mass metric either way.

 'mass'
tag str | None

the snapshot campaign's tag (None for an untagged campaign), so a repeated campaign at the same Reynolds number gets its own cache file.

 None

Returns:

Type Description
qlGalerkin

qlGalerkin compilation of K ROMs (Mw = velocity mass matrix);
qlGalerkin

.pressure_ops[k] the per-chart PODGalerkinODEOperators -- .recover_pressure(a)
qlGalerkin

reconstructs the physical pressure field (solving the PPE for b if not given
qlGalerkin

directly, via .Psi, the chart's pressure POD modes) -- plus .reynolds, .Ntrain.
Source code in qlroms/intrusive_qlroms/pinball/rom.py 
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
def build_local_model(
    Xv,
    Xp,
    *,
    K: int,
    r_velocity: int,
    r_pressure: int,
    reynolds: float,
    dt: float,
    save_dir=None,
    mesh_dir=None,
    pressure_centroid: PressureCentroidMode = "poisson_from_velocity_centroid",
    random_state: int = 1,
    cluster_space: str = "pod_lossless",
    cluster_metric: ClusterMetric = "mass",
    tag: str | None = None,
) -> qlGalerkin:
    """Offline intrusive build (cached): cluster velocity snapshots (qlroms.fit_clusters),
    per cluster do mass-weighted POD (velocity + pressure), supremizer enrichment,
    mass-re-orthonormalization, and UFL Galerkin operator assembly; compile the
    pressure-eliminated members into a qlGalerkin (K=1 is the global case).

    Args:
        Xv, Xp: aligned (Nu, Nt) velocity and (Np, Nt) pressure snapshot matrices --
            chronological, the caller's training window.
        K, r_velocity, r_pressure: clusters and POD ranks. Each chart's final
            velocity dimension is r_velocity + (r_pressure + 1) supremizers.
        reynolds, dt: physics + snapshot spacing (dt sizes the ETDRK4 stepping).
        save_dir: cache directory (default qlroms.utils.paths.PINBALL_MODELS); the
            cache file (see intrusive_model_path) is checked FIRST, so a cached model
            reloads without dolfinx.
        mesh_dir: pinball mesh directory (default PINBALL_MESHES), loaded only on a
            cache miss.
        pressure_centroid: "poisson_from_velocity_centroid" (default: a Poisson solve driven
            by the velocity centroid, consistent with the same PPE the fluctuation obeys) or
            "empirical" (the cluster's plain pressure mean). The pressure POD basis and the
            supremizer enrichment are built around whichever centroid is chosen, not always
            the empirical mean, so this changes the fitted chart materially -- part of the
            cache filename (intrusive_model_path).
        random_state, cluster_space: forwarded to the clustering.
        cluster_metric: "mass" (default) clusters the snapshots and assigns states to charts
            with the FEM velocity mass matrix, the same inner product as the POD and the
            transition maps; for K > 1 the exact atlas is built too. "euclidean" reproduces
            the earlier plain dof-distance partition and affiliation (for comparisons
            only); transitions stay exact in the mass metric either way.
        tag: the snapshot campaign's tag (None for an untagged campaign), so a repeated
            campaign at the same Reynolds number gets its own cache file.

        Every one of the arguments above (plus Xv's own column count as Ntrain) is part of
        the cache filename (see intrusive_model_path) -- no two builds that differ in any of
        them ever collide on the same cached file.

    Returns:
        qlGalerkin compilation of K ROMs (Mw = velocity mass matrix);
        `.pressure_ops[k]` the per-chart PODGalerkinODEOperators -- `.recover_pressure(a)`
        reconstructs the physical pressure field (solving the PPE for b if not given
        directly, via `.Psi`, the chart's pressure POD modes) -- plus `.reynolds`, `.Ntrain`.
    """
    cache = intrusive_model_path(K, r_velocity, r_pressure, reynolds, Ntrain=int(Xv.shape[1]), dt=dt,
                                 save_dir=save_dir, random_state=random_state, cluster_space=cluster_space,
                                 cluster_metric=cluster_metric, pressure_centroid=pressure_centroid,
                                 tag=tag)
    cache.parent.mkdir(parents=True, exist_ok=True)
    if cache.exists():
        try:
            return torch.load(cache, map_location="cpu", weights_only=False)
        except Exception as e:
            print(f"Failed to load cached model (will rebuild): {e}")

    from . import aux_fenics as fx  # dolfinx only to BUILD

    mesh, facet_tags, labels = fx.load_pinball_mesh(mesh_dir)
    V, Q = fx.taylor_hood_spaces(mesh)

    nu = 1.0 / float(reynolds)
    Xv = torch.as_tensor(Xv, dtype=torch.float64)
    Xp = torch.as_tensor(Xp, dtype=torch.float64)
    if Xv.shape[1] != Xp.shape[1]:
        raise ValueError("velocity and pressure snapshot counts differ")

    M_u = fx.assemble_mass_matrix(V)
    M_p = fx.assemble_mass_matrix(Q)
    Mw_u = mass_to_torch_sparse(M_u)
    Mw_p = mass_to_torch_sparse(M_p)
    fx.export_mass_matrices(M_u, M_p, mesh_dir)       # sidecars for the dolfinx-free consumers
    walls = ("inlet", "walls", "obstacle1", "obstacle2", "obstacle3")
    bcs_hom = fx.dirichlet_bcs(V, facet_tags, [labels[k] for k in walls if k in labels])
    bcs_p = fx.dirichlet_bcs(Q, facet_tags, [labels["outlet"]])

    cluster_labels = _cluster_snapshots(Xv, K, r_velocity, cluster_space=cluster_space,
                                        random_state=random_state,
                                        Mw=Mw_u if cluster_metric == "mass" else None)
    sizes = [int((cluster_labels == k).sum()) for k in range(K)]
    print(f"  [rom] intrusive build: K={K} clusters, sizes={sizes}", flush=True)

    r_total = r_velocity + r_pressure + 1          # POD modes + supremizers
    Phi_all = torch.zeros(Xv.shape[0], r_total, K, dtype=torch.float64)
    centroids = torch.zeros(K, Xv.shape[0], dtype=torch.float64)
    b_all = torch.zeros(K, r_total, dtype=torch.float64)
    A_all = torch.zeros(K, r_total, r_total, dtype=torch.float64)
    B_all = torch.zeros(K, r_total, r_total, r_total, dtype=torch.float64)
    pressure_ops: list[PODGalerkinODEOperators] = []

    for k in range(K):
        cols = torch.nonzero(cluster_labels == k, as_tuple=False).flatten()
        Xk, Pk = Xv[:, cols], Xp[:, cols]
        u_mean_col = Xk.mean(dim=1)
        p_mean_col = Pk.mean(dim=1)

        # The pressure centroid is needed BEFORE the pressure POD basis: whichever centroid is
        # chosen is what the PPE expands p around (assemble_galerkin_operators below), so the
        # POD basis and the supremizer enrichment must be centered on that SAME field -- fitting
        # them to the empirical mean while the equations expand around a different centroid
        # (e.g. the Poisson one, which can differ from the mean by a large fraction) makes the
        # basis suboptimal for the fluctuation it is actually asked to represent.
        u_mean_fn = fx.function_from_array(V, u_mean_col.numpy())
        p_mean_fn = fx.function_from_array(Q, p_mean_col.numpy())
        if pressure_centroid == "empirical":
            p_center_fn = p_mean_fn
        elif pressure_centroid == "poisson_from_velocity_centroid":
            p_center_fn = fx.solve_pressure_centroid_from_velocity(Q, u_mean_fn, bcs_p)
        else:
            raise ValueError(f"Unsupported pressure centroid mode: {pressure_centroid}")
        p_center_col = torch.from_numpy(p_center_fn.x.array.copy()).to(dtype=Pk.dtype)

        print(f"  [rom] cluster {k}: POD (rv={r_velocity}, rp={r_pressure}) on {cols.numel()} snapshots",
              flush=True)
        Phi_u = compute_pod_basis(Xk - u_mean_col[:, None], r_velocity, Mw=Mw_u)
        Phi_p = compute_pod_basis(Pk - p_center_col[:, None], r_pressure, Mw=Mw_p)
        if Phi_u.shape[1] < r_velocity or Phi_p.shape[1] < r_pressure:
            raise ValueError(f"cluster {k} has too few snapshots for the requested ranks.")

        Pmode_fns = [fx.function_from_array(Q, Phi_p[:, j].numpy()) for j in range(Phi_p.shape[1])]
        supremizers = fx.compute_supremizers(V, [p_center_fn] + Pmode_fns, bcs_hom)
        sup_cols = torch.from_numpy(np.stack([w.x.array.copy() for w in supremizers], axis=1))

        U_full = mass_orthonormalize(torch.cat([Phi_u, sup_cols], dim=1), Mw=Mw_u)
        if U_full.shape[1] != r_total:
            raise ValueError(
                f"cluster {k}: orthonormalization kept {U_full.shape[1]} of {r_total} velocity "
                "modes (near-dependent supremizers); reduce r_velocity/r_pressure."
            )

        ops = assemble_galerkin_operators(
            [fx.function_from_array(V, U_full[:, j].numpy()) for j in range(U_full.shape[1])],
            Pmode_fns, u_mean_fn, p_center_fn, nu,
            pressure_centroid_mode=pressure_centroid,
        )
        f_eff, A_eff, B_eff = ops.closed_quadratic_operators()

        Phi_all[:, :, k] = U_full
        centroids[k] = u_mean_col
        b_all[k] = torch.from_numpy(f_eff)
        A_all[k] = torch.from_numpy(A_eff)
        B_all[k] = torch.from_numpy(B_eff)
        pressure_ops.append(ops)

    # the quantized-local model is just the compilation of the K single-cluster ROMs
    roms = []
    for k in range(K):
        rom = ROM(Phi_all[:, :, k], centroids[k], dt=dt, Mw=Mw_u)
        rom.set_operators(b_all[k], A_all[k], B_all[k])
        roms.append(rom)
    model = qlGalerkin(roms, dt=dt, qbar_g=Xv.mean(dim=1) if K > 1 else None)    # K > 1: exact atlas in the mass metric
    if K > 1 and cluster_metric == "euclidean":
        model.tree = CentroidIndex(model.centroids, None)
    model.cluster_metric = cluster_metric
    model.b_all, model.A_all, model.B_all = b_all, A_all, B_all
    model.reynolds = float(reynolds)
    model.pressure_ops = pressure_ops       # per-cluster PPE operators (pressure recovery)
    model.Ntrain = int(Xv.shape[1])
    torch.save(model, cache)
    print(f"Saved pinball qlGalerkin cache: K={K}, r={model.r}, file={cache}.")
    return model

free_run_divergence_free(model, x0, n_steps, on_step=None)

Like qlroms.free_run, but after every step projects the reduced velocity coordinate onto the active chart's own reduced-divergence constraint (PODGalerkinODEOperators. correct_divergence): predict (one ETDRK4 step of the closed quadratic ODE, model.step) / correct (minimal-norm projection onto D_k a + d_k = 0). Needs model.pressure_ops (a pinball qlGalerkin from build_local_model); a no-op correction on any chart with r_pressure == 0.

Verified on real data: kills the reduced-divergence drift a plain free_run accumulates over a long closed loop, but does NOT by itself fix instability from an unrelated unstable eigenvalue of A_eff -- it is a kinematic projection, not a stabilizer.

on_step(n, apod), if given, is called after every step with the corrected (r+1, 1) augmented coordinate (chart id in the last row) -- e.g. a pinball.forces.ForceHistory, which records lift and drag on the fly from it.

Returns:

Type Description
(X_rec, cluster_path)

(N, n_steps) recovered trajectory and the (n_steps,)

active-chart id per step, same convention as qlroms.free_run.
Source code in qlroms/intrusive_qlroms/pinball/rom.py 
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
def free_run_divergence_free(model: qlGalerkin, x0, n_steps: int, on_step=None):
    """Like qlroms.free_run, but after every step projects the reduced velocity coordinate
    onto the active chart's own reduced-divergence constraint (PODGalerkinODEOperators.
    correct_divergence): predict (one ETDRK4 step of the closed quadratic ODE, model.step)
    / correct (minimal-norm projection onto D_k a + d_k = 0). Needs model.pressure_ops (a
    pinball qlGalerkin from build_local_model); a no-op correction on any chart with
    r_pressure == 0.

    Verified on real data: kills the reduced-divergence drift a plain free_run accumulates
    over a long closed loop, but does NOT by itself fix instability from an unrelated
    unstable eigenvalue of A_eff -- it is a kinematic projection, not a stabilizer.

    on_step(n, apod), if given, is called after every step with the corrected (r+1, 1)
    augmented coordinate (chart id in the last row) -- e.g. a pinball.forces.ForceHistory,
    which records lift and drag on the fly from it.

    Returns:
        (X_rec, cluster_path): (N, n_steps) recovered trajectory and the (n_steps,)
        active-chart id per step, same convention as qlroms.free_run.
    """
    x0 = torch.as_tensor(x0, dtype=model.rdtype, device=model.device).reshape(-1, 1)
    apod = model.project_state(x0)
    X_rec = torch.zeros(model.N, n_steps, dtype=model.rdtype, device=model.device)
    cluster_path = np.zeros(n_steps, dtype=int)
    for n in range(n_steps):
        apod = model.step(apod)                                  # predict
        k = int(apod[-1, 0].item())
        a_corrected = model.pressure_ops[k].correct_divergence(   # correct
            apod[:-1, 0].detach().cpu().numpy())
        a_col = torch.as_tensor(a_corrected, dtype=model.rdtype, device=model.device)[:, None]
        apod = torch.cat([a_col, apod[-1:, :]], dim=0)
        X_rec[:, n] = model.recover_state(apod).reshape(-1)
        cluster_path[n] = k
        if on_step is not None:
            on_step(n, apod)
    return X_rec, cluster_path

qlroms.intrusive_qlroms.pinball.solver

Full-order incompressible Navier-Stokes solver (dolfinx fractional step).

Two entry points share the same stepping kernel: - run() -- the offline snapshot-generation loop (MPI-capable): warm restart, snapshot/force saving, diagnostics. Needs snapshot_dir. - step() -- one fractional step, the kernel pinball.fom.FOM drives through the dynamodels Model protocol (serial), together with the state accessors get_state / set_state / set_reynolds / num_dofs.

PinballFOMParameters dataclass

Runtime parameters controlling the FOM time integration and snapshot export.

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
@dataclass
class PinballFOMParameters:
    """Runtime parameters controlling the FOM time integration and snapshot export."""

    dt: float
    num_steps: int
    reynolds: float
    snapshot_interval: int = 1
    snapshot_start: int = 0
    store_velocity_star: bool = True
    # Warm restart from saved snapshot arrays (velocity_*.npy / pressure_*.npy);
    # None keeps the historical from-rest start. init_velocity alone collapses the AB2
    # convection history (u_n = u_n1 = the loaded field, see set_state's docstring) --
    # the historical approximation fom_re130's two restarts used. init_velocity_prev
    # closes that gap with the TRUE u_{n-1} (only recoverable when the campaign being
    # continued saved every solver step, snapshot_interval=1): with it set, restarting
    # reproduces the same trajectory an uninterrupted run would have produced from this
    # point on, not an approximation of it.
    init_velocity: Path | None = None
    init_pressure: Path | None = None
    init_velocity_prev: Path | None = None
    # |u| quick-look PNGs (previews/preview_<step>.png in snapshot_dir), this many spread evenly
    # over the run plus the last step; 0 disables. Serial runs only.
    n_previews: int = 10
    # Channel side "walls": "no_slip" (default) are physical rigid walls (zero velocity);
    # "far_field" reuses the inlet's u=(1,0) profile there instead, the large-domain
    # approximation the fom_re130-style production campaigns were actually generated with.
    # Obstacles are always no-slip (physical cylinders) regardless of this choice.
    wall_bc: Literal["no_slip", "far_field"] = "no_slip"

PinballFOM

Fractional step Navier-Stokes solver used as FOM for the bluff-body benchmark cases.

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
class PinballFOM:
    """Fractional step Navier-Stokes solver used as FOM for the bluff-body benchmark cases."""

    def __init__(
        self,
        mesh,
        facet_tags,
        boundary_labels: dict[str, int] | None,
        parameters: PinballFOMParameters,
        snapshot_dir: Path | str | None = None,
    ) -> None:
        self.mesh = mesh
        self.facet_tags = facet_tags
        # Fall back to default labels when metadata is missing from the mesh loader.
        self.boundary_labels = boundary_labels or BOUNDARY_LABELS
        self.obstacle_label_names = _sorted_obstacle_names(self.boundary_labels)
        self.params = parameters
        self.snapshot_dir = Path(snapshot_dir) if snapshot_dir is not None else None
        self.comm = mesh.comm
        self.time = 0.0
        self._setup_spaces()
        self._setup_boundary_data()
        self._setup_forms()
        self._setup_observables()
        self._step_history: list[int] = []
        self._time_history: list[float] = []
        self._drag_history: list[float] = []
        self._lift_history: list[float] = []
        self._drag_obstacle_history: list[np.ndarray] = []
        self._lift_obstacle_history: list[np.ndarray] = []
        self._history_first_step = 0
        self._saved_snapshot_steps: list[int] = []
        self._saved_snapshot_times: list[float] = []
        self._preview_tri = None

    # ------------------------------------------------------------------
    # Discretisation setup
    # ------------------------------------------------------------------
    def _setup_spaces(self) -> None:
        self.V, self.Q = taylor_hood_spaces(self.mesh)

    @staticmethod
    def _inlet_profile(x: np.ndarray) -> np.ndarray:
        # Uniform unit inflow, constant in time. If a time-dependent profile is ever
        # needed, make this depend on self.time and re-interpolate u_inlet in step().
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = 1.0
        return values

    def _setup_boundary_data(self) -> None:
        labels = self.boundary_labels
        ft = self.facet_tags
        if ft is None:
            raise ValueError("Facet tags are required to setup boundary conditions")
        self.u_inlet = Function(self.V)
        self.u_inlet.interpolate(self._inlet_profile)
        # Obstacles are always physical no-slip cylinders (dirichlet_bcs' default value is a
        # zero Function). The channel "walls" are either no-slip too, or far-field (reusing
        # the inlet profile) -- self.params.wall_bc selects it explicitly because it changes
        # the physics materially (see PinballFOMParameters.wall_bc).
        if self.params.wall_bc == "no_slip":
            wall_bcs = dirichlet_bcs(self.V, ft, [labels["walls"]])
        elif self.params.wall_bc == "far_field":
            wall_bcs = dirichlet_bcs(self.V, ft, [labels["walls"]], self.u_inlet)
        else:
            raise ValueError(f"Unsupported wall_bc {self.params.wall_bc!r}; "
                             "choose 'no_slip' or 'far_field'.")
        self.bcu = (dirichlet_bcs(self.V, ft, [labels["inlet"]], self.u_inlet)
                    + wall_bcs
                    + dirichlet_bcs(self.V, ft, [labels[k] for k in self.obstacle_label_names]))
        self.bcp = dirichlet_bcs(self.Q, ft, [labels["outlet"]])

    def _setup_forms(self) -> None:
        params = self.params

        self.k = Constant(self.mesh, PETSc.ScalarType(params.dt))
        self.mu = Constant(self.mesh, PETSc.ScalarType(1.0 / params.reynolds))
        self.rho = Constant(self.mesh, PETSc.ScalarType(1.0))

        u, v = TrialFunction(self.V), TestFunction(self.V)
        p, q = TrialFunction(self.Q), TestFunction(self.Q)

        self.u_ = Function(self.V)
        self.u_star = Function(self.V)
        self.u_n = Function(self.V)
        self.u_n1 = Function(self.V)
        self.p_ = Function(self.Q)
        self.phi = Function(self.Q)

        # Three-step projection scheme:
        # 1. tentative velocity  2. pressure correction  3. divergence-free update
        f = Constant(self.mesh, PETSc.ScalarType((0.0, 0.0)))
        u_n, u_n1 = self.u_n, self.u_n1

        F1 = self.rho / self.k * dot(u - u_n, v) * dx
        F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
        F1 += 0.5 * self.mu * inner(grad(u + u_n), grad(v)) * dx
        F1 += -dot(self.p_, div(v)) * dx
        F1 += dot(f, v) * dx

        self.a1 = fem.form(lhs(F1))
        self.L1 = fem.form(rhs(F1))
        self.A1 = create_matrix(self.a1)
        self.b1 = create_vector(self.L1)

        self.a2 = fem.form(dot(grad(p), grad(q)) * dx)
        self.L2 = fem.form(-self.rho / self.k * dot(div(self.u_star), q) * dx)
        self.A2 = assemble_matrix(self.a2, bcs=self.bcp)
        self.A2.assemble()
        self.b2 = create_vector(self.L2)

        self.a3 = fem.form(self.rho * dot(u, v) * dx)
        self.L3 = fem.form(self.rho * dot(self.u_star, v) * dx - self.k * dot(nabla_grad(self.phi), v) * dx)
        self.A3 = assemble_matrix(self.a3)
        self.A3.assemble()
        self.b3 = create_vector(self.L3)

        self.solver1 = self._make_ksp(self.A1, PETSc.KSP.Type.BCGS, PETSc.PC.Type.JACOBI)
        self.solver2 = self._make_ksp(self.A2, PETSc.KSP.Type.MINRES, PETSc.PC.Type.HYPRE, hypre="boomeramg")
        self.solver3 = self._make_ksp(self.A3, PETSc.KSP.Type.CG, PETSc.PC.Type.SOR)

    def _make_ksp(self, A, ksp_type, pc_type, hypre: str | None = None) -> PETSc.KSP:
        solver = PETSc.KSP().create(self.mesh.comm)
        solver.setOperators(A)
        solver.setType(ksp_type)
        pc = solver.getPC()
        pc.setType(pc_type)
        if hypre is not None:
            pc.setHYPREType(hypre)
        return solver

    def _setup_observables(self) -> None:
        # Build separate force integrals for each obstacle, then sum them during time stepping.
        # (forces.force_integrands is also what forces.assemble_force_functionals differentiates)
        labels = self.boundary_labels
        obs_measures = [
            Measure("ds", domain=self.mesh, subdomain_data=self.facet_tags, subdomain_id=labels[key])
            for key in self.obstacle_label_names
        ]
        drag, lift = force_integrands(self.mesh, self.u_, self.p_, self.mu, self.rho)
        self.drag_forms = [fem.form(drag * measure) for measure in obs_measures]
        self.lift_forms = [fem.form(lift * measure) for measure in obs_measures]

    # ------------------------------------------------------------------
    # State access (the pinball.fom.FOM adapter drives these)
    # ------------------------------------------------------------------
    @property
    def num_dofs(self) -> tuple[int, int]:
        """(Nu, Np): local velocity and pressure dof-array sizes."""
        return self.u_.x.array.size, self.p_.x.array.size

    def get_state(self) -> np.ndarray:
        """Stacked [velocity dofs; pressure dofs] state vector (a copy)."""
        return np.concatenate([self.u_.x.array, self.p_.x.array])

    def set_state(self, u_values, p_values) -> None:
        """Overwrite the solver state; the AB2 convection history collapses onto
        the given field (u_n = u_n1 = u), the same approximation as a warm restart."""
        for fn in (self.u_, self.u_n, self.u_n1):
            fn.x.array[:] = np.asarray(u_values, dtype=float).reshape(-1)
            fn.x.scatter_forward()
        self.p_.x.array[:] = np.asarray(p_values, dtype=float).reshape(-1)
        self.p_.x.scatter_forward()

    def set_reynolds(self, reynolds: float) -> None:
        """Update the viscosity Constant in the assembled forms (mu = 1/Re)."""
        self.params.reynolds = float(reynolds)
        self.mu.value = 1.0 / float(reynolds)

    # ------------------------------------------------------------------
    # Time stepping
    # ------------------------------------------------------------------
    def _solve_substep(self, solver, b, L, out, a_form=None, bcs=None) -> None:
        """Assemble the RHS of one projection substep and solve into `out`
        (with BC lifting when `a_form`/`bcs` are given)."""
        with b.localForm() as loc:
            loc.set(0.0)
        assemble_vector(b, L)
        if bcs:
            apply_lifting(b, [a_form], [bcs])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        if bcs:
            set_bc(b, bcs)
        solver.solve(b, out.x.petsc_vec)
        out.x.scatter_forward()

    def step(self) -> None:
        """Advance one fractional step (tentative velocity, pressure correction,
        divergence-free update) and shift the AB2 history."""
        self.time += self.params.dt

        # Step 1: tentative velocity (the convection matrix depends on u_n: reassemble).
        self.A1.zeroEntries()
        assemble_matrix(self.A1, self.a1, bcs=self.bcu)
        self.A1.assemble()
        self._solve_substep(self.solver1, self.b1, self.L1, self.u_star,
                            a_form=self.a1, bcs=self.bcu)

        # Step 2: Poisson problem for the pressure increment.
        self._solve_substep(self.solver2, self.b2, self.L2, self.phi,
                            a_form=self.a2, bcs=self.bcp)
        self.p_.x.petsc_vec.axpy(1.0, self.phi.x.petsc_vec)
        self.p_.x.scatter_forward()

        # Step 3: correct the velocity so the final field is divergence-free.
        self._solve_substep(self.solver3, self.b3, self.L3, self.u_)

        # Shift the AB2 history: u_n1 <- u_n, u_n <- u_.
        with self.u_.x.petsc_vec.localForm() as u_loc, \
                self.u_n.x.petsc_vec.localForm() as u_n_loc, \
                self.u_n1.x.petsc_vec.localForm() as u_n1_loc:
            u_n_loc.copy(u_n1_loc)
            u_loc.copy(u_n_loc)

    def run(self, overwrite: bool = False) -> None:
        """Offline snapshot-generation loop (needs `snapshot_dir`). Unless `overwrite`, refuses
        to start a from-rest run in a directory that already holds a campaign (even at other
        indices it would overwrite snapshot_metadata.npz and forces_history.png), and a warm
        restart whose snapshot indices are already on disk; a restart that only adds new
        indices goes ahead, first copying the campaign's snapshot_metadata.npz to
        snapshot_metadata_before_<first new index>.npz -- the `previous_metadata` that
        build_pinball.merge_restart_snapshot_metadata folds this leg back into."""
        if self.snapshot_dir is None:
            raise ValueError("PinballFOM.run() needs a snapshot_dir; the Model adapter "
                             "path (pinball.fom.FOM) does not.")
        params = self.params
        step_offset = self._restart_step_offset()
        meta_path = self.snapshot_dir / "snapshot_metadata.npz"
        if not overwrite:
            existing = {int(f.stem.rsplit("_", 1)[1]) for f in self.snapshot_dir.glob("velocity_*.npy")}
            if params.init_velocity is None and (existing or meta_path.exists()):
                raise FileExistsError(
                    f"{self.snapshot_dir} already holds a campaign ({len(existing)} velocity snapshots): a "
                    "from-rest run would overwrite its snapshot_metadata.npz -- use a new campaign tag "
                    "(`tag:` in the run_pinball_fom.py config), or overwrite=True (--overwrite).")
            planned = range(step_offset + params.snapshot_start, step_offset + params.num_steps,
                            params.snapshot_interval)
            clash = sorted(existing.intersection(planned))
            if clash:
                raise FileExistsError(
                    f"{len(clash)} of the {len(planned)} snapshots this run would save already exist in "
                    f"{self.snapshot_dir} (velocity_{clash[0]:05d} ... velocity_{clash[-1]:05d}): use a new "
                    "campaign tag (`tag:` in the run_pinball_fom.py config), or overwrite=True (--overwrite) to replace them.")
        if self.comm.rank == 0 and params.init_velocity is not None and meta_path.exists():
            backup = self.snapshot_dir / f"snapshot_metadata_before_{step_offset:05d}.npz"
            if not backup.exists():   # a retried leg keeps the copy taken before its first attempt
                shutil.copy2(meta_path, backup)

        if self.comm.rank == 0:
            self.snapshot_dir.mkdir(parents=True, exist_ok=True)
            for hist in (self._step_history, self._time_history, self._drag_history, self._lift_history,
                         self._drag_obstacle_history, self._lift_obstacle_history,
                         self._saved_snapshot_steps, self._saved_snapshot_times):
                hist.clear()
        self._history_first_step = step_offset

        progress = None
        if self.comm.rank == 0:
            try:
                import tqdm
                progress = tqdm.tqdm(total=params.num_steps, desc="FOM time stepping")
            except ModuleNotFoundError:
                pass

        # ponytail: serial restart only -- arrays are loaded whole, no MPI scatter.
        if params.init_velocity is not None:
            vals = np.load(params.init_velocity)
            for fn in (self.u_, self.u_n, self.u_n1):
                fn.x.array[:] = vals
                fn.x.scatter_forward()
        if params.init_pressure is not None:
            self.p_.x.array[:] = np.load(params.init_pressure)
            self.p_.x.scatter_forward()
        if params.init_velocity_prev is not None:
            # Overrides only u_n1 (u_/u_n above already hold the true last state): an exact
            # restart, not the u_n1 = u_n collapse -- see PinballFOMParameters' docstring.
            self.u_n1.x.array[:] = np.load(params.init_velocity_prev)
            self.u_n1.x.scatter_forward()

        self.time = step_offset * params.dt   # 0.0 for a from-rest run (step_offset == 0)
        preview_every = max(1, params.num_steps // params.n_previews) if params.n_previews else 0
        for i in range(params.num_steps):
            self.step()
            save_index = step_offset + i

            if preview_every and (i % preview_every == 0 or i == params.num_steps - 1):
                self._save_preview(save_index)

            # Save only the post-transient portion requested by the caller.
            if i >= params.snapshot_start and (i - params.snapshot_start) % params.snapshot_interval == 0:
                self._save_snapshot(save_index, params.store_velocity_star)

            self._update_forces(save_index, i, params.snapshot_start)
            if self.comm.rank == 0 and (i + 1) % FORCE_HISTORY_FLUSH == 0:
                self._save_force_history()

            if progress is not None:
                progress.update(1)

        if progress is not None:
            progress.close()

        if self.comm.rank == 0 and self._time_history:
            self._save_snapshot_metadata()
            self._save_force_history()
            self._save_force_plot()

    # ------------------------------------------------------------------
    # Helpers
    # ------------------------------------------------------------------
    def _restart_step_offset(self) -> int:
        """First new snapshot index, continuing after a warm-restart file's own
        index (velocity_03998.npy -> new snapshots start at 03999) so re-running
        a restart config never overwrites the snapshots it restarted from."""
        if self.params.init_velocity is None:
            return 0
        match = re.search(r"_(\d+)$", self.params.init_velocity.stem)
        if not match:
            raise ValueError(f"init_velocity {self.params.init_velocity} has no trailing "
                             "step index (expected e.g. velocity_03998.npy) to restart from.")
        return int(match.group(1)) + 1

    def _save_snapshot(self, index: int, store_star: bool) -> None:
        # Store reduced-order training data as raw arrays for fast reload by the ROM scripts.
        np.save(self.snapshot_dir / f"velocity_{index:05d}.npy", self.u_.x.petsc_vec.getArray())
        np.save(self.snapshot_dir / f"pressure_{index:05d}.npy", self.p_.x.petsc_vec.getArray())
        if store_star:
            np.save(self.snapshot_dir / f"velocitystar_{index:05d}.npy", self.u_star.x.petsc_vec.getArray())
        if self.comm.rank == 0:
            self._saved_snapshot_steps.append(index)
            self._saved_snapshot_times.append(self.time)

    def _reduced_forces(self, forms) -> tuple[float | None, np.ndarray | None]:
        # Assemble each obstacle's integral locally, then MPI reduces of their sum (the total,
        # summed in the same order as ever) and of the per-obstacle values.
        local = [fem.assemble_scalar(form) for form in forms]
        return (self.comm.reduce(sum(local), op=MPI.SUM, root=0),
                self.comm.reduce(np.asarray(local, dtype=float), op=MPI.SUM, root=0))

    def _update_forces(self, save_index: int, rel_step: int, warmup_steps: int) -> None:
        drag_total, drag_obstacles = self._reduced_forces(self.drag_forms)
        lift_total, lift_obstacles = self._reduced_forces(self.lift_forms)
        if self.comm.rank == 0:
            if rel_step >= warmup_steps:
                np.save(self.snapshot_dir / f"drag_{save_index:05d}.npy", drag_total)
                np.save(self.snapshot_dir / f"lift_{save_index:05d}.npy", lift_total)
            self._step_history.append(save_index)
            self._time_history.append(self.time)
            self._drag_history.append(drag_total)
            self._lift_history.append(lift_total)
            self._drag_obstacle_history.append(drag_obstacles)
            self._lift_obstacle_history.append(lift_obstacles)

    def _save_force_history(self) -> None:
        """Every step's force coefficients so far, the unsaved transient included (drag_/lift_*.npy
        only cover the saved snapshots): forces_history_<first step>.npz, one file per run() leg --
        a restart never overwrites an earlier leg's -- replaced atomically, so a crash leaves the
        last complete one. forces.load_force_history reads the legs back as one series."""
        if not self._step_history:
            return
        path = self.snapshot_dir / f"forces_history_{self._history_first_step:05d}.npz"
        partial = path.with_name("partial_" + path.name)
        np.savez(partial, step=np.asarray(self._step_history, dtype=int),
                 t=np.asarray(self._time_history, dtype=float),
                 Cd=np.asarray(self._drag_history, dtype=float), Cl=np.asarray(self._lift_history, dtype=float),
                 Cd_obstacles=np.stack(self._drag_obstacle_history, axis=1),
                 Cl_obstacles=np.stack(self._lift_obstacle_history, axis=1),
                 obstacles=np.asarray(self.obstacle_label_names), reynolds=float(self.params.reynolds))
        os.replace(partial, path)

    def _save_force_plot(self) -> None:
        # Produce a quick-look diagnostic directly from the FOM run.
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt

        times = np.asarray(self._time_history, dtype=float)
        fig, ax = plt.subplots(figsize=(10, 4.5))
        ax.plot(times, np.asarray(self._drag_history, dtype=float), label="$C_d$")
        ax.plot(times, np.asarray(self._lift_history, dtype=float), label="$C_l$")
        ax.set_xlabel("Time")
        ax.set_ylabel("Force coefficient")
        ax.set_title("Drag and lift history")
        ax.grid(True, alpha=0.3)
        ax.legend(loc="best")
        fig.tight_layout()
        fig.savefig(self.snapshot_dir / "forces_history.png", dpi=200)
        plt.close(fig)

    def _save_preview(self, step: int) -> None:
        # |u| quick-look for eyeballing the run while it goes.
        # ponytail: serial only -- the triangulation covers the local partition; gather to rank 0 if MPI previews are needed
        if self.comm.size != 1:
            return
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt

        if self._preview_tri is None:
            self._preview_tri = scalar_lagrange_triangulation(self.mesh, value_degree=2)
        tri, gather_idx = self._preview_tri
        speed = np.linalg.norm(self.u_.x.array.reshape(-1, 2), axis=1)[gather_idx]
        aspect = np.ptp(tri.y) / np.ptp(tri.x)
        fig, ax = plt.subplots(figsize=(8, 8 * aspect + 0.5), constrained_layout=True)
        im = ax.tricontourf(tri, speed, levels=20, cmap="viridis")
        fig.colorbar(im, ax=ax, shrink=0.9, label="|u|")
        ax.set(aspect="equal", xticks=[], yticks=[],
               title=f"step {step}, t = {self.time:.2f}, Re = {self.params.reynolds:g}")
        path = self.snapshot_dir / "previews" / f"preview_{step:05d}.png"
        path.parent.mkdir(exist_ok=True)
        fig.savefig(path, dpi=150)
        plt.close(fig)

    def _save_snapshot_metadata(self) -> None:
        # Metadata lets the ROM scripts reconstruct the exact saved timeline later.
        np.savez(
            self.snapshot_dir / "snapshot_metadata.npz",
            saved_steps=np.asarray(self._saved_snapshot_steps, dtype=int),
            saved_times=np.asarray(self._saved_snapshot_times, dtype=float),
            dt=float(self.params.dt),
            snapshot_interval=int(self.params.snapshot_interval),
            snapshot_start_step=int(self.params.snapshot_start),
            num_steps=int(self.params.num_steps),
            reynolds=float(self.params.reynolds),
            wall_bc=self.params.wall_bc,
        )

num_dofs property

(Nu, Np): local velocity and pressure dof-array sizes.

get_state()

Stacked [velocity dofs; pressure dofs] state vector (a copy).

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
240
241
242
def get_state(self) -> np.ndarray:
    """Stacked [velocity dofs; pressure dofs] state vector (a copy)."""
    return np.concatenate([self.u_.x.array, self.p_.x.array])

set_state(u_values, p_values)

Overwrite the solver state; the AB2 convection history collapses onto the given field (u_n = u_n1 = u), the same approximation as a warm restart.

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
244
245
246
247
248
249
250
251
def set_state(self, u_values, p_values) -> None:
    """Overwrite the solver state; the AB2 convection history collapses onto
    the given field (u_n = u_n1 = u), the same approximation as a warm restart."""
    for fn in (self.u_, self.u_n, self.u_n1):
        fn.x.array[:] = np.asarray(u_values, dtype=float).reshape(-1)
        fn.x.scatter_forward()
    self.p_.x.array[:] = np.asarray(p_values, dtype=float).reshape(-1)
    self.p_.x.scatter_forward()

set_reynolds(reynolds)

Update the viscosity Constant in the assembled forms (mu = 1/Re).

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
253
254
255
256
def set_reynolds(self, reynolds: float) -> None:
    """Update the viscosity Constant in the assembled forms (mu = 1/Re)."""
    self.params.reynolds = float(reynolds)
    self.mu.value = 1.0 / float(reynolds)

step()

Advance one fractional step (tentative velocity, pressure correction, divergence-free update) and shift the AB2 history.

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
def step(self) -> None:
    """Advance one fractional step (tentative velocity, pressure correction,
    divergence-free update) and shift the AB2 history."""
    self.time += self.params.dt

    # Step 1: tentative velocity (the convection matrix depends on u_n: reassemble).
    self.A1.zeroEntries()
    assemble_matrix(self.A1, self.a1, bcs=self.bcu)
    self.A1.assemble()
    self._solve_substep(self.solver1, self.b1, self.L1, self.u_star,
                        a_form=self.a1, bcs=self.bcu)

    # Step 2: Poisson problem for the pressure increment.
    self._solve_substep(self.solver2, self.b2, self.L2, self.phi,
                        a_form=self.a2, bcs=self.bcp)
    self.p_.x.petsc_vec.axpy(1.0, self.phi.x.petsc_vec)
    self.p_.x.scatter_forward()

    # Step 3: correct the velocity so the final field is divergence-free.
    self._solve_substep(self.solver3, self.b3, self.L3, self.u_)

    # Shift the AB2 history: u_n1 <- u_n, u_n <- u_.
    with self.u_.x.petsc_vec.localForm() as u_loc, \
            self.u_n.x.petsc_vec.localForm() as u_n_loc, \
            self.u_n1.x.petsc_vec.localForm() as u_n1_loc:
        u_n_loc.copy(u_n1_loc)
        u_loc.copy(u_n_loc)

run(overwrite=False)

Offline snapshot-generation loop (needs snapshot_dir). Unless overwrite, refuses to start a from-rest run in a directory that already holds a campaign (even at other indices it would overwrite snapshot_metadata.npz and forces_history.png), and a warm restart whose snapshot indices are already on disk; a restart that only adds new indices goes ahead, first copying the campaign's snapshot_metadata.npz to snapshot_metadata_before_.npz -- the previous_metadata that build_pinball.merge_restart_snapshot_metadata folds this leg back into.

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def run(self, overwrite: bool = False) -> None:
    """Offline snapshot-generation loop (needs `snapshot_dir`). Unless `overwrite`, refuses
    to start a from-rest run in a directory that already holds a campaign (even at other
    indices it would overwrite snapshot_metadata.npz and forces_history.png), and a warm
    restart whose snapshot indices are already on disk; a restart that only adds new
    indices goes ahead, first copying the campaign's snapshot_metadata.npz to
    snapshot_metadata_before_<first new index>.npz -- the `previous_metadata` that
    build_pinball.merge_restart_snapshot_metadata folds this leg back into."""
    if self.snapshot_dir is None:
        raise ValueError("PinballFOM.run() needs a snapshot_dir; the Model adapter "
                         "path (pinball.fom.FOM) does not.")
    params = self.params
    step_offset = self._restart_step_offset()
    meta_path = self.snapshot_dir / "snapshot_metadata.npz"
    if not overwrite:
        existing = {int(f.stem.rsplit("_", 1)[1]) for f in self.snapshot_dir.glob("velocity_*.npy")}
        if params.init_velocity is None and (existing or meta_path.exists()):
            raise FileExistsError(
                f"{self.snapshot_dir} already holds a campaign ({len(existing)} velocity snapshots): a "
                "from-rest run would overwrite its snapshot_metadata.npz -- use a new campaign tag "
                "(`tag:` in the run_pinball_fom.py config), or overwrite=True (--overwrite).")
        planned = range(step_offset + params.snapshot_start, step_offset + params.num_steps,
                        params.snapshot_interval)
        clash = sorted(existing.intersection(planned))
        if clash:
            raise FileExistsError(
                f"{len(clash)} of the {len(planned)} snapshots this run would save already exist in "
                f"{self.snapshot_dir} (velocity_{clash[0]:05d} ... velocity_{clash[-1]:05d}): use a new "
                "campaign tag (`tag:` in the run_pinball_fom.py config), or overwrite=True (--overwrite) to replace them.")
    if self.comm.rank == 0 and params.init_velocity is not None and meta_path.exists():
        backup = self.snapshot_dir / f"snapshot_metadata_before_{step_offset:05d}.npz"
        if not backup.exists():   # a retried leg keeps the copy taken before its first attempt
            shutil.copy2(meta_path, backup)

    if self.comm.rank == 0:
        self.snapshot_dir.mkdir(parents=True, exist_ok=True)
        for hist in (self._step_history, self._time_history, self._drag_history, self._lift_history,
                     self._drag_obstacle_history, self._lift_obstacle_history,
                     self._saved_snapshot_steps, self._saved_snapshot_times):
            hist.clear()
    self._history_first_step = step_offset

    progress = None
    if self.comm.rank == 0:
        try:
            import tqdm
            progress = tqdm.tqdm(total=params.num_steps, desc="FOM time stepping")
        except ModuleNotFoundError:
            pass

    # ponytail: serial restart only -- arrays are loaded whole, no MPI scatter.
    if params.init_velocity is not None:
        vals = np.load(params.init_velocity)
        for fn in (self.u_, self.u_n, self.u_n1):
            fn.x.array[:] = vals
            fn.x.scatter_forward()
    if params.init_pressure is not None:
        self.p_.x.array[:] = np.load(params.init_pressure)
        self.p_.x.scatter_forward()
    if params.init_velocity_prev is not None:
        # Overrides only u_n1 (u_/u_n above already hold the true last state): an exact
        # restart, not the u_n1 = u_n collapse -- see PinballFOMParameters' docstring.
        self.u_n1.x.array[:] = np.load(params.init_velocity_prev)
        self.u_n1.x.scatter_forward()

    self.time = step_offset * params.dt   # 0.0 for a from-rest run (step_offset == 0)
    preview_every = max(1, params.num_steps // params.n_previews) if params.n_previews else 0
    for i in range(params.num_steps):
        self.step()
        save_index = step_offset + i

        if preview_every and (i % preview_every == 0 or i == params.num_steps - 1):
            self._save_preview(save_index)

        # Save only the post-transient portion requested by the caller.
        if i >= params.snapshot_start and (i - params.snapshot_start) % params.snapshot_interval == 0:
            self._save_snapshot(save_index, params.store_velocity_star)

        self._update_forces(save_index, i, params.snapshot_start)
        if self.comm.rank == 0 and (i + 1) % FORCE_HISTORY_FLUSH == 0:
            self._save_force_history()

        if progress is not None:
            progress.update(1)

    if progress is not None:
        progress.close()

    if self.comm.rank == 0 and self._time_history:
        self._save_snapshot_metadata()
        self._save_force_history()
        self._save_force_plot()

qlroms.intrusive_qlroms.pinball.aux_fenics

Pinball dolfinx layer: mesh (boundary labels, loader, gmsh generator), Taylor-Hood spaces and the FEM helpers shared by the intrusive build (pinball.rom, which imports this lazily) and the solver (pinball.solver). Mesh files (pinball_mesh.xdmf / pinball_facet_tags.xdmf / pinball_boundaries.json + .h5) are DATA, default location qlroms.utils.paths.PINBALL_MESHES; generate once with python scripts/pinball/run_pinball_fom.py configs/pinball_fom.yml --mesh-only (needs gmsh) or copy an existing set.

scalar_lagrange_triangulation(mesh, value_degree=2)

(triangulation, gather_idx): a linear (P1) matplotlib.tri.Triangulation of the mesh, plus the index array that gathers -- from a field ordered like a degree-value_degree Lagrange space's own dofs (e.g. velocity_magnitude's output) -- the values at the triangulation's P1 vertices.

A P1 space's dofmap is unambiguous (3 dofs/cell = the cell's 3 vertices, no edge dofs to confuse them with), so building the triangulation from it and matching a higher-degree field's nodes to it by coordinate is more robust than assuming a Pk (k > 1) triangle's local dof order puts vertices first.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def scalar_lagrange_triangulation(mesh, value_degree=2):
    """(triangulation, gather_idx): a linear (P1) `matplotlib.tri.Triangulation` of
    the mesh, plus the index array that gathers -- from a field ordered like a
    degree-`value_degree` Lagrange space's own dofs (e.g. `velocity_magnitude`'s
    output) -- the values at the triangulation's P1 vertices.

    A P1 space's dofmap is unambiguous (3 dofs/cell = the cell's 3 vertices, no
    edge dofs to confuse them with), so building the triangulation from it and
    matching a higher-degree field's nodes to it by coordinate is more robust than
    assuming a Pk (k > 1) triangle's local dof order puts vertices first.
    """
    import matplotlib.tri as mtri
    from dolfinx.fem import functionspace
    from scipy.spatial import cKDTree

    V1 = functionspace(mesh, ("Lagrange", 1))
    coords1 = V1.tabulate_dof_coordinates()[:, :2]
    cells = np.asarray(V1.dofmap.list)
    if cells.shape[1] == 4:
        # the pinball mesh is recombined into quads (gmsh RecombineAll), in a dofmap
        # order that need not run around the perimeter -- sort each quad's corners
        # by angle about its centroid first, so the diagonal split (into 2 triangles
        # for tricontourf) does not self-intersect.
        pts = coords1[cells]                                            # (ncell, 4, 2)
        centroid = pts.mean(axis=1, keepdims=True)
        angle = np.arctan2(pts[..., 1] - centroid[..., 1], pts[..., 0] - centroid[..., 0])
        cells = np.take_along_axis(cells, np.argsort(angle, axis=1), axis=1)
        cells = np.concatenate([cells[:, [0, 1, 2]], cells[:, [0, 2, 3]]], axis=0)
    tri = mtri.Triangulation(coords1[:, 0], coords1[:, 1], cells)

    if value_degree == 1:
        return tri, np.arange(len(coords1))
    Vv = functionspace(mesh, ("Lagrange", value_degree))
    coords_v = Vv.tabulate_dof_coordinates()[:, :2]
    gather_idx = cKDTree(coords_v).query(coords1)[1]
    return tri, gather_idx

taylor_hood_spaces(mesh)

(V, Q): the P2 velocity / P1 pressure pair every pinball consumer builds on.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
77
78
79
80
81
82
83
def taylor_hood_spaces(mesh: Mesh):
    """(V, Q): the P2 velocity / P1 pressure pair every pinball consumer builds on."""
    cell = mesh.topology.cell_name()
    gdim = mesh.geometry.dim
    V = fem.functionspace(mesh, element("Lagrange", cell, 2, shape=(gdim,)))
    Q = fem.functionspace(mesh, element("Lagrange", cell, 1))
    return V, Q

load_pinball_mesh(mesh_dir=None)

Load the mesh and facet tags generated for the pinball benchmark.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def load_pinball_mesh(mesh_dir: Path | str | None = None) -> tuple[Mesh, MeshTags, dict[str, int]]:
    """Load the mesh and facet tags generated for the pinball benchmark."""
    mesh_dir = Path(mesh_dir) if mesh_dir is not None else PINBALL_MESHES

    mesh_path = mesh_dir / "pinball_mesh.xdmf"
    facet_path = mesh_dir / "pinball_facet_tags.xdmf"
    metadata_path = mesh_dir / "pinball_boundaries.json"

    if not mesh_path.exists():
        raise FileNotFoundError(
            f"{mesh_path} -- generate the mesh with "
            "`python scripts/pinball/run_pinball_fom.py configs/pinball_fom.yml --mesh-only` "
            "or copy an existing mesh set there.")

    comm = MPI.COMM_WORLD

    def _read_mesh_any(xdmf: io.XDMFFile) -> Mesh:
        # XDMF naming is not always stable across generators, so try the common candidates.
        candidates = [None, "mesh", "Mesh", "grid", "Grid"]
        last_error: RuntimeError | None = None
        for name in candidates:
            try:
                if name is None:
                    return xdmf.read_mesh()
                return xdmf.read_mesh(name=name)
            except RuntimeError as exc:
                last_error = exc
        raise last_error if last_error is not None else RuntimeError("Unable to read mesh from file")

    with io.XDMFFile(comm, mesh_path.as_posix(), "r") as xdmf:
        mesh = _read_mesh_any(xdmf)

    if not facet_path.exists():
        raise FileNotFoundError(facet_path)

    with io.XDMFFile(comm, facet_path.as_posix(), "r") as xdmf:
        try:
            # Prefer the explicit tag name emitted by the generator, but keep a fallback for older files.
            facet_tags = xdmf.read_meshtags(mesh, name="Facet markers")
        except RuntimeError:
            facet_tags = xdmf.read_meshtags(mesh)

    if facet_tags is None:
        raise RuntimeError("Facet tags could not be read from the mesh file")

    if metadata_path.exists():
        # Load only the boundary ids; the JSON may contain extra geometric metadata as well.
        metadata = json.loads(metadata_path.read_text())
        labels = {key: int(value) for key, value in metadata.items() if key in BOUNDARY_LABELS}
    else:
        labels = BOUNDARY_LABELS.copy()

    return mesh, facet_tags, labels

generate_pinball_mesh(output_dir=None, *, polynomial_order=2)

Generate the gmsh mesh with facet tags and metadata for the pinball setup.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
def generate_pinball_mesh(output_dir: Path | str | None = None, *, polynomial_order: int = 2) -> None:
    """Generate the gmsh mesh with facet tags and metadata for the pinball setup."""

    import gmsh
    try:
        from dolfinx.io.gmsh import model_to_mesh  # dolfinx >= 0.10
    except ModuleNotFoundError:
        from dolfinx.io.gmshio import model_to_mesh  # dolfinx 0.9

    if polynomial_order < 1:
        raise ValueError("polynomial_order must be >= 1")

    output_dir = Path(output_dir) if output_dir is not None else PINBALL_MESHES
    output_dir.mkdir(parents=True, exist_ok=True)

    mesh_path = output_dir / "pinball_mesh.xdmf"
    facet_path = output_dir / "pinball_facet_tags.xdmf"
    metadata_path = output_dir / "pinball_boundaries.json"

    gmsh.initialize()
    try:
        # Geometry in cylinder diameters: an L x H channel with three unit-diameter
        # disks in the equilateral pinball arrangement, `dist` center-to-center.
        r = 0.5
        diameter = 2.0 * r
        L, H = 80.0 * diameter, 60.0 * diameter
        dist = 1.5
        front = (20.0 * diameter, H / 2.0)
        x_off = dist * np.cos(np.radians(30.0))
        centers = [front,
                   (front[0] + x_off, front[1] + 0.5 * dist),
                   (front[0] + x_off, front[1] - 0.5 * dist)]

        comm = MPI.COMM_WORLD
        model_rank = 0
        gdim = 2

        if comm.rank == model_rank:
            # Build the fluid domain by subtracting the three disks from the channel.
            rectangle = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, L, H, tag=1)
            disks = [gmsh.model.occ.addDisk(cx, cy, 0.0, r, r) for cx, cy in centers]
            gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, disk) for disk in disks])
            gmsh.model.occ.synchronize()

            gmsh.model.addPhysicalGroup(gdim, [1], BOUNDARY_LABELS["fluid"])
            gmsh.model.setPhysicalName(gdim, BOUNDARY_LABELS["fluid"], "Fluid")

            # Classify each 1-D entity by its center of mass, then tag one physical
            # group per boundary name so later solvers can apply BCs by tag.
            targets = {"inlet": (0.0, H / 2.0), "outlet": (L, H / 2.0),
                       "obstacle1": centers[0], "obstacle2": centers[1], "obstacle3": centers[2]}
            groups: dict[str, list[int]] = {name: [] for name in ("walls", *targets)}
            for dim, tag in gmsh.model.getBoundary(gmsh.model.getEntities(dim=gdim), oriented=False):
                com = gmsh.model.occ.getCenterOfMass(dim, tag)
                name = next((n for n, (cx, cy) in targets.items() if np.allclose(com, [cx, cy, 0.0])),
                            "walls")
                groups[name].append(tag)
            for name, tags in groups.items():
                gmsh.model.addPhysicalGroup(1, tags, BOUNDARY_LABELS[name])
                gmsh.model.setPhysicalName(1, BOUNDARY_LABELS[name], name.capitalize())

            obstacles = [groups[f"obstacle{i}"] for i in (1, 2, 3)]
            res_min = r / 3.0
            if all(obstacles):
                # Refine the wake region and the obstacle neighborhood more aggressively
                # than the far field: distance threshold + two nested wake boxes.
                distance_field = gmsh.model.mesh.field.add("Distance")
                gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", [obs[0] for obs in obstacles])

                def wake_box(size: float, pad: float) -> int:
                    box = gmsh.model.mesh.field.add("Box")
                    for key, value in (("VIn", size), ("XMin", front[0] - pad * r), ("XMax", L),
                                       ("YMin", centers[2][1] - pad * r), ("YMax", centers[1][1] + pad * r),
                                       ("Thickness", r)):
                        gmsh.model.mesh.field.setNumber(box, key, value)
                    return box

                boxes = [wake_box(1.0 * res_min, 2.0), wake_box(5.0 * res_min, 5.0)]

                threshold_field = gmsh.model.mesh.field.add("Threshold")
                for key, value in (("IField", distance_field), ("LcMin", res_min), ("LcMax", 0.25 * H),
                                   ("DistMin", r), ("DistMax", 2.0 * H)):
                    gmsh.model.mesh.field.setNumber(threshold_field, key, value)

                min_field = gmsh.model.mesh.field.add("Min")
                gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field, *boxes])
                gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

            gmsh.option.setNumber("Mesh.Algorithm", 8)
            gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
            gmsh.option.setNumber("Mesh.RecombineAll", 1)
            gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
            # Generate a quadratic mesh when requested so the FE spaces can use curved geometry consistently.
            gmsh.model.mesh.generate(gdim)
            gmsh.model.mesh.setOrder(polynomial_order)
            gmsh.model.mesh.optimize("Netgen")

        # Support both older tuple returns and newer MeshData returns from dolfinx.
        mesh_data = model_to_mesh(gmsh.model, comm, model_rank, gdim=gdim)
        if hasattr(mesh_data, "mesh"):
            msh = mesh_data.mesh
            cell_tags = getattr(mesh_data, "cell_tags", None)
            facet_tags = getattr(mesh_data, "facet_tags", None)
        else:
            msh, cell_tags, facet_tags = mesh_data
        msh.name = "mesh"
        if cell_tags is not None:
            cell_tags.name = "cell_tags"
        if facet_tags is not None:
            facet_tags.name = "Facet markers"

        with io.XDMFFile(comm, mesh_path.as_posix(), "w") as xdmf:
            xdmf.write_mesh(msh)
            if cell_tags is not None:
                xdmf.write_meshtags(cell_tags, msh.geometry)
            if facet_tags is not None:
                xdmf.write_meshtags(facet_tags, msh.geometry)

        if facet_tags is not None:
            with io.XDMFFile(comm, facet_path.as_posix(), "w") as xdmf:
                xdmf.write_mesh(msh)
                xdmf.write_meshtags(facet_tags, msh.geometry)

        if comm.rank == 0:
            # Persist boundary tag ids so later scripts can recover them without hard-coding.
            metadata = {name: int(tag) for name, tag in BOUNDARY_LABELS.items()}
            metadata.update({"diameter": diameter, "length": L, "height": H})
            metadata_path.write_text(json.dumps(metadata, indent=2))
    finally:
        gmsh.finalize()

function_from_array(space, values)

Function on space whose dof array is the 1-D values (scattered forward).

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
279
280
281
282
283
284
def function_from_array(space, values) -> Function:
    """Function on `space` whose dof array is the 1-D `values` (scattered forward)."""
    fn = Function(space)
    fn.x.array[:] = np.asarray(values, dtype=float).reshape(-1)
    fn.x.scatter_forward()
    return fn

assemble_mass_matrix(V)

Assembled FEM mass matrix (u, v) of the space V.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
287
288
289
290
291
292
def assemble_mass_matrix(V) -> PETSc.Mat:
    """Assembled FEM mass matrix (u, v) of the space V."""
    u, v = TrialFunction(V), TestFunction(V)
    M = assemble_matrix(fem.form(inner(u, v) * dx))
    M.assemble()
    return M

assemble_stiffness_matrix(V)

Assembled FEM stiffness matrix (grad u, grad v) of the space V -- the viscous block: the ROM's linear operator carries -nu * Phi^T K Phi (see pinball.rom).

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
295
296
297
298
299
300
301
def assemble_stiffness_matrix(V) -> PETSc.Mat:
    """Assembled FEM stiffness matrix (grad u, grad v) of the space V -- the viscous
    block: the ROM's linear operator carries -nu * Phi^T K Phi (see pinball.rom)."""
    u, v = TrialFunction(V), TestFunction(V)
    K = assemble_matrix(fem.form(inner(grad(u), grad(v)) * dx))
    K.assemble()
    return K

export_mass_matrices(M_u, M_p, mesh_dir=None)

Write the assembled velocity / pressure mass matrices next to the mesh (pinball_mass_velocity.npz, pinball_mass_pressure.npz, CSR + a weights.json manifest) so everything that needs the inner product -- data-driven models, diagnostics -- can read it without dolfinx (qlroms.intrusive_qlroms.build_pinball.load_mass_matrix). A file that already holds the same matrix (same sha1) is left alone; a read-only directory only warns. Returns {"velocity": sha1, "pressure": sha1}.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
def export_mass_matrices(M_u, M_p, mesh_dir: Path | str | None = None) -> dict[str, str]:
    """Write the assembled velocity / pressure mass matrices next to the mesh
    (pinball_mass_velocity.npz, pinball_mass_pressure.npz, CSR + a weights.json manifest) so
    everything that needs the inner product -- data-driven models, diagnostics -- can read it
    without dolfinx (qlroms.intrusive_qlroms.build_pinball.load_mass_matrix). A file that
    already holds the same matrix (same sha1) is left alone; a read-only directory only warns.
    Returns {"velocity": sha1, "pressure": sha1}."""
    mesh_dir = Path(mesh_dir) if mesh_dir is not None else PINBALL_MESHES
    out = {}
    for kind, M in (("velocity", M_u), ("pressure", M_p)):
        path = mesh_dir / f"pinball_mass_{kind}.npz"
        sha = weight_sha1(M)
        out[kind] = sha
        if _find_sidecar(sha, [mesh_dir]) == str(path.resolve()):
            continue
        try:
            save_weight(M, path)
        except OSError as exc:
            print(f"  [rom] could not write {path} ({exc}); models saved with this mass matrix "
                  "will need a sidecar file (qlroms.metric.save_weight).")
    return out

dirichlet_bcs(space, facet_tags, markers, value=None)

Dirichlet BC with value (a Function; default zero) on the facets tagged with ANY of the integer markers; [] when no dof is tagged (absent tags are fine).

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
327
328
329
330
331
332
333
334
335
336
337
def dirichlet_bcs(space, facet_tags, markers, value=None) -> list:
    """Dirichlet BC with `value` (a Function; default zero) on the facets tagged with
    ANY of the integer `markers`; [] when no dof is tagged (absent tags are fine)."""
    top = space.mesh.topology
    fdim = top.dim - 1
    top.create_connectivity(fdim, top.dim)
    top.create_connectivity(top.dim, fdim)
    facets = (np.unique(np.concatenate([facet_tags.find(m) for m in markers])) if markers
              else np.array([], dtype=np.int32))
    dofs = locate_dofs_topological(space, fdim, facets)
    return [dirichletbc(value if value is not None else Function(space), dofs)] if dofs.size else []

solve_linear_problems(a_form, rhs_forms, bcs)

CG+AMG solve of the UFL bilinear form a_form against each UFL linear form in rhs_forms, BCs applied; one Function per RHS.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
343
344
345
346
347
def solve_linear_problems(a_form, rhs_forms, bcs) -> list[Function]:
    """CG+AMG solve of the UFL bilinear form a_form against each UFL linear form in
    rhs_forms, BCs applied; one Function per RHS."""
    # ponytail: re-assembles a_form per RHS (<= r_pressure + 1 solves); hoist the KSP if it ever shows in a profile
    return [LinearProblem(a_form, L, bcs=bcs, petsc_options=_CG_AMG).solve() for L in rhs_forms]

compute_supremizers(V, pressure_modes, bcs)

Supremizer enrichment: for each pressure mode psi solve (grad w, grad v) = -(psi, div v) with homogeneous velocity BCs -- stabilises the reduced velocity-pressure coupling (inf-sup).

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
350
351
352
353
354
355
def compute_supremizers(V, pressure_modes, bcs) -> list[Function]:
    """Supremizer enrichment: for each pressure mode psi solve (grad w, grad v) = -(psi, div v)
    with homogeneous velocity BCs -- stabilises the reduced velocity-pressure coupling (inf-sup)."""
    u, v = TrialFunction(V), TestFunction(V)
    return solve_linear_problems(inner(grad(u), grad(v)) * dx,
                                 [-psi * div(v) * dx for psi in pressure_modes], bcs)

solve_pressure_centroid_from_velocity(Q, velocity_center, bcs)

Pressure centroid from the velocity centroid via the Poisson solve (grad p, grad q) = -(u_c.grad u_c, grad q) with outlet BCs -- the "poisson_from_velocity_centroid" alternative to the empirical pressure mean.

Source code in qlroms/intrusive_qlroms/pinball/aux_fenics.py 
358
359
360
361
362
363
364
365
def solve_pressure_centroid_from_velocity(Q, velocity_center, bcs) -> Function:
    """Pressure centroid from the velocity centroid via the Poisson solve
    (grad p, grad q) = -(u_c.grad u_c, grad q) with outlet BCs -- the
    "poisson_from_velocity_centroid" alternative to the empirical pressure mean."""
    p, q = TrialFunction(Q), TestFunction(Q)
    return solve_linear_problems(
        inner(grad(p), grad(q)) * dx,
        [-inner(dot(velocity_center, nabla_grad(velocity_center)), grad(q)) * dx], bcs)[0]