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 : // Currently, the RBSCMConstruction should only be available
24 : // if SLEPc support is enabled.
25 : #if defined(LIBMESH_HAVE_SLEPC) && (LIBMESH_HAVE_GLPK)
26 :
27 : #include "libmesh/rb_scm_construction.h"
28 : #include "libmesh/rb_construction.h"
29 : #include "libmesh/rb_scm_evaluation.h"
30 :
31 : #include "libmesh/libmesh_logging.h"
32 : #include "libmesh/numeric_vector.h"
33 : #include "libmesh/sparse_matrix.h"
34 : #include "libmesh/equation_systems.h"
35 : #include "libmesh/getpot.h"
36 : #include "libmesh/dof_map.h"
37 : #include "libmesh/enum_eigen_solver_type.h"
38 :
39 : // For creating a directory
40 : #include <sys/types.h>
41 : #include <sys/stat.h>
42 : #include <errno.h>
43 :
44 : namespace libMesh
45 : {
46 :
47 2 : RBSCMConstruction::RBSCMConstruction (EquationSystems & es,
48 : const std::string & name_in,
49 2 : const unsigned int number_in)
50 : : Parent(es, name_in, number_in),
51 2 : SCM_training_tolerance(0.5),
52 4 : RB_system_name(""),
53 2 : rb_scm_eval(nullptr)
54 : {
55 :
56 : // set assemble_before_solve flag to false
57 : // so that we control matrix assembly.
58 2 : assemble_before_solve = false;
59 :
60 : // We symmetrize all operators hence use symmetric solvers
61 2 : set_eigenproblem_type(GHEP);
62 :
63 2 : }
64 :
65 4 : RBSCMConstruction::~RBSCMConstruction () = default;
66 :
67 0 : void RBSCMConstruction::clear()
68 : {
69 0 : Parent::clear();
70 0 : }
71 :
72 2 : void RBSCMConstruction::set_rb_scm_evaluation(RBSCMEvaluation & rb_scm_eval_in)
73 : {
74 2 : rb_scm_eval = &rb_scm_eval_in;
75 2 : }
76 :
77 113 : RBSCMEvaluation & RBSCMConstruction::get_rb_scm_evaluation()
78 : {
79 113 : libmesh_error_msg_if(!rb_scm_eval, "Error: RBSCMEvaluation object hasn't been initialized yet");
80 :
81 113 : return *rb_scm_eval;
82 : }
83 :
84 112 : RBThetaExpansion & RBSCMConstruction::get_rb_theta_expansion()
85 : {
86 112 : return get_rb_scm_evaluation().get_rb_theta_expansion();
87 : }
88 :
89 1 : void RBSCMConstruction::process_parameters_file(const std::string & parameters_filename)
90 : {
91 : // First read in data from parameters_filename
92 2 : GetPot infile(parameters_filename);
93 1 : const unsigned int n_training_samples = infile("n_training_samples",1);
94 1 : const bool deterministic_training = infile("deterministic_training",false);
95 :
96 : // Read in training_parameters_random_seed value. This is used to
97 : // seed the RNG when picking the training parameters. By default the
98 : // value is -1, which means use std::time to seed the RNG.
99 1 : unsigned int training_parameters_random_seed_in = static_cast<unsigned int>(-1);
100 1 : training_parameters_random_seed_in = infile("training_parameters_random_seed",
101 : training_parameters_random_seed_in);
102 1 : set_training_random_seed(static_cast<int>(training_parameters_random_seed_in));
103 :
104 : // SCM Greedy termination tolerance
105 1 : const Real SCM_training_tolerance_in = infile("SCM_training_tolerance", SCM_training_tolerance);
106 0 : set_SCM_training_tolerance(SCM_training_tolerance_in);
107 :
108 : // Initialize the parameter ranges and the parameters themselves
109 1 : unsigned int n_continuous_parameters = infile.vector_variable_size("parameter_names");
110 1 : RBParameters mu_min_in;
111 1 : RBParameters mu_max_in;
112 4 : for (unsigned int i=0; i<n_continuous_parameters; i++)
113 : {
114 : // Read in the parameter names
115 3 : std::string param_name = infile("parameter_names", "NONE", i);
116 :
117 : {
118 3 : Real min_val = infile(param_name, 0., 0);
119 3 : mu_min_in.set_value(param_name, min_val);
120 : }
121 :
122 : {
123 3 : Real max_val = infile(param_name, 0., 1);
124 3 : mu_max_in.set_value(param_name, max_val);
125 : }
126 : }
127 :
128 0 : std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
129 :
130 1 : unsigned int n_discrete_parameters = infile.vector_variable_size("discrete_parameter_names");
131 1 : for (unsigned int i=0; i<n_discrete_parameters; i++)
132 : {
133 0 : std::string param_name = infile("discrete_parameter_names", "NONE", i);
134 :
135 0 : unsigned int n_vals_for_param = infile.vector_variable_size(param_name);
136 0 : std::vector<Real> vals_for_param(n_vals_for_param);
137 0 : for (unsigned int j=0; j != n_vals_for_param; j++)
138 0 : vals_for_param[j] = infile(param_name, 0., j);
139 :
140 0 : discrete_parameter_values_in[param_name] = vals_for_param;
141 : }
142 :
143 1 : initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_in);
144 :
145 0 : std::map<std::string,bool> log_scaling;
146 1 : const RBParameters & mu = get_parameters();
147 0 : unsigned int i=0;
148 4 : for (const auto & pr : mu)
149 : {
150 3 : const std::string & param_name = pr.first;
151 3 : log_scaling[param_name] = static_cast<bool>(infile("log_scaling", 0, i++));
152 : }
153 :
154 1 : initialize_training_parameters(this->get_parameters_min(),
155 : this->get_parameters_max(),
156 : n_training_samples,
157 : log_scaling,
158 0 : deterministic_training); // use deterministic parameters
159 1 : }
160 :
161 1 : void RBSCMConstruction::print_info()
162 : {
163 : // Print out info that describes the current setup
164 0 : libMesh::out << std::endl << "RBSCMConstruction parameters:" << std::endl;
165 0 : libMesh::out << "system name: " << this->name() << std::endl;
166 0 : libMesh::out << "SCM Greedy tolerance: " << get_SCM_training_tolerance() << std::endl;
167 1 : if (rb_scm_eval)
168 : {
169 1 : libMesh::out << "A_q operators attached: " << get_rb_theta_expansion().get_n_A_terms() << std::endl;
170 : }
171 : else
172 : {
173 0 : libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
174 : }
175 1 : libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
176 4 : for (const auto & pr : get_parameters())
177 : {
178 3 : const std::string & param_name = pr.first;
179 0 : libMesh::out << "Parameter " << param_name
180 3 : << ": Min = " << get_parameter_min(param_name)
181 3 : << ", Max = " << get_parameter_max(param_name) << std::endl;
182 : }
183 1 : print_discrete_parameter_values();
184 1 : libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
185 0 : libMesh::out << std::endl;
186 1 : }
187 :
188 0 : void RBSCMConstruction::resize_SCM_vectors()
189 : {
190 : // Clear SCM data vectors
191 0 : rb_scm_eval->B_min.clear();
192 0 : rb_scm_eval->B_max.clear();
193 0 : rb_scm_eval->C_J.clear();
194 0 : rb_scm_eval->C_J_stability_vector.clear();
195 0 : for (auto & vec : rb_scm_eval->SCM_UB_vectors)
196 0 : vec.clear();
197 0 : rb_scm_eval->SCM_UB_vectors.clear();
198 :
199 : // Resize the bounding box vectors
200 0 : rb_scm_eval->B_min.resize(get_rb_theta_expansion().get_n_A_terms());
201 0 : rb_scm_eval->B_max.resize(get_rb_theta_expansion().get_n_A_terms());
202 0 : }
203 :
204 45 : void RBSCMConstruction::add_scaled_symm_Aq(unsigned int q_a, Number scalar)
205 : {
206 0 : LOG_SCOPE("add_scaled_symm_Aq()", "RBSCMConstruction");
207 : // Load the operators from the RBConstruction
208 0 : EquationSystems & es = this->get_equation_systems();
209 45 : RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);
210 45 : rb_system.add_scaled_Aq(scalar, q_a, &get_matrix_A(), true);
211 45 : }
212 :
213 1 : void RBSCMConstruction::load_matrix_B()
214 : {
215 : // Load the operators from the RBConstruction
216 0 : EquationSystems & es = this->get_equation_systems();
217 1 : RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);
218 :
219 1 : matrix_B->zero();
220 1 : matrix_B->close();
221 1 : matrix_B->add(1.,*rb_system.get_inner_product_matrix());
222 1 : }
223 :
224 1 : void RBSCMConstruction::perform_SCM_greedy()
225 : {
226 0 : LOG_SCOPE("perform_SCM_greedy()", "RBSCMConstruction");
227 :
228 : // initialize rb_scm_eval's parameters
229 1 : rb_scm_eval->initialize_parameters(*this);
230 :
231 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
232 : // Get a list of constrained dofs from rb_system
233 0 : std::set<dof_id_type> constrained_dofs_set;
234 0 : EquationSystems & es = this->get_equation_systems();
235 1 : RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);
236 :
237 677 : for (dof_id_type i=0; i<rb_system.n_dofs(); i++)
238 : {
239 650 : if (rb_system.get_dof_map().is_constrained_dof(i))
240 : {
241 26 : constrained_dofs_set.insert(i);
242 : }
243 : }
244 :
245 : // Use these constrained dofs to identify which dofs we want to "get rid of"
246 : // (i.e. condense) in our eigenproblems.
247 1 : this->initialize_condensed_dofs(constrained_dofs_set);
248 : #endif // LIBMESH_ENABLE_CONSTRAINTS
249 :
250 : // Copy the inner product matrix over from rb_system to be used as matrix_B
251 1 : load_matrix_B();
252 :
253 1 : attach_deflation_space();
254 :
255 1 : compute_SCM_bounding_box();
256 : // This loads the new parameter into current_parameters
257 1 : enrich_C_J(0);
258 :
259 0 : unsigned int SCM_iter=0;
260 : while (true)
261 : {
262 : // matrix_A is reinitialized for the current parameters
263 : // on each call to evaluate_stability_constant
264 7 : evaluate_stability_constant();
265 :
266 7 : std::pair<unsigned int,Real> SCM_error_pair = compute_SCM_bounds_on_training_set();
267 :
268 0 : libMesh::out << "SCM iteration " << SCM_iter
269 0 : << ", max_SCM_error = " << SCM_error_pair.second << std::endl;
270 :
271 7 : if (SCM_error_pair.second < SCM_training_tolerance)
272 : {
273 0 : libMesh::out << std::endl << "SCM tolerance of " << SCM_training_tolerance << " reached."
274 0 : << std::endl << std::endl;
275 1 : break;
276 : }
277 :
278 : // If we need another SCM iteration, then enrich C_J
279 6 : enrich_C_J(SCM_error_pair.first);
280 :
281 0 : libMesh::out << std::endl << "-----------------------------------" << std::endl << std::endl;
282 :
283 6 : SCM_iter++;
284 6 : }
285 1 : }
286 :
287 1 : void RBSCMConstruction::compute_SCM_bounding_box()
288 : {
289 0 : LOG_SCOPE("compute_SCM_bounding_box()", "RBSCMConstruction");
290 :
291 : // Resize the bounding box vectors
292 1 : rb_scm_eval->B_min.resize(get_rb_theta_expansion().get_n_A_terms());
293 1 : rb_scm_eval->B_max.resize(get_rb_theta_expansion().get_n_A_terms());
294 :
295 4 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
296 : {
297 3 : matrix_A->zero();
298 3 : add_scaled_symm_Aq(q, 1.);
299 :
300 : // Compute B_min(q)
301 0 : eigen_solver->set_position_of_spectrum(SMALLEST_REAL);
302 3 : set_eigensolver_properties(q);
303 :
304 3 : solve();
305 : // TODO: Assert convergence for eigensolver
306 :
307 0 : unsigned int nconv = get_n_converged();
308 3 : if (nconv != 0)
309 : {
310 3 : std::pair<Real, Real> eval = get_eigenpair(0);
311 :
312 : // ensure that the eigenvalue is real
313 0 : libmesh_assert_less (eval.second, TOLERANCE);
314 :
315 3 : rb_scm_eval->set_B_min(q, eval.first);
316 3 : libMesh::out << std::endl << "B_min("<<q<<") = " << rb_scm_eval->get_B_min(q) << std::endl;
317 : }
318 : else
319 0 : libmesh_error_msg("Eigen solver for computing B_min did not converge");
320 :
321 : // Compute B_max(q)
322 0 : eigen_solver->set_position_of_spectrum(LARGEST_REAL);
323 3 : set_eigensolver_properties(q);
324 :
325 3 : solve();
326 : // TODO: Assert convergence for eigensolver
327 :
328 0 : nconv = get_n_converged();
329 3 : if (nconv != 0)
330 : {
331 3 : std::pair<Real, Real> eval = get_eigenpair(0);
332 :
333 : // ensure that the eigenvalue is real
334 0 : libmesh_assert_less (eval.second, TOLERANCE);
335 :
336 3 : rb_scm_eval->set_B_max(q,eval.first);
337 3 : libMesh::out << "B_max("<<q<<") = " << rb_scm_eval->get_B_max(q) << std::endl;
338 : }
339 : else
340 0 : libmesh_error_msg("Eigen solver for computing B_max did not converge");
341 : }
342 1 : }
343 :
344 7 : void RBSCMConstruction::evaluate_stability_constant()
345 : {
346 0 : LOG_SCOPE("evaluate_stability_constant()", "RBSCMConstruction");
347 :
348 : // Get current index of C_J
349 7 : const unsigned int j = cast_int<unsigned int>(rb_scm_eval->C_J.size()-1);
350 :
351 0 : eigen_solver->set_position_of_spectrum(SMALLEST_REAL);
352 :
353 : // We assume B is set in system assembly
354 : // For coercive problems, B is set to the inner product matrix
355 : // For non-coercive time-dependent problems, B is set to the mass matrix
356 :
357 : // Set matrix A corresponding to mu_star
358 7 : matrix_A->zero();
359 28 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
360 : {
361 21 : add_scaled_symm_Aq(q, get_rb_theta_expansion().eval_A_theta(q,get_parameters()));
362 : }
363 :
364 7 : set_eigensolver_properties(-1);
365 7 : solve();
366 : // TODO: Assert convergence for eigensolver
367 :
368 0 : unsigned int nconv = get_n_converged();
369 7 : if (nconv != 0)
370 : {
371 7 : std::pair<Real, Real> eval = get_eigenpair(0);
372 :
373 : // ensure that the eigenvalue is real
374 0 : libmesh_assert_less (eval.second, TOLERANCE);
375 :
376 : // Store the coercivity constant corresponding to mu_star
377 7 : rb_scm_eval->set_C_J_stability_constraint(j,eval.first);
378 0 : libMesh::out << std::endl << "Stability constant for C_J("<<j<<") = "
379 7 : << rb_scm_eval->get_C_J_stability_constraint(j) << std::endl << std::endl;
380 :
381 : // Compute and store the vector y = (y_1, \ldots, y_Q) for the
382 : // eigenvector currently stored in eigen_system.solution.
383 : // We use this later to compute the SCM upper bounds.
384 7 : Real norm_B2 = libmesh_real( B_inner_product(*solution, *solution) );
385 :
386 28 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
387 : {
388 21 : Real norm_Aq2 = libmesh_real( Aq_inner_product(q, *solution, *solution) );
389 :
390 21 : rb_scm_eval->set_SCM_UB_vector(j,q,norm_Aq2/norm_B2);
391 : }
392 : }
393 : else
394 0 : libmesh_error_msg("Error: Eigensolver did not converge in evaluate_stability_constant");
395 7 : }
396 :
397 7 : Number RBSCMConstruction::B_inner_product(const NumericVector<Number> & v,
398 : const NumericVector<Number> & w) const
399 : {
400 7 : matrix_B->vector_mult(*inner_product_storage_vector, w);
401 :
402 7 : return v.dot(*inner_product_storage_vector);
403 : }
404 :
405 21 : Number RBSCMConstruction::Aq_inner_product(unsigned int q,
406 : const NumericVector<Number> & v,
407 : const NumericVector<Number> & w)
408 : {
409 21 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
410 : "Error: We must have q < Q_a in Aq_inner_product.");
411 :
412 21 : matrix_A->zero();
413 21 : add_scaled_symm_Aq(q, 1.);
414 21 : matrix_A->vector_mult(*inner_product_storage_vector, w);
415 :
416 21 : return v.dot(*inner_product_storage_vector);
417 : }
418 :
419 7 : std::pair<unsigned int,Real> RBSCMConstruction::compute_SCM_bounds_on_training_set()
420 : {
421 0 : LOG_SCOPE("compute_SCM_bounds_on_training_set()", "RBSCMConstruction");
422 :
423 : // Now compute the maximum bound error over training_parameters
424 0 : unsigned int new_C_J_index = 0;
425 0 : Real max_SCM_error = 0.;
426 :
427 7 : numeric_index_type first_index = get_first_local_training_index();
428 707 : for (unsigned int i=0; i<get_local_n_training_samples(); i++)
429 : {
430 700 : set_params_from_training_set(first_index+i);
431 700 : rb_scm_eval->set_parameters( get_parameters() );
432 700 : Real LB = rb_scm_eval->get_SCM_LB();
433 700 : Real UB = rb_scm_eval->get_SCM_UB();
434 :
435 700 : Real error_i = SCM_greedy_error_indicator(LB, UB);
436 :
437 700 : if (error_i > max_SCM_error)
438 : {
439 0 : max_SCM_error = error_i;
440 0 : new_C_J_index = i;
441 : }
442 : }
443 :
444 7 : numeric_index_type global_index = first_index + new_C_J_index;
445 0 : std::pair<numeric_index_type,Real> error_pair(global_index, max_SCM_error);
446 7 : get_global_max_error_pair(this->comm(),error_pair);
447 :
448 7 : return error_pair;
449 : }
450 :
451 7 : void RBSCMConstruction::enrich_C_J(unsigned int new_C_J_index)
452 : {
453 0 : LOG_SCOPE("enrich_C_J()", "RBSCMConstruction");
454 :
455 7 : set_params_from_training_set_and_broadcast(new_C_J_index);
456 :
457 7 : rb_scm_eval->C_J.push_back(get_parameters());
458 :
459 0 : libMesh::out << std::endl << "SCM: Added mu = (";
460 :
461 0 : bool first = true;
462 28 : for (const auto & pr : get_parameters())
463 : {
464 21 : if (!first)
465 0 : libMesh::out << ",";
466 21 : const std::string & param_name = pr.first;
467 21 : RBParameters C_J_params = rb_scm_eval->C_J[rb_scm_eval->C_J.size()-1];
468 21 : libMesh::out << C_J_params.get_value(param_name);
469 0 : first = false;
470 21 : }
471 0 : libMesh::out << ")" << std::endl;
472 :
473 : // Finally, resize C_J_stability_vector and SCM_UB_vectors
474 7 : rb_scm_eval->C_J_stability_vector.push_back(0.);
475 :
476 7 : std::vector<Real> zero_vector(get_rb_theta_expansion().get_n_A_terms());
477 7 : rb_scm_eval->SCM_UB_vectors.push_back(zero_vector);
478 7 : }
479 :
480 :
481 : } // namespace libMesh
482 :
483 : #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
|
