Skip to content

Commit

Permalink
Snapshot dmd (#289)
Browse files Browse the repository at this point in the history
* 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.
  • Loading branch information
ptranq authored Jun 21, 2024
1 parent 884d5ae commit 1c523cb
Show file tree
Hide file tree
Showing 4 changed files with 279 additions and 8 deletions.
92 changes: 84 additions & 8 deletions examples/dmd/wave_equation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
// =================================================================================
//
Expand All @@ -41,11 +76,13 @@
#include "mfem.hpp"
#include "algo/DMD.h"
#include "algo/NonuniformDMD.h"
#include "algo/SnapshotDMD.h"
#include "linalg/Vector.h"
#include <cmath>
#include <limits>
#include <fstream>
#include <iostream>
#include <sys/stat.h>

using namespace std;
using namespace mfem;
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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())
{
Expand All @@ -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;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -448,8 +503,16 @@ int main(int argc, char *argv[])
fom_timer.Stop();
dmd_training_timer.Start();
int curr_window = 0;

vector<CAROM::DMD*> 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);

Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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++;
}
Expand Down
1 change: 1 addition & 0 deletions lib/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ set(module_list
algo/DMDc
algo/AdaptiveDMD
algo/NonuniformDMD
algo/SnapshotDMD
algo/DifferentialEvolution
algo/greedy/GreedyCustomSampler
algo/greedy/GreedyRandomSampler
Expand Down
63 changes: 63 additions & 0 deletions lib/algo/SnapshotDMD.cpp
Original file line number Diff line number Diff line change
@@ -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 <algorithm>

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<Vector*> new_snapshots;
std::vector<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];
}

}
131 changes: 131 additions & 0 deletions lib/algo/SnapshotDMD.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
#ifndef included_SnapshotDMD_h
#define included_SnapshotDMD_h

#include "DMD.h"
#include <vector>

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<Vector*> getSnapshotVectors()
{
std::vector<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>(SnapshotDMD*& parametric_dmd,
std::vector<Vector*>& parameter_points,
std::vector<SnapshotDMD*>& 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<std::complex<double>> 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

0 comments on commit 1c523cb

Please sign in to comment.