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) distinguishes intrusive compilations from non-intrusive ones (OpInf/ESN) 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 same
        compute as the OpInf family once reduced, so both members subclass the shared
        qlroms.charts.QuadROM; 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

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 same compute as the OpInf family once reduced, so both members subclass the shared qlroms.charts.QuadROM; 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 same
    compute as the OpInf family once reduced, so both members subclass the shared
    qlroms.charts.QuadROM; 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

qlroms.intrusive_qlroms.build_pinball

Pinball snapshot loading -- the build_ks.get_full_trajectory 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).

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

Stack <quantity>_*.npy snapshots into an (Ndof, Nt) matrix, chronological.

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) -> tuple[torch.Tensor, float]:
    """Stack `<quantity>_*.npy` snapshots into an (Ndof, Nt) matrix, chronological.

    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"))
    if not files:
        raise FileNotFoundError(f"No {quantity}_*.npy snapshots in {snapshot_dir}.")
    X = torch.from_numpy(np.stack([np.load(f) for f in files], axis=1)).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