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 : // libMesh includes
21 : #include "libmesh/int_range.h"
22 : #include "libmesh/simple_range.h"
23 : #include "libmesh/xdr_cxx.h"
24 :
25 : // rbOOmit includes
26 : #include "libmesh/rb_parametrized.h"
27 :
28 : // C++ includes
29 : #include <sstream>
30 : #include <fstream>
31 : #include <algorithm> // std::min_element
32 :
33 : namespace libMesh
34 : {
35 :
36 1211 : RBParametrized::RBParametrized()
37 : :
38 1143 : verbose_mode(false),
39 1211 : parameters_initialized(false)
40 : {
41 1211 : }
42 :
43 1143 : RBParametrized::~RBParametrized() = default;
44 :
45 0 : void RBParametrized::clear()
46 : {
47 0 : parameters.clear();
48 0 : parameters_min.clear();
49 0 : parameters_max.clear();
50 0 : parameters_initialized = false;
51 0 : }
52 :
53 1139 : void RBParametrized::initialize_parameters(const RBParameters & mu_min_in,
54 : const RBParameters & mu_max_in,
55 : const std::map<std::string, std::vector<Real>> & discrete_parameter_values)
56 : {
57 : // Check that the min/max vectors have the same size.
58 1210 : libmesh_error_msg_if(mu_min_in.n_parameters() != mu_max_in.n_parameters(),
59 : "Error: Invalid mu_min/mu_max in initialize_parameters(), different number of parameters.");
60 1139 : libmesh_error_msg_if(mu_min_in.n_samples() != 1 ||
61 : mu_max_in.n_samples() != 1,
62 : "Error: Invalid mu_min/mu_max in initialize_parameters(), only 1 sample supported.");
63 :
64 : // Ensure all the values are valid for min and max.
65 1025 : auto pr_min = mu_min_in.begin_serialized();
66 1025 : auto pr_max = mu_max_in.begin_serialized();
67 9757 : for (; pr_min != mu_min_in.end_serialized(); ++pr_min, ++pr_max)
68 4059 : libmesh_error_msg_if((*pr_min).second > (*pr_max).second,
69 : "Error: Invalid mu_min/mu_max in RBParameters constructor.");
70 :
71 926 : parameters_min = mu_min_in;
72 926 : parameters_max = mu_max_in;
73 :
74 : // Add in min/max values due to the discrete parameters
75 926 : for (const auto & [name, vals] : discrete_parameter_values)
76 : {
77 0 : libmesh_error_msg_if(vals.empty(), "Error: List of discrete parameters for " << name << " is empty.");
78 :
79 0 : Real min_val = *std::min_element(vals.begin(), vals.end());
80 0 : Real max_val = *std::max_element(vals.begin(), vals.end());
81 :
82 0 : libmesh_assert_less_equal(min_val, max_val);
83 :
84 0 : parameters_min.set_value(name, min_val);
85 0 : parameters_max.set_value(name, max_val);
86 : }
87 :
88 26 : _discrete_parameter_values = discrete_parameter_values;
89 :
90 926 : parameters_initialized = true;
91 :
92 : // Initialize the current parameters to parameters_min
93 926 : set_parameters(parameters_min);
94 926 : }
95 :
96 285 : void RBParametrized::initialize_parameters(const RBParametrized & rb_parametrized)
97 : {
98 285 : initialize_parameters(rb_parametrized.get_parameters_min(),
99 : rb_parametrized.get_parameters_max(),
100 : rb_parametrized.get_discrete_parameter_values());
101 285 : }
102 :
103 837125 : unsigned int RBParametrized::get_n_params() const
104 : {
105 837125 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_n_params");
106 :
107 68406 : libmesh_assert_equal_to ( parameters_min.n_parameters(), parameters_max.n_parameters() );
108 :
109 837125 : return parameters_min.n_parameters();
110 : }
111 :
112 49 : unsigned int RBParametrized::get_n_continuous_params() const
113 : {
114 49 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_n_continuous_params");
115 :
116 4 : libmesh_assert(get_n_params() >= get_n_discrete_params());
117 :
118 49 : return static_cast<unsigned int>(get_n_params() - get_n_discrete_params());
119 : }
120 :
121 387 : unsigned int RBParametrized::get_n_discrete_params() const
122 : {
123 387 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_n_discrete_params");
124 :
125 : return cast_int<unsigned int>
126 387 : (get_discrete_parameter_values().size());
127 : }
128 :
129 : #ifdef LIBMESH_ENABLE_DEPRECATED
130 0 : std::set<std::string> RBParametrized::get_parameter_names() const
131 : {
132 : libmesh_deprecated();
133 0 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameter_names");
134 :
135 0 : std::set<std::string> parameter_names;
136 0 : for (const auto & pr : parameters_min)
137 0 : parameter_names.insert(pr.first);
138 :
139 0 : return parameter_names;
140 : }
141 : #endif // LIBMESH_ENABLE_DEPRECATED
142 :
143 831830 : bool RBParametrized::set_parameters(const RBParameters & params)
144 : {
145 831830 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::set_parameters");
146 :
147 : // Terminate if params has the wrong number of parameters or samples.
148 : // If the parameters are outside the min/max range, return false.
149 831830 : const bool valid_params = check_if_valid_params(params);
150 :
151 : // Make a copy of params (default assignment operator just does memberwise copy, which is sufficient here)
152 763580 : this->parameters = params;
153 :
154 831830 : return valid_params;
155 : }
156 :
157 11155637 : const RBParameters & RBParametrized::get_parameters() const
158 : {
159 11155637 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters");
160 :
161 11155637 : return parameters;
162 : }
163 :
164 839 : const RBParameters & RBParametrized::get_parameters_min() const
165 : {
166 839 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters_min");
167 :
168 839 : return parameters_min;
169 : }
170 :
171 839 : const RBParameters & RBParametrized::get_parameters_max() const
172 : {
173 839 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters_max");
174 :
175 839 : return parameters_max;
176 : }
177 :
178 5292592 : Real RBParametrized::get_parameter_min(const std::string & param_name) const
179 : {
180 5292592 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameter_min");
181 :
182 5292592 : return parameters_min.get_value(param_name);
183 : }
184 :
185 5292592 : Real RBParametrized::get_parameter_max(const std::string & param_name) const
186 : {
187 5292592 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameter_max");
188 :
189 5292592 : return parameters_max.get_value(param_name);
190 : }
191 :
192 5240 : void RBParametrized::print_parameters() const
193 : {
194 5240 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::print_current_parameters");
195 :
196 5240 : get_parameters().print();
197 5240 : }
198 :
199 49 : void RBParametrized::write_parameter_data_to_files(const std::string & continuous_param_file_name,
200 : const std::string & discrete_param_file_name,
201 : const bool write_binary_data)
202 : {
203 49 : write_parameter_ranges_to_file(continuous_param_file_name, write_binary_data);
204 49 : write_discrete_parameter_values_to_file(discrete_param_file_name, write_binary_data);
205 49 : }
206 :
207 49 : void RBParametrized::write_parameter_ranges_to_file(const std::string & file_name,
208 : const bool write_binary_data)
209 : {
210 : // The writing mode: ENCODE for binary, WRITE for ASCII
211 49 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
212 :
213 : // Write out the parameter ranges
214 57 : Xdr parameter_ranges_out(file_name, mode);
215 49 : unsigned int n_continuous_params = get_n_continuous_params();
216 8 : parameter_ranges_out << n_continuous_params;
217 :
218 : // Note: the following loops are not candidates for structured
219 : // bindings syntax because the Xdr APIs which they call are not
220 : // defined for const references. We must therefore make copies to
221 : // call these functions.
222 269 : for (const auto & pr : get_parameters_min())
223 : {
224 220 : std::string param_name = pr.first;
225 220 : if (!is_discrete_parameter(param_name))
226 : {
227 220 : Real param_value = get_parameters_min().get_value(param_name);
228 36 : parameter_ranges_out << param_name << param_value;
229 : }
230 : }
231 269 : for (const auto & pr : get_parameters_max())
232 : {
233 220 : std::string param_name = pr.first;
234 220 : if (!is_discrete_parameter(param_name))
235 : {
236 220 : Real param_value = get_parameters_max().get_value(param_name);
237 36 : parameter_ranges_out << param_name << param_value;
238 : }
239 : }
240 :
241 49 : parameter_ranges_out.close();
242 49 : }
243 :
244 49 : void RBParametrized::write_discrete_parameter_values_to_file(const std::string & file_name,
245 : const bool write_binary_data)
246 : {
247 : // write out the discrete parameters, if we have any
248 49 : if (get_n_discrete_params() > 0)
249 : {
250 : // The writing mode: ENCODE for binary, WRITE for ASCII
251 0 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
252 :
253 0 : Xdr discrete_parameters_out(file_name, mode);
254 0 : unsigned int n_discrete_params = get_n_discrete_params();
255 0 : discrete_parameters_out << n_discrete_params;
256 :
257 : // Note: the following loops are not candidates for structured
258 : // bindings syntax because the Xdr APIs which they call are not
259 : // defined for const references. We must therefore make copies
260 : // to call these functions.
261 0 : for (const auto & pr : get_discrete_parameter_values())
262 : {
263 0 : std::string param_name = pr.first;
264 0 : auto n_discrete_values = cast_int<unsigned int>(pr.second.size());
265 0 : discrete_parameters_out << param_name << n_discrete_values;
266 :
267 0 : for (unsigned int i=0; i<n_discrete_values; i++)
268 : {
269 0 : Real discrete_value = pr.second[i];
270 0 : discrete_parameters_out << discrete_value;
271 : }
272 : }
273 0 : }
274 49 : }
275 :
276 285 : void RBParametrized::read_parameter_data_from_files(const std::string & continuous_param_file_name,
277 : const std::string & discrete_param_file_name,
278 : const bool read_binary_data)
279 : {
280 301 : RBParameters param_min;
281 301 : RBParameters param_max;
282 285 : read_parameter_ranges_from_file(continuous_param_file_name,
283 : read_binary_data,
284 : param_min,
285 : param_max);
286 :
287 16 : std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
288 285 : read_discrete_parameter_values_from_file(discrete_param_file_name,
289 : read_binary_data,
290 : discrete_parameter_values_in);
291 :
292 285 : initialize_parameters(param_min, param_max, discrete_parameter_values_in);
293 285 : }
294 :
295 285 : void RBParametrized::read_parameter_ranges_from_file(const std::string & file_name,
296 : const bool read_binary_data,
297 : RBParameters & param_min,
298 : RBParameters & param_max)
299 : {
300 : // The reading mode: DECODE for binary, READ for ASCII
301 285 : XdrMODE mode = read_binary_data ? DECODE : READ;
302 :
303 : // Read in the parameter ranges
304 570 : Xdr parameter_ranges_in(file_name, mode);
305 : unsigned int n_continuous_params;
306 16 : parameter_ranges_in >> n_continuous_params;
307 :
308 1567 : for (unsigned int i=0; i<n_continuous_params; i++)
309 : {
310 72 : std::string param_name;
311 : Real param_value;
312 :
313 72 : parameter_ranges_in >> param_name;
314 72 : parameter_ranges_in >> param_value;
315 :
316 1282 : param_min.set_value(param_name, param_value);
317 : }
318 1567 : for (unsigned int i=0; i<n_continuous_params; i++)
319 : {
320 72 : std::string param_name;
321 : Real param_value;
322 :
323 72 : parameter_ranges_in >> param_name;
324 72 : parameter_ranges_in >> param_value;
325 :
326 1282 : param_max.set_value(param_name, param_value);
327 : }
328 :
329 285 : parameter_ranges_in.close();
330 285 : }
331 :
332 285 : void RBParametrized::read_discrete_parameter_values_from_file(const std::string & file_name,
333 : const bool read_binary_data,
334 : std::map<std::string, std::vector<Real>> & discrete_parameter_values)
335 : {
336 : // read in the discrete parameters, if we have any
337 301 : std::ifstream check_if_file_exists(file_name.c_str());
338 285 : if (check_if_file_exists.good())
339 : {
340 : // The reading mode: DECODE for binary, READ for ASCII
341 0 : XdrMODE mode = read_binary_data ? DECODE : READ;
342 :
343 : // Read in the parameter ranges
344 0 : Xdr discrete_parameter_values_in(file_name, mode);
345 : unsigned int n_discrete_params;
346 0 : discrete_parameter_values_in >> n_discrete_params;
347 :
348 0 : for (unsigned int i=0; i<n_discrete_params; i++)
349 : {
350 0 : std::string param_name;
351 0 : discrete_parameter_values_in >> param_name;
352 :
353 : unsigned int n_discrete_values;
354 0 : discrete_parameter_values_in >> n_discrete_values;
355 :
356 0 : std::vector<Real> discrete_values(n_discrete_values);
357 0 : for (auto & val : discrete_values)
358 0 : discrete_parameter_values_in >> val;
359 :
360 0 : discrete_parameter_values[param_name] = discrete_values;
361 : }
362 0 : }
363 285 : }
364 :
365 1719 : bool RBParametrized::is_discrete_parameter(const std::string & mu_name) const
366 : {
367 1719 : libmesh_error_msg_if(!parameters_initialized,
368 : "Error: parameters not initialized in RBParametrized::is_discrete_parameter");
369 :
370 1719 : return _discrete_parameter_values.count(mu_name);
371 : }
372 :
373 10584429 : const std::map<std::string, std::vector<Real>> & RBParametrized::get_discrete_parameter_values() const
374 : {
375 10584429 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_discrete_parameter_values");
376 :
377 10584429 : return _discrete_parameter_values;
378 : }
379 :
380 285 : void RBParametrized::print_discrete_parameter_values() const
381 : {
382 293 : for (const auto & [name, values] : get_discrete_parameter_values())
383 : {
384 0 : libMesh::out << "Discrete parameter " << name << ", values: ";
385 :
386 0 : for (const auto & value : values)
387 0 : libMesh::out << value << " ";
388 0 : libMesh::out << std::endl;
389 : }
390 285 : }
391 :
392 831830 : bool RBParametrized::check_if_valid_params(const RBParameters ¶ms) const
393 : {
394 : // Check if number of parameters are correct.
395 831830 : libmesh_error_msg_if(params.n_parameters() != get_n_params(),
396 : "Error: Number of parameters don't match; found "
397 : << params.n_parameters() << ", expected "
398 : << get_n_params());
399 :
400 68250 : bool is_valid = true;
401 831830 : std::string prev_param_name = "";
402 1705874 : for (const auto & [param_name, sample_vec] : params)
403 : {
404 437022 : std::size_t sample_idx = 0;
405 5291310 : const Real & min_value = get_parameter_min(param_name);
406 5291310 : const Real & max_value = get_parameter_max(param_name);
407 10583046 : for (const auto & value_vec : sample_vec)
408 : {
409 10583472 : for (const auto & value : value_vec)
410 : {
411 : // Check every parameter value (including across samples and vector-values)
412 : // to ensure it's within the min/max range.
413 5291736 : const bool outside_range = ((value < min_value) || (value > max_value));
414 5291736 : is_valid = is_valid && !outside_range;
415 5291736 : if (outside_range && verbose_mode)
416 : {
417 0 : libMesh::out << "Warning: parameter " << param_name << " value="
418 0 : << value << " outside acceptable range: ("
419 0 : << min_value << ", " << max_value << ")";
420 : }
421 :
422 : // For discrete params, make sure params.get_value(param_name) is sufficiently
423 : // close to one of the discrete parameter values.
424 : // Note that vector-values not yet supported in discrete parameters,
425 : // and the .get_sample_value() call will throw an error if the user
426 : // tries to do it.
427 10146438 : if (const auto it = get_discrete_parameter_values().find(param_name);
428 5291736 : it != get_discrete_parameter_values().end())
429 : {
430 : const bool is_value_discrete =
431 0 : is_value_in_list(params.get_sample_value(param_name, sample_idx),
432 0 : it->second,
433 : TOLERANCE);
434 0 : is_valid = is_valid && is_value_discrete;
435 0 : if (!is_value_discrete && verbose_mode)
436 0 : libMesh::out << "Warning: parameter " << param_name << " value="
437 0 : << value << " is not in discrete value list.";
438 : }
439 : }
440 5291736 : ++sample_idx;
441 : }
442 : }
443 900080 : return is_valid;
444 : }
445 :
446 0 : Real RBParametrized::get_closest_value(Real value, const std::vector<Real> & list_of_values)
447 : {
448 0 : libmesh_error_msg_if(list_of_values.empty(), "Error: list_of_values is empty.");
449 :
450 0 : Real min_distance = std::numeric_limits<Real>::max();
451 0 : Real closest_val = 0.;
452 0 : for (const auto & current_value : list_of_values)
453 : {
454 0 : Real distance = std::abs(value - current_value);
455 0 : if (distance < min_distance)
456 : {
457 0 : min_distance = distance;
458 0 : closest_val = current_value;
459 : }
460 : }
461 :
462 0 : return closest_val;
463 : }
464 :
465 0 : bool RBParametrized::is_value_in_list(Real value, const std::vector<Real> & list_of_values, Real tol)
466 : {
467 0 : Real closest_value = get_closest_value(value, list_of_values);
468 :
469 : // Check if relative tolerance is satisfied
470 0 : Real rel_error = std::abs(value - closest_value) / std::abs(value);
471 0 : if (rel_error <= tol)
472 : {
473 0 : return true;
474 : }
475 :
476 : // If relative tolerance isn't satisfied, we should still check an absolute
477 : // error, since relative tolerance can be misleading if value is close to zero
478 0 : Real abs_error = std::abs(value - closest_value);
479 0 : return (abs_error <= tol);
480 : }
481 :
482 : } // namespace libMesh
|
