qlroms.intrusive_qlroms — qlGalerkin and the data pipelines

The intrusive family: qlGalerkin is the compilation of Galerkin-projected members, and the three cases define the members it compiles. Every case package has the same shape — config.py (TEST_CASES + the torch-side grids/operators), fom.py (the full-order model, a dynamodels Model), rom.py (the single-cluster ROM class and build_local_model, which caches itself). The offline data pipelines sit next to them: build_ks (KS FOM + cached trajectory) and build_pinball (snapshot loading).

The cases have their own pages: KS 1-D · KS 2-D · fluidic pinball.

At a glance

Name	One-liner
galerkin.qlGalerkin	Compilation of intrusively-projected members, one GalerkinROM per chart (isinstance tells the family apart); from_operators reloads stacked (b, A, B).
build_ks	KS offline data pipeline: build_fom (a case's FOM + merged run config), get_full_trajectory (the physical trajectory, generated in chunks and cached), build_sensor_indices (sensor placement on the KS grid).
build_pinball	load_snapshots: stack velocity_*.npy / pressure_*.npy into an (Ndof, Nt) matrix plus the snapshot spacing.

qlroms.intrusive_qlroms.galerkin

qlGalerkin: the intrusive model family.

A qlGalerkin is simply the compilation (qlroms.base.qlROM) of single-cluster ROMs whose operators come from an intrusive Galerkin projection of a known governing-equation operator -- the case ROM classes (ks1d.rom.ROM, ks2d.rom.ROM, pinball.rom.ROM) own those operators and step_reduced; compiling K of them yields the quantized-local model. The name records the family, so isinstance(model, qlGalerkin) identifies intrusive compilations in family-agnostic code.

GalerkinROM

Bases: QuadROM

One cluster's intrusive quadratic ROM: the family-neutral QuadROM member (qlroms.charts.QuadROM -- geometry + (b, A, B) + ETDRK4 step), labeled as INTRUSIVE: its operators come from a Galerkin projection of a known governing-equation operator, never from a regression.

Source code in qlroms/intrusive_qlroms/galerkin.py

class GalerkinROM(QuadROM):
    """One cluster's intrusive quadratic ROM: the family-neutral QuadROM member
    (qlroms.charts.QuadROM -- geometry + (b, A, B) + ETDRK4 step), labeled as INTRUSIVE:
    its operators come from a Galerkin projection of a known governing-equation
    operator, never from a regression."""

qlGalerkin

Bases: qlROM

Compilation of intrusively-projected (Galerkin) single-cluster ROMs.

Source code in qlroms/intrusive_qlroms/galerkin.py

class qlGalerkin(qlROM):
    """Compilation of intrusively-projected (Galerkin) single-cluster ROMs."""

    @classmethod
    def from_operators(cls, centroids, Phi_all, b_all, A_all, B_all, dt: float = 1.0,
                       Mw=None, Mcontour: int = 32, Rcontour: float = 15.0) -> "qlGalerkin":
        """Build a qlGalerkin from already-projected reduced operators, skipping any
        fit -- e.g. to reload stacked (b, A, B) from disk (tests/test_qlgalerkin_bridge.py);
        the case builders compile their own ROM classes directly.

        The reduced continuous quadratic ODE  da/dt = b + A a + B(a, a)  is the
        family-neutral qlroms.charts.QuadROM compute; GalerkinROM records the
        intrusive provenance.
        centroids (K, N), Phi_all (N, r, K), b/A/B stacked (K, ...); Mw is the chart
        inner-product weight (e.g. the FEM mass matrix; None for M = I).
        """
        K = int(centroids.shape[0])
        roms = []
        for k in range(K):
            rom = GalerkinROM(Phi_all[:, :, k], centroids[k], dt=dt, Mw=Mw,
                              Mcontour=Mcontour, Rcontour=Rcontour)
            rom.set_operators(b_all[k], A_all[k], B_all[k])
            roms.append(rom)
        model = cls(roms, dt=dt)
        model.b_all, model.A_all, model.B_all = b_all, A_all, B_all
        return model

    def _post_step(self, a, ids):
        """Predict/correct: project each column of a (r, m) onto its chart's reduced-divergence
        constraint D_k a + d_k = 0 (pinball PODGalerkinODEOperators.correct_divergence), after
        the ETDRK4 step and the transition. Only pinball models carry .pressure_ops; a no-op
        for every other qlGalerkin and on charts with r_pressure == 0."""
        ops = getattr(self, "pressure_ops", None)
        if ops is None:
            return a
        a = a.clone()
        for k in ids.unique().tolist():
            Pq = ops[int(k)]._divergence_projection
            if Pq is not None:
                P, q = (torch.as_tensor(v, dtype=a.dtype, device=a.device) for v in Pq)
                m = ids == k
                a[:, m] = P @ a[:, m] + q[:, None]
        return a

from_operators(centroids, Phi_all, b_all, A_all, B_all, dt=1.0, Mw=None, Mcontour=32, Rcontour=15.0) classmethod

Build a qlGalerkin from already-projected reduced operators, skipping any fit -- e.g. to reload stacked (b, A, B) from disk (tests/test_qlgalerkin_bridge.py); the case builders compile their own ROM classes directly.

The reduced continuous quadratic ODE da/dt = b + A a + B(a, a) is the family-neutral qlroms.charts.QuadROM compute; GalerkinROM records the intrusive provenance. centroids (K, N), Phi_all (N, r, K), b/A/B stacked (K, ...); Mw is the chart inner-product weight (e.g. the FEM mass matrix; None for M = I).

Source code in qlroms/intrusive_qlroms/galerkin.py

@classmethod
def from_operators(cls, centroids, Phi_all, b_all, A_all, B_all, dt: float = 1.0,
                   Mw=None, Mcontour: int = 32, Rcontour: float = 15.0) -> "qlGalerkin":
    """Build a qlGalerkin from already-projected reduced operators, skipping any
    fit -- e.g. to reload stacked (b, A, B) from disk (tests/test_qlgalerkin_bridge.py);
    the case builders compile their own ROM classes directly.

    The reduced continuous quadratic ODE  da/dt = b + A a + B(a, a)  is the
    family-neutral qlroms.charts.QuadROM compute; GalerkinROM records the
    intrusive provenance.
    centroids (K, N), Phi_all (N, r, K), b/A/B stacked (K, ...); Mw is the chart
    inner-product weight (e.g. the FEM mass matrix; None for M = I).
    """
    K = int(centroids.shape[0])
    roms = []
    for k in range(K):
        rom = GalerkinROM(Phi_all[:, :, k], centroids[k], dt=dt, Mw=Mw,
                          Mcontour=Mcontour, Rcontour=Rcontour)
        rom.set_operators(b_all[k], A_all[k], B_all[k])
        roms.append(rom)
    model = cls(roms, dt=dt)
    model.b_all, model.A_all, model.B_all = b_all, A_all, B_all
    return model

qlroms.intrusive_qlroms.build_ks

Offline KS build pipeline (ks1d / ks2d): the case FOM (build_fom) and its cached full physical trajectory (get_full_trajectory), which feed the case rom.build_local_model -- see qlroms.utils.builder.build_case for the config-driven orchestration. Also the sensor-placement helper build_sensor_indices.

The full-order models are the cases' own FOM classes (module.FOM: a dynamodels.physical.KS / KS2D subclass); the ROM machinery stays torch. build_fom constructs module.FOM for a case, get_full_trajectory drives model.time_integrate and caches the physical snapshot matrix in the exact same file format/layout as the old torch FOMs, so existing full_trajectory.pth caches keep hitting.

state_dim(model)

Flattened PHYSICAL state length of a KS FOM or config (1D: Nx, 2D: Nx*Ny).

Do not use model.N for dynamodels models: dynamodels.Model.N is the analysis-augmented size (Nphi + Na + Nq), not the physical grid size.

Source code in qlroms/intrusive_qlroms/build_ks.py

def state_dim(model) -> int:
    """Flattened PHYSICAL state length of a KS FOM or config (1D: Nx, 2D: Nx*Ny).

    Do not use model.N for dynamodels models: dynamodels.Model.N is the
    analysis-augmented size (Nphi + Na + Nq), not the physical grid size.
    """
    if hasattr(model, "Ny"):
        return int(model.Nx) * int(model.Ny)
    return int(model.Nx)

build_fom(module, case=None, overrides=None)

Build the case's full-order model and return (fom, cfg).

The FOM is module.FOM(case=...) -- a dynamodels Model subclass (ks1d.FOM(physical.KS) / ks2d.FOM(physical.KS2D)) that reads the case's PHYSICAL parameters straight out of its own torch-side config, which it keeps as fom.cfg. There is no rescaling and nothing is attached after construction: fom.dt IS the case's physical time step, and the case-side quantities the ROM stack reads (visc/Lx/lamb1/device/rdtype/cdtype, 2-D wt) are properties on the class delegating to fom.cfg.

Per-case model params (Nx/Ny/dt/visc/nu/lamb1) live in module.TEST_CASES[case]. Dataset params (Ntrain/Ntest/i0) come from TEST_CASES[case] (1D) or module.TIME_DEFAULTS (2D). overrides (CLI flags; None values ignored) win over both and, because the torch ROM configs read their params from TEST_CASES, are written there before construction.

cfg is the merged dict the scripts run on: TIME_DEFAULTS < TEST_CASES[case] < overrides, plus case/module.

Source code in qlroms/intrusive_qlroms/build_ks.py

def build_fom(module, case=None, overrides=None):
    """Build the case's full-order model and return (fom, cfg).

    The FOM is module.FOM(case=...) -- a dynamodels Model subclass
    (ks1d.FOM(physical.KS) / ks2d.FOM(physical.KS2D)) that reads the case's
    PHYSICAL parameters straight out of its own torch-side config, which it
    keeps as `fom.cfg`. There is no rescaling and nothing is attached after
    construction: `fom.dt` IS the case's physical time step, and the case-side
    quantities the ROM stack reads (visc/Lx/lamb1/device/rdtype/cdtype, 2-D
    `wt`) are properties on the class delegating to `fom.cfg`.

    Per-case model params (Nx/Ny/dt/visc/nu/lamb1) live in module.TEST_CASES[case].
    Dataset params (Ntrain/Ntest/i0) come from TEST_CASES[case] (1D) or
    module.TIME_DEFAULTS (2D). `overrides` (CLI flags; None values ignored)
    win over both and, because the torch ROM configs read their params from
    TEST_CASES, are written there before construction.

    cfg is the merged dict the scripts run on:
        TIME_DEFAULTS < TEST_CASES[case] < overrides, plus case/module.
    """
    TEST_CASES = module.TEST_CASES
    case = case or next(iter(TEST_CASES))
    if case not in TEST_CASES:
        raise ValueError(f"Unknown case '{case}'. Valid: {list(TEST_CASES)}.")

    ov = {k: v for k, v in (overrides or {}).items() if v is not None}
    if "Nx" in ov and "Ny" in TEST_CASES[case]:     # 2D: keep Ny in sync unless given
        ov.setdefault("Ny", ov["Nx"])

    TEST_CASES[case].update(ov)  # ponytail: mutate module global; one FOM per process run, no leak

    cfg = {**getattr(module, "TIME_DEFAULTS", {}), **TEST_CASES[case],
           "case": case, "module": module}
    fom = module.FOM(case=case)
    return fom, cfg

get_full_trajectory(Ntot, model, i0=0)

Compute or load the full physical trajectory of a KS FOM.

Works for both 1D and 2D KS: model is the CONSTRUCTED dynamodels model returned by build_fom (physical.KS or physical.KS2D). Generation runs model.time_integrate in chunks; ks1d spectral history is converted to physical snapshots via the model's own transforms. The cache file format and layout are identical to the old torch-FOM pipeline -- a torch-saved float64 tensor of POST-burn-in physical snapshots, (state_dim, Ncols), at get_simulation_path(model)/full_trajectory.pth -- so existing caches hit. Cache-hit rule unchanged: cached_cols >= Ntot -> last Ntot columns.

Parameters:

Name	Type	Description	Default
Ntot	int	Number of kept snapshots to return after discarding i0 transient steps.	required
model		dynamodels FOM (from build_fom).	required
i0	int	Number of initial transient steps to skip (default 0). Only used when generating from scratch; the cache stores post-burn-in states.	0

Returns:

Type	Description
Tensor	Trajectory tensor with shape (state_dim, Ntot), float64, CPU.

Source code in qlroms/intrusive_qlroms/build_ks.py

def get_full_trajectory(
    Ntot: int,
    model,
    i0: int = 0,
) -> torch.Tensor:
    """Compute or load the full physical trajectory of a KS FOM.

    Works for both 1D and 2D KS: `model` is the CONSTRUCTED dynamodels model
    returned by build_fom (physical.KS or physical.KS2D). Generation runs
    model.time_integrate in chunks; ks1d spectral history is converted to
    physical snapshots via the model's own transforms. The cache file format
    and layout are identical to the old torch-FOM pipeline -- a torch-saved
    float64 tensor of POST-burn-in physical snapshots, (state_dim, Ncols), at
    get_simulation_path(model)/full_trajectory.pth -- so existing caches hit.
    Cache-hit rule unchanged: cached_cols >= Ntot -> last Ntot columns.

    Args:
        Ntot: Number of kept snapshots to return after discarding i0 transient steps.
        model: dynamodels FOM (from build_fom).
        i0: Number of initial transient steps to skip (default 0). Only used
            when generating from scratch; the cache stores post-burn-in states.

    Returns:
        Trajectory tensor with shape (state_dim, Ntot), float64, CPU.
    """
    get_simulation_path, to_columns, to_psi = _fom_conversions(model)
    n_state = state_dim(model)

    # ~8 MB of dynamodels history per chunk keeps memory flat for any grid size.
    chunk = max(1, 1_000_000 // max(n_state, 1))

    def _advance(nsteps: int):
        """Advance the model nsteps without keeping history; returns nothing."""
        done = 0
        while done < nsteps:
            n = min(chunk, nsteps - done)
            psi, t = model.time_integrate(Nt=n)
            model.update_history(psi[[-1]], t=t[[-1]], reset=True)
            done += n

    def _collect(n_keep: int, keep_current: bool) -> torch.Tensor:
        """Collect n_keep physical snapshot columns, optionally starting with the
        model's current state, integrating in chunks (history reset each chunk)."""
        buf = torch.empty(n_state, n_keep, dtype=torch.float64, device="cpu")
        kept = 0
        if keep_current and n_keep > 0:
            cur = model.current_state[np.newaxis, :, :]
            buf[:, 0] = to_columns(cur)[:, 0]
            kept = 1
        t_start = t_last_print = time.time()
        while kept < n_keep:
            n = min(chunk, n_keep - kept)
            psi, t = model.time_integrate(Nt=n)
            model.update_history(psi[[-1]], t=t[[-1]], reset=True)
            buf[:, kept:kept + n] = to_columns(psi)
            kept += n
            now = time.time()
            if now - t_last_print >= 15.0:
                rate = kept / (now - t_start)
                print(f"  ...integrating snapshot {kept}/{n_keep} ({100 * kept / n_keep:.1f}%), "
                      f"{rate:.1f} steps/s, ETA {(n_keep - kept) / rate:.0f}s", flush=True)
                t_last_print = now
        return buf

    # Get save directory and filename
    save_dir = get_simulation_path(model)
    os.makedirs(save_dir, exist_ok=True)
    save_filename = os.path.join(save_dir, "full_trajectory.pth")

    spatial_dim = n_state  # flattened physical state length (1D: Nx, 2D: Nx*Ny)

    # Try to load cached trajectory. Always load to CPU: the trajectory can be very large
    # (hundreds of thousands of snapshots) and callers already move the slices they need to
    # model.device explicitly (build_local_model, sweep_qlrom_diagnosis, etc.) -- loading the
    # whole thing onto the GPU here just risks an OOM for no benefit.
    if os.path.exists(save_filename):
        try:
            print(f"Loading full trajectory from cache: {save_filename}.")
            # Try with safe_globals if available (for compatibility with older PyTorch/NumPy)
            try:
                safe_globals = [np.dtype]
                if hasattr(np, '_core') and hasattr(np._core, 'multiarray'):
                    safe_globals.append(np._core.multiarray.scalar)
                with torch.serialization.safe_globals(safe_globals):
                    full_trajectory = torch.load(
                        save_filename,
                        map_location="cpu",
                        weights_only=False,  # trusted local files
                    )
            except torch.OutOfMemoryError:
                raise
            except (AttributeError, RuntimeError):
                # Fallback: load without safe_globals if not available
                full_trajectory = torch.load(
                    save_filename,
                    map_location="cpu",
                    weights_only=False,
                )

            # Validate spatial dimension
            if spatial_dim is not None and full_trajectory.shape[0] != spatial_dim:
                raise ValueError(
                    f"Loaded trajectory spatial dimension {full_trajectory.shape[0]} "
                    f"does not match FOM spatial dimension {spatial_dim}."
                )

            # Check if we have enough snapshots
            if full_trajectory.shape[1] >= Ntot:
                return full_trajectory[:, -Ntot:]
            else:
                # Need more steps; extend trajectory from the last cached state
                print(
                    f"Loaded trajectory has {full_trajectory.shape[1]} snapshots, "
                    f"need {Ntot}. Extending by {Ntot - full_trajectory.shape[1]} steps."
                )
                seed = full_trajectory[:, -1].numpy()
                model.update_history(psi=to_psi(seed)[np.newaxis, :, :],
                                     t=np.array([0.0]), reset=True)
                steps_to_add = _collect(Ntot - full_trajectory.shape[1], keep_current=False)
                full_trajectory = torch.cat([full_trajectory, steps_to_add], dim=1)
                torch.save(full_trajectory, save_filename)
        except torch.OutOfMemoryError:
            # A resource error, not a corrupt/incompatible cache -- don't silently discard a
            # valid (possibly hours-expensive) cache and trigger a from-scratch rebuild that
            # would hit the exact same memory pressure. Let it surface.
            raise
        except Exception as e:
            print(f"Failed to load cached trajectory (will rebuild): {e}")
            full_trajectory = None
    else:
        full_trajectory = None

    # Build full trajectory if not loaded or incomplete
    if full_trajectory is None:
        print("Building full trajectory from scratch...")
        # From the model's own initial condition (module.FOM seeds the same
        # deterministic IC the old torch FOMs used): burn i0 transient steps,
        # then keep Ntot states starting with the state AT step i0 -- the same
        # kept-state semantics as the old stepper loop.
        model.update_history(psi=model.psi0[np.newaxis, :, :], t=np.array([0.0]), reset=True)
        if i0 > 0:
            _advance(i0)
        full_trajectory = _collect(Ntot, keep_current=True)
        torch.save(full_trajectory, save_filename)

    return full_trajectory

load_case(module, case, Ntrain, Ntest)

Cached trajectory of one test case, split into a training and a test window.

module is ks1d or ks2d. The trajectory comes from get_full_trajectory (generated and cached on first use), and traj keeps all of it for callers that need more than the two windows.

Returns:

Type	Description
dict	dict with name, fom, cfg, dt, shape (the grid (Nx, Ny) that a flattened
dict	2-D snapshot reshapes to, None in 1-D), traj, Xtrain (N, Ntrain), Xtest
dict	(N, Ntest) and save_dir, the cache directory of the ROMs fitted on this window.

Source code in qlroms/intrusive_qlroms/build_ks.py

def load_case(module, case: str, Ntrain: int, Ntest: int) -> dict:
    """Cached trajectory of one test case, split into a training and a test window.

    `module` is ks1d or ks2d. The trajectory comes from `get_full_trajectory` (generated
    and cached on first use), and `traj` keeps all of it for callers that need more than
    the two windows.

    Returns:
        dict with `name`, `fom`, `cfg`, `dt`, `shape` (the grid (Nx, Ny) that a flattened
        2-D snapshot reshapes to, None in 1-D), `traj`, `Xtrain` (N, Ntrain), `Xtest`
        (N, Ntest) and `save_dir`, the cache directory of the ROMs fitted on this window.
    """
    fom, cfg = build_fom(module, case=case, overrides={"Ntrain": Ntrain, "Ntest": Ntest})
    traj = get_full_trajectory(Ntot=cfg["i0"] + Ntrain + Ntest, model=fom, i0=cfg["i0"])
    dt = float(fom.dt)
    return dict(name=case, fom=fom, cfg=cfg, dt=dt,
                shape=(int(fom.Nx), int(fom.Ny)) if hasattr(fom, "Ny") else None,
                traj=traj, Xtrain=traj[:, :Ntrain], Xtest=traj[:, Ntrain:Ntrain + Ntest],
                save_dir=module.get_simulation_path(model=fom, Ntrain=Ntrain, Ntest=Ntest, dt=dt))

qlroms.intrusive_qlroms.build_pinball

Pinball snapshot loading -- the build_ks analog.

Snapshots are written by PinballFOM.run() (scripts/pinball/run_pinball_fom.py) as _*.npy + snapshot_metadata.npz under qlroms.utils.paths.PINBALL_SNAPSHOTS (or any directory passed in). Models are built and cached by qlroms.intrusive_qlroms.pinball.rom.build_local_model (intrusive Galerkin; dolfinx only on a cache miss). No dolfinx anywhere in this module.

load_mass_matrix(kind='velocity', mesh_dir=None, n_dof=None)

The FEM mass matrix of the pinball's Taylor-Hood velocity ("velocity") or pressure ("pressure") space as a torch sparse tensor -- the inner product of the charts, the clustering, the chart assignment and every error of the models built on the snapshots. Read from the sidecar file written next to the mesh (aux_fenics.export_mass_matrices, run by build_local_model, or run_pinball_fom.py --export-mass), so no dolfinx is needed.

Parameters:

Name	Type	Description	Default
kind	str	"velocity" or "pressure".	'velocity'
mesh_dir		directory of the mesh (default qlroms.utils.paths.PINBALL_MESHES).	None
n_dof	int | None	expected dof count (e.g. the snapshot matrix's rows); a mesh that does not match the snapshots fails here instead of producing a wrong metric.	None

Source code in qlroms/intrusive_qlroms/build_pinball.py

def load_mass_matrix(kind: str = "velocity", mesh_dir=None, n_dof: int | None = None) -> torch.Tensor:
    """The FEM mass matrix of the pinball's Taylor-Hood velocity ("velocity") or pressure
    ("pressure") space as a torch sparse tensor -- the inner product of the charts, the
    clustering, the chart assignment and every error of the models built on the snapshots.
    Read from the sidecar file written next to the mesh (aux_fenics.export_mass_matrices, run
    by build_local_model, or `run_pinball_fom.py --export-mass`), so no dolfinx is needed.

    Args:
        kind: "velocity" or "pressure".
        mesh_dir: directory of the mesh (default qlroms.utils.paths.PINBALL_MESHES).
        n_dof: expected dof count (e.g. the snapshot matrix's rows); a mesh that does not
            match the snapshots fails here instead of producing a wrong metric.
    """
    if kind not in ("velocity", "pressure"):
        raise ValueError(f"kind must be 'velocity' or 'pressure', got {kind!r}.")
    mesh_dir = Path(mesh_dir) if mesh_dir is not None else paths.PINBALL_MESHES
    path = mesh_dir / f"pinball_mass_{kind}.npz"
    if not path.is_file():
        raise FileNotFoundError(
            f"{path} not found -- write it once (needs dolfinx) with "
            "`python scripts/pinball/run_pinball_fom.py configs/pinball_fom.yml --export-mass`.")
    M = load_weight(path)
    if n_dof is not None and M.shape[0] != n_dof:
        raise ValueError(f"the {kind} mass matrix in {mesh_dir} has {M.shape[0]} dofs but the "
                         f"snapshots have {n_dof}: the mesh does not belong to this campaign.")
    return M

load_snapshots(snapshot_dir=None, quantity='velocity', device=None, rdtype=torch.float64, start=0, stop=None)

Stack <quantity>_*.npy snapshots into an (Ndof, Nt) matrix, chronological -- only the saved snapshots start .. stop-1 (default: all), and only those files are read, so a window of a long campaign costs its own size, not the campaign's.

Returns:

Type	Description
(X, dt)	snapshot matrix and the saved-snapshot time spacing from
float	snapshot_metadata.npz (dt * snapshot_interval).

Source code in qlroms/intrusive_qlroms/build_pinball.py

def load_snapshots(snapshot_dir=None, quantity: str = "velocity",
                   device=None, rdtype=torch.float64, start: int = 0,
                   stop: int | None = None) -> tuple[torch.Tensor, float]:
    """Stack `<quantity>_*.npy` snapshots into an (Ndof, Nt) matrix, chronological -- only
    the saved snapshots start .. stop-1 (default: all), and only those files are read, so a
    window of a long campaign costs its own size, not the campaign's.

    Returns:
        (X, dt): snapshot matrix and the saved-snapshot time spacing from
        snapshot_metadata.npz (dt * snapshot_interval).
    """
    snapshot_dir = Path(snapshot_dir) if snapshot_dir is not None else PINBALL_SNAPSHOTS
    files = sorted(snapshot_dir.glob(f"{quantity}_*.npy"))[start:stop]
    if not files:
        raise FileNotFoundError(f"No {quantity}_*.npy snapshots in {snapshot_dir} (window [{start}, {stop})).")
    first = np.load(files[0])
    X = np.empty((first.size, len(files)), dtype=first.dtype)   # filled in place: np.stack would hold
    X[:, 0] = first                                              # every file AND the stack, twice the peak
    for j, f in enumerate(files[1:], start=1):
        X[:, j] = np.load(f)
    X = torch.from_numpy(X).to(rdtype)
    if device is not None:
        X = X.to(device)
    meta = np.load(snapshot_dir / "snapshot_metadata.npz")
    dt = float(meta["dt"]) * int(meta["snapshot_interval"])
    return X, dt

snapshot_time_index(snapshot_dir=None, quantity='velocity')

Monotone integer time index for a regression fit's time_index=..., one entry per sorted <quantity>_*.npy file: +1 where consecutive files are exactly snapshot_interval solver steps apart, +2 (a gap -> non-adjacent) anywhere the step sequence is irregular. The fom_re130 warm-restart seam is the motivating case: ..., 3996, 3998, 3999, 4001, ... -- the 3998 -> 3999 pair is half a snapshot interval apart, so it must not become a regression pair.

Source code in qlroms/intrusive_qlroms/build_pinball.py

def snapshot_time_index(snapshot_dir=None, quantity: str = "velocity") -> np.ndarray:
    """Monotone integer time index for a regression fit's time_index=..., one
    entry per sorted `<quantity>_*.npy` file: +1 where consecutive files are exactly
    `snapshot_interval` solver steps apart, +2 (a gap -> non-adjacent) anywhere the
    step sequence is irregular. The fom_re130 warm-restart seam is the motivating
    case: ..., 3996, 3998, 3999, 4001, ... -- the 3998 -> 3999 pair is half a
    snapshot interval apart, so it must not become a regression pair."""
    snapshot_dir = Path(snapshot_dir) if snapshot_dir is not None else PINBALL_SNAPSHOTS
    files = sorted(snapshot_dir.glob(f"{quantity}_*.npy"))
    if not files:
        raise FileNotFoundError(f"No {quantity}_*.npy snapshots in {snapshot_dir}.")
    steps = np.array([int(f.stem.rsplit("_", 1)[1]) for f in files])
    interval = int(np.load(snapshot_dir / "snapshot_metadata.npz")["snapshot_interval"])
    return np.concatenate([[0], np.cumsum(np.where(np.diff(steps) == interval, 1, 2))])

merge_restart_snapshot_metadata(snapshot_dir, previous_metadata)

Fold an exact restart leg's own snapshot_metadata.npz back into the campaign it continued, into one file describing the whole thing.

PinballFOM.run() clears its bookkeeping and overwrites snapshot_metadata.npz on every call (qlroms.intrusive_qlroms.pinball.solver.py's run()), so after init_velocity_prev-based exact restart (see PinballFOMParameters), the file on disk describes only the just-completed leg's own steps -- not the campaign as a whole. This merges it back with a saved copy of the metadata from before that leg ran, producing a single file indistinguishable from what one uninterrupted run of the combined step count would have written: saved_steps/saved_times concatenated and checked for a uniform, gap- and duplicate-free spacing (an exact restart's saved_times are already campaign-relative, not leg-local -- solver.py's run() seeds self.time from the restart's own step offset), everything else taken from the first leg except num_steps (summed).

Parameters:

Name	Type	Description	Default
snapshot_dir		the campaign directory (its CURRENT snapshot_metadata.npz is read as the just-completed leg and then overwritten with the merged result).	required
previous_metadata		path to a copy of snapshot_metadata.npz taken before that leg's run() call (its own overwrite would otherwise destroy the earlier legs' record); run() leaves it as snapshot_metadata_before_.npz.	required

Returns: the (now merged) snapshot_metadata.npz path.

Source code in qlroms/intrusive_qlroms/build_pinball.py

def merge_restart_snapshot_metadata(snapshot_dir, previous_metadata) -> Path:
    """Fold an exact restart leg's own snapshot_metadata.npz back into the campaign it
    continued, into one file describing the whole thing.

    PinballFOM.run() clears its bookkeeping and overwrites snapshot_metadata.npz on every
    call (qlroms.intrusive_qlroms.pinball.solver.py's run()), so after
    init_velocity_prev-based exact restart (see PinballFOMParameters), the file on disk
    describes only the just-completed leg's own steps -- not the campaign as a whole. This
    merges it back with a saved copy of the metadata from before that leg ran, producing a
    single file indistinguishable from what one uninterrupted run of the combined step count
    would have written: saved_steps/saved_times concatenated and checked for a uniform,
    gap- and duplicate-free spacing (an exact restart's saved_times are already
    campaign-relative, not leg-local -- solver.py's run() seeds self.time from the restart's
    own step offset), everything else taken from the first leg except num_steps (summed).

    Args:
        snapshot_dir: the campaign directory (its CURRENT snapshot_metadata.npz is read as
            the just-completed leg and then overwritten with the merged result).
        previous_metadata: path to a copy of snapshot_metadata.npz taken before that leg's
            run() call (its own overwrite would otherwise destroy the earlier legs' record);
            run() leaves it as snapshot_metadata_before_<the leg's first index>.npz.
    Returns:
        the (now merged) snapshot_metadata.npz path.
    """
    snapshot_dir = Path(snapshot_dir)
    meta_path = snapshot_dir / "snapshot_metadata.npz"
    prev = dict(np.load(previous_metadata))
    new = dict(np.load(meta_path))

    for key in ("dt", "snapshot_interval", "reynolds", "wall_bc"):
        if not np.array_equal(prev[key], new[key]):
            raise ValueError(f"{key} differs between the two legs ({prev[key]!r} vs "
                             f"{new[key]!r}) -- these are not two legs of the same campaign.")

    saved_steps = np.concatenate([prev["saved_steps"], new["saved_steps"]])
    saved_times = np.concatenate([prev["saved_times"], new["saved_times"]])
    order = np.argsort(saved_steps)
    saved_steps, saved_times = saved_steps[order], saved_times[order]
    interval = int(prev["snapshot_interval"])
    gaps = np.diff(saved_steps)
    if not np.all(gaps == interval):
        bad = np.nonzero(gaps != interval)[0]
        raise ValueError(f"merged saved_steps is not uniformly spaced by {interval} (a gap, "
                         f"overlap, or duplicate at the join): breaks at "
                         f"{list(zip(saved_steps[bad], saved_steps[bad + 1], strict=True))}.")

    np.savez(
        meta_path,
        saved_steps=saved_steps,
        saved_times=saved_times,
        dt=prev["dt"],
        snapshot_interval=prev["snapshot_interval"],
        snapshot_start_step=prev["snapshot_start_step"],
        num_steps=int(prev["num_steps"]) + int(new["num_steps"]),
        reynolds=prev["reynolds"],
        wall_bc=prev["wall_bc"],
    )
    return meta_path

build_sensor_indices(fom, n_sensors, sampling='random', seed=0)

Sensor dof indices into the stacked velocity state -- the build_ks analog.

The pinball state is an unstructured Taylor-Hood dof vector, so there is no grid to lay sensors on: "random" (the qlROM-DA convention) draws sorted velocity-dof indices, "uniform" strides the dof vector (index-uniform, not space-uniform). fom is any case handle whose .N is the state dimension.

Source code in qlroms/intrusive_qlroms/build_pinball.py

def build_sensor_indices(fom, n_sensors, sampling="random", seed=0):
    """Sensor dof indices into the stacked velocity state -- the build_ks analog.

    The pinball state is an unstructured Taylor-Hood dof vector, so there is no
    grid to lay sensors on: "random" (the qlROM-DA convention) draws sorted
    velocity-dof indices, "uniform" strides the dof vector (index-uniform, not
    space-uniform). `fom` is any case handle whose .N is the state dimension.
    """
    n_state = int(fom.N)
    if not 1 <= n_sensors <= n_state:
        raise ValueError(f"n_sensors must be in [1, {n_state}], got {n_sensors}.")
    if sampling == "random":
        return np.sort(np.random.default_rng(seed).choice(n_state, size=n_sensors,
                                                          replace=False))
    if sampling != "uniform":
        raise ValueError(f"unknown sensor sampling mode: {sampling}")
    return np.linspace(0, n_state - 1, n_sensors, dtype=int)

campaign_name(reynolds, tag=None)

Snapshot-campaign name: fom_re140 untagged (the historical layout), fom_re140__ for a repeated run at the same Reynolds number. The name is the snapshot directory, the prefix of the run results and (via the tag) part of the model cache key, so a tagged campaign shares -- and overwrites -- nothing with the untagged one.

Source code in qlroms/intrusive_qlroms/build_pinball.py

def campaign_name(reynolds: float, tag: str | None = None) -> str:
    """Snapshot-campaign name: fom_re140 untagged (the historical layout), fom_re140__<tag>
    for a repeated run at the same Reynolds number. The name is the snapshot directory, the
    prefix of the run results and (via the tag) part of the model cache key, so a tagged
    campaign shares -- and overwrites -- nothing with the untagged one."""
    if tag and not re.fullmatch(r"[A-Za-z0-9][A-Za-z0-9._-]*", tag):
        raise ValueError(f"campaign tag {tag!r}: use letters, digits and . _ - only")
    return f"fom_re{reynolds:g}" + (f"__{tag}" if tag else "")

config_campaign(config)

Campaign name of a FOM config (configs/pinball_.yml): its fom.reynolds and top-level tag (absent or null -> the untagged fom_re). The one place a campaign's tag is set: run_pinball_fom.py, the build_re_roms.py scripts and the postprocessing all read it here.

Source code in qlroms/intrusive_qlroms/build_pinball.py

def config_campaign(config) -> str:
    """Campaign name of a FOM config (configs/pinball_*.yml): its fom.reynolds and top-level
    `tag` (absent or null -> the untagged fom_re<Re>). The one place a campaign's tag is set:
    run_pinball_fom.py, the build_re*_roms.py scripts and the postprocessing all read it here."""
    cfg = yaml.safe_load(Path(config).read_text()) or {}
    return campaign_name(float((cfg.get("fom") or {}).get("reynolds", 130.0)), cfg.get("tag"))

parse_campaign(dataset)

(reynolds, tag) of a campaign name, the inverse of campaign_name: fom_re140 -> (140.0, None), fom_re140__long -> (140.0, "long").

Source code in qlroms/intrusive_qlroms/build_pinball.py

def parse_campaign(dataset: str) -> tuple[float, str | None]:
    """(reynolds, tag) of a campaign name, the inverse of campaign_name:
    fom_re140 -> (140.0, None), fom_re140__long -> (140.0, "long")."""
    base, _, tag = dataset.partition("__")
    return float(base.split("_re")[-1]), tag or None

load_intrusive_model(dataset, Xtrain, *, K, r_velocity, r_pressure, dt, random_state=1, build=False, cluster_space='pod_lossless', cluster_metric='mass', pressure_centroid='poisson_from_velocity_centroid', i0=0)

Intrusive qlGalerkin for a snapshot campaign, from the build_local_model cache.

The Reynolds number and the tag come from the campaign name (fom_re120 -> 120, fom_re140__long -> 140 with its own cache files, see campaign_name). On a cache miss the model is built only if build=True (dolfinx, about 1 h at K=8 on 6000 snapshots); otherwise this prints the build call and returns None, so notebooks stay cheap. The pressure snapshots are read only when a build is needed.

Parameters:

Name	Type	Description	Default
dataset	str	campaign name, e.g. "fom_re120" or "fom_re140__long".	required
Xtrain		(N, T) velocity training snapshots -- the window [i0, i0+T) of the campaign's own snapshot sequence. Pass i0 whenever Xtrain does not start at the campaign's first saved snapshot (e.g. Xtrain = X[:, i0:i0+Ntrain] from a larger loaded X), or the pressure window read on a cache miss would silently come from the WRONG instants (always the campaign's own first T snapshots), not the ones Xtrain's velocity columns actually are.	required
i0	int	index of Xtrain's first column within the campaign's own snapshot sequence (default 0, i.e. Xtrain starts at the beginning).	0

Returns: the qlGalerkin, or None. Each chart has r_velocity + r_pressure + 1 velocity coordinates (POD modes plus one supremizer per pressure mode and pressure mean).

Source code in qlroms/intrusive_qlroms/build_pinball.py

def load_intrusive_model(dataset: str, Xtrain, *, K: int, r_velocity: int, r_pressure: int,
                         dt: float, random_state: int = 1, build: bool = False,
                         cluster_space: str = "pod_lossless",
                         cluster_metric: str = "mass",
                         pressure_centroid: str = "poisson_from_velocity_centroid",
                         i0: int = 0):
    """Intrusive qlGalerkin for a snapshot campaign, from the `build_local_model` cache.

    The Reynolds number and the tag come from the campaign name (fom_re120 -> 120,
    fom_re140__long -> 140 with its own cache files, see campaign_name). On a cache
    miss the model is built only if `build=True` (dolfinx, about 1 h at K=8 on 6000
    snapshots); otherwise this prints the build call and returns None, so notebooks
    stay cheap. The pressure snapshots are read only when a build is needed.

    Args:
        dataset: campaign name, e.g. "fom_re120" or "fom_re140__long".
        Xtrain: (N, T) velocity training snapshots -- the window [i0, i0+T) of the
            campaign's own snapshot sequence. Pass `i0` whenever Xtrain does not start at
            the campaign's first saved snapshot (e.g. Xtrain = X[:, i0:i0+Ntrain] from a
            larger loaded X), or the pressure window read on a cache miss would silently
            come from the WRONG instants (always the campaign's own first T snapshots),
            not the ones Xtrain's velocity columns actually are.
        K, r_velocity, r_pressure, dt, random_state, cluster_space, cluster_metric,
            pressure_centroid: forwarded to build_local_model AND to intrusive_model_path,
            so the cache key always matches exactly what a build with these settings would
            produce (Ntrain is read off Xtrain itself).
        i0: index of Xtrain's first column within the campaign's own snapshot sequence
            (default 0, i.e. Xtrain starts at the beginning).
    Returns:
        the qlGalerkin, or None. Each chart has r_velocity + r_pressure + 1 velocity
        coordinates (POD modes plus one supremizer per pressure mode and pressure mean).
    """
    from .pinball.rom import build_local_model, intrusive_model_path

    reynolds, tag = parse_campaign(dataset)
    cache = intrusive_model_path(K, r_velocity, r_pressure, reynolds, Ntrain=int(Xtrain.shape[1]), dt=dt,
                                 random_state=random_state, cluster_space=cluster_space,
                                 cluster_metric=cluster_metric, pressure_centroid=pressure_centroid,
                                 tag=tag)
    if not cache.is_file() and not build:
        print(f"no cached intrusive model at {cache} -- build it with build=True (dolfinx, ~1 h at K=8)")
        return None
    Xp = None
    if not cache.is_file():   # only the training window's pressure files are read
        Xp = load_snapshots(_dataset_dir(dataset), "pressure", start=i0, stop=i0 + Xtrain.shape[1])[0]
        if Xp.shape[1] != Xtrain.shape[1]:
            raise ValueError(f"{dataset} has only {Xp.shape[1]} pressure snapshots in [{i0}, "
                             f"{i0 + Xtrain.shape[1]}), the training window of Xtrain")
    model = build_local_model(Xtrain, Xp, K=K, r_velocity=r_velocity, r_pressure=r_pressure,
                              reynolds=reynolds, dt=dt, random_state=random_state,
                              cluster_space=cluster_space, cluster_metric=cluster_metric,
                              pressure_centroid=pressure_centroid, tag=tag)
    print(f"intrusive model (Re {reynolds:g}): K={model.K}, chart dim r={model.r} = "
          f"{r_velocity} velocity POD + {r_pressure + 1} supremizers")
    return model