How do I use Pressio to build and run a ROM?

This tutorial is designed to show the entire, end-to-end process of building and running a reduced-order model with Pressio.

We recommend reading this step-by-step explanation alongside the source code itself, which is fully commented to provide additional context. The source code for this tutorial can be found in $SRC/full-tutorial/burgers.cpp. It is also included at the bottom of this page for reference.

Each step found below (and in the Table of Contents to the right) corresponds to the same-named section of the source code. After Step 0, we recommend you jump to the main function and follow along with the steps as they are explained.

Note

You’ll notice that the Pressio examples, this one included, favor the auto keyword for type deduction. This is because the types in Pressio are often complex template types that would be cumbersome to write out in full. We strongly recommend using auto when working with Pressio to simplify your code and ensure you avoid mistakes with type declarations.

Step 0: Pressio Setup

Pressio is a header-only library that provides building blocks to create reduced-order models (ROMs) for large-scale dynamical systems.

Because it is header-only, we depend on user-set macros to enable or disable features at compile time. These macros may be set via CMake configuration or directly in the source code before including any Pressio headers.

A complete discussion of Pressio’s macros can be found in the pressio-rom documentation.

In this tutorial, we simply define the key macros in the source code, like so:

#define PRESSIO_ENABLE_LOGGING
#define PRESSIO_ENABLE_TPL_EIGEN

Once the macros have been defined, we can include the core Pressio header.

#include "pressio/rom.hpp"

Since rom is the top-level module of the Pressio ecosystem, including this file will pull in all of Pressio–from the built-in logger to the low-level ops to the ODE steppers. One #include to rule them all.

From here, we just include a couple other header files–one with basic helper functions for this tutorial, and another with utilities for performing hyper-reduction later on.

#include "helpers.h"
#include "hyperreduction.h"

Next, we’ll also include the standard C++ vector header to store our snapshots.

#include <vector>

And finally, we’ll define some convenient type aliases for vectors and matrices using the Eigen library, which Pressio supports out of the box.

using vector_t = Eigen::VectorXd;
using matrix_t = Eigen::MatrixXd;

Now, we can jump down to the main function to start the tutorial.

First, we initialize the Pressio logger so that we can see output from the library as we go:

PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::info);

You can change the LogLevel argument to debug for more information, or sparse for less.

Then we’ll just define some types for convenience:

using FomSystem = Burgers1dFom;
using time_t = typename FomSystem::time_type;

This FomSystem type is our implementation of the 1-D Burgers’ equation, which meets the Pressio FOM API as described below.

The purpose of defining this abstract FomSystem type is to indicate how Pressio is agnostic to the actual FOM implementation. As long as the FOM meets the required API, Pressio can work with it.

With that, all set-up is complete and we can move on to building and running the ROM.

Step 1: Define the Full-Order Model (FOM)

The first step in building a ROM is defining the full-order model (FOM) that we want to reduce.

In this tutorial, we will be working with the 1-D Burgers’ equation with periodic boundary conditions, as explained in the Background section.

Pressio expects a FOM to be defined such that it meets a specific API. There are a variety of different APIs depending on the type of problem being solved (e.g., steady vs unsteady, linear vs nonlinear, etc.). These APIs, defined as C++ concepts, can be found in the Pressio documentation.

We’ll use the API for a semi-discrete FOM defined by a system of ODEs:

\[\frac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u}, t; \mu)\]

where \(\mathbf{u}\) is the state vector, \(t\) is time, and \(\mu\) is a parameter.

According to the Pressio API for such a FOM, we need to define a class that meets the following API:

• Three core types: time_type, state_type, and rhs_type

• A method createRhs() that returns a new right-hand side vector

• A method rhs(const state_type & u, time_type t, rhs_type & f) that computes the right-hand side

Note

“RHS” here stands for “right-hand side,” and refers to the function \(\mathbf{f}(\mathbf{u}, t; \mu)\) in the ODE above.

Step 2: Run the FOM to generate snapshots

In order to generate a reduced-order model, we need to collect snapshots of the FOM solution over time. The runFom function in the source code performs this task.

Typically, applications would implement this step (and the prior definition of the FOM) themselves. In this tutorial, we provide a simple implementation of a Forward Euler time integrator that advances the FOM in time and stores snapshots of both the state and the right-hand side at each time step.

for ( int step = 0; step < numSteps; ++step ) {
   stateSnaps.push_back( u );   // Store state snapshot
   auto f = fom.createRhs();    // Create RHS vector
   fom.rhs( u, t, f );          // Compute RHS
   rhsSnaps.push_back( f );     // Store RHS snapshot
   u += dt * f;                 // Update solution
   t += dt;                     // Advance time
 }

We store these snapshots in a SnapshotSet object that will be used later to construct the reduced basis via Proper Orthogonal Decomposition (POD).

Step 3: Build the ROM from the snapshot matrix

Note

The burgers.cpp driver includes code for building both standard and hyper-reduced ROMs. For now, we’ll focus on the standard ROM, and leave a discussion of hyper-reduction for the Hyper-Reduction tutorial that follows this one.

Before we can use Pressio to construct the ROM, we need to find the reduced basis from the collected snapshots. This is done using POD, which identifies the most energetic modes in the snapshot data.

Snapshots to Trial Space

In this tutorial, we compute the singular value decomposition (SVD) of the snapshot matrix and select the first r=5 modes to form the reduced basis. The selection of the number of modes is arbitrary–you may experiment with how changing the value of r affects the results. Heuristically, we find that r=5 typically captures most of the system’s energy.

Once we have the basis, we can use Pressio’s create_trial_column_subspace function to create a trial subspace object that represents the reduced basis:

auto trialSpace = pressio::rom::create_trial_column_subspace<vector_t>(
   basis, affineShift, isAffine
);

Here, basis is a matrix whose columns are the selected POD modes, and affineShift is the mean of the snapshots. The isAffine argument indicates whether we are using an affine trial space.

This is implemented in the helpers.h file by the snapshots_to_trial_space function.

Trial Space to Reduced State

Next, we can define the reduced state by projecting the FOM’s initial condition onto the reduced basis. To do this, we use Pressio to compute the reduced state vector:

auto reducedState = trialSpace.createReducedState();

// reducedState = basis^T * (u0 - affineShift)
reducedState = trialSpace.basisOfTranslatedSpace().transpose() *
               (fomInitialState - trialSpace.translationVector());

This is implemented in the helpers.h file by the trial_space_to_reduced_state function.

Determine the Time Integration Scheme

Since the FOM is typically defined by an application independently of Pressio, we set up our own Forward Euler time integrator in Step 2. However, now that we’re using Pressio to build and run the ROM, we can use Pressio’s built-in time integrators.

The stepScheme defined below will be passed to the ROM constructor and used when we advance the ROM in time.

auto stepScheme = pressio::ode::StepScheme::ForwardEuler;

The step schemes (both implicit and explicit) that are supported by Pressio can be found in the pressio-rom/ode module here.

Build the ROM

With the trial space, reduced state, and time integrator defined, we can now construct the reduced-order model using Pressio’s built-in functions.

We’ll start with a default Galerkin ROM, which projects the FOM’s dynamics onto the trial space using a Galerkin projection. Later on, we’ll revisit this step to implement hyper-reduction.

auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem(
   stepScheme,
   trialSpace,
   fom
);

Step 4: Run the ROM

With the ROM constructed, we can now advance it in time using Pressio’s ode module.

Define the Time Stepping Policy

We create a stepping policy using the same time step size dt and number of steps numSteps as we used for the FOM simulation:

auto policy = pressio::ode::steps_fixed_dt(
    startTime,
    pressio::ode::StepCount( numSteps ),
    dt
);

There are other ways to define stepping policies in Pressio, such as using the start time, end time, and step size instead of the number of steps. These options are documented in the ode module.

Optional: Define the Observer

This step is optional, but it will allow us to monitor the ROM’s progress as it advances in time. Since we want to plot the results at the end, we’ll define an Observer that collects the ROM state at each time step.

Observers in Pressio are user-defined functors that are called at each time step during the time integration process. They can be used to log data, compute quantities of interest, or store states for later analysis.

Observers must meet the API defined by Pressio. Namely, they must implement the operator() method that takes the current time, step number, and state as arguments.

Our RomObserver will create the full state at every time step and store it in a trajectory member variable for later.

// Inside of RomObserver struct...
template <typename StepCountType, typename TimeType, typename StateType>
void operator()(StepCountType step, TimeType /*time*/, const StateType& reducedState) const
{
    // Capture every timestep
    auto fullState = trialSpace.createFullStateFromReducedState( reducedState );
    trajectory.push_back( fullState );
    ++stepCount;
}

The full RomObserver struct is implemented in the helpers.h file.

Advance the ROM in Time

Finally, we can advance the ROM in time using Pressio’s advance function:

pressio::ode::advance( rom, reducedState, policy, observer );

If you defined an observer, it will be called at each time step. Otherwise, simply omit the observer argument.

Reconstruct the Full-Order Solution

After the ROM has been advanced in time, the reducedState we defined earlier contains the final reduced state at the end of the simulation.

To get the solution in the full-order space, we need to reconstruct it using the trial space:

auto romSolution = trialSpace.createFullStateFromReducedState( reducedState );

Step 5: Compare ROM and FOM Solutions

At this point, we have constructed and run our reduced-order model. The remaining steps will compare the results to the FOM and plot the outcomes.

We already have the ROM solution at the final timestep, reconstructed back to the full-order space. Now, to get the FOM solution, we just extract the last snapshot from the FOM SnapshotSet we collected earlier. Specifically, we want the last column of the state snapshot matrix:

auto fomSolution = snapshots.stateSnapshots.col( snapshots.stateSnapshots.cols() - 1 );

With both the FOM and ROM solutions available, we can now compare them to assess the accuracy of the reduced-order model. We compute the L2 norm of the difference between the two solutions:

double error = ( fomSolution - romSolution ).norm() / fomSolution.norm();

This relative error gives us a measure of how well the ROM approximates the FOM solution. You should see something like:

[info] Relative error between FOM and default ROM: 0.0026138604529909602

In other words, the ROM solution is within about 0.26% of the FOM solution, which is quite accurate given that we only used 5 POD modes. Plus, for larger systems, the computational savings from using a ROM can be significant.

Step 6: Analyze the Results

Our Observer in Step 4 collected the ROM states at each time step and stored them in the trajectory vector. We can write these states to CSV files for further analysis and plotting.

After running the main executable:

cd $BUILD/full-tutorial
./burgers

we can generate plots with the output by running the python script in the same directory. It is hard-coded to find the output files in the output/ subdirectory.

cd $BUILD/full-tutorial
python plot.py

The resulting plots show how the ROM solutions compare to the FOM at various points in the simulation. You should see that the ROM closely tracks the FOM solution throughout the time integration.

The plot will be saved to your current directory as burgers_rom_comparison.png.

Note

You’ll find that the hyper-reduced ROM does not perform quite as well as the standard ROM. But the tradeoff is that the hyper-reduced ROM is computationally cheaper to run, especially for large-scale systems.

Go to top

Full Source File

///////////////////////////////////////////////////////////////////////////////
// Step 0: Pressio Setup
//////////////////////////////////////////////////////////////////////////////////

/**
 * Pressio uses macros to enable features like logging and TPLs.
 * Typically, these would be set during configuration, but we will
 * define them explicitly here for demonstration purposes.
 *
 * Importantly, these macros should be defined BEFORE including any Pressio
 * headers.
 */
#define PRESSIO_ENABLE_LOGGING
#define PRESSIO_ENABLE_TPL_EIGEN

/**
 * Due to the hierarchical structure of Pressio, including
 * pressio/rom.hpp will also pull in all necessary dependencies
 * (including Eigen, since we defined the macro above).
 * The logging macros (from pressio-log) are also included here.
 */
#include <pressio/rom.hpp>

/**
 * We'll also need some helper functions along the way.
 */
#include "helpers.h"

/**
 * And some functions for hyper-reduction.
 */
#include "hyperreduction.h"

/**
 * And we'll use std::vectors to store snapshots.
 */
#include <vector>

/**
 * Pressio supports various linear algebra backends (such as
 * Eigen, Tpetra, and Kokkos) and provides a unified interface
 * with the pressio::ops library for common operations. In this
 * tutorial, we will use Eigen.
 */
using vector_t = Eigen::VectorXd;
using matrix_t = Eigen::MatrixXd;

///////////////////////////////////////////////////////////////////////////////
// Step 1: Define the Full-Order Model (FOM)
///////////////////////////////////////////////////////////////////////////////

/**
 * Pressio requires a specific API to be used for the FOM class. This API can
 * vary depending on the type of problem. For the 1D Burgers' equation,
 * we will use the API for a semi-discrete FOM for time-dependent problems of
 * the form:
 *
 *     d/dt y(t) = f(y, t).
 *
 * The various FOM APIs can be found in pressio/rom/concepts.hpp.
 *
 * According to the semi-discrete API, our FOM class must define the following:
 *
 *     1. Three core types
 *         - time_type, state_type, rhs_type
 *
 *     2. A method to create the RHS vector
 *         - createRhs() -> rhs_type
 *
 *     3. A method to compute the RHS
 *         - rhs(const state_type & u, time_type t, rhs_type & f) -> void
 *
 * And that's it!
 */
class Burgers1dFom
{
    public:
        // API Requirement: Core types
        using time_type  = double;
        using state_type = vector_t;
        using rhs_type   = vector_t;

    private:
        std::size_t N_{ };
        double dx_{ };
        double nu_{ };

    public:
        Burgers1dFom( std::size_t N, double domainLength, double viscosity )
        : N_( N ), dx_( domainLength / ( N ) ), nu_( viscosity )
        {
            assert( N_  >= 3 );
            assert( dx_ > 0.0 );
            assert( nu_ > 0.0 );
        }

        std::size_t N() const { return N_; }
        double dx()     const { return dx_; }
        double nu()     const { return nu_; }

        // API Requirement: createRhs() method
        rhs_type createRhs() const
        {
            return rhs_type::Zero( N_ );
        }

        // API Requirement: rhs(...) method
        void rhs(
            const state_type & u,
            time_type /*t*/,
            rhs_type & f
        ) const
        {
            assert( u.size() == N_ );
            assert( f.size() == N_ );

            // periodic indexing helpers
            auto ip = [ this ]( std::size_t i ){ return ( i + 1 ) % N_; };
            auto im = [ this ]( std::size_t i ){ return ( i + N_ - 1 ) % N_; };

            const double inv2dx = 1.0 / ( 2.0 * dx_ );
            const double invdx2 = 1.0 / ( dx_ * dx_ );

            for ( std::size_t i = 0; i < N_; ++i ) {
                const std::size_t iL = im( i );
                const std::size_t iR = ip( i );

                // convective term: -(1/2 d/dx (u^2)) using centered flux difference
                // d/dx (u^2) ~ (u_{i+1}^2 - u_{i-1}^2) / (2 dx)
                const double dudt_conv =
                    -0.5 * ( ( u[ iR ] * u[ iR ] ) - ( u[ iL ] * u[ iL ] ) ) * inv2dx;

                // viscous term: nu * u_xx using second-order centered FD
                const double dudt_diff =
                    nu_ * ( u[ iR ] - 2.0 * u[ i ] + u[ iL ] ) * invdx2;
                f[ i ] = dudt_conv + dudt_diff;
            }
        }
};

///////////////////////////////////////////////////////////////////////////////
// Step 2: Run the FOM to generate snapshots (states and RHS)
////////////////////////////////////////////////////////////////////////////////

struct SnapshotSet
{
    matrix_t stateSnapshots;
    matrix_t rhsSnapshots;
};

/**
 * Building and running the FOM doesn't require any Pressio utilities.
 * For this example, we'll define a simple Forward Euler time integrator
 * to advance the FOM in time and collect the snapshots that will be used
 * to build the ROM with Pressio.
 */
template <typename FomSystem>
SnapshotSet runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt )
{
    // We will store each snapshot vector in a std::vector for now
    // For a typical ROM, we would only need the state snapshots. However,
    // for hyper-reduction, we will also need snapshots of the RHS.
    std::vector< vector_t > stateSnaps;
    std::vector< vector_t > rhsSnaps;

    // Compute the initial condition
    auto u = computeInitialCondition< FomSystem, vector_t >( fom );

    // Time integration loop (using Forward Euler)
    double t = startTime;
    for ( int step = 0; step < numSteps; ++step ) {
        stateSnaps.push_back( u );   // Store state snapshot
        auto f = fom.createRhs();    // Create RHS vector
        fom.rhs( u, t, f );          // Compute RHS
        rhsSnaps.push_back( f );     // Store RHS snapshot
        u += dt * f;                 // Update solution
        t += dt;                     // Advance time
    }

    // Convert our vector of snapshot vectors into matrices
    const std::size_t numSnapshots = stateSnaps.size();
    matrix_t stateMatrix( fom.N(), numSnapshots );
    matrix_t rhsMatrix( fom.N(), numSnapshots );
    for ( std::size_t i = 0; i < numSnapshots; ++i ) {
        stateMatrix.col( i ) = stateSnaps[ i ];
        rhsMatrix.col( i )      = rhsSnaps[ i ];
    }

    // Return both snapshot matrices
    return SnapshotSet{ std::move(stateMatrix), std::move(rhsMatrix) };
}

///////////////////////////////////////////////////////////////////////////////
// Step 3: Build the ROM from the snapshot matrix
///////////////////////////////////////////////////////////////////////////////

/**
 * Now we use the Pressio ecosystem to construct a ROM representation
 * with the snapshot matrix. As before, there are various APIs that we can meet
 * to use different types of ROMs. We'll use a basic Galerkin ROM here.
 */

template <typename FomSystem>
auto buildStandardGalerkinROM( FomSystem& fom, auto& stepScheme, auto& trialSpace )
{
    return pressio::rom::galerkin::create_unsteady_explicit_problem(
        stepScheme, trialSpace, fom
    );
}

/**
 * Hyper-reduced ROMs do not evaluate the full RHS at every time step.
 * Instead, they sample the RHS at selected indices and use a
 * projection to approximate the reduced RHS. This can greatly
 * reduce the computational cost of evaluating the ROM, especially
 * when the FOM is large and the RHS evaluation is expensive.
 *
 * However, there are trade-offs in accuracy, as you'll see.
 */
template <typename FomSystem>
auto buildHyperReducedGalerkinROM( FomSystem& fom, auto& stepScheme,
                                   auto& trialSpace, auto& hyperReducer )
{
    // Build the hyperreduced Galerkin ROM by passing the hyper-reducer
    // functor to the ROM factory function.
    return pressio::rom::galerkin::create_unsteady_explicit_problem(
        stepScheme, trialSpace, fom, hyperReducer
    );
}

///////////////////////////////////////////////////////////////////////////////
// Step 4: Run the ROM with trajectory capture
///////////////////////////////////////////////////////////////////////////////
/**
 * To run the ROM, we will define a simple observer that captures
 * the reduced state at each time step and reconstructs the full-order
 * state using the trial space. This allows us to capture the full
 * trajectory of the ROM solution for later analysis.
 */
template <typename FomSystem>
auto runROM( auto& rom, auto& trialSpace, auto& reducedState,
             typename FomSystem::time_type startTime, int numSteps, double dt,
             std::vector<vector_t>& trajectory
) {
    // Define the ODE time stepping policy
    auto policy = pressio::ode::steps_fixed_dt(
        startTime,
        pressio::ode::StepCount( numSteps ),
        dt
    );

    /**
     * In Pressio, observers are functors that are called at each time step
     * during time integration. Here, we define an observer that captures
     * the full-order state at each time step by reconstructing it from
     * the reduced state using the trial space.
     *
     * Observers must meet a specific API defined by Pressio. Namely, they
     * must implement the operator() with the signature shown in helpers.h.
     */
    auto observer = RomObserver< vector_t, decltype(trialSpace) >( trajectory, trialSpace );

    // Run the ROM time integration with observer
    pressio::ode::advance( rom, reducedState, policy, observer );

    /**
     * At this point, reducedState contains the ROM solution at the
     * final timestep. We just have to reconstruct the full-order
     * solution from the time reduced coordinates via the trial subspace.
     */
    auto romSolution = trialSpace.createFullStateFromReducedState( reducedState );

    return romSolution;
}

///////////////////////////////////////////////////////////////////////////////
// Step 5: Compare FOM and ROM solutions
///////////////////////////////////////////////////////////////////////////////

/**
 * To compare the FOM and ROM solutions, we can compute the relative error
 * between the two solution vectors at the final timestep.
 */
double compareFomAndRom( const matrix_t& fomSolution, const matrix_t& romSolution )
{
    double error = ( fomSolution - romSolution ).norm() / fomSolution.norm();
    return error;
}

///////////////////////////////////////////////////////////////////////////////
// Main driver
///////////////////////////////////////////////////////////////////////////////

int main() {
    // Initialize the logger so we can see Pressio output
    // Change the log level to "debug" for more information
    PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::info);

    using FomSystem = Burgers1dFom;
    using time_t = typename FomSystem::time_type;

    // 1. Create the FOM
    const std::size_t N = 100;          // Number of spatial points
    const double domainLength = 1.0;    // Domain length
    const double viscosity = 0.01;      // Viscosity parameter
    FomSystem fom( N, domainLength, viscosity );
    PRESSIOLOG_INFO( "1. Created FOM with N = {}, dx = {}, nu = {}", fom.N(), fom.dx(), fom.nu() );

    // 2. Run the FOM to generate a snapshot matrix
    const double dt    = 0.001;
    const int numSteps = 1000;
    time_t startTime   = 0.0;
    SnapshotSet snapshots = runFOM< FomSystem >( fom, startTime, numSteps, dt );
    PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshots.stateSnapshots.cols() );

    // 3 and 4. Build and run the ROM(s)

    /**
     * This tutorial covers both standard and hyper-reduced ROMs.
     * These variables below will hold the trajectories and solutions,
     * which we can analyze at the end.
     */
    std::vector<vector_t> standardRomTrajectory;
    std::vector<vector_t> hypredRomTrajectory;
    matrix_t standardSolution;
    matrix_t hypredSolution;

    // Start with the standard ROM
    PRESSIOLOG_INFO( "Standard ROM:" );
    {
        // 3. Build the ROM

        // Build the state trial space (Phi) from state snapshots
        auto standardTrialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom );
        // Create the initial reduced state by projecting the FOM initial condition
        auto standardReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( standardTrialSpace, fom );
        // Select the time integration scheme with Pressio
        auto stepScheme = pressio::ode::StepScheme::ForwardEuler;
        // Build the standard ROM
        auto standardRom = buildStandardGalerkinROM< FomSystem >( fom, stepScheme, standardTrialSpace );
        PRESSIOLOG_INFO( "  3. Built standard Galerkin ROM" );

        // 4. Run the ROM and capture the trajectory
        standardSolution = runROM< FomSystem >( standardRom, standardTrialSpace, standardReducedState, startTime, numSteps, dt, standardRomTrajectory );
        PRESSIOLOG_INFO( "  4. Ran standard ROM and captured trajectory" );
    }

    // Follow a similar process for hyper-reduced ROM
    PRESSIOLOG_INFO( "Hyper-reduced ROM:" );
    {
        // 3. Build the ROM
        auto hypredTrialSpace   = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom );
        auto hypredReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( hypredTrialSpace, fom );
        auto stepScheme = pressio::ode::StepScheme::ForwardEuler;
        /**
         * Here, we build the hyper-reducer using the RHS snapshot matrix
         * and the trial space. The hyper-reducer is a functor that will
         * be used by the ROM to compute the reduced RHS from sampled FOM RHSs.
         */
        auto hyperReducer = buildHyperReducer< vector_t, matrix_t >( snapshots.rhsSnapshots, hypredTrialSpace );
        auto hypredRom    = buildHyperReducedGalerkinROM< FomSystem >( fom, stepScheme, hypredTrialSpace, hyperReducer );
        PRESSIOLOG_INFO( "  3. Built hyper-reduced Galerkin ROM" );

        // 4. Run the ROM and capture the trajectory
        hypredSolution = runROM< FomSystem >( hypredRom, hypredTrialSpace, hypredReducedState, startTime, numSteps, dt, hypredRomTrajectory );
        PRESSIOLOG_INFO( "  4. Ran hyper-reduced ROM and captured trajectory" );
    }

    // 5. Compare ROM solution against FOM solution (the last snapshot)
    auto fomSolution = snapshots.stateSnapshots.col( snapshots.stateSnapshots.cols() - 1 );
    auto romError = compareFomAndRom( fomSolution, standardSolution );
    PRESSIOLOG_INFO( "5. Relative error between FOM and standard ROM: {}", romError );
    auto hypredError = compareFomAndRom( fomSolution, hypredSolution );
    PRESSIOLOG_INFO( "   Relative error between FOM and hyper-reduced ROM: {}", hypredError );

    // 6. Write trajectories to CSV files in output directory
    // Get FOM trajectory from snapshots
    std::vector<vector_t> fomTrajectory;
    for (int i = 0; i < snapshots.stateSnapshots.cols(); ++i) {
        fomTrajectory.push_back(snapshots.stateSnapshots.col(i));
    }
    writeTrajectoryToCSV< vector_t >("output", "fom_trajectory.csv",         fomTrajectory);
    writeTrajectoryToCSV< vector_t >("output", "standard_rom_trajectory.csv", standardRomTrajectory);
    writeTrajectoryToCSV< vector_t >("output", "hypred_rom_trajectory.csv",  hypredRomTrajectory);
    PRESSIOLOG_INFO( "6. Wrote trajectories to output/{fom,default_rom,hypred_rom}_trajectory.csv" );

    // Finalize the logger
    PRESSIOLOG_FINALIZE();
    return 0;
}

Go to top