libROM  v1.0
Data-driven physical simulation library
AdaptiveDMD.cpp
1 
11 // Description: Implementation of the AdaptiveDMD algorithm.
12 
13 #include "AdaptiveDMD.h"
14 #include "manifold_interp/VectorInterpolator.h"
15 
16 #include "linalg/Matrix.h"
17 #include "linalg/Vector.h"
18 #include <algorithm>
19 
20 namespace CAROM {
21 
22 AdaptiveDMD::AdaptiveDMD(int dim, double desired_dt, std::string rbf,
23  std::string interp_method,
24  double closest_rbf_val,
25  bool alt_output_basis,
26  Vector* state_offset) :
27  DMD(dim, alt_output_basis, state_offset)
28 {
29  CAROM_VERIFY(rbf == "G" || rbf == "IQ" || rbf == "IMQ");
30  CAROM_VERIFY(interp_method == "LS" || interp_method == "IDW"
31  || interp_method == "LP");
32  CAROM_VERIFY(closest_rbf_val >= 0.0 && closest_rbf_val <= 1.0);
33  d_dt = desired_dt;
34  d_interp_method = interp_method;
35  d_rbf = rbf;
36  d_closest_rbf_val = closest_rbf_val;
37 }
38 
40 {
41  for (auto interp_snapshot : d_interp_snapshots)
42  {
43  delete interp_snapshot;
44  }
45 }
46 
47 void AdaptiveDMD::train(double energy_fraction, const Matrix* W0,
48  double linearity_tol)
49 {
50  const Matrix* f_snapshots = getInterpolatedSnapshots();
51  CAROM_VERIFY(f_snapshots->numColumns() > 1);
52  CAROM_VERIFY(energy_fraction > 0 && energy_fraction <= 1);
53  d_energy_fraction = energy_fraction;
54  constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol);
55 
56  delete f_snapshots;
57 }
58 
59 void AdaptiveDMD::train(int k, const Matrix* W0, double linearity_tol)
60 {
61  const Matrix* f_snapshots = getInterpolatedSnapshots();
62  CAROM_VERIFY(f_snapshots->numColumns() > 1);
63  CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1);
64  d_energy_fraction = -1.0;
65  d_k = k;
66  constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol);
67 
68  delete f_snapshots;
69 }
70 
71 void AdaptiveDMD::interpolateSnapshots()
72 {
73  CAROM_VERIFY(d_interp_snapshots.size() == 0);
74  CAROM_VERIFY(d_snapshots.size() == d_sampled_times.size());
75  CAROM_VERIFY(d_sampled_times.size() > 1);
76 
77  if (d_rank == 0) std::cout << "Number of snapshots is: " << d_snapshots.size()
78  << std::endl;
79 
80  bool automate_dt = false;
81  if (d_dt <= 0.0)
82  {
83  automate_dt = true;
84  std::vector<double> d_sampled_dts;
85  for (int i = 1; i < d_sampled_times.size(); i++)
86  {
87  d_sampled_dts.push_back(d_sampled_times[i]->item(0) - d_sampled_times[i -
88  1]->item(0));
89  }
90 
91  auto m = d_sampled_dts.begin() + d_sampled_dts.size() / 2;
92  std::nth_element(d_sampled_dts.begin(), m, d_sampled_dts.end());
93  if (d_rank == 0) std::cout <<
94  "Setting desired dt to the median dt between the samples: " <<
95  d_sampled_dts[d_sampled_dts.size() / 2] << std::endl;
96  d_dt = d_sampled_dts[d_sampled_dts.size() / 2];
97  }
98  CAROM_VERIFY(d_sampled_times.back()->item(0) > d_dt);
99 
100  // Find the nearest dt that evenly divides the snapshots.
101  int num_time_steps = std::round(d_sampled_times.back()->item(0) / d_dt);
102  if (automate_dt && num_time_steps < d_sampled_times.size())
103  {
104  num_time_steps = d_sampled_times.size();
105  if (d_rank == 0) std::cout <<
106  "There will be less interpolated snapshots than FOM snapshots. dt will be decreased."
107  << std::endl;
108  }
109  double new_dt = d_sampled_times.back()->item(0) / num_time_steps;
110  if (new_dt != d_dt)
111  {
112  d_dt = new_dt;
113  if (d_rank == 0) std::cout << "Setting desired dt to " << d_dt <<
114  " to ensure a constant dt given the final sampled time." << std::endl;
115  }
116 
117  // Solve the linear system if required.
118  Matrix* f_T = NULL;
119  double epsilon = convertClosestRBFToEpsilon(d_sampled_times, d_rbf,
120  d_closest_rbf_val);
121  if (d_interp_method == "LS")
122  {
123  f_T = solveLinearSystem(d_sampled_times, d_snapshots, d_interp_method, d_rbf,
124  epsilon);
125  }
126 
127  // Create interpolated snapshots using d_dt as the desired dt.
128  for (int i = 0; i <= num_time_steps; i++)
129  {
130  double curr_time = i * d_dt;
131  if (d_rank == 0) std::cout << "Creating new interpolated sample at: " <<
132  d_t_offset + curr_time << std::endl;
133  CAROM::Vector* point = new Vector(&curr_time, 1, false);
134 
135  // Obtain distances from database points to new point
136  std::vector<double> rbf = obtainRBFToTrainingPoints(d_sampled_times,
137  d_interp_method, d_rbf, epsilon, point);
138 
139  // Obtain the interpolated snapshot.
140  CAROM::Vector* curr_interpolated_snapshot = obtainInterpolatedVector(
141  d_snapshots, f_T, d_interp_method, rbf);
142  d_interp_snapshots.push_back(curr_interpolated_snapshot);
143 
144  delete point;
145  }
146 
147  if (d_rank == 0) std::cout << "Number of interpolated snapshots is: " <<
148  d_interp_snapshots.size() << std::endl;
149 }
150 
152 {
153  return d_dt;
154 }
155 
157 {
158  if (d_interp_snapshots.size() == 0) interpolateSnapshots();
159  return createSnapshotMatrix(d_interp_snapshots);
160 }
161 
162 }
~AdaptiveDMD()
Destroy the AdaptiveDMD object.
Definition: AdaptiveDMD.cpp:39
double getTrueDt() const
Get the true dt between interpolated snapshots.
const Matrix * getInterpolatedSnapshots()
Get the interpolated snapshot matrix contained within d_interp_snapshots.
void train(double energy_fraction, const Matrix *W0=NULL, double linearity_tol=0.0)
Definition: AdaptiveDMD.cpp:47
Definition: DMD.h:71
double d_energy_fraction
The energy fraction used to obtain the DMD modes.
Definition: DMD.h:393
const Matrix * createSnapshotMatrix(std::vector< Vector * > snapshots)
Get the snapshot matrix contained within d_snapshots.
Definition: DMD.cpp:734
void constructDMD(const Matrix *f_snapshots, int rank, int num_procs, const Matrix *W0, double linearity_tol)
Construct the DMD object.
Definition: DMD.cpp:286
double d_dt
The time step size between samples.
Definition: DMD.h:348
std::vector< Vector * > d_sampled_times
The stored times of each sample.
Definition: DMD.h:363
double d_t_offset
The time offset of the first sample.
Definition: DMD.h:353
int d_rank
The rank of the process this object belongs to.
Definition: DMD.h:333
std::vector< Vector * > d_snapshots
std::vector holding the snapshots.
Definition: DMD.h:358
int d_num_procs
The number of processors being run on.
Definition: DMD.h:338
int d_k
The number of columns used after obtaining the SVD decomposition.
Definition: DMD.h:398
int numColumns() const
Returns the number of columns in the Matrix. This method will return the same value from each process...
Definition: Matrix.h:221
double convertClosestRBFToEpsilon(std::vector< Vector * > parameter_points, std::string rbf, double closest_rbf_val)
Convert closest RBF value to an epsilon value.
std::vector< double > obtainRBFToTrainingPoints(std::vector< Vector * > parameter_points, std::string interp_method, std::string rbf, double epsilon, Vector *point)
Compute the RBF from the parameter points with the unsampled parameter point.