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 831830 : bool RBParametrized::set_parameters(const RBParameters & params)
130 : {
131 831830 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::set_parameters");
132 :
133 : // Terminate if params has the wrong number of parameters or samples.
134 : // If the parameters are outside the min/max range, return false.
135 831830 : const bool valid_params = check_if_valid_params(params);
136 :
137 : // Make a copy of params (default assignment operator just does memberwise copy, which is sufficient here)
138 763580 : this->parameters = params;
139 :
140 831830 : return valid_params;
141 : }
142 :
143 11155637 : const RBParameters & RBParametrized::get_parameters() const
144 : {
145 11155637 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters");
146 :
147 11155637 : return parameters;
148 : }
149 :
150 839 : const RBParameters & RBParametrized::get_parameters_min() const
151 : {
152 839 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters_min");
153 :
154 839 : return parameters_min;
155 : }
156 :
157 839 : const RBParameters & RBParametrized::get_parameters_max() const
158 : {
159 839 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameters_max");
160 :
161 839 : return parameters_max;
162 : }
163 :
164 5292592 : Real RBParametrized::get_parameter_min(const std::string & param_name) const
165 : {
166 5292592 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameter_min");
167 :
168 5292592 : return parameters_min.get_value(param_name);
169 : }
170 :
171 5292592 : Real RBParametrized::get_parameter_max(const std::string & param_name) const
172 : {
173 5292592 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_parameter_max");
174 :
175 5292592 : return parameters_max.get_value(param_name);
176 : }
177 :
178 5240 : void RBParametrized::print_parameters() const
179 : {
180 5240 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::print_current_parameters");
181 :
182 5240 : get_parameters().print();
183 5240 : }
184 :
185 49 : void RBParametrized::write_parameter_data_to_files(const std::string & continuous_param_file_name,
186 : const std::string & discrete_param_file_name,
187 : const bool write_binary_data)
188 : {
189 49 : write_parameter_ranges_to_file(continuous_param_file_name, write_binary_data);
190 49 : write_discrete_parameter_values_to_file(discrete_param_file_name, write_binary_data);
191 49 : }
192 :
193 49 : void RBParametrized::write_parameter_ranges_to_file(const std::string & file_name,
194 : const bool write_binary_data)
195 : {
196 : // The writing mode: ENCODE for binary, WRITE for ASCII
197 49 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
198 :
199 : // Write out the parameter ranges
200 57 : Xdr parameter_ranges_out(file_name, mode);
201 49 : unsigned int n_continuous_params = get_n_continuous_params();
202 8 : parameter_ranges_out << n_continuous_params;
203 :
204 : // Note: the following loops are not candidates for structured
205 : // bindings syntax because the Xdr APIs which they call are not
206 : // defined for const references. We must therefore make copies to
207 : // call these functions.
208 269 : for (const auto & pr : get_parameters_min())
209 : {
210 220 : std::string param_name = pr.first;
211 220 : if (!is_discrete_parameter(param_name))
212 : {
213 220 : Real param_value = get_parameters_min().get_value(param_name);
214 36 : parameter_ranges_out << param_name << param_value;
215 : }
216 : }
217 269 : for (const auto & pr : get_parameters_max())
218 : {
219 220 : std::string param_name = pr.first;
220 220 : if (!is_discrete_parameter(param_name))
221 : {
222 220 : Real param_value = get_parameters_max().get_value(param_name);
223 36 : parameter_ranges_out << param_name << param_value;
224 : }
225 : }
226 :
227 49 : parameter_ranges_out.close();
228 49 : }
229 :
230 49 : void RBParametrized::write_discrete_parameter_values_to_file(const std::string & file_name,
231 : const bool write_binary_data)
232 : {
233 : // write out the discrete parameters, if we have any
234 49 : if (get_n_discrete_params() > 0)
235 : {
236 : // The writing mode: ENCODE for binary, WRITE for ASCII
237 0 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
238 :
239 0 : Xdr discrete_parameters_out(file_name, mode);
240 0 : unsigned int n_discrete_params = get_n_discrete_params();
241 0 : discrete_parameters_out << n_discrete_params;
242 :
243 : // Note: the following loops are not candidates for structured
244 : // bindings syntax because the Xdr APIs which they call are not
245 : // defined for const references. We must therefore make copies
246 : // to call these functions.
247 0 : for (const auto & pr : get_discrete_parameter_values())
248 : {
249 0 : std::string param_name = pr.first;
250 0 : auto n_discrete_values = cast_int<unsigned int>(pr.second.size());
251 0 : discrete_parameters_out << param_name << n_discrete_values;
252 :
253 0 : for (unsigned int i=0; i<n_discrete_values; i++)
254 : {
255 0 : Real discrete_value = pr.second[i];
256 0 : discrete_parameters_out << discrete_value;
257 : }
258 : }
259 0 : }
260 49 : }
261 :
262 285 : void RBParametrized::read_parameter_data_from_files(const std::string & continuous_param_file_name,
263 : const std::string & discrete_param_file_name,
264 : const bool read_binary_data)
265 : {
266 301 : RBParameters param_min;
267 301 : RBParameters param_max;
268 285 : read_parameter_ranges_from_file(continuous_param_file_name,
269 : read_binary_data,
270 : param_min,
271 : param_max);
272 :
273 16 : std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
274 285 : read_discrete_parameter_values_from_file(discrete_param_file_name,
275 : read_binary_data,
276 : discrete_parameter_values_in);
277 :
278 285 : initialize_parameters(param_min, param_max, discrete_parameter_values_in);
279 285 : }
280 :
281 285 : void RBParametrized::read_parameter_ranges_from_file(const std::string & file_name,
282 : const bool read_binary_data,
283 : RBParameters & param_min,
284 : RBParameters & param_max)
285 : {
286 : // The reading mode: DECODE for binary, READ for ASCII
287 285 : XdrMODE mode = read_binary_data ? DECODE : READ;
288 :
289 : // Read in the parameter ranges
290 570 : Xdr parameter_ranges_in(file_name, mode);
291 : unsigned int n_continuous_params;
292 16 : parameter_ranges_in >> n_continuous_params;
293 :
294 1567 : for (unsigned int i=0; i<n_continuous_params; i++)
295 : {
296 72 : std::string param_name;
297 : Real param_value;
298 :
299 72 : parameter_ranges_in >> param_name;
300 72 : parameter_ranges_in >> param_value;
301 :
302 1282 : param_min.set_value(param_name, param_value);
303 : }
304 1567 : for (unsigned int i=0; i<n_continuous_params; i++)
305 : {
306 72 : std::string param_name;
307 : Real param_value;
308 :
309 72 : parameter_ranges_in >> param_name;
310 72 : parameter_ranges_in >> param_value;
311 :
312 1282 : param_max.set_value(param_name, param_value);
313 : }
314 :
315 285 : parameter_ranges_in.close();
316 285 : }
317 :
318 285 : void RBParametrized::read_discrete_parameter_values_from_file(const std::string & file_name,
319 : const bool read_binary_data,
320 : std::map<std::string, std::vector<Real>> & discrete_parameter_values)
321 : {
322 : // read in the discrete parameters, if we have any
323 301 : std::ifstream check_if_file_exists(file_name.c_str());
324 285 : if (check_if_file_exists.good())
325 : {
326 : // The reading mode: DECODE for binary, READ for ASCII
327 0 : XdrMODE mode = read_binary_data ? DECODE : READ;
328 :
329 : // Read in the parameter ranges
330 0 : Xdr discrete_parameter_values_in(file_name, mode);
331 : unsigned int n_discrete_params;
332 0 : discrete_parameter_values_in >> n_discrete_params;
333 :
334 0 : for (unsigned int i=0; i<n_discrete_params; i++)
335 : {
336 0 : std::string param_name;
337 0 : discrete_parameter_values_in >> param_name;
338 :
339 : unsigned int n_discrete_values;
340 0 : discrete_parameter_values_in >> n_discrete_values;
341 :
342 0 : std::vector<Real> discrete_values(n_discrete_values);
343 0 : for (auto & val : discrete_values)
344 0 : discrete_parameter_values_in >> val;
345 :
346 0 : discrete_parameter_values[param_name] = discrete_values;
347 : }
348 0 : }
349 285 : }
350 :
351 1719 : bool RBParametrized::is_discrete_parameter(const std::string & mu_name) const
352 : {
353 1719 : libmesh_error_msg_if(!parameters_initialized,
354 : "Error: parameters not initialized in RBParametrized::is_discrete_parameter");
355 :
356 1719 : return _discrete_parameter_values.count(mu_name);
357 : }
358 :
359 10584429 : const std::map<std::string, std::vector<Real>> & RBParametrized::get_discrete_parameter_values() const
360 : {
361 10584429 : libmesh_error_msg_if(!parameters_initialized, "Error: parameters not initialized in RBParametrized::get_discrete_parameter_values");
362 :
363 10584429 : return _discrete_parameter_values;
364 : }
365 :
366 285 : void RBParametrized::print_discrete_parameter_values() const
367 : {
368 293 : for (const auto & [name, values] : get_discrete_parameter_values())
369 : {
370 0 : libMesh::out << "Discrete parameter " << name << ", values: ";
371 :
372 0 : for (const auto & value : values)
373 0 : libMesh::out << value << " ";
374 0 : libMesh::out << std::endl;
375 : }
376 285 : }
377 :
378 831830 : bool RBParametrized::check_if_valid_params(const RBParameters ¶ms) const
379 : {
380 : // Check if number of parameters are correct.
381 831830 : libmesh_error_msg_if(params.n_parameters() != get_n_params(),
382 : "Error: Number of parameters don't match; found "
383 : << params.n_parameters() << ", expected "
384 : << get_n_params());
385 :
386 68250 : bool is_valid = true;
387 831830 : std::string prev_param_name = "";
388 6123140 : for (const auto & [param_name, sample_vec] : params)
389 : {
390 437022 : std::size_t sample_idx = 0;
391 5291310 : const Real & min_value = get_parameter_min(param_name);
392 5291310 : const Real & max_value = get_parameter_max(param_name);
393 10583046 : for (const auto & value_vec : sample_vec)
394 : {
395 10583472 : for (const auto & value : value_vec)
396 : {
397 : // Check every parameter value (including across samples and vector-values)
398 : // to ensure it's within the min/max range.
399 5291736 : const bool outside_range = ((value < min_value) || (value > max_value));
400 5291736 : is_valid = is_valid && !outside_range;
401 5291736 : if (outside_range && verbose_mode)
402 : {
403 0 : libMesh::out << "Warning: parameter " << param_name << " value="
404 0 : << value << " outside acceptable range: ("
405 0 : << min_value << ", " << max_value << ")";
406 : }
407 :
408 : // For discrete params, make sure params.get_value(param_name) is sufficiently
409 : // close to one of the discrete parameter values.
410 : // Note that vector-values not yet supported in discrete parameters,
411 : // and the .get_sample_value() call will throw an error if the user
412 : // tries to do it.
413 10146438 : if (const auto it = get_discrete_parameter_values().find(param_name);
414 5291736 : it != get_discrete_parameter_values().end())
415 : {
416 : const bool is_value_discrete =
417 0 : is_value_in_list(params.get_sample_value(param_name, sample_idx),
418 0 : it->second,
419 : TOLERANCE);
420 0 : is_valid = is_valid && is_value_discrete;
421 0 : if (!is_value_discrete && verbose_mode)
422 0 : libMesh::out << "Warning: parameter " << param_name << " value="
423 0 : << value << " is not in discrete value list.";
424 : }
425 : }
426 5291736 : ++sample_idx;
427 : }
428 : }
429 900080 : return is_valid;
430 : }
431 :
432 0 : Real RBParametrized::get_closest_value(Real value, const std::vector<Real> & list_of_values)
433 : {
434 0 : libmesh_error_msg_if(list_of_values.empty(), "Error: list_of_values is empty.");
435 :
436 0 : Real min_distance = std::numeric_limits<Real>::max();
437 0 : Real closest_val = 0.;
438 0 : for (const auto & current_value : list_of_values)
439 : {
440 0 : Real distance = std::abs(value - current_value);
441 0 : if (distance < min_distance)
442 : {
443 0 : min_distance = distance;
444 0 : closest_val = current_value;
445 : }
446 : }
447 :
448 0 : return closest_val;
449 : }
450 :
451 0 : bool RBParametrized::is_value_in_list(Real value, const std::vector<Real> & list_of_values, Real tol)
452 : {
453 0 : Real closest_value = get_closest_value(value, list_of_values);
454 :
455 : // Check if relative tolerance is satisfied
456 0 : Real rel_error = std::abs(value - closest_value) / std::abs(value);
457 0 : if (rel_error <= tol)
458 : {
459 0 : return true;
460 : }
461 :
462 : // If relative tolerance isn't satisfied, we should still check an absolute
463 : // error, since relative tolerance can be misleading if value is close to zero
464 0 : Real abs_error = std::abs(value - closest_value);
465 0 : return (abs_error <= tol);
466 : }
467 :
468 : } // namespace libMesh
|
