From 1c523cb7b28534ea37d9f7ec3dca36cc586bf8f5 Mon Sep 17 00:00:00 2001 From: ptranq Date: Fri, 21 Jun 2024 15:09:32 -0700 Subject: [PATCH] Snapshot dmd (#289) * Implemented snapshot interpolation, and initial pass at integrating with DMD. Pivoting the DMD changes to derived class SnapshotDMD in future commits. Turned off "window overlap" in parametric_dw_csv. * Reduced scope to only the interpolator and a unit test. Will construct SnapshotDMD as a separat PR. * Fixed unit test compilation. * Stylize * Improved description of snapshot interpolator to include references * Stylize * Stylize * Implemented SnapshotDMD a derived class of DMD. Fixed SnapshotInterpolator and its unit test. * Addressed Dylan's Comments from PR#283. * Stylized. * Added energy fraction training function to SnapshotDMD * Added energy fraction training function to the snapshot DMD header. * Added new examples to CMakeLists * Added snapshot interpolator unit test to CMakeLists * Added snapshot_interpolation unit test * Removed old snapshot interpolator test * Added improvements for floating point comparison * Stylize * removed old snapshotinterpolator unit test from CMakeLists * Matching updates on snapshot_interpolator branch. * Finalized version of snapshot DMD. Added snapshot DMD as an option to the wave equation example. Improved wave equation example to include projected initial conditions in later windows, as well as print output to a specified directory to keep the running directory clean. * merged with snapshot_interpolator, and added example command for the wave equation example. * Renamed SnapshotInterpolator to PCHIPInterpolator. Addressed relevant PR Comments. * Removed redundant else statement. * Updated SnapshotDMD to use new names for PCHIPInterpolator, modified example run for wave equation. * Addressed Dylan's and Cole's comments. Added several CAROM_VERIFY statements to guarantee strictly increasing input / output ts, that there are >2 input ts, and that there are >1 output ts. * Removed some empty space. Improved wave_equation examples. * Addressed Dylans comments. * Removed duplicate sampling in wave_equation example. --- examples/dmd/wave_equation.cpp | 92 +++++++++++++++++++++-- lib/CMakeLists.txt | 1 + lib/algo/SnapshotDMD.cpp | 63 ++++++++++++++++ lib/algo/SnapshotDMD.h | 131 +++++++++++++++++++++++++++++++++ 4 files changed, 279 insertions(+), 8 deletions(-) create mode 100644 lib/algo/SnapshotDMD.cpp create mode 100644 lib/algo/SnapshotDMD.h diff --git a/examples/dmd/wave_equation.cpp b/examples/dmd/wave_equation.cpp index 84c98eb03..d90316863 100644 --- a/examples/dmd/wave_equation.cpp +++ b/examples/dmd/wave_equation.cpp @@ -20,10 +20,45 @@ // wave_equation -o 4 -tf 5 -nwinsamp 25 // // Output 1: -// Relative error of DMD solution (u) at t_final: 5 is 3.0483662e-05 -// Elapsed time for solving FOM: 3.122185e+00 second -// Elapsed time for training DMD: 6.904051e-01 second -// Elapsed time for predicting DMD: 2.496171e-03 second +// Relative error of DMD solution (u) at t_final: 5 is 0.00015906082 +// Elapsed time for solving FOM: 2.223148e+00 second +// Elapsed time for training DMD: 1.140992e+00 second +// Elapsed time for predicting DMD: 1.082472e-02 second +// +// Sample run with Time-Windowing DMD with initial condition projection +// +// Command 2: +// wave_equation -o 4 -tf 5 -nwinsamp 25 -proj +// +// Output 2: +// Relative error of DMD solution (u) at t_final: 5 is 0.029780998 +// Elapsed time for solving FOM: 2.348481e+00 second +// Elapsed time for training DMD: 1.141874e+00 second +// Elapsed time for predicting DMD: 1.075082e-02 second +// +// Sample run using snapshotDMD and using projected initial conditions for the +// prediction phase. +// +// Command 3: +// wave_equation -tf 2 -nwinsamp 25 -rdim 27 -visit -snap -proj +// +// Output 3: +// Relative error of DMD solution (u) at t_final: 2 is 0.0029767885 +// Elapsed time for solving FOM: 1.989041e-01 second +// Elapsed time for training DMD: 1.449013e+00 second +// Elapsed time for predicting DMD: 4.555228e-02 second +// +// Sample run using snapshotDMD without projecting initial conditions for the +// prediction phase. +// +// Command 4: +// wave_equation -tf 2 -nwinsamp 25 -rdim 27 -visit -snap -no-proj +// +// Output 4: +// Relative error of DMD solution (u) at t_final: 2 is 0.00049470862 +// Elapsed time for solving FOM: 1.869209e-01 second +// Elapsed time for training DMD: 1.152315e+00 second +// Elapsed time for predicting DMD: 3.621168e-02 second // // ================================================================================= // @@ -41,11 +76,13 @@ #include "mfem.hpp" #include "algo/DMD.h" #include "algo/NonuniformDMD.h" +#include "algo/SnapshotDMD.h" #include "linalg/Vector.h" #include #include #include #include +#include using namespace std; using namespace mfem; @@ -219,6 +256,10 @@ int main(int argc, char *argv[]) bool dirichlet = true; int vis_steps = 5; int precision = 8; + bool snapshotDMD = false; + bool project_initial_condition = false; + const char *temp_io_dir = "./wave_equation_out"; + std::string io_dir; cout.precision(precision); OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", @@ -256,6 +297,17 @@ int main(int argc, char *argv[]) "Reduced dimension for DMD."); args.AddOption(&windowNumSamples, "-nwinsamp", "--numwindowsamples", "Number of samples in DMD windows."); + args.AddOption(&snapshotDMD, "-snap", "--snapshot-dmd", "-no-snap", + "--no-snapshot-dmd", + "Use snapshot DMD instead of standard DMD."); + args.AddOption(&project_initial_condition, "-proj", "--project-initial", + "-no-proj", + "--no-project-initial", + "Project the initial condition in subsequent window as the " + "final value of the current window."); + args.AddOption(&temp_io_dir, "-io", "--io-dir-name", + "Name of the sub-folder to load/dump input/output files " + "within the current directory."); args.Parse(); if (!args.Good()) { @@ -264,6 +316,9 @@ int main(int argc, char *argv[]) } args.PrintOptions(cout); + io_dir = temp_io_dir; + mkdir(io_dir.c_str(), 0777); + if (rdim <= 0 && rdim != -1) { cout << "rdim is set to " << rdim << " , rdim can only be a positive integer or -1" << endl; @@ -406,7 +461,7 @@ int main(int argc, char *argv[]) dudt_gf.Save(osol); } - VisItDataCollection visit_dc("Wave_Equation", mesh); + VisItDataCollection visit_dc(io_dir+"/Wave_Equation", mesh); visit_dc.RegisterField("solution", &u_gf); visit_dc.RegisterField("rate", &dudt_gf); if (visit) @@ -448,8 +503,16 @@ int main(int argc, char *argv[]) fom_timer.Stop(); dmd_training_timer.Start(); int curr_window = 0; + vector dmd_u; - dmd_u.push_back(new CAROM::DMD(u.Size(), dt)); + if(snapshotDMD) + { + dmd_u.push_back(new CAROM::SnapshotDMD(u.Size(), dt)); + } + else + { + dmd_u.push_back(new CAROM::DMD(u.Size(), dt)); + } dmd_u[curr_window]->takeSample(u.GetData(), t); ts.push_back(t); @@ -487,7 +550,14 @@ int main(int argc, char *argv[]) if (!last_step) { curr_window++; - dmd_u.push_back(new CAROM::DMD(u.Size(), dt)); + if(snapshotDMD) + { + dmd_u.push_back(new CAROM::SnapshotDMD(u.Size(), dt)); + } + else + { + dmd_u.push_back(new CAROM::DMD(u.Size(), dt)); + } dmd_u[curr_window]->takeSample(u.GetData(), t); } @@ -528,7 +598,7 @@ int main(int argc, char *argv[]) // 10. Predict using DMD. cout << "Predicting temperature using DMD" << endl; CAROM::Vector* result_u = nullptr; - VisItDataCollection dmd_visit_dc("DMD_Wave_Equation", mesh); + VisItDataCollection dmd_visit_dc(io_dir+"/DMD_Wave_Equation", mesh); dmd_visit_dc.RegisterField("solution", &u_gf); curr_window = 0; if (visit) { @@ -558,6 +628,12 @@ int main(int argc, char *argv[]) if (i % windowNumSamples == 0 && i < ts.size()-1) { + if(project_initial_condition) + { + result_u = dmd_u[curr_window]->predict(ts[i]); + cout << "Projecting solution for new window at " << ts[i] << endl; + dmd_u[curr_window+1]->projectInitialCondition(result_u, ts[i]); + } delete dmd_u[curr_window]; curr_window++; } diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 23343006e..873fd3db2 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -38,6 +38,7 @@ set(module_list algo/DMDc algo/AdaptiveDMD algo/NonuniformDMD + algo/SnapshotDMD algo/DifferentialEvolution algo/greedy/GreedyCustomSampler algo/greedy/GreedyRandomSampler diff --git a/lib/algo/SnapshotDMD.cpp b/lib/algo/SnapshotDMD.cpp new file mode 100644 index 000000000..c76ffd15f --- /dev/null +++ b/lib/algo/SnapshotDMD.cpp @@ -0,0 +1,63 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: Implementation of the AdaptiveDMD algorithm. + +#include "manifold_interp/PCHIPInterpolator.h" +#include "SnapshotDMD.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include "utils/CSVDatabase.h" +#include + +namespace CAROM { + +SnapshotDMD::~SnapshotDMD() +{ +} + +void SnapshotDMD::train(int k, const Matrix* W0, double linearity_tol) +{ + CAROM_VERIFY(d_snapshots.size() > 0); + + if(k >= d_snapshots.size()) + { + interpolateToNSnapshots(k + 1); + } + + const Matrix* f_snapshots = getSnapshotMatrix(); + CAROM_VERIFY(f_snapshots->numColumns() > 1); + CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1); + d_energy_fraction = -1.0; + d_k = k; + constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol); + + delete f_snapshots; +} + +void SnapshotDMD::train(double energy_fraction, const Matrix* W0, + double linearity_tol) +{ + DMD::train(energy_fraction,W0,linearity_tol); +} + +void SnapshotDMD::interpolateToNSnapshots(int n) +{ + PCHIPInterpolator* interp = new PCHIPInterpolator(); + std::vector new_snapshots; + std::vector new_times; + + interp->interpolate(d_sampled_times,d_snapshots,n,new_times,new_snapshots); + d_snapshots = new_snapshots; + d_sampled_times = new_times; + d_dt = d_sampled_times[2]->getData()[0]-d_sampled_times[1]->getData()[0]; +} + +} diff --git a/lib/algo/SnapshotDMD.h b/lib/algo/SnapshotDMD.h new file mode 100644 index 000000000..404147a42 --- /dev/null +++ b/lib/algo/SnapshotDMD.h @@ -0,0 +1,131 @@ +#ifndef included_SnapshotDMD_h +#define included_SnapshotDMD_h + +#include "DMD.h" +#include + +namespace CAROM { + +class Vector; +class SnapshotDMD : public DMD +{ +public: + + /** + * @brief Constructor. Basic DMD with uniform time step size. + * Inherited directly from base DMD class. + * + * @param[in] dim The full-order state dimension. + * @param[in] dt The dt between samples. + * @param[in] alt_output_basis Whether to use the alternative basis for + * output, i.e. phi = U^(+)*V*Omega^(-1)*X. + * @param[in] state_offset The state offset. + */ + SnapshotDMD(int dim, double dt, bool alt_output_basis = false, + Vector* state_offset = NULL) : DMD(dim,dt,alt_output_basis,state_offset) {} + + /** + * @brief Constructor. DMD from saved models. Inherited directly + * from base DMD class. + * + * @param[in] base_file_name The base part of the filename of the + * database to load when restarting from a save. + */ + SnapshotDMD(std::string base_file_name) : DMD(base_file_name) {} + + /** + * @brief Destroy the SnapshotDMD object + */ + ~SnapshotDMD(); + + /** + * @brief Interpolate the current snapshots to n, new snapshots + * distributed uniformly over the currently sampled time + * domain. + * @param[in] n The number of desired snapshots. + */ + void interpolateToNSnapshots(int n); + + /** + * @brief Train the DMD model with specified reduced dimension. If k is + * too large then new snapshots are computed using + * interpolateToNSnapshots(k+1). + * + * @param[in] k The number of modes to keep after doing SVD. + * @param[in] W0 The initial basis to prepend to W. + * @param[in] linearity_tol The tolerance for determining whether a column + * of W is linearly independent with W0. + */ + void train(int k, const Matrix* W0 = NULL, double linearity_tol = 0.0); + + /** + * @brief Train the DMD model with energy fraction criterion. + * + * @param[in] energy_fraction The energy fraction to keep after doing SVD. + * @param[in] W0 The initial basis to prepend to W. + * @param[in] linearity_tol The tolerance for determining whether a column + * of W is linearly independent with W0. + */ + void train(double energy_fraction, const Matrix* W0 = NULL, + double linearity_tol = 0.0); + + /** + * @brief Returns a copy of the current snapshot vector "d_snapshots" + */ + std::vector getSnapshotVectors() + { + std::vector return_snapshots(d_snapshots); + return return_snapshots; + } +protected: + /** + * @brief Obtain DMD model interpolant at desired parameter point by + * interpolation of DMD models from training parameter points. + * + * @param[in] parametric_dmd The interpolant DMD model at the desired point. + * @param[in] parameter_points The training parameter points. + * @param[in] dmds The DMD objects associated with + * each training parameter point. + * @param[in] desired_point The desired point at which to create a parametric DMD. + * @param[in] rbf The RBF type ("G" == gaussian, + * "IQ" == inverse quadratic, + * "IMQ" == inverse multiquadric) + * @param[in] interp_method The interpolation method type + * ("LS" == linear solve, + * "IDW" == inverse distance weighting, + * "LP" == lagrangian polynomials) + * @param[in] closest_rbf_val The RBF parameter determines the width of influence. + * Set the RBF value of the nearest two parameter points to a value between 0.0 to 1.0 + * @param[in] reorthogonalize_W Whether to reorthogonalize the interpolated W (basis) matrix. + */ + friend void getParametricDMD(SnapshotDMD*& parametric_dmd, + std::vector& parameter_points, + std::vector& dmds, + Vector* desired_point, + std::string rbf, + std::string interp_method, + double closest_rbf_val, + bool reorthogonalize_W); + + /** + * @brief Constructor. + * + * @param[in] eigs d_eigs + * @param[in] phi_real d_phi_real + * @param[in] phi_imaginary d_phi_imaginary + * @param[in] k d_k + * @param[in] dt d_dt + * @param[in] t_offset d_t_offset + * @param[in] state_offset d_state_offset + * @param[in] derivative_offset d_derivative_offset + */ + SnapshotDMD(std::vector> eigs, Matrix* phi_real, + Matrix* phi_imaginary, int k, + double dt, double t_offset, + Vector* state_offset) : + DMD(eigs, phi_real, phi_imaginary, k, dt, t_offset, state_offset) {} + +private: +}; +} +#endif \ No newline at end of file