Skip to content

qlroms.intrusive_qlroms.ks1d — 1-D Kuramoto–Sivashinsky

\(u_t + u u_x + u_{xx} + \nu u_{xxxx} = 0\) on a periodic domain, spectral (ETDRK4) FOM. The reduced operators come straight from the FOM's Fourier multipliers, so each ROM member is an exact Galerkin projection (with the nonlinearity lifted through \(\bm{\Phi}_k\) and the centroid at every stage).

ks1d test cases: space-time plots of u(x, t) for the quasi-periodic and chaotic cases

The two TEST_CASES, 1500 steps past the transient: quasi-periodic (\(\nu = 16/71\), \(L_x = 2\pi\), \(N_x = 128\), \(\Delta t = 0.1\), \(\lambda_1 = 0\)) and chaotic (\(\nu = 1\), \(L_x = 20\pi\), \(N_x = 128\), \(\Delta t = 0.05\), \(\lambda_1 \approx 0.062\)).

At a glance

Name One-liner
config.TEST_CASES The cases (quasi-periodic, chaotic): physics, generation windows, usual \((K, r)\), characterization timescales.
config.KSConfig Torch-side grid and spectral operators of a case; base of ROM, carried by FOM as .cfg.
fom.FOM The full-order model: a dynamodels.physical.KS built from a case name (build_ks.build_fom).
rom.ROM One chart: Chart geometry + the projected operators and the ETDRK4 step_reduced.
rom.build_local_model(Xtrain, fom, r, K, save_dir) Cluster, per-cluster POD, atlas, compile the ROMs into a qlGalerkin; cached.
rom.get_simulation_path(model, Ntrain) Cache directory keyed on \((\nu, L_x, N_x, \Delta t)\) and Ntrain.

qlroms.intrusive_qlroms.ks1d.config

ks1d/config.py -- 1D Kuramoto-Sivashinsky case definitions.

TEST_CASES: one dict per case (visc/Lx/Nx/dt/lamb1, the generation windows i0/Ntrain/Ntest, the usual (K, r), characterization timescales). KSConfig: torch-side grids and spectral operators (Lhat, Dhat, Ghat, L, G) of a case; base of rom.ROM and carried by fom.FOM as .cfg.

The full-order model lives in fom.py (FOM); the ROM stack in rom.py (ROM, qlROM, build_local_model, get_simulation_path).

KSConfig dataclass

Configuration class for 1D KS system, used as base for both FOM and ROM dataclasses. Equatuons: u_t + u u_x + u_xx + nu u_xxxx = 0, x in [0, L], t > 0 u(x, 0) = u0(x) u(x+L, t) = u(x, t) (periodic BCs)

Parameters:

Name Type Description Default
visc

Viscosity parameter nu.

required
Lx

Domain length L.

required
Nx

Number of spatial grid points.

required
dt

Time step size for ETDRK4.

required
lamb1

Leading Lyapunov exponent (used for estimating predictability time).

required
M int

Number of contour points for ETDRK4 resolvent integrals (default 32).

32

The domain is multiplied by 2pi for consistency with standard KS literature, so L=102*pi corresponds to the commonly studied chaotic regime with ~40 active modes and lam1 ~= 0.062.

Source code in qlroms/intrusive_qlroms/ks1d/config.py
 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
@dataclass
class KSConfig:
    """ Configuration class for 1D KS system, used as base for both FOM and ROM dataclasses.
        Equatuons:
            u_t + u u_x + u_xx + nu u_xxxx = 0,  x in [0, L], t > 0
            u(x, 0) = u0(x)
            u(x+L, t) = u(x, t) (periodic BCs)

        Parameters:
            visc: Viscosity parameter nu.
            Lx: Domain length L.
            Nx: Number of spatial grid points.
            dt: Time step size for ETDRK4.
            lamb1: Leading Lyapunov exponent (used for estimating predictability time).
            M: Number of contour points for ETDRK4 resolvent integrals (default 32).

        The domain is multiplied by 2*pi for consistency with standard KS literature,
        so L=10*2*pi corresponds to the commonly studied chaotic regime with ~40
        active modes and lam1 ~= 0.062.

    """

    # --- numerics / dtypes ---
    device: torch.device = field(default_factory=_default_device)
    cdtype: torch.dtype = torch.complex128
    rdtype: torch.dtype = torch.float64

    # --- KS parameters ---
    case: str = "TBD"  # "quasi-periodic" or "chaotic"; sets default visc and Lx
    M: int = 32  # contour points for ETDRK4 resolvent integrals

    @property
    def properties(self) -> dict:
        if self.case not in TEST_CASES:
            raise ValueError(f"Unknown case '{self.case}'. Valid cases: {list(TEST_CASES.keys())}.")

        return TEST_CASES[self.case]

    def __post_init__(self):
        props = self.properties
        self.lamb1: float = props["lamb1"]
        self.visc: float = props["visc"]
        self.Lx: float = props["Lx"]
        self.Nx: int = props["Nx"]
        self.dt: float = props["dt"]


        self.t_lyap: float = 1.0 / self.lamb1 if self.lamb1 > 0. else float('inf')

        self.x: torch.Tensor = (self.Lx * torch.arange(self.Nx, dtype=self.rdtype, device=self.device) / self.Nx).reshape(self.Nx, 1)

        eye = torch.eye(self.Nx, dtype=self.cdtype, device=self.device)
        self.FFTmat: torch.Tensor = torch.fft.fft(eye, dim=0) / torch.sqrt(torch.tensor(self.Nx, dtype=self.rdtype, device=self.device))
        self.iFFTmat_H: torch.Tensor = self.FFTmat.conj().T

        k = torch.cat((torch.arange(0, self.Nx // 2, device=self.device),
                       torch.arange(-self.Nx // 2, 0, device=self.device)), dim=0)
        self.kx: torch.Tensor = k.reshape(self.Nx, 1).to(self.rdtype)

        self.alpha: torch.Tensor = (2 * torch.pi * self.kx / self.Lx).to(self.rdtype)
        self.Lhat: torch.Tensor = (self.alpha**2 - self.visc * self.alpha**4).to(self.cdtype)
        self.Dhat: torch.Tensor = (1j * self.alpha).to(self.cdtype)
        self.Ghat: torch.Tensor = (-0.5 * self.Dhat).to(self.cdtype)

    @property
    def N(self) -> int:
        """Flattened state-vector length. 1D: N == Nx (matches 2D's Nx*Ny)."""
        return self.Nx

    @property
    def base_cfg(self) -> dict:
        return {f.name: getattr(self, f.name) for f in fields(KSConfig)}

    # ----- physical-space matrices -----
    @cached_property
    def L(self) -> torch.Tensor:
        return torch.real(self.iFFTmat_H @ torch.diag(self.Lhat.flatten()) @ self.FFTmat)

    @cached_property
    def G(self) -> torch.Tensor:
        return torch.real(self.iFFTmat_H @ torch.diag(self.Ghat.flatten()) @ self.FFTmat)

N property

Flattened state-vector length. 1D: N == Nx (matches 2D's Nx*Ny).

qlroms.intrusive_qlroms.ks1d.fom

Full-order model of a ks1d case: dynamodels.physical.KS built from a case name; carries its KSConfig as .cfg -- see build_ks.build_fom.

FOM

Bases: KS

Full-order model of a ks1d case: a dynamodels.physical.KS built from a case name.

The case's PHYSICAL parametrization is passed straight through -- KS(Nx=Nx, dt=dt, L=Lx, nu=visc) -- because dynamodels >= 0.3.3 honours an independent (nu, L) pair. So there is no rescaling anywhere: self.dt IS the case's physical time step, self.L IS Lx, self.nu IS visc, and the spectral state is the physical field's rfft (no sqrt(visc) amplitude factor).

The torch-side case config is carried whole as self.cfg (a KSConfig); pass THAT to qlroms.utils.diagnosis, whose FOM argument reads N/dt in case units -- dynamodels' Model.N is the analysis-augmented size Nphi + Na + Nq, not Nx. Call sites wanting the state size use fom.Nx or fom.cfg.N.

The handful of case-side quantities the torch ROM stack reads off the FOM are real properties delegating to self.cfg: visc, Lx, lamb1, device, rdtype, cdtype. Nothing else is aliased -- in particular N, M, m, x and dt keep their dynamodels Model meaning.

The initial condition is the deterministic u0 = cos(x) on the case grid, spectrally seeded, so trajectories (and the caches keyed on visc/Lx/Nx/dt) are unchanged.

Source code in qlroms/intrusive_qlroms/ks1d/fom.py
12
13
14
15
16
17
18
19
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
class FOM(KS):
    """Full-order model of a ks1d case: a dynamodels.physical.KS built from a case name.

    The case's PHYSICAL parametrization is passed straight through --
    ``KS(Nx=Nx, dt=dt, L=Lx, nu=visc)`` -- because dynamodels >= 0.3.3 honours an
    independent (nu, L) pair. So there is no rescaling anywhere: ``self.dt`` IS the
    case's physical time step, ``self.L`` IS Lx, ``self.nu`` IS visc, and the
    spectral state is the physical field's rfft (no sqrt(visc) amplitude factor).

    The torch-side case config is carried whole as ``self.cfg`` (a `KSConfig`); pass
    THAT to `qlroms.utils.diagnosis`, whose FOM argument reads ``N``/``dt`` in case units --
    dynamodels' ``Model.N`` is the analysis-augmented size Nphi + Na + Nq, not Nx.
    Call sites wanting the state size use ``fom.Nx`` or ``fom.cfg.N``.

    The handful of case-side quantities the torch ROM stack reads off the FOM are
    real properties delegating to ``self.cfg``: `visc`, `Lx`, `lamb1`, `device`,
    `rdtype`, `cdtype`. Nothing else is aliased -- in particular `N`, `M`, `m`, `x`
    and `dt` keep their dynamodels `Model` meaning.

    The initial condition is the deterministic u0 = cos(x) on the case grid,
    spectrally seeded, so trajectories (and the caches keyed on visc/Lx/Nx/dt) are
    unchanged.
    """

    # `case` joins KS's structural params so ntsa.respawn rebuilds the SAME case
    # (it is not a scalar class default, so Model.filename ignores it).
    fixed_params = [*KS.fixed_params, 'case']

    def __init__(self, case: str = "chaotic", **model_dict):
        # cfg BEFORE super().__init__: Model's kwarg-setattr loop probes hasattr(),
        # and the properties below read self.cfg.
        self.cfg = KSConfig(case=case)
        c = self.cfg
        u0 = np.cos(c.Lx * np.arange(c.Nx) / c.Nx)
        defaults = dict(Nx=int(c.Nx), dt=float(c.dt), L=float(c.Lx), nu=float(c.visc),
                        psi0=np.fft.rfft(u0)[:, None])
        super().__init__(**{**defaults, **model_dict})

    # ---- case-side surface, delegating to self.cfg (no attributes bolted on) ----

    @property
    def case(self) -> str:
        """Key into TEST_CASES this FOM was built from."""
        return self.cfg.case

    @property
    def visc(self) -> float:
        """Physical viscosity nu of the case (== self.nu, in case notation)."""
        return float(self.cfg.visc)

    @property
    def Lx(self) -> float:
        """Physical domain length (== self.L, in case notation)."""
        return float(self.cfg.Lx)

    @property
    def lamb1(self) -> float:
        """Leading Lyapunov exponent of the case."""
        return float(self.cfg.lamb1)

    @property
    def device(self) -> torch.device:
        return self.cfg.device

    @property
    def rdtype(self) -> torch.dtype:
        return self.cfg.rdtype

    @property
    def cdtype(self) -> torch.dtype:
        return self.cfg.cdtype

case property

Key into TEST_CASES this FOM was built from.

visc property

Physical viscosity nu of the case (== self.nu, in case notation).

Lx property

Physical domain length (== self.L, in case notation).

lamb1 property

Leading Lyapunov exponent of the case.

qlroms.intrusive_qlroms.ks1d.rom

Reduced-order models for the 1D KS system, on the unified qlroms stack.

ROM (a qlroms.charts.Chart) is the ONLY class this case defines: one cluster's geometry plus its intrusive Galerkin operators and step_reduced. The quantized-local model is the generic compilation qlroms.base.qlROM -- a compilation of Galerkin members IS a qlGalerkin -- assembled by build_local_model. The non-intrusive OpInf family is entirely equation-free and lives in qlroms.data_driven_qlroms.opinf (qlOpinf.from_snapshots builds it from the same charts/atlas). Stepping is model.step(apod); there is no case timestepping module.

ROM dataclass

Bases: Chart, KSConfig, ABC

Single-cluster reduced-order model: a qlroms.charts.Chart (local POD geometry) plus this case's intrusive Galerkin operators. Phi and centroid are fixed after construction, so all cached_property operators are computed once and never invalidated.

Source code in qlroms/intrusive_qlroms/ks1d/rom.py
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
@dataclass
class ROM(Chart, KSConfig, ABC):
    """Single-cluster reduced-order model: a qlroms.charts.Chart (local POD geometry)
    plus this case's intrusive Galerkin operators. Phi and centroid are fixed after
    construction, so all cached_property operators are computed once and never
    invalidated."""


    Phi: torch.Tensor = field(kw_only=True)       # (Nx, r) -- required, set at construction
    centroid: torch.Tensor = field(kw_only=True)  # (Nx, 1) -- required, set at construction

    @property
    def r(self) -> int:
        return self.Phi.shape[1]

    @cached_property
    def L_rom(self) -> torch.Tensor:
        Phic = self.Phi.clone().to(self.cdtype)
        return torch.real(Phic.conj().T @ self.L.to(self.cdtype) @ Phic)

    @cached_property
    def galerkin_rom(self) -> tuple[torch.Tensor, torch.Tensor]:
        """(Lc, G) -- constant centroid linear forcing and projected nonlinear operator."""
        Phic = self.Phi.clone().to(self.cdtype).conj().T
        Lc = torch.real(Phic @ self.L.to(self.cdtype) @ self.centroid.clone().to(self.cdtype))
        G  = torch.real(Phic @ self.G.to(self.cdtype))
        return Lc.to(self.rdtype), G.to(self.rdtype)

    # ------------------------------------------------------------------ #
    # ETDRK4 coefficients for the reduced linear operator                #
    # ------------------------------------------------------------------ #

    @cached_property
    def etdrk4_rom(self) -> tuple:
        """(E, E2, Q, f1, f2, f3) ETDRK4 coefficients for reduced model."""
        L_np = self.L_rom.detach().cpu().numpy()
        E, E2, Q, f1, f2, f3 = build_etdrk4_coeffs(L_np, self.dt, M=self.M, R=15.0)

        def _t(arr: np.ndarray) -> torch.Tensor:
            return torch.from_numpy(arr).to(dtype=self.rdtype, device=self.device)

        return tuple(_t(x) for x in (E, E2, Q, f1, f2, f3))

    # ------------------------------------------------------------------ #
    # One reduced ETDRK4 step (pure (r, 1) state, this cluster's operators)
    # ------------------------------------------------------------------ #

    def step_reduced(self, a: torch.Tensor) -> torch.Tensor:
        """One ETDRK4 step of the pure (r, 1) reduced state with this cluster's
        intrusive Galerkin dynamics (nonlinearity lifted through Phi/centroid, as
        the KS equation requires). The non-intrusive family lives entirely in
        qlroms.data_driven_qlroms.opinf (qlOpinf.from_snapshots)."""
        E, E2, Q, f1, f2, f3 = self.etdrk4_rom
        Lc, G = self.galerkin_rom
        Phi, cent = self.Phi, self.centroid

        Nu = Lc + G @ ((Phi @ a + cent).pow(2))
        aa = E2 @ a + Q @ Nu
        Na = Lc + G @ ((Phi @ aa + cent).pow(2))
        bb = E2 @ a + Q @ Na
        Nb = Lc + G @ ((Phi @ bb + cent).pow(2))
        cc = E2 @ aa + Q @ (2 * Nb - Nu)
        Nc = Lc + G @ ((Phi @ cc + cent).pow(2))
        return E @ a + f1 @ Nu + 2 * f2 @ (Na + Nb) + f3 @ Nc

galerkin_rom cached property

(Lc, G) -- constant centroid linear forcing and projected nonlinear operator.

etdrk4_rom cached property

(E, E2, Q, f1, f2, f3) ETDRK4 coefficients for reduced model.

step_reduced(a)

One ETDRK4 step of the pure (r, 1) reduced state with this cluster's intrusive Galerkin dynamics (nonlinearity lifted through Phi/centroid, as the KS equation requires). The non-intrusive family lives entirely in qlroms.data_driven_qlroms.opinf (qlOpinf.from_snapshots).

Source code in qlroms/intrusive_qlroms/ks1d/rom.py
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def step_reduced(self, a: torch.Tensor) -> torch.Tensor:
    """One ETDRK4 step of the pure (r, 1) reduced state with this cluster's
    intrusive Galerkin dynamics (nonlinearity lifted through Phi/centroid, as
    the KS equation requires). The non-intrusive family lives entirely in
    qlroms.data_driven_qlroms.opinf (qlOpinf.from_snapshots)."""
    E, E2, Q, f1, f2, f3 = self.etdrk4_rom
    Lc, G = self.galerkin_rom
    Phi, cent = self.Phi, self.centroid

    Nu = Lc + G @ ((Phi @ a + cent).pow(2))
    aa = E2 @ a + Q @ Nu
    Na = Lc + G @ ((Phi @ aa + cent).pow(2))
    bb = E2 @ a + Q @ Na
    Nb = Lc + G @ ((Phi @ bb + cent).pow(2))
    cc = E2 @ aa + Q @ (2 * Nb - Nu)
    Nc = Lc + G @ ((Phi @ cc + cent).pow(2))
    return E @ a + f1 @ Nu + 2 * f2 @ (Na + Nb) + f3 @ Nc

build_local_model(Xtrain, FOM, r=30, K=5, save_dir='.', clustering_kwargs=None)

Train a local POD-ROM from training snapshots and return model data. Parameters: Xtrain: Training snapshot matrix, shape (Nx, Ntrain). FOM: FOM instance providing KS parameters and operators. r: Number of POD modes for each local ROM. K: Number of clusters/local ROMs. save_dir: Directory path to save the trained model. clustering_kwargs: Forwarded to fit_clusters (random_state, kmeans_method, kmeans_n_init, kmeans_max_iter, assign_overlapping, overlap_tolerance). Returns: qlGalerkin instance with trained local ROM data.

Source code in qlroms/intrusive_qlroms/ks1d/rom.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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def build_local_model(Xtrain,
                      FOM,
                      r: int = 30,
                      K: int = 5,
                      save_dir: str = '.',
                      clustering_kwargs: dict | None = None) -> qlGalerkin:
    """Train a local POD-ROM from training snapshots and return model data.
    Parameters:
        Xtrain: Training snapshot matrix, shape (Nx, Ntrain).
        FOM: FOM instance providing KS parameters and operators.
        r: Number of POD modes for each local ROM.
        K: Number of clusters/local ROMs.
        save_dir: Directory path to save the trained model.
        clustering_kwargs: Forwarded to fit_clusters (random_state, kmeans_method,
            kmeans_n_init, kmeans_max_iter, assign_overlapping, overlap_tolerance).
    Returns:
        qlGalerkin instance with trained local ROM data.
    """

    os.makedirs(save_dir, exist_ok=True)
    local_model_path = os.path.join(save_dir, f"local_model_K{K}.pth")

    assert FOM.Nx == Xtrain.shape[0], f"Snapshot spatial dimension {Xtrain.shape[0]} does not match FOM Nx={FOM.Nx}."

    # FOM may be the ks1d.FOM (a dynamodels model carrying its KSConfig as .cfg) or a
    # torch-side KSConfig directly; the member ROMs are built from the KSConfig kwargs.
    base_cfg = FOM.base_cfg if isinstance(FOM, KSConfig) else FOM.cfg.base_cfg


    # Internal function to load a cached model and check its compatibility with the requested K and r.

    def _load_local_model(path):
        loaded_model = torch.load(path, map_location=FOM.device,
                                  weights_only=False)  # trusted local files

        #  Basic sanity checks on the loaded model
        if loaded_model.K != K:
            raise ValueError(f"Cached model at {local_model_path} has K={loaded_model.K}, expected K={K}.")
        if loaded_model.Phi_all is None or loaded_model.centroids is None:
            raise ValueError(f"Cached model at {local_model_path} is incomplete.")
        if r > loaded_model.Phi_all.shape[1]:
            print(f"Cached model has insufficient r={loaded_model.Phi_all.shape[1]} for requested r={r}. Rebuilding model.")
            return None
        else:
            return loaded_model



    def _build_and_save_model() -> qlGalerkin:
        clusters, visited_clusters, cluster_sizes, Xtrain_augment, _ = fit_clusters(
            Xtrain, K, **(clustering_kwargs or {}))

        # assign_overlapping duplicates boundary snapshots into every chart that claims
        # them, so the labels index the AUGMENTED matrix, not Xtrain (ks2d does the same).
        Xsrc = Xtrain
        if len(visited_clusters) > Xtrain.shape[1]:
            print(f"Augmented snapshot matrix has {Xtrain_augment.shape[1]} snapshots "
                  f"(original {Xtrain.shape[1]}).")
            Xsrc = Xtrain_augment
        print(f"Fitted K={K} clusters with sizes: {cluster_sizes}.")
        r_store = min(max(r, 30), FOM.Nx, min(cluster_sizes))

        if r > r_store:
            raise ValueError(
                f"Requested r={r} exceeds r_store={r_store} for K={K}. "
                f"The snapshot matrix per cluster is ({FOM.Nx} x min_cluster={min(cluster_sizes)}), "
                f"so at most min(Nx={FOM.Nx}, min_cluster={min(cluster_sizes)})={r_store} modes are available. "
                f"Reduce r or Nx."
            )

        # local_rom_full = _build_model_from_assignments(clusters, visited_clusters, r_store)
        Phimat = torch.zeros((FOM.Nx, r_store, K), dtype=FOM.rdtype, device=FOM.device)

        for k in range(K):
            cluster_indices = (visited_clusters == k)
            if K == 1:
                # Single cluster spans the whole array: reuse Xtrain in place instead of a
                # boolean-mask copy, which would otherwise double peak memory for large Ntrain.
                Xk = Xsrc
            else:
                Xk = Xsrc[:, cluster_indices]  # boolean-mask indexing always copies; safe to mutate in place
            Xk -= clusters[[k]].T  # in-place center

            Uk = torch.linalg.svd(Xk, full_matrices=False)[0]
            Phimat[:, :, k] = Uk[:, :r_store]
            del Xk

        # Global assimilation basis + local<->global maps (M = I in 1D); only meaningful for K>1.
        Phi_g = qbar_g = Tgk = dgk = Tkg = dkg = None
        if K > 1:
            qbar_g = Xtrain.mean(dim=1)
            Phi_g, Tgk, dgk, Tkg, dkg = build_global_assimilation_basis(
                Phimat, clusters, qbar_g, Mw=None)

        # the quantized-local model is just the compilation of the K single-cluster ROMs
        roms = [ROM(Phi=Phimat[:, :, k], centroid=clusters[[k]].T, **base_cfg)
                for k in range(K)]
        local_rom_full = qlGalerkin(roms, cfg=base_cfg)
        local_rom_full.Ntrain = Xtrain.shape[1]
        local_rom_full.Phi_g, local_rom_full.qbar_g = Phi_g, qbar_g
        local_rom_full.Tgk, local_rom_full.dgk = Tgk, dgk
        local_rom_full.Tkg, local_rom_full.dkg = Tkg, dkg
        local_rom_full.compute_transition_matrix(visited_clusters)

        torch.save(local_rom_full, local_model_path)
        print(f"Saved qlROM cache: K={K}, r_store={r_store}, file={'/'.join(local_model_path.split('/')[-3:])}.")
        return local_rom_full

    def _truncate_model(model_in: qlGalerkin, r_target: int) -> qlGalerkin:
        """Re-compile the model with every member ROM truncated to r_target modes."""
        if model_in.Phi_all is None:
            raise ValueError("Cached model has no Phi_all basis.")
        r_available = model_in.Phi_all.shape[1]
        if r_target > r_available:
            raise ValueError(
                f"Requested r={r_target} exceeds cached basis size r={r_available} for K={model_in.K}."
            )
        # Old caches predate the compilation API and lack base_cfg; rebuild it from
        # the KSConfig fields they do carry, so big cached bases stay usable.
        cfg = getattr(model_in, "base_cfg", None) or             {f.name: getattr(model_in, f.name) for f in fields(KSConfig)}
        roms = [ROM(Phi=model_in.Phi_all[:, :r_target, k].clone(),
                    centroid=model_in.centroids[[k]].T.clone(), **cfg)
                for k in range(model_in.K)]
        model_out = qlGalerkin(roms, cfg=cfg)
        model_out.Ntrain = getattr(model_in, "Ntrain", 0)
        model_out.T = getattr(model_in, "T", None)
        # Phi_g/qbar_g/dgk are independent of the local rank; the maps slice with the
        # local basis. Phi_g still contains the truncated local spaces, so the maps
        # stay exact inverses.
        model_out.Phi_g = getattr(model_in, "Phi_g", None)
        model_out.qbar_g = getattr(model_in, "qbar_g", None)
        if getattr(model_in, "Tgk", None) is not None:
            model_out.Tgk = model_in.Tgk[:, :, :r_target].clone()
            model_out.Tkg = model_in.Tkg[:, :r_target, :].clone()
            model_out.dgk = model_in.dgk.clone()
            model_out.dkg = model_in.dkg[:, :r_target].clone()
        return model_out

    # Load cached model if available and compatible, otherwise build a new model and cache it.
    if os.path.exists(local_model_path):
        try:
            cached_rom = _load_local_model(local_model_path)
            if cached_rom is not None:
                return _truncate_model(cached_rom, r)
        except Exception as e:
            print(f"Failed to load cached model (will rebuild): {e}")

    # Build new model
    cached_rom = _build_and_save_model()
    return _truncate_model(cached_rom, r)

get_simulation_path(model, Ntrain=None, Ntest=None, dt=None)

Cache directory of the models fitted on this case at this training length.

model is either the ks1d FOM or a torch-side KSConfig/ROM; both carry the case's PHYSICAL (visc, Lx, Nx, dt). Only Ntrain enters the path: a model is defined by what it was fitted on, and the test window is chosen afterwards, out of the same cached trajectory. Ntest is accepted and ignored so older call sites keep working.

Source code in qlroms/intrusive_qlroms/ks1d/rom.py
258
259
260
261
262
263
264
265
266
267
268
269
270
def get_simulation_path(model, Ntrain=None, Ntest=None, dt=None):
    """Cache directory of the models fitted on this case at this training length.

    `model` is either the ks1d FOM or a torch-side KSConfig/ROM; both carry the case's
    PHYSICAL (visc, Lx, Nx, dt). Only `Ntrain` enters the path: a model is defined by
    what it was fitted on, and the test window is chosen afterwards, out of the same
    cached trajectory. `Ntest` is accepted and ignored so older call sites keep working.
    """
    _dt = dt if dt is not None else model.dt
    base = os.path.join(str(data_dir() / "ks1d"), f"nu_{model.visc:.4f}_Lx{model.Lx:.4f}_Nx{model.Nx}", f"dt{_dt:.4f}")
    save_dir = os.path.join(base, f"Ntrain{Ntrain}") if Ntrain is not None else base
    os.makedirs(save_dir, exist_ok=True)
    return save_dir