From bbbed05803a691ad538239162d7a378d42a20753 Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Wed, 14 Jan 2026 11:17:29 -0500 Subject: [PATCH 1/7] #63: write end-to-end demo of burgers' eq --- CMakeLists.txt | 3 + full-tutorial/CMakeLists.txt | 25 +++ full-tutorial/burgers.cpp | 372 +++++++++++++++++++++++++++++++++++ full-tutorial/helpers.h | 68 +++++++ 4 files changed, 468 insertions(+) create mode 100644 full-tutorial/CMakeLists.txt create mode 100644 full-tutorial/burgers.cpp create mode 100644 full-tutorial/helpers.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 770c555..653067b 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -62,3 +62,6 @@ set(CMAKE_CXX_EXTENSIONS OFF) enable_testing() add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/end-to-end-roms) + +# Build the standalone full tutorial (burgers example) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/full-tutorial) diff --git a/full-tutorial/CMakeLists.txt b/full-tutorial/CMakeLists.txt new file mode 100644 index 0000000..9752504 --- /dev/null +++ b/full-tutorial/CMakeLists.txt @@ -0,0 +1,25 @@ +cmake_minimum_required(VERSION 3.18) + +project(burgers_tutorial CXX) + +# Ensure C++17 (root also enforces this, but keep local for safety) +include(CheckCXXCompilerFlag) +check_cxx_compiler_flag("-std=c++17" COMPILER_SUPPORT_CPP17) +if(NOT COMPILER_SUPPORT_CPP17) + message(FATAL_ERROR "Compiler does not support -std=c++17. This is required.") +endif() +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +set(exename burgers) +add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/burgers.cpp) + +# Pull in the headers +target_include_directories(${exename} PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/eigen-3.4.0 + ${CMAKE_CURRENT_SOURCE_DIR}/../tpls # pressio_cmake_config.h + ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/pressio-rom/include + ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/pressio-ops/include + ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/pressio-log/include +) diff --git a/full-tutorial/burgers.cpp b/full-tutorial/burgers.cpp new file mode 100644 index 0000000..937a5e9 --- /dev/null +++ b/full-tutorial/burgers.cpp @@ -0,0 +1,372 @@ +/** + * End to end demo for solving 1D Burgers' equation using reduced-order + * modeling with Pressio. + * + * MATHEMATICAL FORMULATION + * ===================== + * + * The one-dimensional Burgers' equation is a classic nonlinear PDE + * combining advection and diffusion: + * + * ∂u/∂t + u * ∂u/∂x = ν * ∂²u/∂x² [1] + * + * where: + * - u(x, t) is the velocity field (scalar in 1D) + * - ν (nu) is the kinematic viscosity parameter + * - The domain is periodic: x ∈ [0, L] with u(0,t) = u(L,t) + * - Initial condition: u(x, 0) = u₀(x) + * + * PHYSICAL INTERPRETATION + * ======================= + * + * Burgers' equation models the balance between: + * + * 1. CONVECTIVE TRANSPORT (left side): u*∂u/∂x + * A fluid parcel carries its momentum forward. Steep gradients + * can steepen further (nonlinear shock formation). This term is + * NONLINEAR and generates fine-scale structures. + * + * 2. VISCOUS DIFFUSION (right side): ν*∂²u/∂x² + * Friction smooths velocity gradients, dissipating kinetic + * energy. High viscosity => smooth, diffuse solution. Low + * viscosity => sharp fronts, potential shock layers. + * + * The competing dynamics (steepening vs smoothing) produce rich + * behavior: traveling waves, shock-like structures, and energy + * cascade from large to small scales before dissipation. + * + * SPATIAL DISCRETIZATION + * ====================== + * + * We discretize on a uniform grid x_i = i*Δx, i = 0,...,N-1: + * + * du_i/dt = -1/2 * (u²_{i+1} - u²_{i-1}) / (2Δx) + * + ν * (u_{i+1} - 2u_i + u_{i-1}) / Δx² [2] + * + * Using centered finite differences: + * - Convective term: -(1/2) * d(u²)/dx via flux form + * - Diffusive term: ν * ∂²u/∂x² via standard 3-point stencil + * - Periodicity: indices wrap (i+1 mod N, i-1 mod N) + * + * This discretization produces a system of N ODEs: + * du/dt = f(u, t) [3] + * which we call the Full-Order Model (FOM). + * + * REDUCED-ORDER MODELING WITH PRESSIO + * ==================================== + * + * Motivation: The FOM [3] has N equations (thousands/millions + * in practical 2D/3D). We want to approximate u(t) ≈ ũ(t) + * using M << N basis vectors, reducing computational cost by + * orders of magnitude. + * + * POD Basis Construction: + * 1. Simulate the FOM for representative parameters/ICs, + * collecting M_snapshots state vectors: {u^(1), u^(2), ..., u^(p)} + * 2. Compute SVD: U ≈ U_r Σ_r V_r^T + * Columns of U_r are the r dominant POD modes (r << N) + * 3. Basis truncation: keep modes explaining 99.999% of energy + * + * Reduced State Representation: + * ũ(t) = U_r * a(t) + u_mean [4] + * + * where a(t) ∈ ℝ^r are reduced coordinates (r coefficients) + * and u_mean is the snapshot mean. + * + * Galerkin Projection (default ROM approach): + * Project equation [1] onto the r-dimensional subspace: + * + * da/dt = U_r^T * f(U_r * a + u_mean, t) [5] + * + * This reduces N equations to r equations. Since r << N, + * the ROM is dramatically cheaper than the FOM for + * evaluation and sensitivity analysis. + * + * ACCURACY & LIMITATIONS + * ====================== + * + * The ROM accuracy depends on: + * - Basis quality: POD basis should capture dominant dynamics + * - Extrapolation risk: ROM is trained on training parameter + * ranges; testing outside these ranges may be inaccurate + * - Nonlinear effects: Strong nonlinearity (Burgers convection) + * can cause Galerkin ROM errors to accumulate in time + * - Viscosity regime: Low-ν (high Reynolds number) regimes + * with shock-like structures are harder to capture with + * a fixed POD basis (manifold non-linearity) + * + * WORKFLOW IN THIS CODE + * ===================== + * + * This tutorial demonstrates: + * 1. Instantiate the FOM (class Burgers1dFom) + * 2. Run FOM simulations for training parameter values + * 3. Collect snapshots and compute POD basis + * 4. Construct Galerkin ROM [5] with Pressio + * 5. Run ROM on test parameters and compare against FOM + * + * Expected Result: ROM should deliver ~100-1000× speedup + * with minimal error (< 1%) on test cases within training range. + */ + +//////////////////////////////////////////////////////////////////////////////// + +/** + * 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 + +/** + * We'll also need some helper functions along the way. + */ +#include "helpers.h" + +/** + * 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 +//////////////////////////////////////////////////////////////////////////////// + +/** + * 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 +matrix_t runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt ) +{ + // We will store each snapshot vector in a std::vector for now + std::vector< vector_t > snapshots; + + // 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 ) { + snapshots.push_back( u ); // Store snapshot + auto f = fom.createRhs(); // Create RHS vector + fom.rhs( u, t, f ); // Compute RHS + u += dt * f; // Update solution + t += dt; // Advance time + } + + // Convert our vector of snapshot vectors into a matrix + const std::size_t numSnapshots = snapshots.size(); + matrix_t snapshotMatrix( fom.N(), numSnapshots ); + for ( std::size_t i = 0; i < numSnapshots; ++i ) { + snapshotMatrix.col( i ) = snapshots[ i ]; + } + + return snapshotMatrix; +} + +/////////////////////////////////////////////////////////////////////////////// +// Step 3: Build the ROM from the snapshot matrix +/////////////////////////////////////////////////////////////////////////////// + +/** + * Now we use the Pressio ecosystem to construct a Galerkin ROM representation + * with the snapshot matrix. As before, there are various APIs that we can meet + * to use different types of ROMs. We'll start with a default Galerkin here, + * but there are other options for features like hyper-reduction. + */ +template +auto buildAndRunROM( const matrix_t& snapshots, FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt ) +{ + // Use helper functions (helpers.h) to compute the trial subspace from snapshots + // and initialize the reduced state by projecting the FOM initial condition. + // Typically this would be done offline by the app using Pressio. + auto trialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots, fom ); + auto reducedState = trial_space_to_reduced_state< FomSystem, vector_t >( trialSpace, fom ); + + // Now we can set up the time integration scheme and policy with Pressio + auto stepScheme = pressio::ode::StepScheme::ForwardEuler; + auto policy = pressio::ode::steps_fixed_dt( + startTime, + pressio::ode::StepCount( numSteps ), + dt + ); + + // And finally build the default Galerkin ROM + auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom + ); + + // Then we use Pressio to run the ROM time integration + pressio::ode::advance( rom, reducedState, policy ); /* might need observer */ + + // At this point, reducedState contains the ROM solution at final time. + // We just have to reconstruct the full-order solution from the + // reduced coordinates via the trial subspace. + auto romSolution = trialSpace.createFullStateFromReducedState( reducedState ); + + return romSolution; +} + +/////////////////////////////////////////////////////////////////////////////// +// Step 4: 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; + matrix_t snapshotMatrix = runFOM< FomSystem >( fom, startTime, numSteps, dt ); + PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshotMatrix.cols() ); + + // 3. Build the ROM from the snapshot matrix and run it + auto romSolution = buildAndRunROM< FomSystem >( snapshotMatrix, fom, startTime, numSteps, dt ); + PRESSIOLOG_INFO( "3. Built and ran the ROM" ); + + // 4. Compare ROM solution against FOM solution (the last snapshot) + auto fomSolution = snapshotMatrix.col( snapshotMatrix.cols() - 1 ); + auto error = compareFomAndRom( fomSolution, romSolution ); + PRESSIOLOG_INFO( "4. Relative error between FOM and ROM: {}", error ); + + // Finalize the logger + PRESSIOLOG_FINALIZE(); + return 0; +} diff --git a/full-tutorial/helpers.h b/full-tutorial/helpers.h new file mode 100644 index 0000000..21f76ff --- /dev/null +++ b/full-tutorial/helpers.h @@ -0,0 +1,68 @@ +#include + +/** + * Determine the FOM initial condition + */ + +template +vector_t computeInitialCondition( const FomSystem& fom ) +{ + // Compuite the initial condition: e.g., u(x,0) = sin(2*pi*x/L) + vector_t u0( fom.N() ); + for ( std::size_t i = 0; i < fom.N(); ++i ) { + double x = i * fom.dx(); + u0( i ) = std::sin( 2.0 * M_PI * x / ( fom.dx() * fom.N() ) ); + } + return u0; +} + +/** + * Computes the POD basis from snapshots and constructs the trial + * subspace with the basis and affine shift. + */ +template +auto snapshots_to_trial_space( const matrix_t& snapshots, const FomSystemType& fom ) +{ + // Compute the SVD of the snapshot matrix + Eigen::JacobiSVD< matrix_t > svd( + snapshots, + Eigen::ComputeThinU | Eigen::ComputeThinV + ); + + // Extract the first r left singular vectors as the POD basis + int r = 5; + matrix_t basis = svd.matrixU().leftCols( r ).eval(); + + // Determine the affine shift (mean of snapshots) + vector_t affineShift = snapshots.rowwise().mean().eval(); + + // Construct the trial subspace + auto trialSpace = pressio::rom::create_trial_column_subspace( + basis, affineShift, false + ); + + return trialSpace; +} + +/** + * Projects the FOM initial condition onto the reduced basis + * to initialize the reduced state. + */ +template +auto trial_space_to_reduced_state( const TrialSpaceType& trialSpace, const FomSystemType& fom ) +{ + // Create the reduced state vector + auto reducedState = trialSpace.createReducedState(); + + // Reconstruct the same initial condition that the FOM used + // For Burgers: u(x, 0) = sin(2*pi*x/L) + auto u0 = computeInitialCondition< FomSystemType, vector_t >( fom ); + + // Project the initial condition onto the reduced basis: + // reducedState = basis^T * (u0 - affineShift) + vector_t centered = u0 - trialSpace.translationVector(); + reducedState = trialSpace.basisOfTranslatedSpace().transpose() * centered; + + // Return the reduced state + return reducedState; +} From 68609f7911fa32dabc4f373cdcf660f92f78a024 Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Tue, 27 Jan 2026 09:54:36 -0500 Subject: [PATCH 2/7] #63: add hyperreduction --- full-tutorial/CMakeLists.txt | 17 ++--- full-tutorial/burgers.cpp | 136 ++++++++++++++++++++++++--------- full-tutorial/helpers.h | 32 ++++++++ full-tutorial/hyperreduction.h | 76 ++++++++++++++++++ 4 files changed, 213 insertions(+), 48 deletions(-) create mode 100644 full-tutorial/hyperreduction.h diff --git a/full-tutorial/CMakeLists.txt b/full-tutorial/CMakeLists.txt index 9752504..4d63151 100644 --- a/full-tutorial/CMakeLists.txt +++ b/full-tutorial/CMakeLists.txt @@ -2,19 +2,16 @@ cmake_minimum_required(VERSION 3.18) project(burgers_tutorial CXX) -# Ensure C++17 (root also enforces this, but keep local for safety) -include(CheckCXXCompilerFlag) -check_cxx_compiler_flag("-std=c++17" COMPILER_SUPPORT_CPP17) -if(NOT COMPILER_SUPPORT_CPP17) - message(FATAL_ERROR "Compiler does not support -std=c++17. This is required.") -endif() -set(CMAKE_CXX_STANDARD 17) -set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_CXX_EXTENSIONS OFF) - set(exename burgers) add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/burgers.cpp) +# Require C++20 for this example (needed for auto in template parameters) +set_target_properties(${exename} PROPERTIES + CXX_STANDARD 20 + CXX_STANDARD_REQUIRED ON + CXX_EXTENSIONS OFF +) + # Pull in the headers target_include_directories(${exename} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/eigen-3.4.0 diff --git a/full-tutorial/burgers.cpp b/full-tutorial/burgers.cpp index 937a5e9..4ea3b04 100644 --- a/full-tutorial/burgers.cpp +++ b/full-tutorial/burgers.cpp @@ -103,7 +103,8 @@ * 2. Run FOM simulations for training parameter values * 3. Collect snapshots and compute POD basis * 4. Construct Galerkin ROM [5] with Pressio - * 5. Run ROM on test parameters and compare against FOM + * 5. (Optional) Implement hyper-reduction for efficiency + * 6. Run ROM on test parameters and compare against FOM * * Expected Result: ROM should deliver ~100-1000× speedup * with minimal error (< 1%) on test cases within training range. @@ -129,12 +130,18 @@ * The logging macros (from pressio-log) are also included here. */ #include +#include /** * We'll also need some helper functions along the way. */ #include "helpers.h" +/** + * And some functions for hyper-reduction. + */ +#include "hyperreduction.h" + /** * Pressio supports various linear algebra backends (such as * Eigen, Tpetra, and Kokkos) and provides a unified interface @@ -238,9 +245,15 @@ class Burgers1dFom }; /////////////////////////////////////////////////////////////////////////////// -// Step 2: Run the FOM to generate snapshots +// 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 @@ -248,10 +261,13 @@ class Burgers1dFom * to build the ROM with Pressio. */ template -matrix_t runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt ) +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 - std::vector< vector_t > snapshots; + // 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 ); @@ -259,21 +275,25 @@ matrix_t runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int nu // Time integration loop (using Forward Euler) double t = startTime; for ( int step = 0; step < numSteps; ++step ) { - snapshots.push_back( u ); // Store snapshot - auto f = fom.createRhs(); // Create RHS vector - fom.rhs( u, t, f ); // Compute RHS - u += dt * f; // Update solution - t += dt; // Advance time + 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 a matrix - const std::size_t numSnapshots = snapshots.size(); - matrix_t snapshotMatrix( fom.N(), numSnapshots ); + // 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 ) { - snapshotMatrix.col( i ) = snapshots[ i ]; + stateMatrix.col( i ) = stateSnaps[ i ]; + rhsMatrix.col( i ) = rhsSnaps[ i ]; } - return snapshotMatrix; + // Return both snapshot matrices + return SnapshotSet{ std::move(stateMatrix), std::move(rhsMatrix) }; } /////////////////////////////////////////////////////////////////////////////// @@ -281,21 +301,25 @@ matrix_t runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int nu /////////////////////////////////////////////////////////////////////////////// /** - * Now we use the Pressio ecosystem to construct a Galerkin ROM representation + * 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 start with a default Galerkin here, - * but there are other options for features like hyper-reduction. + * to use different types of ROMs. We'll use a Galerkin ROM here, with the option + * for hyper-reduction. Typically, you would choose which ROM is best suited for + * your application (as opposed to trying both as we do here for demonstration). */ template -auto buildAndRunROM( const matrix_t& snapshots, FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt ) +auto buildAndRunROM(const SnapshotSet& snaps, FomSystem& fom, + typename FomSystem::time_type startTime, int numSteps, double dt, + bool useHyperReduction = false ) { - // Use helper functions (helpers.h) to compute the trial subspace from snapshots - // and initialize the reduced state by projecting the FOM initial condition. - // Typically this would be done offline by the app using Pressio. - auto trialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots, fom ); + const auto& stateSnaps = snaps.stateSnapshots; + const auto& rhsSnaps = snaps.rhsSnapshots; + + // Build state trial space (Phi) and reduced state init + auto trialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( stateSnaps, fom ); auto reducedState = trial_space_to_reduced_state< FomSystem, vector_t >( trialSpace, fom ); - // Now we can set up the time integration scheme and policy with Pressio + // Now we can set up the time integration scheme and ODE policy with Pressio auto stepScheme = pressio::ode::StepScheme::ForwardEuler; auto policy = pressio::ode::steps_fixed_dt( startTime, @@ -303,13 +327,43 @@ auto buildAndRunROM( const matrix_t& snapshots, FomSystem& fom, typename FomSyst dt ); - // And finally build the default Galerkin ROM - auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( - stepScheme, trialSpace, fom - ); - - // Then we use Pressio to run the ROM time integration - pressio::ode::advance( rom, reducedState, policy ); /* might need observer */ + // Check if hyperreduction is enabled + if ( useHyperReduction ) + { + /** + * 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. + * + * 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 >( rhsSnaps, trialSpace ); + + // Build the hyperreduced Galerkin ROM by passing the hyper-reducer + // functor to the ROM factory function. + auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom, hyperReducer + ); + + // Run the hyperreduced ROM time integration + pressio::ode::advance( rom, reducedState, policy ); + } + else + { + // Build the default Galerkin ROM (no hyperreduction) + auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom + ); + + // Run the ROM time integration + pressio::ode::advance( rom, reducedState, policy ); + } // At this point, reducedState contains the ROM solution at final time. // We just have to reconstruct the full-order solution from the @@ -333,8 +387,10 @@ double compareFomAndRom( const matrix_t& fomSolution, const matrix_t& romSolutio return error; } - +/////////////////////////////////////////////////////////////////////////////// // Main driver +/////////////////////////////////////////////////////////////////////////////// + int main() { // Initialize the logger so we can see Pressio output // Change the log level to "debug" for more information @@ -354,17 +410,21 @@ int main() { const double dt = 0.001; const int numSteps = 1000; time_t startTime = 0.0; - matrix_t snapshotMatrix = runFOM< FomSystem >( fom, startTime, numSteps, dt ); - PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshotMatrix.cols() ); + SnapshotSet snapshots = runFOM< FomSystem >( fom, startTime, numSteps, dt ); + PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshots.stateSnapshots.cols() ); // 3. Build the ROM from the snapshot matrix and run it - auto romSolution = buildAndRunROM< FomSystem >( snapshotMatrix, fom, startTime, numSteps, dt ); - PRESSIOLOG_INFO( "3. Built and ran the ROM" ); + auto hypredSolution = buildAndRunROM< FomSystem >( snapshots, fom, startTime, numSteps, dt, true ); + auto romSolution = buildAndRunROM< FomSystem >( snapshots, fom, startTime, numSteps, dt ); + + PRESSIOLOG_INFO( "3. Built and ran the ROM (both default and hyper-reduced)" ); // 4. Compare ROM solution against FOM solution (the last snapshot) - auto fomSolution = snapshotMatrix.col( snapshotMatrix.cols() - 1 ); - auto error = compareFomAndRom( fomSolution, romSolution ); - PRESSIOLOG_INFO( "4. Relative error between FOM and ROM: {}", error ); + auto fomSolution = snapshots.stateSnapshots.col( snapshots.stateSnapshots.cols() - 1 ); + auto romError = compareFomAndRom( fomSolution, romSolution ); + PRESSIOLOG_INFO( "4. Relative error between FOM and default ROM: {}", romError ); + auto hypredError = compareFomAndRom( fomSolution, hypredSolution ); + PRESSIOLOG_INFO( " Relative error between FOM and hyper-reduced ROM: {}", hypredError ); // Finalize the logger PRESSIOLOG_FINALIZE(); diff --git a/full-tutorial/helpers.h b/full-tutorial/helpers.h index 21f76ff..b6ac146 100644 --- a/full-tutorial/helpers.h +++ b/full-tutorial/helpers.h @@ -66,3 +66,35 @@ auto trial_space_to_reduced_state( const TrialSpaceType& trialSpace, const FomSy // Return the reduced state return reducedState; } + +// Simple pseudoinverse using SVD +template +matrix_t pinv(const matrix_t& A, double tol = 1e-10) +{ + Eigen::JacobiSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); + const auto& S = svd.singularValues(); + const int r = static_cast(S.size()); + + // Build S^{-1} as r x r + matrix_t Sinv = matrix_t::Zero(r, r); + for (int i = 0; i < r; ++i) + if (S(i) > tol) Sinv(i,i) = 1.0 / S(i); + + // Using thin SVD: A = U_r S_r V_r^T => A^+ = V_r S_r^{-1} U_r^T + matrix_t V_r = svd.matrixV().leftCols(r); // n x r + matrix_t U_r = svd.matrixU().leftCols(r); // m x r + return (V_r * Sinv * U_r.transpose()).eval(); +} + +std::vector make_stride_samples(std::size_t N, std::size_t m) +{ + m = std::max(1, std::min(N, m)); + std::vector idx; + idx.reserve(m); + const double stride = static_cast(N) / static_cast(m); + for (std::size_t j = 0; j < m; ++j) { + int i = static_cast(std::floor(j * stride)) % static_cast(N); + idx.push_back(i); + } + return idx; +} diff --git a/full-tutorial/hyperreduction.h b/full-tutorial/hyperreduction.h new file mode 100644 index 0000000..215bb08 --- /dev/null +++ b/full-tutorial/hyperreduction.h @@ -0,0 +1,76 @@ +#include +#include +#include + +#include + +/** + * To build a hyperreduced ROM, we need to define a hyperreduction + * functor that will be used to project the sampled RHS to the + * reduced RHS. In this example, we will use an explicit Galerkin + * hyperreduction approach similar to the one described in + * + * "Hyper-reduction of nonlinear operators and the discrete empirical interpolation method + * for large nonlinear structural dynamics models" + * Chaturantabut and Sorensen, 2010, SIAM J. Sci. Comput. + * + * The hyperreducer functor will take as input the full RHS vector, + * sample it at selected indices, and then project it using a + * precomputed matrix H. + */ + +// Hyperreducer functor: projects sampled RHS to reduced RHS +template +class ExplicitGalerkinHyperReducer +{ + matrix_t H_; // r x m + std::vector samp_; // size m + +public: + ExplicitGalerkinHyperReducer(matrix_t H, std::vector samp) + : H_(std::move(H)), samp_(std::move(samp)) {} + + template + void operator()(const RhsType& rhsFull, const double&, ResultType& result) const + { + const int m = static_cast(samp_.size()); + vector_t sampled(m); + for (int j = 0; j < m; ++j) sampled(j) = rhsFull(samp_[j]); + result = H_ * sampled; + } +}; + +template +auto buildHyperReducer(const matrix_t& rhsSnaps, auto& trialSpace) { + // POD on RHS snapshots + Eigen::JacobiSVD svdRhs(rhsSnaps, Eigen::ComputeThinU | Eigen::ComputeThinV); + const int r = static_cast(trialSpace.basisOfTranslatedSpace().cols()); + int K = std::max(r+1, std::min(rhsSnaps.rows(), 3*r + 1)); + matrix_t Theta = svdRhs.matrixU().leftCols(K); + + // Sample selection + const int N = static_cast(rhsSnaps.rows()); + const int m = std::min(N, std::max(K, std::max(r+1, 20))); // at least r+1, K, 20 + auto samp = make_stride_samples(N, m); + + // Slice Theta at sample points -> m x K + matrix_t ThetaS(m, K); + for (int j = 0; j < m; ++j) ThetaS.row(j) = Theta.row(samp[j]); + + // Phi + matrix_t Phi = trialSpace.basisOfTranslatedSpace(); // N x r + + // Compute H = (Phi^T Theta) * pinv(ThetaS) + matrix_t cross = Phi.transpose() * Theta; // r x K + matrix_t ThetaS_pinv = pinv(ThetaS); + + PRESSIOLOG_DEBUG("hyper dims: r={} K={} m={} cross=({},{}) pinv=({},{})", + r, K, m, cross.rows(), cross.cols(), ThetaS_pinv.rows(), ThetaS_pinv.cols()); + + matrix_t H = cross * ThetaS_pinv; // r x m + + // Hyperreducer functor + ExplicitGalerkinHyperReducer hyperreducer(std::move(H), std::move(samp)); + + return hyperreducer; +} From bd97a1243d36951b7604964e91d6814f75450a74 Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Thu, 5 Feb 2026 09:18:23 -0500 Subject: [PATCH 3/7] #63: separate build and run; add plotting utility --- full-tutorial/CMakeLists.txt | 7 ++ full-tutorial/burgers.cpp | 167 +++++++++++++++++++++------------ full-tutorial/helpers.h | 68 ++++++++++++-- full-tutorial/hyperreduction.h | 25 ++++- full-tutorial/plot.py | 135 ++++++++++++++++++++++++++ 5 files changed, 330 insertions(+), 72 deletions(-) create mode 100644 full-tutorial/plot.py diff --git a/full-tutorial/CMakeLists.txt b/full-tutorial/CMakeLists.txt index 4d63151..41f1712 100644 --- a/full-tutorial/CMakeLists.txt +++ b/full-tutorial/CMakeLists.txt @@ -20,3 +20,10 @@ target_include_directories(${exename} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/pressio-ops/include ${CMAKE_CURRENT_SOURCE_DIR}/../tpls/pressio-log/include ) + +# Copy plot.py to the build directory +add_custom_command(TARGET ${exename} POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/plot.py + ${CMAKE_CURRENT_BINARY_DIR}/plot.py +) diff --git a/full-tutorial/burgers.cpp b/full-tutorial/burgers.cpp index 4ea3b04..f1b5f96 100644 --- a/full-tutorial/burgers.cpp +++ b/full-tutorial/burgers.cpp @@ -307,74 +307,82 @@ SnapshotSet runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int * for hyper-reduction. Typically, you would choose which ROM is best suited for * your application (as opposed to trying both as we do here for demonstration). */ + template -auto buildAndRunROM(const SnapshotSet& snaps, FomSystem& fom, - typename FomSystem::time_type startTime, int numSteps, double dt, - bool useHyperReduction = false ) +auto buildDefaultGalerkinROM( FomSystem& fom, auto& stepScheme, auto& trialSpace ) { - const auto& stateSnaps = snaps.stateSnapshots; - const auto& rhsSnaps = snaps.rhsSnapshots; + return pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom + ); +} - // Build state trial space (Phi) and reduced state init - auto trialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( stateSnaps, fom ); - auto reducedState = trial_space_to_reduced_state< FomSystem, vector_t >( 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 +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 + ); +} - // Now we can set up the time integration scheme and ODE policy with Pressio - auto stepScheme = pressio::ode::StepScheme::ForwardEuler; +/////////////////////////////////////////////////////////////////////////////// +// 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 +auto runROM( auto& rom, auto& trialSpace, auto& reducedState, + typename FomSystem::time_type startTime, int numSteps, double dt, + std::vector& trajectory +) { + // Define the ODE time stepping policy auto policy = pressio::ode::steps_fixed_dt( startTime, pressio::ode::StepCount( numSteps ), dt ); - // Check if hyperreduction is enabled - if ( useHyperReduction ) - { - /** - * 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. - * - * 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 >( rhsSnaps, trialSpace ); - - // Build the hyperreduced Galerkin ROM by passing the hyper-reducer - // functor to the ROM factory function. - auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( - stepScheme, trialSpace, fom, hyperReducer - ); - - // Run the hyperreduced ROM time integration - pressio::ode::advance( rom, reducedState, policy ); - } - else - { - // Build the default Galerkin ROM (no hyperreduction) - auto rom = pressio::rom::galerkin::create_unsteady_explicit_problem( - stepScheme, trialSpace, fom - ); - - // Run the ROM time integration - pressio::ode::advance( rom, reducedState, policy ); - } - - // At this point, reducedState contains the ROM solution at final time. - // We just have to reconstruct the full-order solution from the - // reduced coordinates via the trial subspace. + /** + * 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 4: Compare FOM and ROM solutions +// Step 5: Compare FOM and ROM solutions /////////////////////////////////////////////////////////////////////////////// /** @@ -413,19 +421,60 @@ int main() { SnapshotSet snapshots = runFOM< FomSystem >( fom, startTime, numSteps, dt ); PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshots.stateSnapshots.cols() ); - // 3. Build the ROM from the snapshot matrix and run it - auto hypredSolution = buildAndRunROM< FomSystem >( snapshots, fom, startTime, numSteps, dt, true ); - auto romSolution = buildAndRunROM< FomSystem >( snapshots, fom, startTime, numSteps, dt ); - - PRESSIOLOG_INFO( "3. Built and ran the ROM (both default and hyper-reduced)" ); + // 3. Build the ROMs from the snapshot matrix - // 4. Compare ROM solution against FOM solution (the last snapshot) + ////////////////////////////////////////////////////////////////////// + // Start with the default Galerkin ROM + ////////////////////////////////////////////////////////////////////// + // Build the state trial space (Phi) from state snapshots + auto defaultTrialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom ); + // Create the initial reduced state by projecting the FOM initial condition + auto defaultReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( defaultTrialSpace, fom ); + // Select the time integration scheme with Pressio + auto stepScheme = pressio::ode::StepScheme::ForwardEuler; + // Build the default ROM + auto defaultRom = buildDefaultGalerkinROM< FomSystem >( fom, stepScheme, defaultTrialSpace ); + + ////////////////////////////////////////////////////////////////////// + // Repeat the process for the hyper-reduced 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 ); + /** + * 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 both default and hyper-reduced ROMs" ); + + // 4. Run the ROMs and capture trajectories + std::vector defaultRomTrajectory; + std::vector hypredRomTrajectory; + auto romSolution = runROM< FomSystem >( defaultRom, defaultTrialSpace, defaultReducedState, startTime, numSteps, dt, defaultRomTrajectory ); + auto hypredSolution = runROM< FomSystem >( hypredRom, hypredTrialSpace, hypredReducedState, startTime, numSteps, dt, hypredRomTrajectory ); + PRESSIOLOG_INFO( "4. Ran the ROMs and captured trajectories" ); + + // 5. Compare ROM solution against FOM solution (the last snapshot) auto fomSolution = snapshots.stateSnapshots.col( snapshots.stateSnapshots.cols() - 1 ); auto romError = compareFomAndRom( fomSolution, romSolution ); - PRESSIOLOG_INFO( "4. Relative error between FOM and default ROM: {}", romError ); + PRESSIOLOG_INFO( "5. Relative error between FOM and default 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 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", "default_rom_trajectory.csv", defaultRomTrajectory); + 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; diff --git a/full-tutorial/helpers.h b/full-tutorial/helpers.h index b6ac146..be310a8 100644 --- a/full-tutorial/helpers.h +++ b/full-tutorial/helpers.h @@ -1,9 +1,37 @@ +#include +#include +#include + #include /** - * Determine the FOM initial condition + * Observer to capture ROM solutions at each timestep + * Signature: void operator()(StepCount, Time, State const&) const */ +template +struct RomObserver +{ + std::vector& trajectory; + const TrialSpaceType& trialSpace; + mutable int stepCount; + + RomObserver(std::vector& traj, const TrialSpaceType& ts) + : trajectory(traj), trialSpace(ts), stepCount(0) {} + // Observer call signature required by Pressio + template + void operator()(StepCountType step, TimeType /*time*/, const StateType& reducedState) const + { + // Capture every timestep + auto fullState = trialSpace.createFullStateFromReducedState(reducedState); + trajectory.push_back(fullState); + ++stepCount; + } +}; + +/** + * Determine the FOM initial condition + */ template vector_t computeInitialCondition( const FomSystem& fom ) { @@ -86,15 +114,35 @@ matrix_t pinv(const matrix_t& A, double tol = 1e-10) return (V_r * Sinv * U_r.transpose()).eval(); } -std::vector make_stride_samples(std::size_t N, std::size_t m) +/** + * Write solutions to CSV file in the specified output directory + */ +template +void writeTrajectoryToCSV( const std::string& outputDir, + const std::string& filename, + const std::vector< vector_t >& trajectory ) { - m = std::max(1, std::min(N, m)); - std::vector idx; - idx.reserve(m); - const double stride = static_cast(N) / static_cast(m); - for (std::size_t j = 0; j < m; ++j) { - int i = static_cast(std::floor(j * stride)) % static_cast(N); - idx.push_back(i); + // Ensure directory exists + std::filesystem::create_directories(outputDir); + + // Construct full file path + std::filesystem::path filepath = std::filesystem::path(outputDir) / filename; + + std::ofstream file(filepath); + for (int i = 0; i < trajectory[0].size(); ++i) { + file << "x" << i; + if (i < trajectory[0].size() - 1) file << ","; + } + file << "\n"; + + // Write data: each row is one timestep + for (const auto& state : trajectory) { + for (int i = 0; i < state.size(); ++i) { + file << std::scientific << std::setprecision(12) << state(i); + if (i < state.size() - 1) file << ","; + } + file << "\n"; } - return idx; + + file.close(); } diff --git a/full-tutorial/hyperreduction.h b/full-tutorial/hyperreduction.h index 215bb08..06eed69 100644 --- a/full-tutorial/hyperreduction.h +++ b/full-tutorial/hyperreduction.h @@ -40,6 +40,25 @@ class ExplicitGalerkinHyperReducer } }; +/** + * Typical hyperreducers would use a more refined sampling strategy, + * such as the Discrete Empirical Interpolation Method (DEIM) or + * Q-DEIM. Here, for simplicity, we use uniform striding to select + * sample indices. + */ +std::vector make_stride_samples(std::size_t N, std::size_t m) +{ + m = std::max(1, std::min(N, m)); + std::vector idx; + idx.reserve(m); + const double stride = static_cast(N) / static_cast(m); + for (std::size_t j = 0; j < m; ++j) { + int i = static_cast(std::floor(j * stride)) % static_cast(N); + idx.push_back(i); + } + return idx; +} + template auto buildHyperReducer(const matrix_t& rhsSnaps, auto& trialSpace) { // POD on RHS snapshots @@ -48,7 +67,7 @@ auto buildHyperReducer(const matrix_t& rhsSnaps, auto& trialSpace) { int K = std::max(r+1, std::min(rhsSnaps.rows(), 3*r + 1)); matrix_t Theta = svdRhs.matrixU().leftCols(K); - // Sample selection + // Sample selection (via striding here for simplicity) const int N = static_cast(rhsSnaps.rows()); const int m = std::min(N, std::max(K, std::max(r+1, 20))); // at least r+1, K, 20 auto samp = make_stride_samples(N, m); @@ -61,13 +80,13 @@ auto buildHyperReducer(const matrix_t& rhsSnaps, auto& trialSpace) { matrix_t Phi = trialSpace.basisOfTranslatedSpace(); // N x r // Compute H = (Phi^T Theta) * pinv(ThetaS) - matrix_t cross = Phi.transpose() * Theta; // r x K + matrix_t cross = Phi.transpose() * Theta; // r x K matrix_t ThetaS_pinv = pinv(ThetaS); PRESSIOLOG_DEBUG("hyper dims: r={} K={} m={} cross=({},{}) pinv=({},{})", r, K, m, cross.rows(), cross.cols(), ThetaS_pinv.rows(), ThetaS_pinv.cols()); - matrix_t H = cross * ThetaS_pinv; // r x m + matrix_t H = cross * ThetaS_pinv; // r x m // Hyperreducer functor ExplicitGalerkinHyperReducer hyperreducer(std::move(H), std::move(samp)); diff --git a/full-tutorial/plot.py b/full-tutorial/plot.py new file mode 100644 index 0000000..e1c3dfb --- /dev/null +++ b/full-tutorial/plot.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python3 +""" +Plotting script for Burgers' equation ROM comparison. + +This script reads the CSV trajectories generated by the C++ burgers simulation +and creates a visualization comparing FOM vs ROM solutions at four time snapshots. +""" + +import numpy as np +import matplotlib.pyplot as plt +import sys +import os + +def load_trajectory(filepath): + """Load trajectory from CSV file.""" + if not os.path.exists(filepath): + print(f"Warning: {filepath} not found") + return None + + data = np.genfromtxt(filepath, delimiter=',') + return data + +def main(): + # Define output directory + output_dir = "output" + + # Load trajectories from output directory + fom_traj = load_trajectory(os.path.join(output_dir, "fom_trajectory.csv")) + default_rom_traj = load_trajectory(os.path.join(output_dir, "default_rom_trajectory.csv")) + hypred_rom_traj = load_trajectory(os.path.join(output_dir, "hypred_rom_trajectory.csv")) + + if fom_traj is None: + print("Error: Could not load FOM trajectory from output directory") + sys.exit(1) + + # Select time indices for snapshots: early, mid, final + num_timesteps = fom_traj.shape[0] + time_indices = [ + num_timesteps // 3, # early + 2 * num_timesteps // 3, # mid + num_timesteps - 1 # final + ] + time_labels = ["Early (t≈T/3)", "Mid (t≈2T/3)", "Final (t=T)"] + + # Create figure with 1x3 subplots (one for each time snapshot) + fig, axes = plt.subplots(1, 3, figsize=(15, 5)) + fig.suptitle("1D Burgers' Equation: FOM vs ROM Solutions Over Time", fontsize=16, fontweight='bold') + axes = axes.flatten() + + # Spatial coordinate array + x_space = np.arange(fom_traj.shape[1]) + + # Plot each time snapshot + for plot_idx, (time_idx, time_label) in enumerate(zip(time_indices, time_labels)): + ax = axes[plot_idx] + + # Extract states at this time index + fom_state = fom_traj[time_idx, :] + + # Plot FOM solution + ax.plot(x_space, fom_state, 'b-', linewidth=2.5, label='FOM', marker='o', + markersize=3, markevery=8, zorder=3) + + # Plot default ROM if available and has this timestep + if default_rom_traj is not None and time_idx < default_rom_traj.shape[0]: + default_state = default_rom_traj[time_idx, :] + ax.plot(x_space, default_state, 'r--', linewidth=2, label='Default Galerkin ROM', + marker='s', markersize=3, markevery=8, zorder=2) + + # Compute error + default_error = np.linalg.norm(fom_state - default_state) / np.linalg.norm(fom_state) + else: + default_error = None + + # Plot hyperreduced ROM if available and has this timestep + if hypred_rom_traj is not None and time_idx < hypred_rom_traj.shape[0]: + hypred_state = hypred_rom_traj[time_idx, :] + ax.plot(x_space, hypred_state, 'g-.', linewidth=2, label='Hyper-Reduced ROM', + marker='^', markersize=3, markevery=8, zorder=1) + + # Compute error + hypred_error = np.linalg.norm(fom_state - hypred_state) / np.linalg.norm(fom_state) + else: + hypred_error = None + + # Format subplot + ax.set_xlabel('Spatial Index', fontsize=10) + ax.set_ylabel('Solution Value', fontsize=10) + ax.set_title(f'{time_label}', fontsize=11, fontweight='bold') + ax.legend(loc='best', fontsize=9) + ax.grid(True, alpha=0.3) + + # Add error text if available + error_text = "" + if default_error is not None: + error_text += f"ROM err: {default_error:.3e}\n" + if hypred_error is not None: + error_text += f"HR-ROM err: {hypred_error:.3e}" + + if error_text: + ax.text(0.98, 0.05, error_text, transform=ax.transAxes, + fontsize=8, verticalalignment='bottom', horizontalalignment='right', + bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) + + plt.tight_layout() + plt.savefig('burgers_rom_comparison.png', dpi=150, bbox_inches='tight') + print("✓ Saved plot to burgers_rom_comparison.png") + + # Print summary statistics + print("\n" + "="*60) + print("TRAJECTORY SUMMARY") + print("="*60) + print(f"FOM trajectory shape: {fom_traj.shape}") + if default_rom_traj is not None: + print(f"Default ROM trajectory shape: {default_rom_traj.shape}") + if hypred_rom_traj is not None: + print(f"Hyper-Reduced ROM trajectory shape: {hypred_rom_traj.shape}") + + # Compute errors at final time + fom_final = fom_traj[-1, :] + + if default_rom_traj is not None: + default_final = default_rom_traj[-1, :] + default_error = np.linalg.norm(fom_final - default_final) / np.linalg.norm(fom_final) + print(f"\nFinal-time error (Default ROM vs FOM): {default_error:.6e}") + + if hypred_rom_traj is not None: + hypred_final = hypred_rom_traj[-1, :] + hypred_error = np.linalg.norm(fom_final - hypred_final) / np.linalg.norm(fom_final) + print(f"Final-time error (Hyper-Reduced ROM vs FOM): {hypred_error:.6e}") + + print("="*60) + +if __name__ == "__main__": + main() From 8203118db49c6ee17892afac328a663ea958022c Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Thu, 5 Feb 2026 16:37:02 -0500 Subject: [PATCH 4/7] #63: add full write-up to docs --- docs/source/build.rst | 2 +- docs/source/index.rst | 12 +- docs/source/tutorial/background.rst | 25 ++ docs/source/tutorial/hyperreduction.rst | 116 +++++++ docs/source/tutorial/overview.rst | 24 ++ docs/source/tutorial/start.rst | 60 ++++ docs/source/tutorial/understand.rst | 405 ++++++++++++++++++++++++ full-tutorial/burgers.cpp | 14 +- 8 files changed, 648 insertions(+), 10 deletions(-) create mode 100644 docs/source/tutorial/background.rst create mode 100644 docs/source/tutorial/hyperreduction.rst create mode 100644 docs/source/tutorial/overview.rst create mode 100644 docs/source/tutorial/start.rst create mode 100644 docs/source/tutorial/understand.rst diff --git a/docs/source/build.rst b/docs/source/build.rst index 47af8b5..8a5be72 100644 --- a/docs/source/build.rst +++ b/docs/source/build.rst @@ -36,7 +36,7 @@ Verify the build .. code-block:: cpp - cd $BUIlDDIR + cd $BUILDDIR ctest Then what? diff --git a/docs/source/index.rst b/docs/source/index.rst index 75f2c9f..3e3697f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -39,19 +39,21 @@ Find us on `Slack `_. .. toctree:: :maxdepth: 1 - :hidden: build - Join our Slack - GitHub Repo - Open an issue/feature req. license +.. toctree:: + :maxdepth: 0 + :hidden: + :caption: 1. Full Tutorial: 1-D Burgers' Equation + + ./tutorial/overview .. toctree:: :maxdepth: 0 :hidden: - :caption: 1. End-to-end ROMs using pressio-demoapps + :caption: 2. End-to-end ROMs Using pressio-demoapps ./endtoend/readthisfirst ./endtoend/swe_galerkin_default diff --git a/docs/source/tutorial/background.rst b/docs/source/tutorial/background.rst new file mode 100644 index 0000000..d861256 --- /dev/null +++ b/docs/source/tutorial/background.rst @@ -0,0 +1,25 @@ +Background +========== + +1-D Burgers' Equation +---------------------- + +The 1-D Burgers' equation is a fundamental partial differential equation +that models various physical phenomena, including fluid dynamics and traffic flow. +It is given by: + +.. math:: + + \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2} + +where :math:`u` is the velocity field, :math:`t` is time, :math:`x` is the spatial coordinate, +and :math:`\nu` is the viscosity coefficient. + +While a more comprehensive introduction to the Burgers' equation is beyond the scope of this tutorial, +numerical methods for solving it typically involve discretizing the spatial domain +using finite difference or finite volume methods, and then integrating the resulting +system of ordinary differential equations (ODEs) in time using methods such as +Forward Euler, Runge-Kutta, or implicit schemes. + +In this tutorial, we will focus on setting up a reduced-order model (ROM) +for the 1-D Burgers' equation using the Pressio library. diff --git a/docs/source/tutorial/hyperreduction.rst b/docs/source/tutorial/hyperreduction.rst new file mode 100644 index 0000000..c9eed4b --- /dev/null +++ b/docs/source/tutorial/hyperreduction.rst @@ -0,0 +1,116 @@ +Hyper-Reduction +=============== + +This tutorial builds off of the information presented in +:doc:`understand` to demonstrate how to set up and run a hyper-reduced +reduced-order model (ROM) for the 1-D Burgers' equation using the Pressio library. + +If you have not completed that initial tutorial, please do so before proceeding. + +This page will follow the same structure as the previous tutorial, +going through the enumerated steps but only highlighting the differences +needed to implement hyper-reduction. + +Step 2: Run the FOM to generate snapshots +----------------------------------------- + +In the original tutorial, we defined a ``SnapshotSet`` +to collect snapshots for both the ``state`` and the +right-hand side ``rhs`` of the full-order model (FOM). + +You may have noticed that the default ROM did not use +the RHS snapshots at all. However, for hyper-reduction, +we will need these RHS snapshots to build the +hyper-reduction basis. + +We will explain this further in subsequent steps. + +If you are not implementing hyper-reduction, there is no +need to store the RHS snapshots. You can simply return +the state snapshots directly. + +Step 3: Build the ROM from the snapshot matrix +---------------------------------------------- + +A typical ROM uses only the state snapshots to build +the reduced basis via Proper Orthogonal Decomposition (POD). + +This reduces the dimensionality of the system, but +we still evaluate the full-order RHS at each time step, +which can be computationally expensive. **Hyper-reduction** +addresses this issue by approximating the RHS using a +reduced basis constructed from the RHS snapshots. + +So instead of a reduced state basis alone, we will +also build a reduced basis for the RHS using the +RHS snapshots we collected earlier. + +The process begins identically to before: we use the +state snapshots to compute the trial space, and then use +that to get the reduced state. We will also use the same +``ForwardEuler`` time stepper as before. + +The key difference is in the definition of a **hyper-reducer**. + +The hyper-reducer is a functor that takes in a sampled +FOM RHS vector and projects it into a reduced space. + +The process is similar to that of the trial space: we +compute a POD basis from the RHS snapshots using SVD, +and then use that basis to project the FOM RHS onto +a lower-dimensional subspace. + +There are refined techniques for selecting which rows +of the RHS to sample, such as the Discrete Empirical +Interpolation Method (DEIM). For simplicity in this tutorial, +we will use a basic sampling approach that selects rows at +a uniform interval (or stride). + +Here is a general outline of how the hyper-reducer functor +is defined: + +.. code-block:: cpp + + auto svdRhs = /* compute SVD of RHS snapshot matrix */; + auto rhsBasis = /* extract leading left singular vectors from SVD */; + + auto sampleIndices = /* select sample indices uniformly */; + auto sampledBasis = /* extract rows of rhsBasis at sampleIndices */; + + auto trialBasis /* aka Phi */ = trialSpace.basisOfTranslatedSpace(); + auto hypredMatrix = /* (Phi^T rhsBasis) * pinv(sampledBasis) */; + + // Construct hyperreducer functor + ExplicitGalerkinHyperReducer hyperReducer(hypredMatrix, sampleIndices); + +This code is implemented fully in the ``hyperreduction.h`` file +by the ``buildHyperReducer`` function. + +The ``ExplicitGalerkinHyperReducer`` above is a +user-defined functor that implements an ``operator()`` method +that takes in a sampled FOM RHS vector and outputs the reduced RHS vector. + +Like with the other Pressio interfaces, you can define your +hyperreducer in any way you like, as long as it meets the API requirements. + +Now that we have defined the hyper-reducer, we can +build the ROM almost identically to before, but passing +the hyperreducer functor as an additional argument: + +.. code-block:: cpp + + auto rom = pressio::rom::galerkin::create_default_problem_with_hyperreduction( + stepScheme, + trialSpace, + fom, + hyperReducer + ); + +Step 4-6: As Before +------------------- + +The rest of the steps are identical regardless of whether +you are using hyper-reduction or not. You can refer to the +:doc:`understand` page for details on how to run the ROM, +reconstruct the full-order solution, compare to the FOM, +and examine the output files. diff --git a/docs/source/tutorial/overview.rst b/docs/source/tutorial/overview.rst new file mode 100644 index 0000000..8410b5a --- /dev/null +++ b/docs/source/tutorial/overview.rst @@ -0,0 +1,24 @@ +End-to-end Tutorial: 1-D Burgers' Equation +========================================== + +This tutorial demonstrates how to set up and run a reduced order model (ROM) +for the 1-D Burgers' equation using the Pressio library. The tutorial covers +the following key steps: + +1. Setting up the full order model (FOM) and generating snapshots. +2. Constructing the reduced basis using Proper Orthogonal Decomposition (POD). +3. Building both default and hyper-reduced ROMs. +4. Running the ROMs and capturing trajectories. +5. Comparing the ROM solutions with the FOM solution. +6. Writing trajectories to CSV files for further analysis. + +The tutorial is implemented in C++ and leverages Pressio's capabilities for +model reduction and time integration. + +.. toctree:: + :maxdepth: 1 + + background + start + understand + hyperreduction diff --git a/docs/source/tutorial/start.rst b/docs/source/tutorial/start.rst new file mode 100644 index 0000000..d2336c3 --- /dev/null +++ b/docs/source/tutorial/start.rst @@ -0,0 +1,60 @@ +Quick Start +=========== + +These instructions assume the following env variables +are set: + +.. code-block:: bash + + export SRC= + export BUILD= + +The main driver for the tutorial is defined in +``$SRC/full-tutorial/burgers.cpp``. +It is fully commented to help you understand each +step of the implementation. + + +Additionally, a full discussion of the tutorial can be found +in the :doc:`understand` page. + +Build +----- + + +The full tutorial will build with the rest of the repository. +See :doc:`../build` for instructions. Be sure to install the Python dependencies +listed in ``py_requirements.txt`` in order to use the plotting +utility at the end of the tutorial. + +Once you have built the repository, the tutorial will be in +``$BUILD/full-tutorial/burgers``. + +.. note:: + + This tutorial requires C++20 support in your compiler. + +Run +--- + +Run the tutorial with: + +.. code-block:: bash + + cd $BUILD/full-tutorial + ./burgers + +This will run the tutorial and generate output files +in the ``output/`` subdirectory. + +You 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. + +.. code-block:: bash + + cd $BUILD/full-tutorial + python plot.py + +For a detailed breakdown of what you just ran, see +the :doc:`understand` page. \ No newline at end of file diff --git a/docs/source/tutorial/understand.rst b/docs/source/tutorial/understand.rst new file mode 100644 index 0000000..6702e40 --- /dev/null +++ b/docs/source/tutorial/understand.rst @@ -0,0 +1,405 @@ +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`` + +Each step found below 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 depends 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 dicussion 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: + +.. code-block:: cpp + + #define PRESSIO_ENABLE_LOGGING + #define PRESSIO_ENABLE_TPL_EIGEN + +Once the macros have been defined, we can include the core Pressio +header. + +.. code-block:: cpp + + #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. + +.. code-block:: cpp + + #include "helpers.h" + #include "hyperreduction.h" + +Next, we'll also include the standard C++ vector header +to store our snapshots. + +.. code-block:: cpp + + #include + +And finally, we'll define some convenient type aliases for vectors and matrices +using the Eigen library, which Pressio supports out of the box. + +.. code-block:: cpp + + 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: + +.. code-block:: cpp + + pressio::log::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: + +.. code-block:: cpp + + 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 :doc:`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: + +.. math:: + + \frac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u}, t; \mu) + +where :math:`\mathbf{u}` is the state vector, :math:`t` is time, and :math:`\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 + +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. + +.. code-block:: cpp + + 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 +---------------------------------------------- + +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: + +.. code-block:: cpp + + auto trialSpace = pressio::rom::create_trial_column_subspace( + 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: + +.. code-block:: cpp + + 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. + +.. code-block:: cpp + + 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. + +.. code-block:: cpp + + 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: + +.. code-block:: cpp + + 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. + +.. code-block:: cpp + + // Inside of RomObserver struct... + template + 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: + +.. code-block:: cpp + + 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: + +.. code-block:: cpp + + 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: + +.. code-block:: cpp + + 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: + +.. code-block:: cpp + + 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: + +.. code-block:: console + + [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: + +.. code-block:: bash + + 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. + +.. code-block:: bash + + 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. diff --git a/full-tutorial/burgers.cpp b/full-tutorial/burgers.cpp index f1b5f96..b7ed347 100644 --- a/full-tutorial/burgers.cpp +++ b/full-tutorial/burgers.cpp @@ -112,6 +112,10 @@ //////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// Step 0: Pressio Setup +////////////////////////////////////////////////////////////////////////////////// + /** * Pressio uses macros to enable features like logging and TPLs. * Typically, these would be set during configuration, but we will @@ -130,7 +134,6 @@ * The logging macros (from pressio-log) are also included here. */ #include -#include /** * We'll also need some helper functions along the way. @@ -142,6 +145,11 @@ */ #include "hyperreduction.h" +/** + * And we'll use std::vectors to store snapshots. + */ +#include + /** * Pressio supports various linear algebra backends (such as * Eigen, Tpetra, and Kokkos) and provides a unified interface @@ -303,9 +311,7 @@ SnapshotSet runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int /** * 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 Galerkin ROM here, with the option - * for hyper-reduction. Typically, you would choose which ROM is best suited for - * your application (as opposed to trying both as we do here for demonstration). + * to use different types of ROMs. We'll use a basic Galerkin ROM here. */ template From ff05267d277eef97d773840908d63c4ebfdc45cb Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Fri, 6 Feb 2026 10:08:28 -0500 Subject: [PATCH 5/7] #63: small fixes --- docs/source/tutorial/overview.rst | 6 +++--- docs/source/tutorial/understand.rst | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/source/tutorial/overview.rst b/docs/source/tutorial/overview.rst index 8410b5a..d3ad089 100644 --- a/docs/source/tutorial/overview.rst +++ b/docs/source/tutorial/overview.rst @@ -5,9 +5,9 @@ This tutorial demonstrates how to set up and run a reduced order model (ROM) for the 1-D Burgers' equation using the Pressio library. The tutorial covers the following key steps: -1. Setting up the full order model (FOM) and generating snapshots. -2. Constructing the reduced basis using Proper Orthogonal Decomposition (POD). -3. Building both default and hyper-reduced ROMs. +1. Setting up the full order model (FOM). +2. Running the FOM to generate a snapshot matrix. +3. Building both basic and hyper-reduced ROMs. 4. Running the ROMs and capturing trajectories. 5. Comparing the ROM solutions with the FOM solution. 6. Writing trajectories to CSV files for further analysis. diff --git a/docs/source/tutorial/understand.rst b/docs/source/tutorial/understand.rst index 6702e40..52a0b93 100644 --- a/docs/source/tutorial/understand.rst +++ b/docs/source/tutorial/understand.rst @@ -31,7 +31,7 @@ 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 depends on user-set macros to +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. @@ -91,7 +91,7 @@ output from the library as we go: .. code-block:: cpp - pressio::log::initialize(pressiolog::LogLevel::info); + PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::info); You can change the ``LogLevel`` argument to ``debug`` for more information, or ``sparse`` for less. From 41651b5289858e84969c7b1c7449e5f140f3dfe9 Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Mon, 9 Feb 2026 12:08:02 -0500 Subject: [PATCH 6/7] #63: clarify hyper-reduction workflow in tutorial --- docs/source/tutorial/background.rst | 4 ++ docs/source/tutorial/hyperreduction.rst | 37 ++++++---- docs/source/tutorial/understand.rst | 14 +++- full-tutorial/burgers.cpp | 89 +++++++++++++++---------- 4 files changed, 94 insertions(+), 50 deletions(-) diff --git a/docs/source/tutorial/background.rst b/docs/source/tutorial/background.rst index d861256..f5d5136 100644 --- a/docs/source/tutorial/background.rst +++ b/docs/source/tutorial/background.rst @@ -1,6 +1,10 @@ Background ========== +Reduced-order models (ROMs) can be used to simulate a variety of physical phenomena +modeled by partial differential equations (PDEs). In this tutorial, we will focus +on the 1-D Burgers' equation as a representative example. + 1-D Burgers' Equation ---------------------- diff --git a/docs/source/tutorial/hyperreduction.rst b/docs/source/tutorial/hyperreduction.rst index c9eed4b..4d06010 100644 --- a/docs/source/tutorial/hyperreduction.rst +++ b/docs/source/tutorial/hyperreduction.rst @@ -11,6 +11,30 @@ This page will follow the same structure as the previous tutorial, going through the enumerated steps but only highlighting the differences needed to implement hyper-reduction. +First: Why use hyper-reduction? +------------------------------- + +As discussed above, a typical ROM uses only the state snapshots +to build the reduced basis. + +This reduces the dimensionality of the system, but +we still evaluate the full-order RHS at each time step, +which can be computationally expensive. **Hyper-reduction** +addresses this issue by approximating the RHS using a +reduced basis constructed from the RHS snapshots. + +So instead of a reduced state basis alone, we will +also build a reduced basis for the RHS using the +RHS snapshots we collected earlier. + +With this, we can begin the tutorial. + +Steps 0-1: As Before +-------------------- + +Setting up Pressio and constructing your FOM is done exactly +the same, regardless of the type of ROM you will build. + Step 2: Run the FOM to generate snapshots ----------------------------------------- @@ -32,19 +56,6 @@ the state snapshots directly. Step 3: Build the ROM from the snapshot matrix ---------------------------------------------- -A typical ROM uses only the state snapshots to build -the reduced basis via Proper Orthogonal Decomposition (POD). - -This reduces the dimensionality of the system, but -we still evaluate the full-order RHS at each time step, -which can be computationally expensive. **Hyper-reduction** -addresses this issue by approximating the RHS using a -reduced basis constructed from the RHS snapshots. - -So instead of a reduced state basis alone, we will -also build a reduced basis for the RHS using the -RHS snapshots we collected earlier. - The process begins identically to before: we use the state snapshots to compute the trial space, and then use that to get the reduced state. We will also use the same diff --git a/docs/source/tutorial/understand.rst b/docs/source/tutorial/understand.rst index 52a0b93..4664494 100644 --- a/docs/source/tutorial/understand.rst +++ b/docs/source/tutorial/understand.rst @@ -36,7 +36,7 @@ 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 dicussion of Pressio's macros can be found in the +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, @@ -144,6 +144,11 @@ that meets the following API: - 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 + :math:`\mathbf{f}(\mathbf{u}, t; \mu)` in the ODE above. + Step 2: Run the FOM to generate snapshots ----------------------------------------- @@ -173,6 +178,13 @@ 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 :doc:`hyperreduction` 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. diff --git a/full-tutorial/burgers.cpp b/full-tutorial/burgers.cpp index b7ed347..8797829 100644 --- a/full-tutorial/burgers.cpp +++ b/full-tutorial/burgers.cpp @@ -315,7 +315,7 @@ SnapshotSet runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int */ template -auto buildDefaultGalerkinROM( FomSystem& fom, auto& stepScheme, auto& trialSpace ) +auto buildStandardGalerkinROM( FomSystem& fom, auto& stepScheme, auto& trialSpace ) { return pressio::rom::galerkin::create_unsteady_explicit_problem( stepScheme, trialSpace, fom @@ -427,46 +427,63 @@ int main() { SnapshotSet snapshots = runFOM< FomSystem >( fom, startTime, numSteps, dt ); PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshots.stateSnapshots.cols() ); - // 3. Build the ROMs from the snapshot matrix - - ////////////////////////////////////////////////////////////////////// - // Start with the default Galerkin ROM - ////////////////////////////////////////////////////////////////////// - // Build the state trial space (Phi) from state snapshots - auto defaultTrialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom ); - // Create the initial reduced state by projecting the FOM initial condition - auto defaultReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( defaultTrialSpace, fom ); - // Select the time integration scheme with Pressio - auto stepScheme = pressio::ode::StepScheme::ForwardEuler; - // Build the default ROM - auto defaultRom = buildDefaultGalerkinROM< FomSystem >( fom, stepScheme, defaultTrialSpace ); - - ////////////////////////////////////////////////////////////////////// - // Repeat the process for the hyper-reduced 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 ); + // 3 and 4. Build and run the ROM(s) + /** - * 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. + * 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. */ - auto hyperReducer = buildHyperReducer< vector_t, matrix_t >( snapshots.rhsSnapshots, hypredTrialSpace ); - auto hypredRom = buildHyperReducedGalerkinROM< FomSystem >( fom, stepScheme, hypredTrialSpace, hyperReducer ); - - PRESSIOLOG_INFO( "3. Built both default and hyper-reduced ROMs" ); - - // 4. Run the ROMs and capture trajectories - std::vector defaultRomTrajectory; + std::vector standardRomTrajectory; std::vector hypredRomTrajectory; - auto romSolution = runROM< FomSystem >( defaultRom, defaultTrialSpace, defaultReducedState, startTime, numSteps, dt, defaultRomTrajectory ); - auto hypredSolution = runROM< FomSystem >( hypredRom, hypredTrialSpace, hypredReducedState, startTime, numSteps, dt, hypredRomTrajectory ); - PRESSIOLOG_INFO( "4. Ran the ROMs and captured trajectories" ); + 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, romSolution ); - PRESSIOLOG_INFO( "5. Relative error between FOM and default ROM: {}", romError ); + 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 ); @@ -477,7 +494,7 @@ int main() { fomTrajectory.push_back(snapshots.stateSnapshots.col(i)); } writeTrajectoryToCSV< vector_t >("output", "fom_trajectory.csv", fomTrajectory); - writeTrajectoryToCSV< vector_t >("output", "default_rom_trajectory.csv", defaultRomTrajectory); + 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" ); From bb1a8fa17f024da67b33e5c7d2baf9282c2d60fe Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Mon, 9 Feb 2026 12:34:36 -0500 Subject: [PATCH 7/7] #65: driveby fix to versioning --- docs/build_requirements.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/build_requirements.txt b/docs/build_requirements.txt index 09de70f..7b462ba 100644 --- a/docs/build_requirements.txt +++ b/docs/build_requirements.txt @@ -1,6 +1,6 @@ -Sphinx==5.0.1 -furo==2022.6.4.1 -myst-parser==0.18.0 +Sphinx>=7.2.0 +furo>=2024.1.29 +myst-parser>=2.0.0 sphinx-copybutton==0.5.0 m2r2==0.3.2 sphinx-design