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/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 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..f5d5136 --- /dev/null +++ b/docs/source/tutorial/background.rst @@ -0,0 +1,29 @@ +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 +---------------------- + +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..4d06010 --- /dev/null +++ b/docs/source/tutorial/hyperreduction.rst @@ -0,0 +1,127 @@ +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. + +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 +----------------------------------------- + +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 +---------------------------------------------- + +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..d3ad089 --- /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). +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. + +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..4664494 --- /dev/null +++ b/docs/source/tutorial/understand.rst @@ -0,0 +1,417 @@ +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 depend on user-set macros to +enable or disable features at compile time. These macros may +be set via CMake configuration or directly in the source code +before including any Pressio headers. + +A complete discussion of Pressio's macros can be found in the +`pressio-rom documentation `_. + +In this tutorial, we simply define the key macros in the source code, +like so: + +.. 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 + + PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::info); + +You can change the ``LogLevel`` argument to ``debug`` for more information, +or ``sparse`` for less. + +Then we'll just define some types for convenience: + +.. 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 + +.. 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 +----------------------------------------- + +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 +---------------------------------------------- + +.. 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. + +**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/CMakeLists.txt b/full-tutorial/CMakeLists.txt new file mode 100644 index 0000000..41f1712 --- /dev/null +++ b/full-tutorial/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required(VERSION 3.18) + +project(burgers_tutorial CXX) + +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 + ${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 +) + +# 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 new file mode 100644 index 0000000..8797829 --- /dev/null +++ b/full-tutorial/burgers.cpp @@ -0,0 +1,504 @@ +/** + * 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. (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. + */ + +//////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// Step 0: Pressio Setup +////////////////////////////////////////////////////////////////////////////////// + +/** + * Pressio uses macros to enable features like logging and TPLs. + * Typically, these would be set during configuration, but we will + * define them explicitly here for demonstration purposes. + * + * Importantly, these macros should be defined BEFORE including any Pressio + * headers. + */ +#define PRESSIO_ENABLE_LOGGING +#define PRESSIO_ENABLE_TPL_EIGEN + +/** + * Due to the hierarchical structure of Pressio, including + * pressio/rom.hpp will also pull in all necessary dependencies + * (including Eigen, since we defined the macro above). + * The logging macros (from pressio-log) are also included here. + */ +#include + +/** + * We'll also need some helper functions along the way. + */ +#include "helpers.h" + +/** + * And some functions for hyper-reduction. + */ +#include "hyperreduction.h" + +/** + * And we'll use std::vectors to store snapshots. + */ +#include + +/** + * Pressio supports various linear algebra backends (such as + * Eigen, Tpetra, and Kokkos) and provides a unified interface + * with the pressio::ops library for common operations. In this + * tutorial, we will use Eigen. + */ +using vector_t = Eigen::VectorXd; +using matrix_t = Eigen::MatrixXd; + +/////////////////////////////////////////////////////////////////////////////// +// Step 1: Define the Full-Order Model (FOM) +/////////////////////////////////////////////////////////////////////////////// + +/** + * Pressio requires a specific API to be used for the FOM class. This API can + * vary depending on the type of problem. For the 1D Burgers' equation, + * we will use the API for a semi-discrete FOM for time-dependent problems of + * the form: + * + * d/dt y(t) = f(y, t). + * + * The various FOM APIs can be found in pressio/rom/concepts.hpp. + * + * According to the semi-discrete API, our FOM class must define the following: + * + * 1. Three core types + * - time_type, state_type, rhs_type + * + * 2. A method to create the RHS vector + * - createRhs() -> rhs_type + * + * 3. A method to compute the RHS + * - rhs(const state_type & u, time_type t, rhs_type & f) -> void + * + * And that's it! + */ +class Burgers1dFom +{ + public: + // API Requirement: Core types + using time_type = double; + using state_type = vector_t; + using rhs_type = vector_t; + + private: + std::size_t N_{ }; + double dx_{ }; + double nu_{ }; + + public: + Burgers1dFom( std::size_t N, double domainLength, double viscosity ) + : N_( N ), dx_( domainLength / ( N ) ), nu_( viscosity ) + { + assert( N_ >= 3 ); + assert( dx_ > 0.0 ); + assert( nu_ > 0.0 ); + } + + std::size_t N() const { return N_; } + double dx() const { return dx_; } + double nu() const { return nu_; } + + // API Requirement: createRhs() method + rhs_type createRhs() const + { + return rhs_type::Zero( N_ ); + } + + // API Requirement: rhs(...) method + void rhs( + const state_type & u, + time_type /*t*/, + rhs_type & f + ) const + { + assert( u.size() == N_ ); + assert( f.size() == N_ ); + + // periodic indexing helpers + auto ip = [ this ]( std::size_t i ){ return ( i + 1 ) % N_; }; + auto im = [ this ]( std::size_t i ){ return ( i + N_ - 1 ) % N_; }; + + const double inv2dx = 1.0 / ( 2.0 * dx_ ); + const double invdx2 = 1.0 / ( dx_ * dx_ ); + + for ( std::size_t i = 0; i < N_; ++i ) { + const std::size_t iL = im( i ); + const std::size_t iR = ip( i ); + + // convective term: -(1/2 d/dx (u^2)) using centered flux difference + // d/dx (u^2) ~ (u_{i+1}^2 - u_{i-1}^2) / (2 dx) + const double dudt_conv = + -0.5 * ( ( u[ iR ] * u[ iR ] ) - ( u[ iL ] * u[ iL ] ) ) * inv2dx; + + // viscous term: nu * u_xx using second-order centered FD + const double dudt_diff = + nu_ * ( u[ iR ] - 2.0 * u[ i ] + u[ iL ] ) * invdx2; + f[ i ] = dudt_conv + dudt_diff; + } + } +}; + +/////////////////////////////////////////////////////////////////////////////// +// Step 2: Run the FOM to generate snapshots (states and RHS) +//////////////////////////////////////////////////////////////////////////////// + +struct SnapshotSet +{ + matrix_t stateSnapshots; + matrix_t rhsSnapshots; +}; + +/** + * Building and running the FOM doesn't require any Pressio utilities. + * For this example, we'll define a simple Forward Euler time integrator + * to advance the FOM in time and collect the snapshots that will be used + * to build the ROM with Pressio. + */ +template +SnapshotSet runFOM( FomSystem& fom, typename FomSystem::time_type startTime, int numSteps, double dt ) +{ + // We will store each snapshot vector in a std::vector for now + // For a typical ROM, we would only need the state snapshots. However, + // for hyper-reduction, we will also need snapshots of the RHS. + std::vector< vector_t > stateSnaps; + std::vector< vector_t > rhsSnaps; + + // Compute the initial condition + auto u = computeInitialCondition< FomSystem, vector_t >( fom ); + + // Time integration loop (using Forward Euler) + double t = startTime; + for ( int step = 0; step < numSteps; ++step ) { + stateSnaps.push_back( u ); // Store state snapshot + auto f = fom.createRhs(); // Create RHS vector + fom.rhs( u, t, f ); // Compute RHS + rhsSnaps.push_back( f ); // Store RHS snapshot + u += dt * f; // Update solution + t += dt; // Advance time + } + + // Convert our vector of snapshot vectors into matrices + const std::size_t numSnapshots = stateSnaps.size(); + matrix_t stateMatrix( fom.N(), numSnapshots ); + matrix_t rhsMatrix( fom.N(), numSnapshots ); + for ( std::size_t i = 0; i < numSnapshots; ++i ) { + stateMatrix.col( i ) = stateSnaps[ i ]; + rhsMatrix.col( i ) = rhsSnaps[ i ]; + } + + // Return both snapshot matrices + return SnapshotSet{ std::move(stateMatrix), std::move(rhsMatrix) }; +} + +/////////////////////////////////////////////////////////////////////////////// +// Step 3: Build the ROM from the snapshot matrix +/////////////////////////////////////////////////////////////////////////////// + +/** + * Now we use the Pressio ecosystem to construct a ROM representation + * with the snapshot matrix. As before, there are various APIs that we can meet + * to use different types of ROMs. We'll use a basic Galerkin ROM here. + */ + +template +auto buildStandardGalerkinROM( FomSystem& fom, auto& stepScheme, auto& trialSpace ) +{ + return pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom + ); +} + +/** + * Hyper-reduced ROMs do not evaluate the full RHS at every time step. + * Instead, they sample the RHS at selected indices and use a + * projection to approximate the reduced RHS. This can greatly + * reduce the computational cost of evaluating the ROM, especially + * when the FOM is large and the RHS evaluation is expensive. + * + * However, there are trade-offs in accuracy, as you'll see. + */ +template +auto buildHyperReducedGalerkinROM( FomSystem& fom, auto& stepScheme, + auto& trialSpace, auto& hyperReducer ) +{ + // Build the hyperreduced Galerkin ROM by passing the hyper-reducer + // functor to the ROM factory function. + return pressio::rom::galerkin::create_unsteady_explicit_problem( + stepScheme, trialSpace, fom, hyperReducer + ); +} + +/////////////////////////////////////////////////////////////////////////////// +// Step 4: Run the ROM with trajectory capture +/////////////////////////////////////////////////////////////////////////////// +/** + * To run the ROM, we will define a simple observer that captures + * the reduced state at each time step and reconstructs the full-order + * state using the trial space. This allows us to capture the full + * trajectory of the ROM solution for later analysis. + */ +template +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 + ); + + /** + * In Pressio, observers are functors that are called at each time step + * during time integration. Here, we define an observer that captures + * the full-order state at each time step by reconstructing it from + * the reduced state using the trial space. + * + * Observers must meet a specific API defined by Pressio. Namely, they + * must implement the operator() with the signature shown in helpers.h. + */ + auto observer = RomObserver< vector_t, decltype(trialSpace) >( trajectory, trialSpace ); + + // Run the ROM time integration with observer + pressio::ode::advance( rom, reducedState, policy, observer ); + + /** + * At this point, reducedState contains the ROM solution at the + * final timestep. We just have to reconstruct the full-order + * solution from the time reduced coordinates via the trial subspace. + */ + auto romSolution = trialSpace.createFullStateFromReducedState( reducedState ); + + return romSolution; +} + +/////////////////////////////////////////////////////////////////////////////// +// Step 5: Compare FOM and ROM solutions +/////////////////////////////////////////////////////////////////////////////// + +/** + * To compare the FOM and ROM solutions, we can compute the relative error + * between the two solution vectors at the final timestep. + */ +double compareFomAndRom( const matrix_t& fomSolution, const matrix_t& romSolution ) +{ + double error = ( fomSolution - romSolution ).norm() / fomSolution.norm(); + return error; +} + +/////////////////////////////////////////////////////////////////////////////// +// Main driver +/////////////////////////////////////////////////////////////////////////////// + +int main() { + // Initialize the logger so we can see Pressio output + // Change the log level to "debug" for more information + PRESSIOLOG_INITIALIZE(pressiolog::LogLevel::info); + + using FomSystem = Burgers1dFom; + using time_t = typename FomSystem::time_type; + + // 1. Create the FOM + const std::size_t N = 100; // Number of spatial points + const double domainLength = 1.0; // Domain length + const double viscosity = 0.01; // Viscosity parameter + FomSystem fom( N, domainLength, viscosity ); + PRESSIOLOG_INFO( "1. Created FOM with N = {}, dx = {}, nu = {}", fom.N(), fom.dx(), fom.nu() ); + + // 2. Run the FOM to generate a snapshot matrix + const double dt = 0.001; + const int numSteps = 1000; + time_t startTime = 0.0; + SnapshotSet snapshots = runFOM< FomSystem >( fom, startTime, numSteps, dt ); + PRESSIOLOG_INFO( "2. Generated snapshot matrix with {} snapshots", snapshots.stateSnapshots.cols() ); + + // 3 and 4. Build and run the ROM(s) + + /** + * This tutorial covers both standard and hyper-reduced ROMs. + * These variables below will hold the trajectories and solutions, + * which we can analyze at the end. + */ + std::vector standardRomTrajectory; + std::vector hypredRomTrajectory; + matrix_t standardSolution; + matrix_t hypredSolution; + + // Start with the standard ROM + PRESSIOLOG_INFO( "Standard ROM:" ); + { + // 3. Build the ROM + + // Build the state trial space (Phi) from state snapshots + auto standardTrialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom ); + // Create the initial reduced state by projecting the FOM initial condition + auto standardReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( standardTrialSpace, fom ); + // Select the time integration scheme with Pressio + auto stepScheme = pressio::ode::StepScheme::ForwardEuler; + // Build the standard ROM + auto standardRom = buildStandardGalerkinROM< FomSystem >( fom, stepScheme, standardTrialSpace ); + PRESSIOLOG_INFO( " 3. Built standard Galerkin ROM" ); + + // 4. Run the ROM and capture the trajectory + standardSolution = runROM< FomSystem >( standardRom, standardTrialSpace, standardReducedState, startTime, numSteps, dt, standardRomTrajectory ); + PRESSIOLOG_INFO( " 4. Ran standard ROM and captured trajectory" ); + } + + // Follow a similar process for hyper-reduced ROM + PRESSIOLOG_INFO( "Hyper-reduced ROM:" ); + { + // 3. Build the ROM + auto hypredTrialSpace = snapshots_to_trial_space< FomSystem, matrix_t, vector_t >( snapshots.stateSnapshots, fom ); + auto hypredReducedState = trial_space_to_reduced_state< FomSystem, vector_t >( hypredTrialSpace, fom ); + auto stepScheme = pressio::ode::StepScheme::ForwardEuler; + /** + * Here, we build the hyper-reducer using the RHS snapshot matrix + * and the trial space. The hyper-reducer is a functor that will + * be used by the ROM to compute the reduced RHS from sampled FOM RHSs. + */ + auto hyperReducer = buildHyperReducer< vector_t, matrix_t >( snapshots.rhsSnapshots, hypredTrialSpace ); + auto hypredRom = buildHyperReducedGalerkinROM< FomSystem >( fom, stepScheme, hypredTrialSpace, hyperReducer ); + PRESSIOLOG_INFO( " 3. Built hyper-reduced Galerkin ROM" ); + + // 4. Run the ROM and capture the trajectory + hypredSolution = runROM< FomSystem >( hypredRom, hypredTrialSpace, hypredReducedState, startTime, numSteps, dt, hypredRomTrajectory ); + PRESSIOLOG_INFO( " 4. Ran hyper-reduced ROM and captured trajectory" ); + } + + // 5. Compare ROM solution against FOM solution (the last snapshot) + auto fomSolution = snapshots.stateSnapshots.col( snapshots.stateSnapshots.cols() - 1 ); + auto romError = compareFomAndRom( fomSolution, standardSolution ); + PRESSIOLOG_INFO( "5. Relative error between FOM and standard ROM: {}", romError ); + auto hypredError = compareFomAndRom( fomSolution, hypredSolution ); + PRESSIOLOG_INFO( " Relative error between FOM and hyper-reduced ROM: {}", hypredError ); + + // 6. Write trajectories to CSV files in output directory + // Get FOM trajectory from snapshots + std::vector fomTrajectory; + for (int i = 0; i < snapshots.stateSnapshots.cols(); ++i) { + fomTrajectory.push_back(snapshots.stateSnapshots.col(i)); + } + writeTrajectoryToCSV< vector_t >("output", "fom_trajectory.csv", fomTrajectory); + writeTrajectoryToCSV< vector_t >("output", "standard_rom_trajectory.csv", standardRomTrajectory); + writeTrajectoryToCSV< vector_t >("output", "hypred_rom_trajectory.csv", hypredRomTrajectory); + PRESSIOLOG_INFO( "6. Wrote trajectories to output/{fom,default_rom,hypred_rom}_trajectory.csv" ); + + // Finalize the logger + PRESSIOLOG_FINALIZE(); + return 0; +} diff --git a/full-tutorial/helpers.h b/full-tutorial/helpers.h new file mode 100644 index 0000000..be310a8 --- /dev/null +++ b/full-tutorial/helpers.h @@ -0,0 +1,148 @@ +#include +#include +#include + +#include + +/** + * 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 ) +{ + // 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; +} + +// 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(); +} + +/** + * 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 ) +{ + // 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"; + } + + file.close(); +} diff --git a/full-tutorial/hyperreduction.h b/full-tutorial/hyperreduction.h new file mode 100644 index 0000000..06eed69 --- /dev/null +++ b/full-tutorial/hyperreduction.h @@ -0,0 +1,95 @@ +#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; + } +}; + +/** + * 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 + 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 (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); + + // 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; +} 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()