HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
DMDROMObject.cpp
Go to the documentation of this file.
1 #ifdef with_librom
2 
8 #include <basic.h>
9 #include <rom_object_dmd.h>
10 
35 DMDROMObject::DMDROMObject( const int a_vec_size,
36  const double a_dt,
37  const int a_rdim,
38  const int a_rank,
39  const int a_nproc,
40  const int a_sim_idx,
41  const int a_var_idx )
42 {
43  m_rank = a_rank;
44  m_nproc = a_nproc;
45 
46  m_dmd.clear();
47  m_dmd_is_trained.clear();
48  m_intervals.clear();
49  m_vec_size = a_vec_size;
50  m_dt = a_dt;
51  m_t_final = -1;
52  m_rdim = a_rdim;
53 
54  m_sim_idx = a_sim_idx;
55  m_var_idx = a_var_idx;
56 
57  m_num_window_samples = INT_MAX;
58 
59  char dirname_c_str[_MAX_STRING_SIZE_] = "DMD";
60  char write_snapshot_mat[_MAX_STRING_SIZE_] = "false";
61 
62  if (!m_rank) {
63 
64  FILE *in;
65  in = fopen(_LIBROM_INP_FNAME_,"r");
66 
67  if (in) {
68 
69  int ferr;
70  char word[_MAX_STRING_SIZE_];
71  ferr = fscanf(in,"%s", word); if (ferr != 1) return;
72 
73  if (std::string(word) == "begin") {
74  while (std::string(word) != "end") {
75  ferr = fscanf(in,"%s",word); if (ferr != 1) return;
76  if (std::string(word) == "dmd_num_win_samples") {
77  ferr = fscanf(in,"%d", &m_num_window_samples); if (ferr != 1) return;
78  } else if (std::string(word) == "dmd_dirname") {
79  ferr = fscanf(in,"%s", dirname_c_str); if (ferr != 1) return;
80  } else if (std::string(word) == "dmd_write_snapshot_mat") {
81  ferr = fscanf(in,"%s", write_snapshot_mat); if (ferr != 1) return;
82  }
83  if (ferr != 1) return;
84  }
85  } else {
86  fprintf( stderr, "Error: Illegal format in file \"%s\". Word read is: %s\n",
87  _LIBROM_INP_FNAME_, word);
88  return;
89  }
90 
91  fclose(in);
92 
93  }
94 
95  /* print useful stuff to screen */
96  printf("DMDROMObject details:\n");
97  printf(" number of samples per window: %d\n", m_num_window_samples);
98  printf(" directory name for DMD onjects: %s\n", dirname_c_str);
99  printf(" write snapshot matrix to file: %s\n", write_snapshot_mat);
100  if (m_sim_idx >= 0) {
101  printf(" simulation domain: %d\n", m_sim_idx);
102  }
103  if (m_var_idx >= 0) {
104  printf(" Vector component: %d\n", m_var_idx);
105  }
106  }
107 
108 #ifndef serial
109  MPI_Bcast(&m_num_window_samples,1,MPI_INT,0,MPI_COMM_WORLD);
110  MPI_Bcast(dirname_c_str,_MAX_STRING_SIZE_,MPI_CHAR,0,MPI_COMM_WORLD);
111  MPI_Bcast(write_snapshot_mat,_MAX_STRING_SIZE_,MPI_CHAR,0,MPI_COMM_WORLD);
112 #endif
113 
114  m_dirname = std::string( dirname_c_str );
115  m_write_snapshot_mat = (std::string(write_snapshot_mat) == "true");
116 
117  if (m_num_window_samples <= m_rdim) {
118  printf("ERROR:DMDROMObject::DMDROMObject() - m_num_window_samples <= m_rdim!!");
119  }
120 
121  m_tic = 0;
122  m_curr_win = 0;
123 }
124 
126 void DMDROMObject::takeSample( const CAROM::Vector& a_U,
127  const double a_time )
128 {
129  if (m_tic == 0) {
130 
131  m_dmd.push_back( new CAROM::DMD(m_vec_size, m_dt) );
132  m_dmd_is_trained.push_back(false);
133  m_intervals.push_back( Interval(a_time, m_t_final) );
134 
135  if (!m_rank) {
136  printf( "DMDROMObject::takeSample() - creating new DMD object for sim. domain %d, var %d, t=%f (total: %d).\n",
137  m_sim_idx, m_var_idx, m_intervals[m_curr_win].first, m_dmd.size());
138  }
139  m_dmd[m_curr_win]->takeSample( a_U.getData(), a_time );
140 
141  } else {
142 
143  m_dmd[m_curr_win]->takeSample( a_U.getData(), a_time );
144 
145  if (m_tic%m_num_window_samples == 0) {
146 
147  m_intervals[m_curr_win].second = a_time;
148  int ncol = m_dmd[m_curr_win]->getSnapshotMatrix()->numColumns();
149  if (!m_rank) {
150  printf( "DMDROMObject::train() - training DMD object %d for sim. domain %d, var %d with %d samples.\n",
151  m_curr_win, m_sim_idx, m_var_idx, ncol );
152  }
153 
154  if (m_write_snapshot_mat) {
155  char idx_string[_MAX_STRING_SIZE_];
156  sprintf(idx_string, "%04d", m_curr_win);
157  std::string fname_root(m_dirname + "/snapshot_mat_"+std::string(idx_string));
158  m_dmd[m_curr_win]->getSnapshotMatrix()->write(fname_root);
159  }
160  m_dmd[m_curr_win]->train(m_rdim);
162 
163  m_curr_win++;
164 
165  m_dmd.push_back( new CAROM::DMD(m_vec_size, m_dt) );
166  m_dmd_is_trained.push_back(false);
167  m_intervals.push_back( Interval(a_time, m_t_final) );
168  m_dmd[m_curr_win]->takeSample( a_U.getData(), a_time );
169  if (!m_rank) {
170  printf("DMDROMObject::takeSample() - creating new DMD object for sim. domain %d, var %d, t=%f (total: %d).\n",
171  m_sim_idx, m_var_idx, m_intervals[m_curr_win].first, m_dmd.size());
172  }
173  }
174 
175  }
176 
177  m_tic++;
178  return;
179 }
180 
183 {
184  /* make sure the number of columns for the last DMD isn't less than m_rdim */
185  {
186  int last_win = m_dmd.size() - 1;
187  int num_columns( m_dmd[last_win]->getSnapshotMatrix()->numColumns() );
188  if (num_columns <= m_rdim) {
189  m_dmd.pop_back();
190  m_intervals.pop_back();
191  m_intervals[m_intervals.size()-1].second = m_t_final;
192  if (!m_rank) {
193  printf("DMDROMObject::train() - last window DMD for sim. domain %d, var %d has %d sample(s) only; deleted it ",
194  m_sim_idx,
195  m_var_idx,
196  num_columns );
197  printf("(total: %d).\n", m_dmd.size());
198  }
199  }
200  }
201 
202  if (m_dmd.size() > 0) {
203  for (int i = 0; i < m_dmd.size(); i++) {
204  if (!m_dmd_is_trained[i]) {
205  int ncol = m_dmd[i]->getSnapshotMatrix()->numColumns();
206  if (!m_rank) {
207  printf( "DMDROMObject::train() - training DMD object %d for sim. domain %d, var %d with %d samples.\n",
208  m_curr_win, m_sim_idx, m_var_idx, ncol );
209  }
210  if (m_write_snapshot_mat) {
211  char idx_string[_MAX_STRING_SIZE_];
212  sprintf(idx_string, "%04d", i);
213  std::string fname_root(m_dirname + "/snapshot_mat_"+std::string(idx_string));
214  m_dmd[i]->getSnapshotMatrix()->write(fname_root);
215  }
216  m_dmd[i]->train(m_rdim);
217  }
218  }
219  } else {
220  printf("ERROR in DMDROMObject::train(): m_dmd is of size zero!");
221  }
222 
223  return;
224 }
225 
235 void DMDROMObject::save(const std::string& a_fname_root ) const
236 {
237  std::string fname_root = m_dirname + "/";
238  std::string summary_fname_root = m_dirname + "/";
239  std::string header_fname = m_dirname + "/";
240  if (a_fname_root == "") {
241  fname_root += "dmdobj_";
242  summary_fname_root += "dmd_summary_";
243  header_fname += "dmd_header.dat";
244  } else {
245  fname_root += (a_fname_root+"_dmdobj_");
246  summary_fname_root += (a_fname_root+"_dmd_summary_");
247  header_fname += (a_fname_root+"_dmd_header.dat");
248  }
249 
250  if (!m_rank) {
251  FILE* out;
252  out = fopen(header_fname.c_str(), "w");
253  fprintf(out, "%d\n", m_dmd.size());
254  for (int i = 0; i < m_dmd.size(); i++) {
255  fprintf(out, "%1.16e %1.16e\n", m_intervals[i].first, m_intervals[i].second);
256  }
257  fclose(out);
258  }
259 
260  for (int i = 0; i < m_dmd.size(); i++) {
261  char idx_string[_MAX_STRING_SIZE_];
262  sprintf(idx_string, "%04d", i);
263  std::string fname = fname_root + std::string(idx_string);
264  std::string summary_fname = summary_fname_root + std::string(idx_string);
265  if (!m_rank) {
266  printf( " Saving DMD object and summary (%s, %s).\n",
267  fname.c_str(), summary_fname.c_str() );
268  }
269  m_dmd[i]->save(fname);
270  //m_dmd[i]->summary(summary_fname);
271  }
272 
273  return;
274 }
275 
281 void DMDROMObject::load(const std::string& a_fname_root )
282 {
283  std::string fname_root = m_dirname + "/";
284  std::string header_fname = m_dirname + "/";
285  if (a_fname_root == "") {
286  fname_root += "dmdobj_";
287  header_fname += "dmd_header.dat";
288  } else {
289  fname_root += (a_fname_root+"_dmdobj_");
290  header_fname += (a_fname_root+"_dmd_header.dat");
291  }
292 
293  int num_dmds = 0;
294  std::vector<double> intervals_start(0), intervals_end(0);
295  if (!m_rank) {
296  FILE* in;
297  in = fopen(header_fname.c_str(), "r");
298  fscanf(in, "%d\n", &num_dmds);
299  intervals_start.resize(num_dmds);
300  intervals_end.resize(num_dmds);
301  for (int i = 0; i < num_dmds; i++) {
302  fscanf(in, "%lf", &intervals_start[i]);
303  fscanf(in, "%lf", &intervals_end[i]);
304  }
305  fclose(in);
306  }
307 #ifndef serial
308  MPI_Bcast(&num_dmds,1,MPI_INT,0,MPI_COMM_WORLD);
309  if (m_rank) {
310  intervals_start.resize(num_dmds);
311  intervals_end.resize(num_dmds);
312  }
313  MPI_Bcast(intervals_start.data(),num_dmds,MPI_DOUBLE,0,MPI_COMM_WORLD);
314  MPI_Bcast(intervals_end.data(),num_dmds,MPI_DOUBLE,0,MPI_COMM_WORLD);
315 #endif
316 
317  for (int i = 0; i < num_dmds; i++) {
318 
319  m_intervals.push_back( Interval(intervals_start[i],intervals_end[i]) );
320 
321  char idx_string[_MAX_STRING_SIZE_];
322  sprintf(idx_string, "%04d", i);
323  std::string fname = fname_root + std::string(idx_string);
324  if (!m_rank) {
325  printf( " Loading DMD object (%s), time window=[%1.2e,%1.2e].\n",
326  fname.c_str(),
327  m_intervals[i].first, m_intervals[i].second );
328  }
329  m_dmd.push_back( new CAROM::DMD(fname) );
330  }
331 
332  return;
333 }
334 
335 #endif
std::vector< bool > m_dmd_is_trained
virtual void save(const std::string &a_fname_root) const
Dynamic Mode Decomposition ROM object.
#define _LIBROM_INP_FNAME_
Definition: rom_object.h:15
virtual void load(const std::string &a_fname_root)
#define _MAX_STRING_SIZE_
Definition: basic.h:14
bool m_write_snapshot_mat
virtual void takeSample(const CAROM::Vector &, const double)
int m_num_window_samples
virtual void train()
std::pair< double, double > Interval
std::vector< Interval > m_intervals
Some basic definitions and macros.
std::string m_dirname
std::vector< CAROM::DMD * > m_dmd
DMDROMObject(const int, const double, const int, const int, const int, const int, const int a_var_idx=-1)