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