[RomApp] ROMs on top of the Kratos Future strategies: architecture and open questions

[Rbravo555]

This document describes how the projection-based ROMs are organised on top of the
KratosMultiphysics.Future solving-strategy stack: which components exist, what each one is responsible
for, and how they talk to each other. The implementation is in #14850.

1. The idea in one paragraph

A projection-based ROM needs three things: a decoder u = D(q) that turns reduced coordinates into a
full state, a full order operator that builds the tangent matrix A(u) and the residual b(u), and a
projection that reduces them to a small system for dq. In the legacy code the three live inside one
BuilderAndSolver class per ROM variant, so every new decoder or projection means another copy of the
class. Here they are three separate components with narrow interfaces. The decoder is Python, the full
order operator is the Kratos core, and the projection is a small C++ class of this application.

2. Layers

flowchart TB
    subgraph L4["Workflow layer (Python)"]
        RM["RomManager<br/>training, basis, rom_data folder"]
        AN["FutureRomAnalysis<br/>wraps any AnalysisStage"]
    end
    subgraph L3["ROM solver layer (Python)"]
        SOL["FutureRomSolver<br/>Newton loop on q"]
        DEC["RomDecoder<br/>Decode(q), Jacobian(q)"]
    end
    subgraph L2["Projection layer (C++, RomApplication)"]
        PRJ["RomProjector<br/>build, project, solve, set solution"]
    end
    subgraph L1["Full order layer (C++, Kratos core, Future)"]
        SCH["StaticScheme + BlockBuilder"]
        DAT["ImplicitStrategyData<br/>linear system, DOF sets"]
    end
    subgraph L0["Model (C++, Kratos core)"]
        MP["ModelPart<br/>elements, conditions, DOFs"]
    end

    RM -. "rom_data: basis, node ids, settings" .-> AN
    AN --> SOL
    SOL --> DEC
    SOL --> PRJ
    PRJ --> SCH
    PRJ --> DAT
    SCH --> DAT
    SCH --> MP

Loading

Each layer only calls the one below it. Nothing in the core knows about ROMs, and the projector knows
nothing about decoders.

Layer	Component	File	Responsibility
Workflow	RomManager	python_scripts/rom_manager.py	Unchanged. Runs the FOM, computes the basis, writes rom_data.
Workflow	FutureRomAnalysis	python_scripts/future/future_rom_analysis.py	Reuses a standard analysis stage for the model, DOFs and processes, and replaces its solution step by the ROM solver. Builds the decoder from rom_data.
ROM solver	FutureRomSolver	python_scripts/future/rom_future_solver.py	Owns the scheme and the strategy data. Newton iteration on q and convergence check.
ROM solver	RomDecoder, LinearDecoder	python_scripts/future/rom_decoders.py	Decode(q) -> u and Jacobian(q) -> Phi.
Projection	RomProjector	future/rom_projector.h	Builds the constrained full order system, projects it, solves the reduced system, writes a state to the DOFs.
Projection	NumpyRomProjector	python_scripts/future/rom_future_solver.py	Same interface in NumPy/SciPy. Reference to validate the C++ one.
Full order	StaticScheme, BlockBuilder, ImplicitStrategyData	kratos/future/solving_strategies/	Core. DOF sets, sparse graph, assembly, constraints, update. Not modified.

3. The two interfaces that matter

Decoder (Python, one class per kind of ROM):

class RomDecoder:
    def NumberOfRomDofs(self): ...
    def Decode(self, q): ...      # u, one entry per effective DOF
    def Jacobian(self, q): ...    # Phi = dD/dq, dense, one row per effective DOF

Projector (C++ RomProjector, or its NumPy twin):

projector.BuildEffectiveSystem()   # A(u), b(u) with constraints applied
projector.Project(Phi)             # returns (Phi^T A Phi, Phi^T b)
projector.SolveReduced()           # returns dq
projector.SetSolution(u)           # writes u to the free DOFs through the scheme

The only convention shared by the two is the row ordering: row i belongs to the DOF with effective
equation id i. The ordering exists once the scheme has been initialised, which is why the solver receives
its decoder after Initialize().

4. One Newton iteration

sequenceDiagram
    participant S as FutureRomSolver
    participant D as RomDecoder
    participant P as RomProjector
    participant K as StaticScheme
    participant M as ModelPart

    S->>D: Decode(q)
    D-->>S: u
    S->>P: SetSolution(u)
    P->>K: Update (Dx = u - u_current)
    K->>M: write free DOFs, move mesh
    S->>P: BuildEffectiveSystem()
    P->>K: Build, BuildLinearSystemConstraints, ApplyLinearSystemConstraints
    K->>M: CalculateLocalSystem on elements and conditions
    S->>D: Jacobian(q)
    D-->>S: Phi
    S->>P: Project(Phi)
    P-->>S: A_r, b_r
    S->>P: SolveReduced()
    P-->>S: dq
    Note over S: q = q + dq, check |dq| / |q|

Loading

Per iteration Python makes one decoder evaluation and four calls into C++. The dense basis crosses the
boundary without being copied. All loops over elements, DOFs and matrix entries are in C++.

Design decisions visible in this loop:

• Absolute update. The state is always set as u = D(q), not incremented by Phi dq. For a linear
  decoder both are the same; for a nonlinear decoder only the absolute form is correct.
• The scheme does the update. SetSolution goes through StaticScheme::Update, so fixed DOFs are
  skipped, constraints are honoured and the mesh is moved by the same code as in a FOM.
• Fixed DOFs. Their rows of Phi are taken as zero inside Project, as the legacy
  GlobalROMBuilderAndSolver does. The decoder does not need to know the boundary conditions.

5. Legacy and Future side by side

flowchart LR
    subgraph Legacy
        direction TB
        LS["Legacy strategy<br/>Newton-Raphson"]
        LB["XxxROMBuilderAndSolver<br/>DOF set + build + projection<br/>+ reduced solve + decoder"]
        LS --> LB
    end
    subgraph Future
        direction TB
        FS["FutureRomSolver"]
        FD["RomDecoder"]
        FP["RomProjector"]
        FK["Future StaticScheme"]
        FS --> FD
        FS --> FP
        FP --> FK
    end

Loading

Concern	Legacy	Future
DOF set and equation ids	Re-implemented in each ROM builder and solver	ImplicitScheme::Initialize
Assembly	Build of each class	ImplicitScheme::Build
Projection	ProjectROM of each class	RomProjector::Project
Reduced solve	SolveROM of each class	RomProjector::SolveReduced
Basis storage	ROM_BASIS matrix on every node	One dense array owned by the decoder
Decoder	Hard-coded: linear in five classes, ANN in two more with a hand-written network and gradient	A Python object; a new decoder is a new Python class
Nonlinear loop and convergence	Legacy strategy and criteria	Python loop in FutureRomSolver
Classes to add a ROM variant	A new builder and solver (400 to 900 lines)	A decoder, or a projection method

6. Where each extension goes

Feature	Component that changes	Other components
ANN, RBF or quadratic decoder	New RomDecoder subclass, Jacobian by automatic differentiation	None
LSPG, Petrov-Galerkin	New projection methods in RomProjector	Solver selects the method
HROM (weights, selected elements)	A RomApplication scheme HRomScheme<TScheme> overriding the virtual ImplicitScheme::Build	None (no core change)
Dynamic problems	A dynamic scheme in the core Future	None here: the projector takes any ImplicitScheme
Other convergence criteria	FutureRomSolver	None
RomManager integration	Workflow layer: choose FutureRomAnalysis instead of RomAnalysis	None

7. What is implemented and checked

Implemented: Galerkin projection, linear decoder, static scheme, serial and shared-memory runs.

Case	Type	Outcome
Three-node Laplacian (C++ test)	Linear	Same reduced solution as the legacy test
Structural static test of the application	Linear, 2 steps	Differs from legacy by 3e-10
Thermal static test of the application	Nonlinear (radiation)	Equal to legacy after one iteration; the legacy ROM stops there
Stanford bunny, 26,679 DOFs, 6 modes	Large displacements, follower pressure, 30 steps	Differs from legacy by 1e-7; 11.6 s against 9.3 s
Flow past a cylinder, 9,216 DOFs, 200 modes	Transient Navier-Stokes, 500 steps	About 1 % from the FOM; differs from legacy by 0.5 %, cause not yet identified

Notes on the last two, which are run from scripts outside the repository:

• The bunny needs move_mesh on in the scheme settings, as the pressure follows the deformed surface.
• The cylinder only runs with BDF2 integrated by the element (element_manages_time_integration), because
  Future has no dynamic scheme. Bossak is not available.

8. Limits that come from the core

• No Newton-Raphson strategy or convergence criteria in Future; the nonlinear loop is Python.
• Only a static scheme.
• No MPI.
• No per-entity hook in ImplicitScheme::Build; the HROM can do without it (Build is virtual and
  CalculateLocalSystemContribution is protected), see question 3.
• No application solver uses Future yet. FutureRomAnalysis therefore still lets the parent analysis
  stage create its legacy solver to read the model and add the DOFs, and bypasses it in the solution steps.

9. Questions for discussion

1. Is the boundary right: decoder in Python, projector in C++?
2. Should the nonlinear loop stay in Python until the core has a Newton strategy, or should the
   RomApplication provide a C++ strategy that calls a decoder object?
3. HROM: the plan is a RomApplication scheme HRomScheme<TScheme> that overrides the virtual Build and assembles only the
   HROM entities, scaled by their weights (no core change). Would the core rather offer a hook for this (entity weights, a list
   of entities to assemble, or public local contribution methods)?
4. Should FutureRomAnalysis remain a wrapper around existing analysis stages, or wait for application
   solvers based on Future?

[Khalil-Work]

Thank you for this write-up — the layering is clean and the decoder/projector
split reads well. I work on verified reduced-order tooling (POD/DMD pipelines
checked against exactly known spectra), so I'll answer the four questions from
a verification-first angle.

1. Decoder in Python, projector in C++ — right boundary?
Yes. The decoder is the part that changes most often (linear, ANN, RBF,
quadratic…), and it is the part that benefits most from fast iteration and
automatic differentiation — Python is the right home for it. The projector is
the part whose correctness must be bit-trusted (Phi^T A Phi against the
effective-equation row ordering, constraint application, fixed-DOF zeroing),
and that is exactly the part you want in C++ with a NumPy twin. One
suggestion to harden the boundary: pin the row-ordering contract with a test
that permutes the effective equation ids and checks the projected system
permutes accordingly — that contract is currently implicit and it is the one
thing both sides must agree on.

2. Nonlinear loop in Python or a C++ strategy calling the decoder?
Keep the loop in Python until the core has its own Newton strategy, but make
the Python loop swappable from day one (a ConvergenceCriterion object now,
even if only one exists). The per-iteration cost is four C++ calls plus one
decoder evaluation — the decoder evaluation will dominate for nonlinear
decoders anyway, so moving the loop to C++ buys little until decoders are
also C++. What I would not postpone: recording the reduced residual history
per step, because for nonlinear cases the reduced residual is the primary
diagnostic (a slowly converging reduced Newton usually means the basis is
missing a direction, not that the loop is slow).

3. HROM via an application scheme overriding Build, or a core hook?
The HRomScheme<TScheme> override is the right first move — it keeps the
core untouched and the weights logic inside the application where the hyper-
reduction tooling lives. I would only promote a hook to the core once a
second consumer appears. A verification note from experience with element-
subset methods: add a test where the HROM weights are all exactly 1 and the
selected-entity list is the complete list — the HROM result must then match
the non-hyper-reduced ROM to roundoff. That single test catches the classic
class of weight/indexing bugs that otherwise only shows up as a mysterious
few-percent error on large cases.

4. FutureRomAnalysis as a wrapper, or wait for Future-based solvers?
Keep the wrapper. It gives you a migration path where the ROM solver is the
first real consumer of Future, which is exactly the kind of pressure test
the Future interfaces need before application solvers adopt them. Waiting
for application solvers would couple two timelines for no gain.

One more observation — the 0.5% cylinder gap. Since NumpyRomProjector
exists as the reference implementation, the fastest localization is to run
the cylinder case with the legacy projection and the Future projection on the
identical basis and identical snapshots, then compare the reduced operators
directly (Phi^T A Phi and Phi^T b at the first step) rather than the end
states. If the operators agree to roundoff, the gap is in the time loop or
the basis; if they differ, it is in the projection/constraint path — that
separates the two in one experiment.

Congrats again on the structure — the "decoder as a Python object, projector
with a NumPy twin" pattern is one I would happily steal.

[Rbravo555]

Thanks for the detailed answers. We will take up:

• a test that permutes the effective equation ids and checks that the projected system permutes accordingly;
• an HROM test with all weights equal to 1 on the complete entity list, against the ROM to roundoff: added to [RomApp] HROM for the ROMs on top of the Kratos Future strategies #14903 (1e-15). Also checked for the quadratic, ANN and RBF decoders (1e-11);
• the reduced residual history per step in the Python nonlinear loop;
• for the cylinder, comparing Phi^T A Phi and Phi^T b at the first step between the legacy and the Future projection on the same basis. Related: [RomApp] HROM training for the ROMs on top of the Kratos Future strategies #14953 (HROM training) tests that the residual of each element projected onto the basis matches the legacy one (1e-13).

[rubenzorrilla]

Thanks for the current status and perspectives wrap up. I'll go thought it in deep when I'm back. In the meantime, I'd like to say that I'd consider switching the RomManager to a SequentialOrchestrator, namely RomOrchestrator if you want.

[Rbravo555]

Agreed: a RomOrchestrator derived from SequentialOrchestrator replacing the RomManager workflow. I'll do it in a separate PR after #14850.