libMesh
rb_evaluation.C
Go to the documentation of this file.
1 // rbOOmit: An implementation of the Certified Reduced Basis method.
2 // Copyright (C) 2009, 2010 David J. Knezevic
3 
4 // This file is part of rbOOmit.
5 
6 // rbOOmit is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 
11 // rbOOmit is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
15 
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20 // rbOOmit includes
21 #include "libmesh/rb_evaluation.h"
22 #include "libmesh/rb_theta_expansion.h"
23 
24 // libMesh includes
25 #include "libmesh/libmesh_version.h"
26 #include "libmesh/system.h"
27 #include "libmesh/numeric_vector.h"
28 #include "libmesh/parallel.h"
29 #include "libmesh/libmesh_logging.h"
30 #include "libmesh/xdr_cxx.h"
31 #include "libmesh/mesh_tools.h"
32 #include "libmesh/utility.h"
33 
34 // C/C++ includes
35 #include <sys/types.h>
36 #include <sys/stat.h>
37 #include <errno.h>
38 #include <fstream>
39 #include <sstream>
40 #include <iomanip>
41 
42 namespace libMesh
43 {
44 
46  :
47  ParallelObject(comm_in),
48  evaluate_RB_error_bound(true),
49  compute_RB_inner_product(false),
50  rb_theta_expansion(libmesh_nullptr)
51 {
52 
53 }
54 
56 {
57  this->clear();
58 }
59 
61 {
62  LOG_SCOPE("clear()", "RBEvaluation");
63 
64  // Clear the basis functions
65  for (std::size_t i=0; i<basis_functions.size(); i++)
66  {
67  if (basis_functions[i])
68  {
69  basis_functions[i]->clear();
70  delete basis_functions[i];
72  }
73  }
75 
77 
78  // Clear the Greedy param list
79  for (std::size_t i=0; i<greedy_param_list.size(); i++)
81  greedy_param_list.clear();
82 }
83 
85 {
86  rb_theta_expansion = &rb_theta_expansion_in;
87 }
88 
90 {
92  libmesh_error_msg("Error: rb_theta_expansion hasn't been initialized yet");
93 
94  return *rb_theta_expansion;
95 }
96 
98 {
100  {
101  return true;
102  }
103  else
104  {
105  return false;
106  }
107 }
108 
109 void RBEvaluation::resize_data_structures(const unsigned int Nmax,
110  bool resize_error_bound_data)
111 {
112  LOG_SCOPE("resize_data_structures()", "RBEvaluation");
113 
114  if (Nmax < this->get_n_basis_functions())
115  libmesh_error_msg("Error: Cannot set Nmax to be less than the current number of basis functions.");
116 
117  // Resize/clear inner product matrix
119  RB_inner_product_matrix.resize(Nmax,Nmax);
120 
121  // Allocate dense matrices for RB solves
123 
124  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
125  {
126  // Initialize the memory for the RB matrices
127  RB_Aq_vector[q].resize(Nmax,Nmax);
128  }
129 
131 
132  for (unsigned int q=0; q<rb_theta_expansion->get_n_F_terms(); q++)
133  {
134  // Initialize the memory for the RB vectors
135  RB_Fq_vector[q].resize(Nmax);
136  }
137 
138 
139  // Initialize the RB output vectors
141  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
142  {
144  for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
145  {
146  RB_output_vectors[n][q_l].resize(Nmax);
147  }
148  }
149 
150  // Initialize vectors storing output data
152 
153 
154  if (resize_error_bound_data)
155  {
156  // Initialize vectors for the norms of the Fq representors
157  unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
158  Fq_representor_innerprods.resize(Q_f_hat);
159 
160  // Initialize vectors for the norms of the representors
162  for (unsigned int i=0; i<rb_theta_expansion->get_n_F_terms(); i++)
163  {
165  for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
166  {
167  Fq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
168  }
169  }
170 
171  unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
172  Aq_Aq_representor_innerprods.resize(Q_a_hat);
173  for (unsigned int i=0; i<Q_a_hat; i++)
174  {
175  Aq_Aq_representor_innerprods[i].resize(Nmax);
176  for (unsigned int j=0; j<Nmax; j++)
177  {
178  Aq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
179  }
180  }
181 
183 
184  // Resize the output dual norm vectors
186  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
187  {
188  unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
189  output_dual_innerprods[n].resize(Q_l_hat);
190  }
191 
192  // Clear and resize the vector of Aq_representors
194 
196  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
197  {
198  Aq_representor[q_a].resize(Nmax);
199  }
200  }
201 }
202 
204 {
205  libmesh_assert_less (i, basis_functions.size());
206 
207  return *(basis_functions[i]);
208 }
209 
211 {
212  LOG_SCOPE("rb_solve()", "RBEvaluation");
213 
214  if (N > get_n_basis_functions())
215  libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");
216 
217  const RBParameters & mu = get_parameters();
218 
219  // Resize (and clear) the solution vector
220  RB_solution.resize(N);
221 
222  // Assemble the RB system
223  DenseMatrix<Number> RB_system_matrix(N,N);
224  RB_system_matrix.zero();
225 
226  DenseMatrix<Number> RB_Aq_a;
227  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
228  {
229  RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);
230 
231  RB_system_matrix.add(rb_theta_expansion->eval_A_theta(q_a, mu), RB_Aq_a);
232  }
233 
234  // Assemble the RB rhs
235  DenseVector<Number> RB_rhs(N);
236  RB_rhs.zero();
237 
238  DenseVector<Number> RB_Fq_f;
239  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
240  {
241  RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);
242 
243  RB_rhs.add(rb_theta_expansion->eval_F_theta(q_f, mu), RB_Fq_f);
244  }
245 
246  // Solve the linear system
247  if (N > 0)
248  {
249  RB_system_matrix.lu_solve(RB_rhs, RB_solution);
250  }
251 
252  // Evaluate RB outputs
253  DenseVector<Number> RB_output_vector_N;
254  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
255  {
256  RB_outputs[n] = 0.;
257  for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
258  {
259  RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
260  RB_outputs[n] += rb_theta_expansion->eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
261  }
262  }
263 
264  if (evaluate_RB_error_bound) // Calculate the error bounds
265  {
266  // Evaluate the dual norm of the residual for RB_solution_vector
267  Real epsilon_N = compute_residual_dual_norm(N);
268 
269  // Get lower bound for coercivity constant
270  const Real alpha_LB = get_stability_lower_bound();
271  // alpha_LB needs to be positive to get a valid error bound
272  libmesh_assert_greater ( alpha_LB, 0. );
273 
274  // Evaluate the (absolute) error bound
275  Real abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB);
276 
277  // Now compute the output error bounds
278  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
279  {
280  RB_output_error_bounds[n] = abs_error_bound * eval_output_dual_norm(n, mu);
281  }
282 
283  return abs_error_bound;
284  }
285  else // Don't calculate the error bounds
286  {
287  // Just return -1. if we did not compute the error bound
288  return -1.;
289  }
290 }
291 
293 {
294  // Normalize the error based on the error bound in the
295  // case of an empty reduced basis. The error bound is based
296  // on the residual F - AU, so with an empty basis this gives
297  // a value based on the norm of F at the current parameters.
298 
299  Real normalization = rb_solve(0);
300  return normalization;
301 }
302 
304 {
305  LOG_SCOPE("compute_residual_dual_norm()", "RBEvaluation");
306 
307  const RBParameters & mu = get_parameters();
308 
309  // Use the stored representor inner product values
310  // to evaluate the residual norm
311  Number residual_norm_sq = 0.;
312 
313  unsigned int q=0;
314  for (unsigned int q_f1=0; q_f1<rb_theta_expansion->get_n_F_terms(); q_f1++)
315  {
316  for (unsigned int q_f2=q_f1; q_f2<rb_theta_expansion->get_n_F_terms(); q_f2++)
317  {
318  Real delta = (q_f1==q_f2) ? 1. : 2.;
319  residual_norm_sq += delta * libmesh_real(
322 
323  q++;
324  }
325  }
326 
327  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
328  {
329  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
330  {
331  for (unsigned int i=0; i<N; i++)
332  {
333  Real delta = 2.;
334  residual_norm_sq +=
335  delta * libmesh_real( rb_theta_expansion->eval_F_theta(q_f, mu) *
338  }
339  }
340  }
341 
342  q=0;
343  for (unsigned int q_a1=0; q_a1<rb_theta_expansion->get_n_A_terms(); q_a1++)
344  {
345  for (unsigned int q_a2=q_a1; q_a2<rb_theta_expansion->get_n_A_terms(); q_a2++)
346  {
347  Real delta = (q_a1==q_a2) ? 1. : 2.;
348 
349  for (unsigned int i=0; i<N; i++)
350  {
351  for (unsigned int j=0; j<N; j++)
352  {
353  residual_norm_sq +=
355  rb_theta_expansion->eval_A_theta(q_a2, mu) *
357  }
358  }
359 
360  q++;
361  }
362  }
363 
364  if (libmesh_real(residual_norm_sq) < 0.)
365  {
366  // libMesh::out << "Warning: Square of residual norm is negative "
367  // << "in RBSystem::compute_residual_dual_norm()" << std::endl;
368 
369  // Sometimes this is negative due to rounding error,
370  // but when this occurs the error is on the order of 1.e-10,
371  // so shouldn't affect error bound much...
372  residual_norm_sq = std::abs(residual_norm_sq);
373  }
374 
375  return std::sqrt( libmesh_real(residual_norm_sq) );
376 }
377 
379 {
380  // Return a default value of 1, this function should
381  // be overloaded to specify a problem-dependent stability
382  // factor lower bound
383  return 1.;
384 }
385 
387 {
388  // Here we implement the residual scaling for a coercive
389  // problem.
390  return alpha_LB;
391 }
392 
394 {
395  Number output_bound_sq = 0.;
396  unsigned int q=0;
397  for (unsigned int q_l1=0; q_l1<rb_theta_expansion->get_n_output_terms(n); q_l1++)
398  {
399  for (unsigned int q_l2=q_l1; q_l2<rb_theta_expansion->get_n_output_terms(n); q_l2++)
400  {
401  Real delta = (q_l1==q_l2) ? 1. : 2.;
402  output_bound_sq += delta * libmesh_real(
405  q++;
406  }
407  }
408 
409  return libmesh_real(std::sqrt( output_bound_sq ));
410 }
411 
413 {
414 
415  // Clear the Aq_representors
416  for (std::size_t q_a=0; q_a<Aq_representor.size(); q_a++)
417  for (std::size_t i=0; i<Aq_representor[q_a].size(); i++)
418  {
419  delete Aq_representor[q_a][i];
421  }
422 }
423 
424 void RBEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
425  const bool write_binary_data)
426 {
427  LOG_SCOPE("legacy_write_offline_data_to_files()", "RBEvaluation");
428 
429  // Get the number of basis functions
430  unsigned int n_bfs = get_n_basis_functions();
431 
432  // The writing mode: ENCODE for binary, WRITE for ASCII
433  XdrMODE mode = write_binary_data ? ENCODE : WRITE;
434 
435  // The suffix to use for all the files that are written out
436  const std::string suffix = write_binary_data ? ".xdr" : ".dat";
437 
438  if (this->processor_id() == 0)
439  {
440 
441  // Make a directory to store all the data files
442  Utility::mkdir(directory_name.c_str());
443  // if (mkdir(directory_name.c_str(), 0777) == -1)
444  // {
445  // libMesh::out << "In RBEvaluation::write_offline_data_to_files, directory "
446  // << directory_name << " already exists, overwriting contents." << std::endl;
447  // }
448 
449  // First, write out how many basis functions we have generated
450  std::ostringstream file_name;
451  {
452  file_name << directory_name << "/n_bfs" << suffix;
453  Xdr n_bfs_out(file_name.str(), mode);
454  n_bfs_out << n_bfs;
455  n_bfs_out.close();
456  }
457 
458  // Write out the parameter ranges
459  file_name.str("");
460  file_name << directory_name << "/parameter_ranges" << suffix;
461  std::string continuous_param_file_name = file_name.str();
462 
463  // Write out the discrete parameter values
464  file_name.str("");
465  file_name << directory_name << "/discrete_parameter_values" << suffix;
466  std::string discrete_param_file_name = file_name.str();
467 
468  write_parameter_data_to_files(continuous_param_file_name,
469  discrete_param_file_name,
470  write_binary_data);
471 
472  // Write out Fq representor norm data
473  file_name.str("");
474  file_name << directory_name << "/Fq_innerprods" << suffix;
475  Xdr RB_Fq_innerprods_out(file_name.str(), mode);
476  unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
477  for (unsigned int i=0; i<Q_f_hat; i++)
478  {
479  RB_Fq_innerprods_out << Fq_representor_innerprods[i];
480  }
481  RB_Fq_innerprods_out.close();
482 
483  // Write out output data
484  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
485  {
486  file_name.str("");
487  file_name << directory_name << "/output_";
488  file_name << std::setw(3)
489  << std::setprecision(0)
490  << std::setfill('0')
491  << std::right
492  << n;
493 
494  file_name << "_dual_innerprods" << suffix;
495  Xdr output_dual_innerprods_out(file_name.str(), mode);
496 
497  unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
498  for (unsigned int q=0; q<Q_l_hat; q++)
499  {
500  output_dual_innerprods_out << output_dual_innerprods[n][q];
501  }
502  output_dual_innerprods_out.close();
503  }
504 
505 
506  // Write out output data to multiple files
507  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
508  {
509  for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
510  {
511  file_name.str("");
512  file_name << directory_name << "/output_";
513  file_name << std::setw(3)
514  << std::setprecision(0)
515  << std::setfill('0')
516  << std::right
517  << n;
518  file_name << "_";
519  file_name << std::setw(3)
520  << std::setprecision(0)
521  << std::setfill('0')
522  << std::right
523  << q_l;
524  file_name << suffix;
525  Xdr output_n_out(file_name.str(), mode);
526 
527  for (unsigned int j=0; j<n_bfs; j++)
528  {
529  output_n_out << RB_output_vectors[n][q_l](j);
530  }
531  output_n_out.close();
532  }
533  }
534 
536  {
537  // Next write out the inner product matrix
538  file_name.str("");
539  file_name << directory_name << "/RB_inner_product_matrix" << suffix;
540  Xdr RB_inner_product_matrix_out(file_name.str(), mode);
541  for (unsigned int i=0; i<n_bfs; i++)
542  {
543  for (unsigned int j=0; j<n_bfs; j++)
544  {
545  RB_inner_product_matrix_out << RB_inner_product_matrix(i,j);
546  }
547  }
548  RB_inner_product_matrix_out.close();
549  }
550 
551  // Next write out the Fq vectors
552  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
553  {
554  file_name.str("");
555  file_name << directory_name << "/RB_F_";
556  file_name << std::setw(3)
557  << std::setprecision(0)
558  << std::setfill('0')
559  << std::right
560  << q_f;
561  file_name << suffix;
562  Xdr RB_Fq_f_out(file_name.str(), mode);
563 
564  for (unsigned int i=0; i<n_bfs; i++)
565  {
566  RB_Fq_f_out << RB_Fq_vector[q_f](i);
567  }
568  RB_Fq_f_out.close();
569  }
570 
571  // Next write out the Aq matrices
572  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
573  {
574  file_name.str("");
575  file_name << directory_name << "/RB_A_";
576  file_name << std::setw(3)
577  << std::setprecision(0)
578  << std::setfill('0')
579  << std::right
580  << q_a;
581  file_name << suffix;
582  Xdr RB_Aq_a_out(file_name.str(), mode);
583 
584  for (unsigned int i=0; i<n_bfs; i++)
585  {
586  for (unsigned int j=0; j<n_bfs; j++)
587  {
588  RB_Aq_a_out << RB_Aq_vector[q_a](i,j);
589  }
590  }
591  RB_Aq_a_out.close();
592  }
593 
594  // Next write out Fq_Aq representor norm data
595  file_name.str("");
596  file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
597  Xdr RB_Fq_Aq_innerprods_out(file_name.str(), mode);
598 
599  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
600  {
601  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
602  {
603  for (unsigned int i=0; i<n_bfs; i++)
604  {
605  RB_Fq_Aq_innerprods_out << Fq_Aq_representor_innerprods[q_f][q_a][i];
606  }
607  }
608  }
609  RB_Fq_Aq_innerprods_out.close();
610 
611  // Next write out Aq_Aq representor norm data
612  file_name.str("");
613  file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
614  Xdr RB_Aq_Aq_innerprods_out(file_name.str(), mode);
615 
616  unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
617  for (unsigned int i=0; i<Q_a_hat; i++)
618  {
619  for (unsigned int j=0; j<n_bfs; j++)
620  {
621  for (unsigned int l=0; l<n_bfs; l++)
622  {
623  RB_Aq_Aq_innerprods_out << Aq_Aq_representor_innerprods[i][j][l];
624  }
625  }
626  }
627  RB_Aq_Aq_innerprods_out.close();
628 
629  // Also, write out the greedily selected parameters
630  {
631  file_name.str("");
632  file_name << directory_name << "/greedy_params" << suffix;
633  Xdr greedy_params_out(file_name.str(), mode);
634 
635  for (std::size_t i=0; i<greedy_param_list.size(); i++)
636  {
639  for ( ; it != it_end; ++it)
640  {
641  // Need to make a copy of the value so that it's not const
642  // Xdr is not templated on const's
643  Real param_value = it->second;
644  greedy_params_out << param_value;
645  }
646  }
647  greedy_params_out.close();
648  }
649 
650  }
651 }
652 
653 void RBEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
654  bool read_error_bound_data,
655  const bool read_binary_data)
656 {
657  LOG_SCOPE("legacy_read_offline_data_from_files()", "RBEvaluation");
658 
659  // The reading mode: DECODE for binary, READ for ASCII
660  XdrMODE mode = read_binary_data ? DECODE : READ;
661 
662  // The suffix to use for all the files that are written out
663  const std::string suffix = read_binary_data ? ".xdr" : ".dat";
664 
665  // The string stream we'll use to make the file names
666  std::ostringstream file_name;
667 
668  // First, find out how many basis functions we had when Greedy terminated
669  unsigned int n_bfs;
670  {
671  file_name << directory_name << "/n_bfs" << suffix;
672  assert_file_exists(file_name.str());
673 
674  Xdr n_bfs_in(file_name.str(), mode);
675  n_bfs_in >> n_bfs;
676  n_bfs_in.close();
677  }
678  resize_data_structures(n_bfs, read_error_bound_data);
679 
680  // Read in the parameter ranges
681  file_name.str("");
682  file_name << directory_name << "/parameter_ranges" << suffix;
683  std::string continuous_param_file_name = file_name.str();
684 
685  // Read in the discrete parameter values
686  file_name.str("");
687  file_name << directory_name << "/discrete_parameter_values" << suffix;
688  std::string discrete_param_file_name = file_name.str();
689  read_parameter_data_from_files(continuous_param_file_name,
690  discrete_param_file_name,
691  read_binary_data);
692 
693  // Read in output data in multiple files
694  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
695  {
696  for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
697  {
698  file_name.str("");
699  file_name << directory_name << "/output_";
700  file_name << std::setw(3)
701  << std::setprecision(0)
702  << std::setfill('0')
703  << std::right
704  << n;
705  file_name << "_";
706  file_name << std::setw(3)
707  << std::setprecision(0)
708  << std::setfill('0')
709  << std::right
710  << q_l;
711  file_name << suffix;
712  assert_file_exists(file_name.str());
713 
714  Xdr output_n_in(file_name.str(), mode);
715 
716  for (unsigned int j=0; j<n_bfs; j++)
717  {
718  Number value;
719  output_n_in >> value;
720  RB_output_vectors[n][q_l](j) = value;
721  }
722  output_n_in.close();
723  }
724  }
725 
727  {
728  // Next read in the inner product matrix
729  file_name.str("");
730  file_name << directory_name << "/RB_inner_product_matrix" << suffix;
731  assert_file_exists(file_name.str());
732 
733  Xdr RB_inner_product_matrix_in(file_name.str(), mode);
734 
735  for (unsigned int i=0; i<n_bfs; i++)
736  {
737  for (unsigned int j=0; j<n_bfs; j++)
738  {
739  Number value;
740  RB_inner_product_matrix_in >> value;
742  }
743  }
744  RB_inner_product_matrix_in.close();
745  }
746 
747  // Next read in the Fq vectors
748  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
749  {
750  file_name.str("");
751  file_name << directory_name << "/RB_F_";
752  file_name << std::setw(3)
753  << std::setprecision(0)
754  << std::setfill('0')
755  << std::right
756  << q_f;
757  file_name << suffix;
758  assert_file_exists(file_name.str());
759 
760  Xdr RB_Fq_f_in(file_name.str(), mode);
761 
762  for (unsigned int i=0; i<n_bfs; i++)
763  {
764  Number value;
765  RB_Fq_f_in >> value;
766  RB_Fq_vector[q_f](i) = value;
767  }
768  RB_Fq_f_in.close();
769  }
770 
771  // Next read in the Aq matrices
772  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
773  {
774  file_name.str("");
775  file_name << directory_name << "/RB_A_";
776  file_name << std::setw(3)
777  << std::setprecision(0)
778  << std::setfill('0')
779  << std::right
780  << q_a;
781  file_name << suffix;
782  assert_file_exists(file_name.str());
783 
784  Xdr RB_Aq_a_in(file_name.str(), mode);
785 
786  for (unsigned int i=0; i<n_bfs; i++)
787  {
788  for (unsigned int j=0; j<n_bfs; j++)
789  {
790  Number value;
791  RB_Aq_a_in >> value;
792  RB_Aq_vector[q_a](i,j) = value;
793  }
794  }
795  RB_Aq_a_in.close();
796  }
797 
798 
799  if (read_error_bound_data)
800  {
801  // Next read in Fq representor norm data
802  file_name.str("");
803  file_name << directory_name << "/Fq_innerprods" << suffix;
804  assert_file_exists(file_name.str());
805 
806  Xdr RB_Fq_innerprods_in(file_name.str(), mode);
807 
808  unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
809  for (unsigned int i=0; i<Q_f_hat; i++)
810  {
811  RB_Fq_innerprods_in >> Fq_representor_innerprods[i];
812  }
813  RB_Fq_innerprods_in.close();
814 
815  // Read in output data
816  for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
817  {
818  file_name.str("");
819  file_name << directory_name << "/output_";
820  file_name << std::setw(3)
821  << std::setprecision(0)
822  << std::setfill('0')
823  << std::right
824  << n;
825  file_name << "_dual_innerprods" << suffix;
826  assert_file_exists(file_name.str());
827 
828  Xdr output_dual_innerprods_in(file_name.str(), mode);
829 
830  unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
831  for (unsigned int q=0; q<Q_l_hat; q++)
832  {
833  output_dual_innerprods_in >> output_dual_innerprods[n][q];
834  }
835  output_dual_innerprods_in.close();
836  }
837 
838 
839  // Next read in Fq_Aq representor norm data
840  file_name.str("");
841  file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
842  assert_file_exists(file_name.str());
843 
844  Xdr RB_Fq_Aq_innerprods_in(file_name.str(), mode);
845 
846  for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
847  {
848  for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
849  {
850  for (unsigned int i=0; i<n_bfs; i++)
851  {
852  RB_Fq_Aq_innerprods_in >> Fq_Aq_representor_innerprods[q_f][q_a][i];
853  }
854  }
855  }
856  RB_Fq_Aq_innerprods_in.close();
857 
858  // Next read in Aq_Aq representor norm data
859  file_name.str("");
860  file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
861  assert_file_exists(file_name.str());
862 
863  Xdr RB_Aq_Aq_innerprods_in(file_name.str(), mode);
864 
865  unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
866  for (unsigned int i=0; i<Q_a_hat; i++)
867  {
868  for (unsigned int j=0; j<n_bfs; j++)
869  {
870  for (unsigned int l=0; l<n_bfs; l++)
871  {
872  RB_Aq_Aq_innerprods_in >> Aq_Aq_representor_innerprods[i][j][l];
873  }
874  }
875  }
876  RB_Aq_Aq_innerprods_in.close();
877  }
878 
879  // Resize basis_functions even if we don't read them in so that
880  // get_n_bfs() returns the correct value. Initialize the pointers
881  // to NULL
882  set_n_basis_functions(n_bfs);
883  for (std::size_t i=0; i<basis_functions.size(); i++)
884  {
885  if (basis_functions[i])
886  {
887  basis_functions[i]->clear();
888  delete basis_functions[i];
889  }
891  }
892 }
893 
894 void RBEvaluation::assert_file_exists(const std::string & file_name)
895 {
896  if (!std::ifstream(file_name.c_str()))
897  libmesh_error_msg("File missing: " + file_name);
898 }
899 
901  const std::string & directory_name,
902  const bool write_binary_basis_functions)
903 {
904  LOG_SCOPE("write_out_basis_functions()", "RBEvaluation");
905 
906  write_out_vectors(sys,
908  directory_name,
909  "bf",
910  write_binary_basis_functions);
911 }
912 
914  std::vector<NumericVector<Number> *> & vectors,
915  const std::string & directory_name,
916  const std::string & data_name,
917  const bool write_binary_vectors)
918 {
919  LOG_SCOPE("write_out_vectors()", "RBEvaluation");
920 
921  if (this->processor_id() == 0)
922  {
923  // Make a directory to store all the data files
924  Utility::mkdir(directory_name.c_str());
925  }
926 
927  // Make sure processors are synced up before we begin
928  this->comm().barrier();
929 
930  std::ostringstream file_name;
931  const std::string basis_function_suffix = (write_binary_vectors ? ".xdr" : ".dat");
932 
933  file_name << directory_name << "/" << data_name << "_data" << basis_function_suffix;
934  Xdr bf_data(file_name.str(),
935  write_binary_vectors ? ENCODE : WRITE);
936 
937  std::string version("libMesh-" + libMesh::get_io_compatibility_version());
938 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
939  version += " with infinite elements";
940 #endif
941  bf_data.data(version ,"# File Format Identifier");
942 
943  sys.write_header(bf_data, /*(unused arg)*/ version, /*write_additional_data=*/false);
944 
945  // Following EquationSystemsIO::write, we use a temporary numbering (node major)
946  // before writing out the data
948 
949  // Write all vectors at once.
950  {
951  // Note the API wants pointers to constant vectors, hence this...
952  std::vector<const NumericVector<Number> *> bf_out(vectors.begin(),
953  vectors.end());
954  // for (std::size_t i=0; i<vectors.size(); i++)
955  // bf_out.push_back(vectors[i]);
956  sys.write_serialized_vectors (bf_data, bf_out);
957  }
958 
959 
960  // set the current version
961  bf_data.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
962  LIBMESH_MINOR_VERSION,
963  LIBMESH_MICRO_VERSION));
964 
965 
966  // Undo the temporary renumbering
968 }
969 
971  const std::string & directory_name,
972  const bool read_binary_basis_functions)
973 {
974  LOG_SCOPE("read_in_basis_functions()", "RBEvaluation");
975 
976  read_in_vectors(sys,
978  directory_name,
979  "bf",
980  read_binary_basis_functions);
981 }
982 
984  std::vector<NumericVector<Number> *> & vectors,
985  const std::string & directory_name,
986  const std::string & data_name,
987  const bool read_binary_vectors)
988 {
989  std::vector<std::vector<NumericVector<Number> *> *> vectors_vec;
990  vectors_vec.push_back(&vectors);
991 
992  std::vector<std::string> directory_name_vec;
993  directory_name_vec.push_back(directory_name);
994 
995  std::vector<std::string> data_name_vec;
996  data_name_vec.push_back(data_name);
997 
999  vectors_vec,
1000  directory_name_vec,
1001  data_name_vec,
1002  read_binary_vectors);
1003 }
1004 
1006  std::vector<std::vector<NumericVector<Number> *> *> multiple_vectors,
1007  const std::vector<std::string> & multiple_directory_names,
1008  const std::vector<std::string> & multiple_data_names,
1009  const bool read_binary_vectors)
1010 {
1011  LOG_SCOPE("read_in_vectors_from_multiple_files()", "RBEvaluation");
1012 
1013  unsigned int n_files = multiple_vectors.size();
1014  unsigned int n_directories = multiple_directory_names.size();
1015  libmesh_assert((n_files == n_directories) && (n_files == multiple_data_names.size()));
1016 
1017  if (n_files == 0)
1018  return;
1019 
1020  // Make sure processors are synced up before we begin
1021  this->comm().barrier();
1022 
1023  std::ostringstream file_name;
1024  const std::string basis_function_suffix = (read_binary_vectors ? ".xdr" : ".dat");
1025 
1026  // Following EquationSystemsIO::read, we use a temporary numbering (node major)
1027  // before writing out the data. For the sake of efficiency, we do this once for
1028  // all the vectors that we read in.
1030 
1031  for (unsigned int data_index=0; data_index<n_directories; data_index++)
1032  {
1033  std::vector<NumericVector<Number> *> & vectors = *multiple_vectors[data_index];
1034 
1035  // Allocate storage for each vector
1036  for (std::size_t i=0; i<vectors.size(); i++)
1037  {
1038  // vectors should all be NULL, otherwise we get a memory leak when
1039  // we create the new vectors in RBEvaluation::read_in_vectors.
1040  if (vectors[i])
1041  libmesh_error_msg("Non-NULL vector passed to read_in_vectors_from_multiple_files");
1042 
1043  vectors[i] = NumericVector<Number>::build(sys.comm()).release();
1044 
1045  vectors[i]->init (sys.n_dofs(),
1046  sys.n_local_dofs(),
1047  false,
1048  PARALLEL);
1049  }
1050 
1051  file_name.str("");
1052  file_name << multiple_directory_names[data_index]
1053  << "/" << multiple_data_names[data_index]
1054  << "_data" << basis_function_suffix;
1055 
1056  // On processor zero check to be sure the file exists
1057  if (this->processor_id() == 0)
1058  {
1059  struct stat stat_info;
1060  int stat_result = stat(file_name.str().c_str(), &stat_info);
1061 
1062  if (stat_result != 0)
1063  libmesh_error_msg("File does not exist: " << file_name.str());
1064  }
1065 
1066  assert_file_exists(file_name.str());
1067  Xdr vector_data(file_name.str(),
1068  read_binary_vectors ? DECODE : READ);
1069 
1070  // Read the header data. This block of code is based on EquationSystems::_read_impl.
1071  {
1072  std::string version;
1073  vector_data.data(version);
1074 
1075  const std::string libMesh_label = "libMesh-";
1076  std::string::size_type lm_pos = version.find(libMesh_label);
1077  if (lm_pos==std::string::npos)
1078  {
1079  libmesh_error_msg("version info missing in Xdr header");
1080  }
1081 
1082  std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
1083  int ver_major = 0, ver_minor = 0, ver_patch = 0;
1084  char dot;
1085  iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
1086  vector_data.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));
1087 
1088  // Actually read the header data. When we do this, set read_header=false
1089  // so taht we do not reinit sys, since we assume that it has already been
1090  // set up properly (e.g. the appropriate variables have already been added).
1091  sys.read_header(vector_data, version, /*read_header=*/false, /*read_additional_data=*/false);
1092  }
1093 
1094  sys.read_serialized_vectors (vector_data, vectors);
1095  }
1096 
1097  // Undo the temporary renumbering
1099 }
1100 
1101 } // namespace libMesh
T libmesh_real(T a)
std::vector< NumericVector< Number > * > basis_functions
The libMesh vectors storing the finite element coefficients of the RB basis functions.
void data(T &a, const char *comment="")
Inputs or outputs a single value.
Definition: xdr_cxx.C:761
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool evaluate_RB_error_bound
Boolean to indicate whether we evaluate a posteriori error bounds when rb_solve is called...
virtual void write_out_basis_functions(System &sys, const std::string &directory_name="offline_data", const bool write_binary_basis_functions=true)
Write out all the basis functions to file.
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
Adds factor times vec to this vector.
Definition: dense_vector.h:431
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
virtual Number eval_A_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_a at the current parameter.
T libmesh_conj(T a)
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true)
Resize and clear the data vectors corresponding to the value of Nmax.
DenseVector< Number > RB_solution
The RB solution vector.
void write_parameter_data_to_files(const std::string &continuous_param_file_name, const std::string &discrete_param_file_name, const bool write_binary_data)
Write out the parameter ranges to files.
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
std::vector< std::vector< NumericVector< Number > * > > Aq_representor
Vector storing the residual representors associated with the left-hand side.
ImplicitSystem & sys
DenseMatrix< Number > RB_inner_product_matrix
The inner product matrix.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
CompareTypes< T, T2 >::supertype dot(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:446
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
void close()
Closes the file if it is open.
Definition: xdr_cxx.C:271
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
virtual void set_n_basis_functions(unsigned int n_bfs)
Set the number of basis functions.
virtual Real compute_residual_dual_norm(const unsigned int N)
Compute the dual norm of the residual for the solution saved in RB_solution_vector.
const class libmesh_nullptr_t libmesh_nullptr
virtual void fix_broken_node_and_element_numbering()=0
There is no reason for a user to ever call this function.
int mkdir(const char *pathname)
Create a directory.
Definition: utility.C:140
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
virtual Real get_error_bound_normalization()
void read_parameter_data_from_files(const std::string &continuous_param_file_name, const std::string &discrete_param_file_name, const bool read_binary_data)
Read in the parameter ranges from files.
The libMesh namespace provides an interface to certain functionality in the library.
std::vector< RBParameters > greedy_param_list
The list of parameters selected by the Greedy algorithm in generating the Reduced Basis associated wi...
std::size_t write_serialized_vectors(Xdr &io, const std::vector< const NumericVector< Number > * > &vectors) const
Serialize & write a number of identically distributed vectors.
Definition: system_io.C:2330
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual Number eval_output_theta(unsigned int output_index, unsigned int q_l, const RBParameters &mu)
Evaluate theta_q_l at the current parameter.
std::vector< DenseVector< Number > > RB_Fq_vector
Dense vector for the RHS.
const MeshBase & get_mesh() const
Definition: system.h:2014
void write_header(Xdr &io, const std::string &version, const bool write_additional_data) const
Writes the basic data header for this System.
Definition: system_io.C:1313
std::string get_io_compatibility_version()
Specifier for I/O file compatibility features.
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:32
void read_header(Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
Reads the basic data header for this System.
Definition: system_io.C:115
virtual Real get_stability_lower_bound()
Get a lower bound for the stability constant (e.g.
virtual void legacy_write_offline_data_to_files(const std::string &directory_name="offline_data", const bool write_binary_data=true)
Write out all the data to text files in order to segregate the Offline stage from the Online stage...
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", bool read_error_bound_data=true, const bool read_binary_data=true)
Read in the saved Offline reduced basis data to initialize the system for Online solves.
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
NumericVector< Number > & get_basis_function(unsigned int i)
Get a reference to the i^th basis function.
void read_in_vectors(System &sys, std::vector< NumericVector< Number > * > &vectors, const std::string &directory_name, const std::string &data_name, const bool read_binary_vectors)
Same as read_in_basis_functions, except in this case we pass in the vectors to be written...
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
void barrier() const
Pause execution until all processors reach a certain point.
std::vector< Real > RB_output_error_bounds
virtual void read_in_basis_functions(System &sys, const std::string &directory_name="offline_data", const bool read_binary_basis_functions=true)
Read in all the basis functions from file.
virtual ~RBEvaluation()
Destructor.
Definition: rb_evaluation.C:55
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
void globally_renumber_nodes_and_elements(MeshBase &)
There is no reason for a user to ever call this function.
Definition: mesh_tools.C:1968
This class forms the base class for all other classes that are expected to be implemented in parallel...
virtual Real residual_scaling_denom(Real alpha_LB)
Specifies the residual scaling on the denominator to be used in the a posteriori error bound...
void set_rb_theta_expansion(RBThetaExpansion &rb_theta_expansion_in)
Set the RBThetaExpansion object.
Definition: rb_evaluation.C:84
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
std::vector< std::vector< std::vector< Number > > > Fq_Aq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
virtual Real rb_solve(unsigned int N)
Perform online solve with the N RB basis functions, for the set of parameters in current_params, where 0 <= N <= RB_size.
dof_id_type n_local_dofs() const
Definition: system.C:185
RBEvaluation(const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
Definition: rb_evaluation.C:45
unsigned int get_n_outputs() const
Get n_outputs, the number output functionals.
virtual Number eval_F_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_f at the current parameter.
virtual void clear_riesz_representors()
Clear all the Riesz representors that are used to compute the RB residual (and hence error bound)...
virtual void write_out_vectors(System &sys, std::vector< NumericVector< Number > * > &vectors, const std::string &directory_name="offline_data", const std::string &data_name="bf", const bool write_binary_basis_functions=true)
Same as write_out_basis_functions, except in this case we pass in the vectors to be written...
std::vector< DenseMatrix< Number > > RB_Aq_vector
Dense matrices for the RB computations.
virtual void clear() libmesh_override
Clear this RBEvaluation object.
Definition: rb_evaluation.C:60
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real eval_output_dual_norm(unsigned int n, const RBParameters &mu)
Evaluate the dual norm of output n for the current parameters.
void assert_file_exists(const std::string &file_name)
Helper function that checks if file_name exists.
const Parallel::Communicator & comm() const
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
std::size_t read_serialized_vectors(Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
Read a number of identically distributed vectors.
Definition: system_io.C:2236
dof_id_type n_dofs() const
Definition: system.C:148
bool is_rb_theta_expansion_initialized() const
Definition: rb_evaluation.C:97
std::vector< Number > RB_outputs
The vectors storing the RB output values and corresponding error bounds.
const RBParameters & get_parameters() const
Get the current parameters.
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
std::map< std::string, Real >::const_iterator const_iterator
Definition: rb_parameters.h:57
RBThetaExpansion * rb_theta_expansion
A pointer to to the object that stores the theta expansion.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the rb_theta_expansion.
Definition: rb_evaluation.C:89
void read_in_vectors_from_multiple_files(System &sys, std::vector< std::vector< NumericVector< Number > * > * > multiple_vectors, const std::vector< std::string > &multiple_directory_names, const std::vector< std::string > &multiple_data_names, const bool read_binary_vectors)
Performs read_in_vectors for a list of directory names and data names.
std::vector< std::vector< std::vector< Number > > > Aq_Aq_representor_innerprods
processor_id_type processor_id() const
std::vector< std::vector< DenseVector< Number > > > RB_output_vectors
The vectors storing the RB output vectors.