libMesh
rb_scm_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 // Configuration data
21 #include "libmesh/libmesh_config.h"
22 
23 // RBSCMEvaluation should only be available
24 // if SLEPc and GLPK support is enabled.
25 #if defined(LIBMESH_HAVE_SLEPC) && (LIBMESH_HAVE_GLPK)
26 
27 // rbOOmit includes
28 #include "libmesh/rb_scm_evaluation.h"
29 #include "libmesh/rb_theta_expansion.h"
30 
31 // libMesh includes
32 #include "libmesh/libmesh_logging.h"
33 #include "libmesh/numeric_vector.h"
34 #include "libmesh/sparse_matrix.h"
35 #include "libmesh/equation_systems.h"
36 #include "libmesh/getpot.h"
37 #include "libmesh/parallel.h"
38 #include "libmesh/dof_map.h"
39 #include "libmesh/xdr_cxx.h"
40 
41 // For creating a directory
42 #include <sys/types.h>
43 #include <sys/stat.h>
44 #include <errno.h>
45 
46 // glpk includes
47 #include <glpk.h>
48 
49 namespace libMesh
50 {
51 
53  ParallelObject(comm_in)
54 {
55  // Clear SCM data vectors
56  B_min.clear();
57  B_max.clear();
58  C_J.clear();
59  C_J_stability_vector.clear();
60  SCM_UB_vectors.clear();
61 }
62 
64 {
65 }
66 
68 {
69  rb_theta_expansion = &rb_theta_expansion_in;
70 }
71 
73 {
74  if (!rb_theta_expansion)
75  libmesh_error_msg("Error: rb_theta_expansion hasn't been initialized yet");
76 
77  return *rb_theta_expansion;
78 }
79 
80 void RBSCMEvaluation::set_C_J_stability_constraint(unsigned int j, Real stability_const_in)
81 {
82  if (j >= C_J_stability_vector.size())
83  libmesh_error_msg("Error: Input parameter j is too large in set_C_J_stability_constraint.");
84 
85  // we assume that C_J_stability_vector is resized elsewhere
86  // to be the same size as C_J.
87  libmesh_assert_equal_to (C_J_stability_vector.size(), C_J.size());
88 
89  C_J_stability_vector[j] = stability_const_in;
90 }
91 
93 {
94  if (j >= C_J_stability_vector.size())
95  libmesh_error_msg("Error: Input parameter j is too large in get_C_J_stability_constraint.");
96 
97  return C_J_stability_vector[j];
98 }
99 
100 void RBSCMEvaluation::set_SCM_UB_vector(unsigned int j, unsigned int q, Real y_q)
101 {
102  // First make sure that j <= J
103  if (j >= SCM_UB_vectors.size())
104  libmesh_error_msg("Error: We must have j < J in set_SCM_UB_vector.");
105 
106  // Next make sure that q <= Q_a or Q_a_hat
107  if (q >= SCM_UB_vectors[0].size())
108  libmesh_error_msg("Error: q is too large in set_SCM_UB_vector.");
109 
110  SCM_UB_vectors[j][q] = y_q;
111 }
112 
113 Real RBSCMEvaluation::get_SCM_UB_vector(unsigned int j, unsigned int q)
114 {
115  // First make sure that j <= J
116  if (j >= SCM_UB_vectors.size())
117  libmesh_error_msg("Error: We must have j < J in get_SCM_UB_vector.");
118 
119  if (q >= SCM_UB_vectors[0].size())
120  libmesh_error_msg("Error: q is too large in get_SCM_UB_vector.");
121 
122  return SCM_UB_vectors[j][q];
123 }
124 
126 {
127  if (j >= C_J.size())
128  libmesh_error_msg("Error: Input parameter j is too large in get_C_J.");
129 
130  return C_J[j];
131 }
132 
133 Real RBSCMEvaluation::get_B_min(unsigned int q) const
134 {
135  if (q >= B_min.size())
136  libmesh_error_msg("Error: q is too large in get_B_min.");
137 
138  return B_min[q];
139 }
140 
141 
142 Real RBSCMEvaluation::get_B_max(unsigned int q) const
143 {
144  if (q >= B_max.size())
145  libmesh_error_msg("Error: q is too large in get_B_max.");
146 
147  return B_max[q];
148 }
149 
150 void RBSCMEvaluation::set_B_min(unsigned int q, Real B_min_val)
151 {
152  if (q >= B_min.size())
153  libmesh_error_msg("Error: q is too large in set_B_min.");
154 
155  B_min[q] = B_min_val;
156 }
157 
158 void RBSCMEvaluation::set_B_max(unsigned int q, Real B_max_val)
159 {
160  if (q >= B_max.size())
161  libmesh_error_msg("Error: q is too large in set_B_max.");
162 
163  B_max[q] = B_max_val;
164 }
165 
167 {
168  LOG_SCOPE("get_SCM_LB()", "RBSCMEvaluation");
169 
170  // Initialize the LP
171  glp_prob * lp;
172  lp = glp_create_prob();
173  glp_set_obj_dir(lp,GLP_MIN);
174 
175  // Add columns to the LP: corresponds to
176  // the variables y_1,...y_Q_a.
177  // These are the same for each \mu in the SCM
178  // training set, hence can do this up front.
179  glp_add_cols(lp,rb_theta_expansion->get_n_A_terms());
180 
181  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
182  {
183  if (B_max[q] < B_min[q]) // Invalid bound, set as free variable
184  {
185  // GLPK indexing is not zero based!
186  glp_set_col_bnds(lp, q+1, GLP_FR, 0., 0.);
187  }
188  else
189  {
190  // GLPK indexing is not zero based!
191  glp_set_col_bnds(lp, q+1, GLP_DB, B_min[q], B_max[q]);
192  }
193 
194  // If B_max is not defined, just set lower bounds...
195  // glp_set_col_bnds(lp, q+1, GLP_LO, B_min[q], 0.);
196  }
197 
198 
199  // Add rows to the LP: corresponds to the auxiliary
200  // variables that define the constraints at each
201  // mu \in C_J_M
202  unsigned int n_rows = C_J.size();
203  glp_add_rows(lp, n_rows);
204 
205  // Now put current_parameters in saved_parameters
207 
208  unsigned int matrix_size = n_rows*rb_theta_expansion->get_n_A_terms();
209  std::vector<int> ia(matrix_size+1);
210  std::vector<int> ja(matrix_size+1);
211  std::vector<double> ar(matrix_size+1);
212  unsigned int count=0;
213  for (unsigned int m=0; m<n_rows; m++)
214  {
216 
217  // Set the lower bound on the auxiliary variable
218  // due to the stability constant at mu_index
219  glp_set_row_bnds(lp, m+1, GLP_LO, C_J_stability_vector[m], 0.);
220 
221  // Now define the matrix that relates the y's
222  // to the auxiliary variables at the current
223  // value of mu.
224  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
225  {
226  count++;
227 
228  ia[count] = m+1;
229  ja[count] = q+1;
230 
231  // This can only handle Reals right now
233  }
234  }
235 
236  // Now load the original parameters back into current_parameters
237  // in order to set the coefficients of the objective function
239 
240  glp_load_matrix(lp, matrix_size, &ia[0], &ja[0], &ar[0]);
241 
242  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
243  {
244  glp_set_obj_coef(lp,q+1, libmesh_real( rb_theta_expansion->eval_A_theta(q,get_parameters()) ) );
245  }
246 
247  // Use this command to initialize the basis for the LP
248  // since default behavior is to use the basis from
249  // the previous solve, but that might become singular
250  // if we switch the order of constraints (as can
251  // happen when we generate a new C_J_M)
252  //lpx_cpx_basis(lp); //glp_cpx_basis(lp);
253 
254  glp_smcp parm;
255  glp_init_smcp(&parm);
256  parm.msg_lev = GLP_MSG_ERR;
257  parm.meth = GLP_DUAL;
258 
259 
260  // use the simplex method and solve the LP
261  glp_simplex(lp, &parm);
262 
263  Real min_J_obj = glp_get_obj_val(lp);
264 
265  // int simplex_status = glp_get_status(lp);
266  // if (simplex_status == GLP_UNBND)
267  // {
268  // libMesh::out << "Simplex method gave unbounded solution." << std::endl;
269  // min_J_obj = std::numeric_limits<Real>::quiet_NaN();
270  // }
271  // else
272  // {
273  // min_J_obj = glp_get_obj_val(lp);
274  // }
275 
276  // Destroy the LP
277  glp_delete_prob(lp);
278 
279  return min_J_obj;
280 }
281 
283 {
284  LOG_SCOPE("get_SCM_UB()", "RBSCMEvaluation");
285 
286  // Add rows to the LP: corresponds to the auxiliary
287  // variables that define the constraints at each
288  // mu \in C_J
289  unsigned int n_rows = C_J.size();
290 
291  // For each mu, we just find the minimum of J_obj over
292  // the subset of vectors in SCM_UB_vectors corresponding
293  // to C_J_M (SCM_UB_vectors contains vectors for all of
294  // C_J).
295  Real min_J_obj = 0.;
296  for (unsigned int m=0; m<n_rows; m++)
297  {
298  const std::vector<Real> UB_vector = SCM_UB_vectors[m];
299 
300  Real J_obj = 0.;
301  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
302  {
303  J_obj += libmesh_real( rb_theta_expansion->eval_A_theta(q,get_parameters()) )*UB_vector[q];
304  }
305 
306  if ((m==0) || (J_obj < min_J_obj))
307  {
308  min_J_obj = J_obj;
309  }
310  }
311 
312  return min_J_obj;
313 }
314 
316 {
317  set_parameters(C_J[C_J_index]);
318 }
319 
321 {
323 }
324 
326 {
328 }
329 
330 void RBSCMEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
331  const bool write_binary_data)
332 {
333  LOG_SCOPE("legacy_write_offline_data_to_files()", "RBSCMEvaluation");
334 
335  if (this->processor_id() == 0)
336  {
337  // Make a directory to store all the data files
338  if (mkdir(directory_name.c_str(), 0777) == -1)
339  {
340  libMesh::out << "In RBSCMEvaluation::write_offline_data_to_files, directory "
341  << directory_name << " already exists, overwriting contents." << std::endl;
342  }
343 
344  // The writing mode: ENCODE for binary, WRITE for ASCII
345  XdrMODE mode = write_binary_data ? ENCODE : WRITE;
346 
347  // The suffix to use for all the files that are written out
348  const std::string suffix = write_binary_data ? ".xdr" : ".dat";
349 
350  // Stream for building the file names
351  std::ostringstream file_name;
352 
353  // Write out the parameter ranges
354  file_name.str("");
355  file_name << directory_name << "/parameter_ranges" << suffix;
356  std::string continuous_param_file_name = file_name.str();
357 
358  // Write out the discrete parameter values
359  file_name.str("");
360  file_name << directory_name << "/discrete_parameter_values" << suffix;
361  std::string discrete_param_file_name = file_name.str();
362 
363  write_parameter_data_to_files(continuous_param_file_name,
364  discrete_param_file_name,
365  write_binary_data);
366 
367  // Write out the bounding box min values
368  file_name.str("");
369  file_name << directory_name << "/B_min" << suffix;
370  Xdr B_min_out(file_name.str(), mode);
371 
372  for (std::size_t i=0; i<B_min.size(); i++)
373  {
374  Real B_min_i = get_B_min(i);
375  B_min_out << B_min_i;
376  }
377  B_min_out.close();
378 
379 
380  // Write out the bounding box max values
381  file_name.str("");
382  file_name << directory_name << "/B_max" << suffix;
383  Xdr B_max_out(file_name.str(), mode);
384 
385  for (std::size_t i=0; i<B_max.size(); i++)
386  {
387  Real B_max_i = get_B_max(i);
388  B_max_out << B_max_i;
389  }
390  B_max_out.close();
391 
392  // Write out the length of the C_J data
393  file_name.str("");
394  file_name << directory_name << "/C_J_length" << suffix;
395  Xdr C_J_length_out(file_name.str(), mode);
396 
397  unsigned int C_J_length = C_J.size();
398  C_J_length_out << C_J_length;
399  C_J_length_out.close();
400 
401  // Write out C_J_stability_vector
402  file_name.str("");
403  file_name << directory_name << "/C_J_stability_vector" << suffix;
404  Xdr C_J_stability_vector_out(file_name.str(), mode);
405 
406  for (std::size_t i=0; i<C_J_stability_vector.size(); i++)
407  {
408  Real C_J_stability_constraint_i = get_C_J_stability_constraint(i);
409  C_J_stability_vector_out << C_J_stability_constraint_i;
410  }
411  C_J_stability_vector_out.close();
412 
413  // Write out C_J
414  file_name.str("");
415  file_name << directory_name << "/C_J" << suffix;
416  Xdr C_J_out(file_name.str(), mode);
417 
418  for (std::size_t i=0; i<C_J.size(); i++)
419  {
420  RBParameters::const_iterator it = C_J[i].begin();
421  RBParameters::const_iterator it_end = C_J[i].end();
422  for ( ; it != it_end; ++it)
423  {
424  // Need to make a copy of the value so that it's not const
425  // Xdr is not templated on const's
426  Real param_value = it->second;
427  C_J_out << param_value;
428  }
429  }
430  C_J_out.close();
431 
432  // Write out SCM_UB_vectors get_SCM_UB_vector
433  file_name.str("");
434  file_name << directory_name << "/SCM_UB_vectors" << suffix;
435  Xdr SCM_UB_vectors_out(file_name.str(), mode);
436 
437  for (std::size_t i=0; i<SCM_UB_vectors.size(); i++)
438  for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
439  {
440  Real SCM_UB_vector_ij = get_SCM_UB_vector(i,j);
441  SCM_UB_vectors_out << SCM_UB_vector_ij;
442  }
443  SCM_UB_vectors_out.close();
444  }
445 }
446 
447 
448 void RBSCMEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
449  const bool read_binary_data)
450 {
451  LOG_SCOPE("legacy_read_offline_data_from_files()", "RBSCMEvaluation");
452 
453  // The reading mode: DECODE for binary, READ for ASCII
454  XdrMODE mode = read_binary_data ? DECODE : READ;
455 
456  // The suffix to use for all the files that are written out
457  const std::string suffix = read_binary_data ? ".xdr" : ".dat";
458 
459  // The string stream we'll use to make the file names
460  std::ostringstream file_name;
461 
462  // Read in the parameter ranges
463  file_name.str("");
464  file_name << directory_name << "/parameter_ranges" << suffix;
465  std::string continuous_param_file_name = file_name.str();
466 
467  // Read in the discrete parameter values
468  file_name.str("");
469  file_name << directory_name << "/discrete_parameter_values" << suffix;
470  std::string discrete_param_file_name = file_name.str();
471  read_parameter_data_from_files(continuous_param_file_name,
472  discrete_param_file_name,
473  read_binary_data);
474 
475  // Read in the bounding box min values
476  // Note that there are Q_a values
477  file_name.str("");
478  file_name << directory_name << "/B_min" << suffix;
479  Xdr B_min_in(file_name.str(), mode);
480 
481  B_min.clear();
482  for (unsigned int i=0; i<rb_theta_expansion->get_n_A_terms(); i++)
483  {
484  Real B_min_val;
485  B_min_in >> B_min_val;
486  B_min.push_back(B_min_val);
487  }
488  B_min_in.close();
489 
490 
491  // Read in the bounding box max values
492  // Note that there are Q_a values
493  file_name.str("");
494  file_name << directory_name << "/B_max" << suffix;
495  Xdr B_max_in(file_name.str(), mode);
496 
497  B_max.clear();
498  for (unsigned int i=0; i<rb_theta_expansion->get_n_A_terms(); i++)
499  {
500  Real B_max_val;
501  B_max_in >> B_max_val;
502  B_max.push_back(B_max_val);
503  }
504 
505  // Read in the length of the C_J data
506  file_name.str("");
507  file_name << directory_name << "/C_J_length" << suffix;
508  Xdr C_J_length_in(file_name.str(), mode);
509 
510  unsigned int C_J_length;
511  C_J_length_in >> C_J_length;
512  C_J_length_in.close();
513 
514  // Read in C_J_stability_vector
515  file_name.str("");
516  file_name << directory_name << "/C_J_stability_vector" << suffix;
517  Xdr C_J_stability_vector_in(file_name.str(), mode);
518 
519  C_J_stability_vector.clear();
520  for (unsigned int i=0; i<C_J_length; i++)
521  {
522  Real C_J_stability_val;
523  C_J_stability_vector_in >> C_J_stability_val;
524  C_J_stability_vector.push_back(C_J_stability_val);
525  }
526  C_J_stability_vector_in.close();
527 
528  // Read in C_J
529  file_name.str("");
530  file_name << directory_name << "/C_J" << suffix;
531  Xdr C_J_in(file_name.str(), mode);
532 
533  // Resize C_J based on C_J_stability_vector and Q_a
534  C_J.resize( C_J_length );
535  for (std::size_t i=0; i<C_J.size(); i++)
536  {
539  for ( ; it != it_end; ++it)
540  {
541  std::string param_name = it->first;
542  Real param_value;
543  C_J_in >> param_value;
544  C_J[i].set_value(param_name, param_value);
545  }
546  }
547  C_J_in.close();
548 
549 
550  // Read in SCM_UB_vectors get_SCM_UB_vector
551  file_name.str("");
552  file_name << directory_name << "/SCM_UB_vectors" << suffix;
553  Xdr SCM_UB_vectors_in(file_name.str(), mode);
554 
555  // Resize SCM_UB_vectors based on C_J_stability_vector and Q_a
556  SCM_UB_vectors.resize( C_J_stability_vector.size() );
557  for (std::size_t i=0; i<SCM_UB_vectors.size(); i++)
558  {
560  for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
561  {
562  SCM_UB_vectors_in >> SCM_UB_vectors[i][j];
563  }
564  }
565  SCM_UB_vectors_in.close();
566 }
567 
568 } // namespace libMesh
569 
570 #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
T libmesh_real(T a)
std::vector< Real > B_min
B_min, B_max define the bounding box.
virtual Real get_SCM_UB()
Evaluate single SCM upper bound.
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
Real get_B_min(unsigned int i) const
Get B_min and B_max.
void set_rb_theta_expansion(RBThetaExpansion &rb_theta_expansion_in)
Set the RBThetaExpansion object.
virtual Number eval_A_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_a at the current parameter.
void set_SCM_UB_vector(unsigned int j, unsigned int q, Real y_q)
Set entries of SCM_UB_vector, which stores the vector y, corresponding to the minimizing eigenvectors...
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.
void close()
Closes the file if it is open.
Definition: xdr_cxx.C:271
int mkdir(const char *pathname)
Create a directory.
Definition: utility.C:140
virtual void save_current_parameters()
Helper function to save current_parameters in saved_parameters.
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.
Real get_B_max(unsigned int i) const
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the rb_theta_expansion.
std::vector< Real > B_max
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", const bool read_binary_data=true)
Read in the saved Offline reduced basis data to initialize the system for Online solves.
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
std::vector< std::vector< Real > > SCM_UB_vectors
This matrix stores the infimizing vectors y_1( ),...,y_Q_a( ), for each in C_J, which are used in co...
std::vector< Real > C_J_stability_vector
Vector storing the (truth) stability values at the parameters in C_J.
virtual Real get_SCM_LB()
Evaluate single SCM lower bound.
const RBParameters & get_C_J_entry(unsigned int j)
Get entry of C_J.
Real get_C_J_stability_constraint(unsigned int j) const
Get stability constraints (i.e.
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:32
void set_parameters(const RBParameters &params)
Set the current parameters to params.
virtual void reload_current_parameters()
Helper function to (re)load current_parameters from saved_parameters.
std::vector< RBParameters > C_J
Vector storing the greedily selected parameters during SCM training.
RBThetaExpansion * rb_theta_expansion
A pointer to to the object that stores the theta expansion.
void set_B_min(unsigned int i, Real B_min_val)
Set B_min and B_max.
void set_C_J_stability_constraint(unsigned int j, Real stability_constraint_in)
Set stability constraints (i.e.
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...
const_iterator end() const
Get a constant iterator to the end of this RBParameters object.
Definition: rb_parameters.C:85
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
This class forms the base class for all other classes that are expected to be implemented in parallel...
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
virtual void set_current_parameters_from_C_J(unsigned int C_J_index)
Set parameters based on values saved in "C_J".
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_B_max(unsigned int i, Real B_max_val)
OStreamProxy out
virtual ~RBSCMEvaluation()
Destructor.
const_iterator begin() const
Get a constant iterator to beginning of this RBParameters object.
Definition: rb_parameters.C:80
const RBParameters & get_parameters() const
Get the current parameters.
RBSCMEvaluation(const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
std::map< std::string, Real >::const_iterator const_iterator
Definition: rb_parameters.h:57
RBParameters saved_parameters
Vector in which to save a parameter set.
Real get_SCM_UB_vector(unsigned int j, unsigned int q)
Get entries of SCM_UB_vector, which stores the vector y, corresponding to the minimizing eigenvectors...
processor_id_type processor_id() const