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.
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.
rom.assemble_galerkin_operators(...) The UFL assembly of one chart's operators.
solver.PinballFOM, solver.PinballFOMParameters The fractional-step solver: run() for snapshot campaigns (MPI-capable), step() / get_state() / set_state() for the FOM adapter.
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.

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 
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
@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

    @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

    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

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 
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
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 same Model protocol the fitted qlROMs get through qlroms.model.QLModel, so romda twin experiments can pair them: truth = FOM, forecast = QLModel, observing the SAME velocity-dof indices (FOM reads psi[obs_idx], QLModel reads H z + c).

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, matching QLModel's H = Phi[obs_idx] rows).

 (0,)
dt_solver float

Time step of the fractional-step scheme itself.

 0.01
Source code in qlroms/intrusive_qlroms/pinball/fom.py 
 20
 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
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, matching QLModel's
        H = Phi[obs_idx] rows).
    dt_solver : float
        Time step of the fractional-step scheme itself.
    """

    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 needs 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']

    reynolds = 60.0
    dt_solver = 0.01

    def __init__(self, psi0=None, dt=None, reynolds=60.0, mesh_dir=None,
                 obs_idx=(0,), dt_solver=0.01, **model_kwargs):
        self.reynolds = float(reynolds)
        self.mesh_dir = mesh_dir
        self.dt_solver = float(dt_solver)
        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)
        return PinballFOM(mesh, facet_tags, labels, parameters, snapshot_dir=None)

    # ---- discrete-map timing (mirrors QLModel) ----
    @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 DA ensembles on the fitted "
                "qlROM (qlroms.model.QLModel) 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 
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
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 
145
146
147
148
149
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 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
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):

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 
55
56
57
58
59
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='empirical')

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 
 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
def assemble_galerkin_operators(
    Umodes,
    Pmodes,
    u_mean,
    pressure_center,
    nu: float,
    pressure_centroid_mode: str = "empirical",
) -> 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)

    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,
    )

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

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 intrusive_model_K{K}_rv{r_velocity}_rp{r_pressure}_Re{reynolds:g}.pth 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

"empirical" (cluster pressure mean) or "poisson_from_velocity_centroid".

 'empirical'
random_state, cluster_space

forwarded to the clustering.

 required

Returns:

Type Description
qlGalerkin

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

.pressure_ops[k] the per-chart PODGalerkinODEOperators (pressure recovery
qlGalerkin

diagnostics), .reynolds, .Ntrain.
Source code in qlroms/intrusive_qlroms/pinball/rom.py 
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
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 = "empirical",
    random_state: int = 1,
    cluster_space: str = "pod_lossless",
) -> 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 intrusive_model_K{K}_rv{r_velocity}_rp{r_pressure}_Re{reynolds:g}.pth
            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: "empirical" (cluster pressure mean) or
            "poisson_from_velocity_centroid".
        random_state, cluster_space: forwarded to the clustering.

    Returns:
        qlGalerkin compilation of K ROMs (Mw = velocity mass matrix);
        `.pressure_ops[k]` the per-chart PODGalerkinODEOperators (pressure recovery
        diagnostics), `.reynolds`, `.Ntrain`.
    """
    save_dir = Path(save_dir) if save_dir is not None else PINBALL_MODELS
    save_dir.mkdir(parents=True, exist_ok=True)
    cache = save_dir / f"intrusive_model_K{K}_rv{r_velocity}_rp{r_pressure}_Re{reynolds:g}.pth"
    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)
    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"]])

    if K == 1:
        cluster_labels = torch.zeros(Xv.shape[1], dtype=torch.long)
    else:
        feats = clustering_features(Xv, r=r_velocity, cluster_space=cluster_space, random_state=random_state)
        cluster_labels = fit_clusters(feats, K, random_state=random_state)[1]
    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)

        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_mean_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.")

        p_mean_fn = fx.function_from_array(Q, p_mean_col.numpy())
        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_mean_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."
            )

        u_mean_fn = fx.function_from_array(V, u_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}")

        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)
    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

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 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
@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: Path | None = None
    init_pressure: Path | None = None

PinballFOM

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

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
 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
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
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._time_history: list[float] = []
        self._drag_history: list[float] = []
        self._lift_history: list[float] = []
        self._saved_snapshot_steps: list[int] = []
        self._saved_snapshot_times: list[float] = []

    # ------------------------------------------------------------------
    # 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)
        # Walls and obstacles are no-slip for the cropped domain (dirichlet_bcs' default
        # value is a zero Function); the big-domain runs reused the inlet profile on the
        # walls -- keep this explicit because it changes the physics materially.
        self.bcu = (dirichlet_bcs(self.V, ft, [labels["inlet"]], self.u_inlet)
                    + dirichlet_bcs(self.V, ft, [labels["walls"],
                                                 *[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.
        labels = self.boundary_labels
        n = -FacetNormal(self.mesh)
        obs_measures = [
            Measure("ds", domain=self.mesh, subdomain_data=self.facet_tags, subdomain_id=labels[key])
            for key in self.obstacle_label_names
        ]
        u_t = inner(as_vector((n[1], -n[0])), self.u_)
        mu, rho = self.mu, self.rho
        self.drag_forms = [fem.form(2.0 * (mu / rho * inner(grad(u_t), n) * n[1] - self.p_ * n[0]) * measure)
                           for measure in obs_measures]
        self.lift_forms = [fem.form(-2.0 * (mu / rho * inner(grad(u_t), n) * n[0] + self.p_ * n[1]) * 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) -> None:
        """Offline snapshot-generation loop (needs `snapshot_dir`)."""
        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

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

        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()

        self.time = 0.0
        for i in range(params.num_steps):
            self.step()

            # 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(i, params.store_velocity_star)

            self._update_forces(i, params.snapshot_start)

            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_plot()

    # ------------------------------------------------------------------
    # Helpers
    # ------------------------------------------------------------------
    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_sum(self, forms) -> float | None:
        # Assemble each obstacle's integral locally, then one MPI reduce of the sum.
        local = sum(fem.assemble_scalar(form) for form in forms)
        return self.comm.reduce(local, op=MPI.SUM, root=0)

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

    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_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),
        )

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 
209
210
211
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 
213
214
215
216
217
218
219
220
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 
222
223
224
225
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 
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
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()

Offline snapshot-generation loop (needs snapshot_dir).

Source code in qlroms/intrusive_qlroms/pinball/solver.py 
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
def run(self) -> None:
    """Offline snapshot-generation loop (needs `snapshot_dir`)."""
    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

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

    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()

    self.time = 0.0
    for i in range(params.num_steps):
        self.step()

        # 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(i, params.store_velocity_star)

        self._update_forces(i, params.snapshot_start)

        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_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.

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 
38
39
40
41
42
43
44
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 
 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
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 
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
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 
240
241
242
243
244
245
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 
248
249
250
251
252
253
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

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 
256
257
258
259
260
261
262
263
264
265
266
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 
272
273
274
275
276
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 
279
280
281
282
283
284
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 
287
288
289
290
291
292
293
294
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]