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