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 : // C++ includes
21 : #include <fstream>
22 : #include <sstream>
23 :
24 : // LibMesh includes
25 : #include "libmesh/sparse_matrix.h"
26 : #include "libmesh/numeric_vector.h"
27 : #include "libmesh/dense_matrix.h"
28 : #include "libmesh/dense_vector.h"
29 : #include "libmesh/dof_map.h"
30 : #include "libmesh/libmesh_logging.h"
31 : #include "libmesh/equation_systems.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/parallel_algebra.h"
34 : #include "libmesh/fe.h"
35 : #include "libmesh/quadrature.h"
36 : #include "libmesh/utility.h"
37 : #include "libmesh/fe_interface.h"
38 : #include "libmesh/fe_compute_data.h"
39 : #include "libmesh/getpot.h"
40 : #include "libmesh/exodusII_io.h"
41 : #include "libmesh/fem_context.h"
42 : #include "libmesh/elem.h"
43 : #include "libmesh/int_range.h"
44 :
45 : // rbOOmit includes
46 : #include "libmesh/rb_construction_base.h"
47 : #include "libmesh/rb_eim_construction.h"
48 : #include "libmesh/rb_eim_evaluation.h"
49 : #include "libmesh/rb_parameters.h"
50 : #include "libmesh/rb_parametrized_function.h"
51 :
52 : // C++ include
53 : #include <limits>
54 : #include <memory>
55 : #include <random>
56 :
57 : namespace libMesh
58 : {
59 :
60 : namespace
61 : {
62 : // We make an anonymous namespace for local helper functions. The helper functions below
63 : // are used to do some basic operators for QpDataMap and SideQpDataMap.
64 :
65 : // Recall that we have:
66 : // typedef std::map<dof_id_type, std::vector<std::vector<Number>>> QpDataMap;
67 : // typedef std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Number>>> SideQpDataMap;
68 :
69 : // Implement u <- u + k*v
70 : template <typename DataMap>
71 0 : void add(DataMap & u, const Number k, const DataMap & v)
72 : {
73 0 : for (auto & [key, vec_vec_u] : u)
74 : {
75 0 : const std::vector<std::vector<Number>> & vec_vec_v = libmesh_map_find(v, key);
76 :
77 0 : libmesh_error_msg_if (vec_vec_u.size() != vec_vec_v.size(), "Size mismatch");
78 :
79 0 : for (auto i : index_range(vec_vec_u))
80 : {
81 0 : libmesh_error_msg_if (vec_vec_u[i].size() != vec_vec_v[i].size(), "Size mismatch");
82 0 : for (auto j : index_range(vec_vec_u[i]))
83 : {
84 0 : vec_vec_u[i][j] += k*vec_vec_v[i][j];
85 : }
86 : }
87 : }
88 0 : }
89 :
90 0 : void add_node_data_map(RBEIMConstruction::NodeDataMap & u, const Number k, const RBEIMConstruction::NodeDataMap & v)
91 : {
92 0 : for (auto & [key, vec_u] : u)
93 : {
94 0 : const std::vector<Number> & vec_v = libmesh_map_find(v, key);
95 :
96 0 : libmesh_error_msg_if (vec_u.size() != vec_v.size(), "Size mismatch");
97 :
98 0 : for (auto i : index_range(vec_u))
99 : {
100 0 : vec_u[i] += k*vec_v[i];
101 : }
102 : }
103 0 : }
104 :
105 : }
106 :
107 0 : RBEIMConstruction::RBEIMConstruction (EquationSystems & es,
108 : const std::string & name_in,
109 0 : const unsigned int number_in)
110 : : RBConstructionBase(es, name_in, number_in),
111 0 : best_fit_type_flag(PROJECTION_BEST_FIT),
112 0 : _Nmax(0),
113 0 : _set_Nmax_from_n_snapshots(false),
114 0 : _Nmax_from_n_snapshots_increment(0),
115 0 : _rel_training_tolerance(1.e-4),
116 0 : _abs_training_tolerance(1.e-12),
117 0 : _max_abs_value_in_training_set(0.),
118 0 : _max_abs_value_in_training_set_index(0)
119 : {
120 : // The training set should be the same on all processors in the
121 : // case of EIM training.
122 0 : serial_training_set = true;
123 0 : }
124 :
125 0 : RBEIMConstruction::~RBEIMConstruction () = default;
126 :
127 0 : void RBEIMConstruction::clear()
128 : {
129 0 : RBConstructionBase::clear();
130 :
131 0 : _rb_eim_assembly_objects.clear();
132 :
133 0 : _local_parametrized_functions_for_training.clear();
134 0 : _local_quad_point_locations.clear();
135 0 : _local_quad_point_JxW.clear();
136 0 : _local_quad_point_subdomain_ids.clear();
137 :
138 0 : _local_side_parametrized_functions_for_training.clear();
139 0 : _local_side_quad_point_locations.clear();
140 0 : _local_side_quad_point_JxW.clear();
141 0 : _local_side_quad_point_subdomain_ids.clear();
142 0 : _local_side_quad_point_boundary_ids.clear();
143 0 : _local_side_quad_point_side_types.clear();
144 :
145 0 : _local_node_parametrized_functions_for_training.clear();
146 0 : _local_node_locations.clear();
147 0 : _local_node_boundary_ids.clear();
148 :
149 0 : _eim_projection_matrix.resize(0,0);
150 0 : }
151 :
152 0 : void RBEIMConstruction::set_rb_eim_evaluation(RBEIMEvaluation & rb_eim_eval_in)
153 : {
154 0 : _rb_eim_eval = &rb_eim_eval_in;
155 0 : }
156 :
157 0 : RBEIMEvaluation & RBEIMConstruction::get_rb_eim_evaluation()
158 : {
159 0 : libmesh_error_msg_if(!_rb_eim_eval, "Error: RBEIMEvaluation object hasn't been initialized yet");
160 0 : return *_rb_eim_eval;
161 : }
162 :
163 0 : const RBEIMEvaluation & RBEIMConstruction::get_rb_eim_evaluation() const
164 : {
165 0 : libmesh_error_msg_if(!_rb_eim_eval, "Error: RBEIMEvaluation object hasn't been initialized yet");
166 0 : return *_rb_eim_eval;
167 : }
168 :
169 0 : void RBEIMConstruction::set_best_fit_type_flag (const std::string & best_fit_type_string)
170 : {
171 0 : if (best_fit_type_string == "projection")
172 : {
173 0 : best_fit_type_flag = PROJECTION_BEST_FIT;
174 : }
175 0 : else if (best_fit_type_string == "eim")
176 : {
177 0 : best_fit_type_flag = EIM_BEST_FIT;
178 : }
179 0 : else if (best_fit_type_string == "pod")
180 : {
181 0 : best_fit_type_flag = POD_BEST_FIT;
182 : }
183 : else
184 0 : libmesh_error_msg("Error: invalid best_fit_type in input file");
185 0 : }
186 :
187 0 : void RBEIMConstruction::print_info()
188 : {
189 : // Print out info that describes the current setup
190 0 : libMesh::out << std::endl << "RBEIMConstruction parameters:" << std::endl;
191 0 : libMesh::out << "system name: " << this->name() << std::endl;
192 0 : libMesh::out << "Nmax: " << get_Nmax() << std::endl;
193 0 : if (_set_Nmax_from_n_snapshots)
194 : {
195 0 : libMesh::out << "Overruling Nmax based on number of snapshots, with increment set to "
196 0 : << _Nmax_from_n_snapshots_increment
197 0 : << std::endl;
198 : }
199 0 : libMesh::out << "Greedy relative error tolerance: " << get_rel_training_tolerance() << std::endl;
200 0 : libMesh::out << "Greedy absolute error tolerance: " << get_abs_training_tolerance() << std::endl;
201 0 : libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
202 0 : for (const auto & pr : get_parameters())
203 0 : if (!is_discrete_parameter(pr.first))
204 : {
205 0 : libMesh::out << "Parameter " << pr.first
206 0 : << ": Min = " << get_parameter_min(pr.first)
207 0 : << ", Max = " << get_parameter_max(pr.first) << std::endl;
208 : }
209 :
210 0 : print_discrete_parameter_values();
211 0 : libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
212 0 : libMesh::out << "quiet mode? " << is_quiet() << std::endl;
213 :
214 0 : if (best_fit_type_flag == PROJECTION_BEST_FIT)
215 : {
216 0 : libMesh::out << "EIM best fit type: projection" << std::endl;
217 : }
218 : else
219 0 : if (best_fit_type_flag == EIM_BEST_FIT)
220 : {
221 0 : libMesh::out << "EIM best fit type: eim" << std::endl;
222 : }
223 0 : libMesh::out << std::endl;
224 0 : }
225 :
226 0 : void RBEIMConstruction::initialize_eim_construction()
227 : {
228 0 : initialize_parametrized_functions_in_training_set();
229 0 : }
230 :
231 0 : void RBEIMConstruction::process_parameters_file (const std::string & parameters_filename)
232 : {
233 : // First read in data from input_filename
234 0 : GetPot infile(parameters_filename);
235 :
236 0 : std::string best_fit_type_string = infile("best_fit_type","projection");
237 0 : set_best_fit_type_flag(best_fit_type_string);
238 :
239 0 : const unsigned int n_training_samples = infile("n_training_samples",0);
240 0 : const bool deterministic_training = infile("deterministic_training",false);
241 0 : unsigned int training_parameters_random_seed_in =
242 : static_cast<unsigned int>(-1);
243 0 : training_parameters_random_seed_in = infile("training_parameters_random_seed",
244 : training_parameters_random_seed_in);
245 0 : const bool quiet_mode_in = infile("quiet_mode", quiet_mode);
246 0 : const unsigned int Nmax_in = infile("Nmax", _Nmax);
247 0 : const Real rel_training_tolerance_in = infile("rel_training_tolerance",
248 0 : _rel_training_tolerance);
249 0 : const Real abs_training_tolerance_in = infile("abs_training_tolerance",
250 0 : _abs_training_tolerance);
251 :
252 : // Read in the parameters from the input file too
253 0 : unsigned int n_continuous_parameters = infile.vector_variable_size("parameter_names");
254 0 : RBParameters mu_min_in;
255 0 : RBParameters mu_max_in;
256 0 : for (unsigned int i=0; i<n_continuous_parameters; i++)
257 : {
258 : // Read in the parameter names
259 0 : std::string param_name = infile("parameter_names", "NONE", i);
260 :
261 : {
262 0 : Real min_val = infile(param_name, 0., 0);
263 0 : mu_min_in.set_value(param_name, min_val);
264 : }
265 :
266 : {
267 0 : Real max_val = infile(param_name, 0., 1);
268 0 : mu_max_in.set_value(param_name, max_val);
269 : }
270 : }
271 :
272 0 : std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
273 :
274 0 : unsigned int n_discrete_parameters = infile.vector_variable_size("discrete_parameter_names");
275 0 : for (unsigned int i=0; i<n_discrete_parameters; i++)
276 : {
277 0 : std::string param_name = infile("discrete_parameter_names", "NONE", i);
278 :
279 0 : unsigned int n_vals_for_param = infile.vector_variable_size(param_name);
280 0 : std::vector<Real> vals_for_param(n_vals_for_param);
281 0 : for (auto j : make_range(vals_for_param.size()))
282 0 : vals_for_param[j] = infile(param_name, 0., j);
283 :
284 0 : discrete_parameter_values_in[param_name] = vals_for_param;
285 : }
286 :
287 0 : std::map<std::string,bool> log_scaling_in;
288 : // For now, just set all entries to false.
289 : // TODO: Implement a decent way to specify log-scaling true/false
290 : // in the input text file
291 0 : for (const auto & pr : mu_min_in)
292 0 : log_scaling_in[pr.first] = false;
293 :
294 : // Set the parameters that have been read in
295 0 : set_rb_construction_parameters(n_training_samples,
296 : deterministic_training,
297 : training_parameters_random_seed_in,
298 : quiet_mode_in,
299 : Nmax_in,
300 : rel_training_tolerance_in,
301 : abs_training_tolerance_in,
302 : mu_min_in,
303 : mu_max_in,
304 : discrete_parameter_values_in,
305 : log_scaling_in);
306 0 : }
307 :
308 0 : void RBEIMConstruction::set_rb_construction_parameters(unsigned int n_training_samples_in,
309 : bool deterministic_training_in,
310 : int training_parameters_random_seed_in,
311 : bool quiet_mode_in,
312 : unsigned int Nmax_in,
313 : Real rel_training_tolerance_in,
314 : Real abs_training_tolerance_in,
315 : const RBParameters & mu_min_in,
316 : const RBParameters & mu_max_in,
317 : const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
318 : const std::map<std::string,bool> & log_scaling_in,
319 : std::map<std::string, std::vector<RBParameter>> * training_sample_list)
320 : {
321 : // Read in training_parameters_random_seed value. This is used to
322 : // seed the RNG when picking the training parameters. By default the
323 : // value is -1, which means use std::time to seed the RNG.
324 0 : set_training_random_seed(training_parameters_random_seed_in);
325 :
326 : // Set quiet mode
327 0 : set_quiet_mode(quiet_mode_in);
328 :
329 : // Initialize RB parameters
330 0 : set_Nmax(Nmax_in);
331 :
332 0 : set_rel_training_tolerance(rel_training_tolerance_in);
333 0 : set_abs_training_tolerance(abs_training_tolerance_in);
334 :
335 0 : if (get_rb_eim_evaluation().get_parametrized_function().is_lookup_table)
336 : {
337 : const std::string & lookup_table_param_name =
338 0 : get_rb_eim_evaluation().get_parametrized_function().lookup_table_param_name;
339 :
340 0 : libmesh_error_msg_if(!discrete_parameter_values_in.count(lookup_table_param_name),
341 : "Lookup table parameter should be discrete");
342 :
343 : // Make an editable copy of discrete_parameters_values_in.
344 : std::map<std::string, std::vector<Real>> discrete_parameter_values_final(
345 0 : discrete_parameter_values_in);
346 :
347 : std::vector<Real> & lookup_table_param_values =
348 0 : libmesh_map_find(discrete_parameter_values_final, lookup_table_param_name);
349 :
350 : // Overwrite the discrete values for lookup_table_param to make sure that
351 : // it is: 0, 1, 2, ..., size-1.
352 0 : std::iota(lookup_table_param_values.begin(), lookup_table_param_values.end(), 0);
353 :
354 : // Also, overwrite n_training_samples_in to make sure it matches
355 : // lookup_table_size so that we will get full coverage of the
356 : // lookup table in our training set.
357 0 : n_training_samples_in = lookup_table_param_values.size();
358 :
359 : // Initialize the parameter ranges and the parameters themselves
360 0 : initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_final);
361 : }
362 : else
363 : {
364 : // Initialize the parameter ranges and the parameters themselves
365 0 : initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_in);
366 : }
367 :
368 0 : bool updated_deterministic_training = deterministic_training_in;
369 0 : if (training_sample_list && (this->get_parameters_min().n_parameters() > 3))
370 : {
371 : // In this case we force deterministic_training to be false because
372 : // a) deterministic training samples are not currrently supported with
373 : // more than 3 parameters, and
374 : // b) we will overwrite the training samples anyway in the call to
375 : // load_training_set() below, so we do not want to generate an
376 : // error due to deterministic training sample generation when
377 : // the samples will be overwritten anyway.
378 0 : updated_deterministic_training = false;
379 : }
380 :
381 0 : initialize_training_parameters(this->get_parameters_min(),
382 : this->get_parameters_max(),
383 : n_training_samples_in,
384 : log_scaling_in,
385 0 : updated_deterministic_training);
386 :
387 0 : if (training_sample_list)
388 : {
389 : // Note that we must call initialize_training_parameters() before
390 : // load_training_set() in order to initialize the parameter vectors.
391 0 : load_training_set(*training_sample_list);
392 : }
393 :
394 :
395 0 : if (get_rb_eim_evaluation().get_parametrized_function().is_lookup_table)
396 : {
397 : // Also, now that we've initialized the training set, overwrite the training
398 : // samples to ensure that we have full coverage of the lookup tbale.
399 : const std::string & lookup_table_param_name =
400 0 : get_rb_eim_evaluation().get_parametrized_function().lookup_table_param_name;
401 :
402 : // Fill the lookup_table_training_samples with sequential single-entry vectors,
403 : // i.e. {{0.0}, {1.0}, {2.0}, ...}
404 0 : Real val = 0.0;
405 0 : std::vector<RBParameter> lookup_table_training_samples(n_training_samples_in, {val});
406 0 : for (auto & vec : lookup_table_training_samples)
407 : {
408 0 : vec[0] = val;
409 0 : val += 1.0; // Could use val++, but better to be explicit for doubles.
410 : }
411 :
412 0 : set_training_parameter_values(lookup_table_param_name, lookup_table_training_samples);
413 0 : }
414 0 : }
415 :
416 0 : Real RBEIMConstruction::train_eim_approximation()
417 : {
418 0 : if (_normalize_solution_snapshots)
419 0 : apply_normalization_to_solution_snapshots();
420 :
421 0 : if(best_fit_type_flag == POD_BEST_FIT)
422 : {
423 0 : train_eim_approximation_with_POD();
424 0 : return 0.;
425 : }
426 : else
427 : {
428 0 : return train_eim_approximation_with_greedy();
429 : }
430 : }
431 :
432 0 : Real RBEIMConstruction::train_eim_approximation_with_greedy()
433 : {
434 0 : LOG_SCOPE("train_eim_approximation_with_greedy()", "RBEIMConstruction");
435 :
436 0 : RBEIMEvaluation & rbe = get_rb_eim_evaluation();
437 :
438 : // We need space for one extra interpolation point if we're using the
439 : // EIM error indicator.
440 0 : unsigned int max_matrix_size = rbe.use_eim_error_indicator() ? get_Nmax()+1 : get_Nmax();
441 0 : _eim_projection_matrix.resize(max_matrix_size,max_matrix_size);
442 :
443 0 : rbe.initialize_parameters(*this);
444 0 : rbe.resize_data_structures(max_matrix_size);
445 :
446 : // If we are continuing from a previous training run,
447 : // we might already be at the max number of basis functions.
448 : // If so, we can just return.
449 0 : libmesh_error_msg_if(rbe.get_n_basis_functions() > 0,
450 : "Error: We currently only support EIM training starting from an empty basis");
451 :
452 0 : libMesh::out << std::endl << "---- Performing Greedy EIM basis enrichment ----" << std::endl;
453 :
454 0 : if (_max_abs_value_in_training_set <= get_abs_training_tolerance())
455 : {
456 0 : std::cout << "Maximum absolute value in the training set is "
457 0 : << _max_abs_value_in_training_set << ", which is less than or equal to the abs. training tolerance of "
458 0 : << get_abs_training_tolerance() << ", hence exiting Greedy basis enrichment with empty basis" << std::endl;
459 0 : return 0.0;
460 : }
461 :
462 : // Initialize greedy_error so that we do not incorrectly set is_zero_bf=true on
463 : // the first iteration.
464 0 : Real greedy_error = -1.;
465 0 : std::vector<RBParameters> greedy_param_list;
466 :
467 : // Initialize the current training index to the index that corresponds
468 : // to the largest (in terms of infinity norm) function in the training set.
469 : // We do this to ensure that the first EIM basis function is not zero.
470 0 : unsigned int current_training_index = _max_abs_value_in_training_set_index;
471 0 : set_params_from_training_set(current_training_index);
472 :
473 : // We use this boolean to indicate if we will run one more iteration
474 : // before exiting the loop below. We use this when computing the EIM
475 : // error indicator, which requires one extra EIM iteration.
476 0 : bool exit_on_next_iteration = false;
477 :
478 : // We also initialize a boolean to keep track of whether we have
479 : // reached "n_samples" EIM basis functions, since we need to
480 : // handle the EIM error indicator in a special way in this case.
481 0 : bool bfs_equals_n_samples = false;
482 :
483 : while (true)
484 : {
485 0 : if (rbe.get_n_basis_functions() >= get_n_training_samples())
486 : {
487 0 : libMesh::out << "Number of basis functions (" << rbe.get_n_basis_functions()
488 0 : << ") equals number of training samples." << std::endl;
489 :
490 0 : bfs_equals_n_samples = true;
491 :
492 : // If exit_on_next_iteration==true then we don't exit yet, since
493 : // we still need to add data for the error indicator before exiting.
494 0 : if (!exit_on_next_iteration)
495 0 : break;
496 : }
497 :
498 0 : libMesh::out << "Greedily selected parameter vector:" << std::endl;
499 0 : print_parameters();
500 0 : greedy_param_list.emplace_back(get_parameters());
501 :
502 0 : libMesh::out << "Enriching the EIM approximation" << std::endl;
503 : libmesh_try
504 : {
505 0 : bool is_zero_bf = bfs_equals_n_samples || (greedy_error == 0.);
506 :
507 : // If is_zero_bf==true then we add an "extra point" because we
508 : // cannot add a usual EIM interpolation point in that case since
509 : // the full EIM space is already covered. This is necessary when we
510 : // want to add an extra point for error indicator purposes in the
511 : // is_zero_bf==true case, for example.
512 0 : std::unique_ptr<EimPointData> eim_point_data;
513 0 : if (is_zero_bf)
514 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point_from_training_sample());
515 :
516 : // If exit_on_next_iteration==true then we do not add a basis function in
517 : // that case since in that case we only need to add data for the EIM error
518 : // indicator.
519 0 : enrich_eim_approximation(current_training_index,
520 0 : /*add_basis_function*/ !exit_on_next_iteration,
521 : eim_point_data.get());
522 0 : update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
523 :
524 0 : libMesh::out << std::endl << "---- Basis dimension: "
525 0 : << rbe.get_n_basis_functions() << " ----" << std::endl;
526 :
527 0 : if (get_rb_eim_evaluation().get_parametrized_function().is_lookup_table &&
528 0 : best_fit_type_flag == EIM_BEST_FIT)
529 : {
530 : // If this is a lookup table and we're using "EIM best fit" then we
531 : // need to update the eim_solutions after each EIM enrichment so that
532 : // we can call rb_eim_eval.rb_eim_solve() from within compute_max_eim_error().
533 0 : store_eim_solutions_for_training_set();
534 : }
535 :
536 0 : libMesh::out << "Computing EIM error on training set" << std::endl;
537 0 : std::tie(greedy_error, current_training_index) = compute_max_eim_error();
538 0 : set_params_from_training_set(current_training_index);
539 :
540 0 : libMesh::out << "Maximum EIM error is " << greedy_error << std::endl << std::endl;
541 0 : }
542 0 : libmesh_catch (const std::exception & e)
543 : {
544 : // If we hit an exception when performing the enrichment for the error indicator, then
545 : // we just continue and skip the error indicator. Otherwise we rethrow the exception.
546 0 : if (exit_on_next_iteration)
547 : {
548 0 : std::cout << "Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
549 0 : break;
550 : }
551 : else
552 0 : throw;
553 0 : }
554 :
555 0 : if (exit_on_next_iteration)
556 : {
557 0 : libMesh::out << "Extra EIM iteration for error indicator is complete, hence exiting EIM training now" << std::endl;
558 0 : break;
559 : }
560 :
561 : // Convergence and/or termination tests
562 : {
563 0 : bool exit_condition_satisfied = false;
564 :
565 0 : if (rbe.get_n_basis_functions() >= this->get_Nmax())
566 : {
567 0 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
568 0 : << get_Nmax() << std::endl;
569 0 : exit_condition_satisfied = true;
570 : }
571 :
572 : // We consider the relative tolerance as relative to the maximum value in the training
573 : // set, since we assume that this maximum value provides a relevant scaling.
574 0 : if (!exit_condition_satisfied)
575 0 : if (greedy_error < (get_rel_training_tolerance() * get_max_abs_value_in_training_set()))
576 : {
577 0 : libMesh::out << "Relative error tolerance reached." << std::endl;
578 0 : exit_condition_satisfied = true;
579 : }
580 :
581 0 : if (!exit_condition_satisfied)
582 0 : if (greedy_error < get_abs_training_tolerance())
583 : {
584 0 : libMesh::out << "Absolute error tolerance reached." << std::endl;
585 0 : exit_condition_satisfied = true;
586 : }
587 :
588 0 : bool has_parameters = (get_parameters().n_parameters() > 0);
589 0 : if (!exit_condition_satisfied)
590 : {
591 0 : bool do_exit = false;
592 : // In the check for repeated parameters we have to make sure this isn't a case
593 : // with no parameters, since in that case we would always report repeated
594 : // parameters.
595 0 : for (auto & param : greedy_param_list)
596 0 : if (param == get_parameters() && has_parameters)
597 : {
598 0 : libMesh::out << "Exiting greedy because the same parameters were selected twice"
599 0 : << std::endl;
600 0 : do_exit = true;
601 0 : break;
602 : }
603 :
604 0 : if (do_exit)
605 0 : exit_condition_satisfied = true;
606 : }
607 :
608 0 : if (exit_condition_satisfied)
609 : {
610 : // If we're using the EIM error indicator then we need to run
611 : // one extra EIM iteration since we use the extra EIM point
612 : // to obtain our error indicator. If we're not using the EIM
613 : // error indicator, then we just exit now.
614 0 : if (get_rb_eim_evaluation().use_eim_error_indicator() && has_parameters)
615 : {
616 0 : exit_on_next_iteration = true;
617 0 : libMesh::out << "EIM error indicator is active, hence we will run one extra EIM iteration before exiting"
618 0 : << std::endl;
619 : }
620 : else
621 0 : break;
622 : }
623 : }
624 0 : } // end while(true)
625 :
626 0 : if (rbe.get_parametrized_function().is_lookup_table &&
627 0 : best_fit_type_flag != EIM_BEST_FIT)
628 : {
629 : // We only enter here if best_fit_type_flag != EIM_BEST_FIT because we
630 : // already called this above in the EIM_BEST_FIT case.
631 0 : store_eim_solutions_for_training_set();
632 : }
633 :
634 0 : return greedy_error;
635 0 : }
636 :
637 0 : void RBEIMConstruction::apply_normalization_to_solution_snapshots()
638 : {
639 0 : LOG_SCOPE("apply_normalization_to_solution_snapshots()", "RBEIMConstruction");
640 :
641 0 : libMesh::out << "Normalizing solution snapshots" << std::endl;
642 :
643 0 : bool apply_comp_scaling = !get_rb_eim_evaluation().scale_components_in_enrichment().empty();
644 0 : unsigned int n_snapshots = get_n_training_samples();
645 0 : RBEIMEvaluation & rbe = get_rb_eim_evaluation();
646 :
647 0 : for (unsigned int i=0; i<n_snapshots; i++)
648 : {
649 0 : Real norm_val = 0.;
650 0 : if (rbe.get_parametrized_function().on_mesh_sides())
651 : {
652 0 : norm_val = std::sqrt(std::real(side_inner_product(
653 0 : _local_side_parametrized_functions_for_training[i],
654 0 : _local_side_parametrized_functions_for_training[i],
655 : apply_comp_scaling)));
656 :
657 0 : if (norm_val > 0.)
658 0 : scale_parametrized_function(_local_side_parametrized_functions_for_training[i], 1./norm_val);
659 : }
660 0 : else if (rbe.get_parametrized_function().on_mesh_nodes())
661 : {
662 0 : norm_val = std::sqrt(std::real(node_inner_product(
663 0 : _local_node_parametrized_functions_for_training[i],
664 0 : _local_node_parametrized_functions_for_training[i],
665 : apply_comp_scaling)));
666 :
667 0 : if (norm_val > 0.)
668 0 : scale_node_parametrized_function(_local_node_parametrized_functions_for_training[i], 1./norm_val);
669 : }
670 : else
671 : {
672 0 : norm_val = std::sqrt(std::real(inner_product(
673 0 : _local_parametrized_functions_for_training[i],
674 0 : _local_parametrized_functions_for_training[i],
675 : apply_comp_scaling)));
676 :
677 0 : if (norm_val > 0.)
678 0 : scale_parametrized_function(_local_parametrized_functions_for_training[i], 1./norm_val);
679 : }
680 :
681 : // Since we're rescaling the training samples, we should also rescale
682 : // _max_abs_value_in_training_set in the same way. We obtained the
683 : // _max_abs_value_in_training_set value from the sample with index
684 : // _max_abs_value_in_training_set_index, so we rescale using the norm
685 : // of this sample here. Note that the value we obtain may not be exactly
686 : // equal to the max value we would obtain after looping over all the
687 : // normalized samples, because we normalized in the L2 norm rather than
688 : // in the max norm, but it should be close to this "true max" value. The
689 : // main use-case of _max_abs_value_in_training_set is to scale the relative
690 : // error tolerance in the Greedy training, and the L2-norm-scaled value of
691 : // _max_abs_value_in_training_set that we obtain here should be sufficient
692 : // for that purpose.
693 0 : if ((i == _max_abs_value_in_training_set_index) && (norm_val > 0.))
694 0 : _max_abs_value_in_training_set /= norm_val;
695 : }
696 :
697 0 : libMesh::out << "Maximum absolute value in the training set after normalization: "
698 0 : << _max_abs_value_in_training_set << std::endl << std::endl;
699 0 : }
700 :
701 0 : Real RBEIMConstruction::train_eim_approximation_with_POD()
702 : {
703 0 : LOG_SCOPE("train_eim_approximation_with_POD()", "RBEIMConstruction");
704 :
705 0 : RBEIMEvaluation & rbe = get_rb_eim_evaluation();
706 :
707 0 : unsigned int n_snapshots = get_n_training_samples();
708 :
709 : // If _set_Nmax_from_n_snapshots=true, then we overrule Nmax.
710 0 : if (_set_Nmax_from_n_snapshots)
711 : {
712 0 : int updated_Nmax = (static_cast<int>(n_snapshots) + _Nmax_from_n_snapshots_increment);
713 :
714 : // We only overrule _Nmax if updated_Nmax is positive, since if Nmax=0 then we'll skip
715 : // training here entirely, which is typically not what we want.
716 0 : if (updated_Nmax > 0)
717 0 : _Nmax = static_cast<unsigned int>(updated_Nmax);
718 : }
719 :
720 : // _eim_projection_matrix is not used in the POD case, but we resize it here in any case
721 : // to be consistent with what we do in train_eim_approximation_with_greedy().
722 : // We need space for one extra interpolation point if we're using the
723 : // EIM error indicator.
724 0 : unsigned int max_matrix_size = rbe.use_eim_error_indicator() ? get_Nmax()+1 : get_Nmax();
725 0 : _eim_projection_matrix.resize(max_matrix_size,max_matrix_size);
726 :
727 0 : rbe.initialize_parameters(*this);
728 0 : rbe.resize_data_structures(get_Nmax());
729 :
730 0 : libmesh_error_msg_if(rbe.get_n_basis_functions() > 0,
731 : "Error: We currently only support EIM training starting from an empty basis");
732 :
733 0 : libMesh::out << std::endl << "---- Performing POD EIM basis enrichment ----" << std::endl;
734 :
735 : // Set up the POD "correlation matrix". This enables us to compute the POD via the
736 : // "method of snapshots", in which we compute a low rank representation of the
737 : // n_snapshots x n_snapshots matrix.
738 0 : DenseMatrix<Number> correlation_matrix(n_snapshots,n_snapshots);
739 :
740 0 : std::cout << "Start computing correlation matrix" << std::endl;
741 :
742 0 : bool apply_comp_scaling = !get_rb_eim_evaluation().scale_components_in_enrichment().empty();
743 0 : for (unsigned int i=0; i<n_snapshots; i++)
744 : {
745 0 : for (unsigned int j=0; j<=i; j++)
746 : {
747 0 : Number inner_prod = 0.;
748 0 : if (rbe.get_parametrized_function().on_mesh_sides())
749 : {
750 0 : inner_prod = side_inner_product(
751 0 : _local_side_parametrized_functions_for_training[i],
752 0 : _local_side_parametrized_functions_for_training[j],
753 : apply_comp_scaling);
754 : }
755 0 : else if (rbe.get_parametrized_function().on_mesh_nodes())
756 : {
757 0 : inner_prod = node_inner_product(
758 0 : _local_node_parametrized_functions_for_training[i],
759 0 : _local_node_parametrized_functions_for_training[j],
760 : apply_comp_scaling);
761 : }
762 : else
763 : {
764 0 : inner_prod = inner_product(
765 0 : _local_parametrized_functions_for_training[i],
766 0 : _local_parametrized_functions_for_training[j],
767 : apply_comp_scaling);
768 : }
769 :
770 :
771 0 : correlation_matrix(i,j) = inner_prod;
772 0 : if(i != j)
773 : {
774 0 : correlation_matrix(j,i) = libmesh_conj(inner_prod);
775 : }
776 : }
777 :
778 : // Print out every 10th row so that we can see the progress
779 0 : if ( (i+1) % 10 == 0)
780 0 : std::cout << "Finished row " << (i+1) << " of " << n_snapshots << std::endl;
781 : }
782 0 : std::cout << "Finished computing correlation matrix" << std::endl;
783 :
784 : // Compute SVD of correlation matrix.
785 : // Let Y = U S V^T, then the SVD below corresponds
786 : // to Y^T Y = V S U^T U S V^T = V S^2 V^T.
787 : // The POD basis we use is then given by U, which
788 : // we can compute via U = Y V S^{-1}, which is what
789 : // we compute below.
790 : //
791 : // Note that the formulation remains the same in the
792 : // case that we use a weighted inner product (as
793 : // in the case that we used apply_comp_scaling=true
794 : // when computing the correlation matrix), see (1.28)
795 : // from the lecture notes from Volkwein on POD for more
796 : // details.
797 0 : DenseVector<Real> sigma( n_snapshots );
798 0 : DenseMatrix<Number> U( n_snapshots, n_snapshots );
799 0 : DenseMatrix<Number> VT( n_snapshots, n_snapshots );
800 0 : correlation_matrix.svd(sigma, U, VT );
801 :
802 : // We use this boolean to indicate if we will run one more iteration
803 : // before exiting the loop below. We use this when computing the EIM
804 : // error indicator, which requires one extra EIM iteration.
805 0 : bool exit_on_next_iteration = false;
806 :
807 : // Add dominant vectors from the POD as basis functions.
808 0 : unsigned int j = 0;
809 0 : Real rel_err = 0.;
810 :
811 : // We also initialize a boolean to keep track of whether we have
812 : // reached "n_snapshots" EIM basis functions, since we need to
813 : // handle the EIM error indicator in a special way in this case.
814 0 : bool j_equals_n_snapshots = false;
815 : while (true)
816 : {
817 0 : bool exit_condition_satisfied = false;
818 :
819 0 : if ((j == 0) && (sigma(0) == 0.))
820 : {
821 0 : libMesh::out << "Terminating EIM POD with empty basis because first singular value is zero" << std::endl;
822 0 : exit_condition_satisfied = true;
823 0 : rel_err = 0.;
824 : }
825 0 : else if (j >= n_snapshots)
826 : {
827 0 : libMesh::out << "Number of basis functions equals number of training samples." << std::endl;
828 0 : exit_condition_satisfied = true;
829 0 : j_equals_n_snapshots = true;
830 :
831 : // In this case we set the rel. error to be zero since we've filled up the
832 : // entire space. We cannot use the formula below for rel_err since
833 : // sigma(n_snapshots) is not defined.
834 0 : rel_err = 0.;
835 : }
836 : else
837 : {
838 : // The "energy" error in the POD approximation is determined by the first omitted
839 : // singular value, i.e. sigma(j). We normalize by sigma(0), which gives the total
840 : // "energy", in order to obtain a relative error.
841 0 : rel_err = std::sqrt(sigma(j)) / std::sqrt(sigma(0));
842 : }
843 :
844 0 : if (exit_on_next_iteration)
845 : {
846 0 : libMesh::out << "Extra EIM iteration for error indicator is complete, POD error norm for extra iteration: " << rel_err << std::endl;
847 0 : break;
848 : }
849 :
850 0 : libMesh::out << "Number of basis functions: " << j
851 0 : << ", POD error norm: " << rel_err << std::endl;
852 :
853 0 : if (!exit_condition_satisfied)
854 0 : if (j >= get_Nmax())
855 : {
856 0 : libMesh::out << "Maximum number of basis functions (" << j << ") reached." << std::endl;
857 0 : exit_condition_satisfied = true;
858 : }
859 :
860 0 : if (!exit_condition_satisfied)
861 0 : if (rel_err < get_rel_training_tolerance())
862 : {
863 0 : libMesh::out << "Training tolerance reached." << std::endl;
864 0 : exit_condition_satisfied = true;
865 : }
866 :
867 0 : if (exit_condition_satisfied)
868 : {
869 : // If we're using the EIM error indicator then we need to run
870 : // one extra EIM iteration since we use the extra EIM point
871 : // to obtain our error indicator. If we're not using the EIM
872 : // error indicator, then we just exit now.
873 0 : bool has_parameters = (get_parameters().n_parameters() > 0);
874 0 : if (get_rb_eim_evaluation().use_eim_error_indicator() && has_parameters)
875 : {
876 0 : exit_on_next_iteration = true;
877 0 : libMesh::out << "EIM error indicator is active, hence we will run one extra EIM iteration before exiting"
878 0 : << std::endl;
879 : }
880 : else
881 0 : break;
882 : }
883 :
884 0 : bool is_zero_bf = j_equals_n_snapshots || (rel_err == 0.);
885 0 : if (rbe.get_parametrized_function().on_mesh_sides())
886 : {
887 : // Make a "zero clone" by copying to get the same data layout, and then scaling by zero
888 0 : SideQpDataMap v = _local_side_parametrized_functions_for_training[0];
889 :
890 0 : if (!is_zero_bf)
891 : {
892 0 : scale_parametrized_function(v, 0.);
893 :
894 0 : for ( unsigned int i=0; i<n_snapshots; ++i )
895 0 : add(v, U.el(i, j), _local_side_parametrized_functions_for_training[i] );
896 :
897 0 : Real norm_v = std::sqrt(sigma(j));
898 0 : scale_parametrized_function(v, 1./norm_v);
899 : }
900 :
901 : libmesh_try
902 : {
903 : // If is_zero_bf==true then we add an "extra point" because we cannot
904 : // add a usual EIM interpolation point in that case since the full EIM
905 : // space is already covered. This is necessary when we want to add an
906 : // extra point for error indicator purposes in the is_zero_bf==true
907 : // case, for example.
908 0 : std::unique_ptr<EimPointData> eim_point_data;
909 0 : if (is_zero_bf)
910 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
911 :
912 : // If exit_on_next_iteration==true then we do not add a basis function in
913 : // that case since in that case we only need to add data for the EIM error
914 : // indicator.
915 0 : bool is_linearly_dependent = enrich_eim_approximation_on_sides(v,
916 0 : /*add_basis_function*/ !exit_on_next_iteration,
917 : eim_point_data.get());
918 :
919 : // If we encountered linearly dependent data, then we treat this the same as
920 : // when we have exit_condition_satisfied==true because the EIM training cannot
921 : // proceed any further. As a result, we set exit_on_next_iteration=true here,
922 : // as we do above in the case that exit_condition_satisfied==true.
923 0 : if (is_linearly_dependent)
924 : {
925 0 : bool has_parameters = (get_parameters().n_parameters() > 0);
926 0 : if (get_rb_eim_evaluation().use_eim_error_indicator() && has_parameters)
927 : {
928 0 : exit_on_next_iteration = true;
929 0 : libMesh::out << "Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting"
930 0 : << std::endl;
931 : }
932 : else
933 0 : break;
934 : }
935 :
936 0 : if (is_linearly_dependent && !is_zero_bf)
937 : {
938 : // In this case we detected that v is actually linearly dependent and that is_zero_bf
939 : // was previously not correct --- it should have been true. We typically
940 : // catch this earlier (e.g. by checking rel_err) but in some cases we do not catch
941 : // this until we call the enrichment method. In this situation we update is_zero_bf
942 : // to true and call the enrichment again.
943 0 : is_zero_bf = true;
944 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
945 :
946 0 : enrich_eim_approximation_on_sides(v,
947 0 : /*add_basis_function*/ !exit_on_next_iteration,
948 : eim_point_data.get());
949 : }
950 :
951 0 : update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
952 0 : }
953 0 : libmesh_catch (const std::exception & e)
954 : {
955 : // If we hit an exception when performing the enrichment for the error indicator, then
956 : // we just continue and skip the error indicator. Otherwise we rethrow the exception.
957 0 : if (exit_on_next_iteration)
958 : {
959 0 : std::cout << "Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
960 0 : break;
961 : }
962 : else
963 0 : throw;
964 0 : }
965 : }
966 0 : else if (rbe.get_parametrized_function().on_mesh_nodes())
967 : {
968 : // Make a "zero clone" by copying to get the same data layout, and then scaling by zero
969 0 : NodeDataMap v = _local_node_parametrized_functions_for_training[0];
970 :
971 0 : if (!is_zero_bf)
972 : {
973 0 : scale_node_parametrized_function(v, 0.);
974 :
975 0 : for ( unsigned int i=0; i<n_snapshots; ++i )
976 0 : add_node_data_map(v, U.el(i, j), _local_node_parametrized_functions_for_training[i] );
977 :
978 0 : Real norm_v = std::sqrt(sigma(j));
979 0 : scale_node_parametrized_function(v, 1./norm_v);
980 : }
981 :
982 : libmesh_try
983 : {
984 : // If is_zero_bf==true then we add an "extra point" because we cannot
985 : // add a usual EIM interpolation point in that case since the full EIM
986 : // space is already covered. This is necessary when we want to add an
987 : // extra point for error indicator purposes in the is_zero_bf==true
988 : // case, for example.
989 0 : std::unique_ptr<EimPointData> eim_point_data;
990 0 : if (is_zero_bf)
991 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
992 :
993 : // If exit_on_next_iteration==true then we do not add a basis function in
994 : // that case since in that case we only need to add data for the EIM error
995 : // indicator.
996 0 : bool is_linearly_dependent = enrich_eim_approximation_on_nodes(v,
997 0 : /*add_basis_function*/ !exit_on_next_iteration,
998 : eim_point_data.get());
999 :
1000 : // If we encountered linearly dependent data, then we treat this the same as
1001 : // when we have exit_condition_satisfied==true because the EIM training cannot
1002 : // proceed any further. As a result, we set exit_on_next_iteration=true here,
1003 : // as we do above in the case that exit_condition_satisfied==true.
1004 0 : if (is_linearly_dependent)
1005 : {
1006 0 : bool has_parameters = (get_parameters().n_parameters() > 0);
1007 0 : if (get_rb_eim_evaluation().use_eim_error_indicator() && has_parameters)
1008 : {
1009 0 : exit_on_next_iteration = true;
1010 0 : libMesh::out << "Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting"
1011 0 : << std::endl;
1012 : }
1013 : else
1014 0 : break;
1015 : }
1016 :
1017 0 : if (is_linearly_dependent && !is_zero_bf)
1018 : {
1019 : // In this case we detected that v is actually linearly dependent and that is_zero_bf
1020 : // was previously not correct --- it should have been true. We typically
1021 : // catch this earlier (e.g. by checking rel_err) but in some cases we do not catch
1022 : // this until we call the enrichment method. In this situation we update is_zero_bf
1023 : // to true and call the enrichment again.
1024 0 : is_zero_bf = true;
1025 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
1026 :
1027 0 : enrich_eim_approximation_on_nodes(v,
1028 0 : /*add_basis_function*/ !exit_on_next_iteration,
1029 : eim_point_data.get());
1030 : }
1031 :
1032 0 : update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
1033 0 : }
1034 0 : libmesh_catch (const std::exception & e)
1035 : {
1036 : // If we hit an exception when performing the enrichment for the error indicator, then
1037 : // we just continue and skip the error indicator. Otherwise we rethrow the exception.
1038 0 : if (exit_on_next_iteration)
1039 : {
1040 0 : std::cout << "Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
1041 0 : break;
1042 : }
1043 : else
1044 0 : throw;
1045 0 : }
1046 : }
1047 : else
1048 : {
1049 : // Make a "zero clone" by copying to get the same data layout, and then scaling by zero
1050 0 : QpDataMap v = _local_parametrized_functions_for_training[0];
1051 :
1052 0 : if (!is_zero_bf)
1053 : {
1054 0 : scale_parametrized_function(v, 0.);
1055 :
1056 0 : for ( unsigned int i=0; i<n_snapshots; ++i )
1057 0 : add(v, U.el(i, j), _local_parametrized_functions_for_training[i] );
1058 :
1059 0 : Real norm_v = std::sqrt(sigma(j));
1060 0 : scale_parametrized_function(v, 1./norm_v);
1061 : }
1062 :
1063 : libmesh_try
1064 : {
1065 : // If is_zero_bf==true then we add an "extra point" because we cannot
1066 : // add a usual EIM interpolation point in that case since the full EIM
1067 : // space is already covered. This is necessary when we want to add an
1068 : // extra point for error indicator purposes in the is_zero_bf==true
1069 : // case, for example.
1070 0 : std::unique_ptr<EimPointData> eim_point_data;
1071 0 : if (is_zero_bf)
1072 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
1073 :
1074 : // If exit_on_next_iteration==true then we do not add a basis function in
1075 : // that case since in that case we only need to add data for the EIM error
1076 : // indicator.
1077 0 : bool is_linearly_dependent = enrich_eim_approximation_on_interiors(v,
1078 0 : /*add_basis_function*/ !exit_on_next_iteration,
1079 : eim_point_data.get());
1080 :
1081 : // If we encountered linearly dependent data, then we treat this the same as
1082 : // when we have exit_condition_satisfied==true because the EIM training cannot
1083 : // proceed any further. As a result, we set exit_on_next_iteration=true here,
1084 : // as we do above in the case that exit_condition_satisfied==true.
1085 0 : if (is_linearly_dependent)
1086 : {
1087 0 : bool has_parameters = (get_parameters().n_parameters() > 0);
1088 0 : if (get_rb_eim_evaluation().use_eim_error_indicator() && has_parameters)
1089 : {
1090 0 : exit_on_next_iteration = true;
1091 0 : libMesh::out << "Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting"
1092 0 : << std::endl;
1093 : }
1094 : else
1095 0 : break;
1096 : }
1097 :
1098 0 : if (is_linearly_dependent && !is_zero_bf)
1099 : {
1100 : // In this case we detected that v is actually linearly dependent and that is_zero_bf
1101 : // was previously not correct --- it should have been true. We typically
1102 : // catch this earlier (e.g. by checking rel_err) but in some cases we do not catch
1103 : // this until we call the enrichment method. In this situation we update is_zero_bf
1104 : // to true and call the enrichment again.
1105 0 : is_zero_bf = true;
1106 0 : eim_point_data = std::make_unique<EimPointData>(get_random_point(v));
1107 :
1108 0 : enrich_eim_approximation_on_interiors(v,
1109 0 : /*add_basis_function*/ !exit_on_next_iteration,
1110 : eim_point_data.get());
1111 : }
1112 :
1113 0 : update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
1114 0 : }
1115 0 : libmesh_catch (const std::exception & e)
1116 : {
1117 : // If we hit an exception when performing the enrichment for the error indicator, then
1118 : // we just continue and skip the error indicator. Otherwise we rethrow the exception.
1119 0 : if (exit_on_next_iteration)
1120 : {
1121 0 : std::cout << "Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
1122 0 : break;
1123 : }
1124 : else
1125 0 : throw;
1126 0 : }
1127 : }
1128 :
1129 0 : if (is_zero_bf)
1130 : {
1131 : // In this case we exit here instead of increment j and continuing because
1132 : // if we've encountered a zero EIM basis function then we must not have
1133 : // any more valid data to add.
1134 0 : std::cout << "Zero basis function encountered, hence exiting." << std::endl;
1135 0 : break;
1136 : }
1137 :
1138 0 : j++;
1139 0 : }
1140 0 : libMesh::out << std::endl;
1141 :
1142 0 : return rel_err;
1143 0 : }
1144 :
1145 0 : void RBEIMConstruction::initialize_eim_assembly_objects()
1146 : {
1147 0 : _rb_eim_assembly_objects.clear();
1148 0 : for (auto i : make_range(get_rb_eim_evaluation().get_n_basis_functions()))
1149 0 : _rb_eim_assembly_objects.push_back(build_eim_assembly(i));
1150 0 : }
1151 :
1152 0 : std::vector<std::unique_ptr<ElemAssembly>> & RBEIMConstruction::get_eim_assembly_objects()
1153 : {
1154 0 : return _rb_eim_assembly_objects;
1155 : }
1156 :
1157 0 : void RBEIMConstruction::init_context(FEMContext & c)
1158 : {
1159 : // Pre-request FE data for all element dimensions present in the
1160 : // mesh. Note: we currently pre-request FE data for all variables
1161 : // in the current system but in some cases that may be overkill, for
1162 : // example if only variable 0 is used.
1163 0 : const System & sys = c.get_system();
1164 0 : const MeshBase & mesh = sys.get_mesh();
1165 :
1166 0 : for (unsigned int dim=1; dim<=3; ++dim)
1167 0 : if (mesh.elem_dimensions().count(dim))
1168 0 : for (auto var : make_range(sys.n_vars()))
1169 : {
1170 0 : auto fe = c.get_element_fe(var, dim);
1171 0 : fe->get_JxW();
1172 0 : fe->get_xyz();
1173 0 : fe->get_phi();
1174 :
1175 0 : auto side_fe = c.get_side_fe(var, dim);
1176 0 : side_fe->get_JxW();
1177 0 : side_fe->get_xyz();
1178 0 : side_fe->get_phi();
1179 : }
1180 0 : }
1181 :
1182 0 : void RBEIMConstruction::set_rel_training_tolerance(Real new_training_tolerance)
1183 : {
1184 0 : _rel_training_tolerance = new_training_tolerance;
1185 0 : }
1186 :
1187 0 : Real RBEIMConstruction::get_rel_training_tolerance()
1188 : {
1189 0 : return _rel_training_tolerance;
1190 : }
1191 :
1192 0 : void RBEIMConstruction::set_abs_training_tolerance(Real new_training_tolerance)
1193 : {
1194 0 : _abs_training_tolerance = new_training_tolerance;
1195 0 : }
1196 :
1197 0 : Real RBEIMConstruction::get_abs_training_tolerance()
1198 : {
1199 0 : return _abs_training_tolerance;
1200 : }
1201 :
1202 0 : unsigned int RBEIMConstruction::get_Nmax() const
1203 : {
1204 0 : return _Nmax;
1205 : }
1206 :
1207 0 : void RBEIMConstruction::set_Nmax(unsigned int Nmax)
1208 : {
1209 0 : _Nmax = Nmax;
1210 0 : }
1211 :
1212 0 : void RBEIMConstruction::enable_set_Nmax_from_n_snapshots(int increment)
1213 : {
1214 0 : _set_Nmax_from_n_snapshots = true;
1215 0 : _Nmax_from_n_snapshots_increment = increment;
1216 0 : }
1217 :
1218 0 : void RBEIMConstruction::disable_set_Nmax_from_n_snapshots()
1219 : {
1220 0 : _set_Nmax_from_n_snapshots = false;
1221 0 : _Nmax_from_n_snapshots_increment = 0;
1222 0 : }
1223 :
1224 0 : Real RBEIMConstruction::get_max_abs_value_in_training_set() const
1225 : {
1226 0 : return _max_abs_value_in_training_set;
1227 : }
1228 :
1229 0 : void RBEIMConstruction::store_eim_solutions_for_training_set()
1230 : {
1231 0 : LOG_SCOPE("store_eim_solutions_for_training_set()", "RBEIMConstruction");
1232 :
1233 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
1234 :
1235 0 : std::vector<DenseVector<Number>> & eim_solutions = get_rb_eim_evaluation().get_eim_solutions_for_training_set();
1236 0 : eim_solutions.clear();
1237 0 : eim_solutions.resize(get_n_training_samples());
1238 :
1239 0 : unsigned int RB_size = get_rb_eim_evaluation().get_n_basis_functions();
1240 :
1241 0 : for (auto i : make_range(get_n_training_samples()))
1242 : {
1243 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
1244 : {
1245 0 : const auto & local_side_pf = _local_side_parametrized_functions_for_training[i];
1246 :
1247 0 : if (RB_size > 0)
1248 : {
1249 : // Get the right-hand side vector for the EIM approximation
1250 : // by sampling the parametrized function (stored in solution)
1251 : // at the interpolation points.
1252 0 : DenseVector<Number> EIM_rhs(RB_size);
1253 0 : for (unsigned int j=0; j<RB_size; j++)
1254 : {
1255 0 : EIM_rhs(j) =
1256 0 : RBEIMEvaluation::get_parametrized_function_side_value(comm(),
1257 : local_side_pf,
1258 : eim_eval.get_interpolation_points_elem_id(j),
1259 : eim_eval.get_interpolation_points_side_index(j),
1260 : eim_eval.get_interpolation_points_comp(j),
1261 : eim_eval.get_interpolation_points_qp(j));
1262 : }
1263 0 : eim_solutions[i] = eim_eval.rb_eim_solve(EIM_rhs);
1264 : }
1265 : }
1266 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
1267 : {
1268 0 : const auto & local_node_pf = _local_node_parametrized_functions_for_training[i];
1269 :
1270 0 : if (RB_size > 0)
1271 : {
1272 : // Get the right-hand side vector for the EIM approximation
1273 : // by sampling the parametrized function (stored in solution)
1274 : // at the interpolation points.
1275 0 : DenseVector<Number> EIM_rhs(RB_size);
1276 0 : for (unsigned int j=0; j<RB_size; j++)
1277 : {
1278 0 : EIM_rhs(j) =
1279 0 : RBEIMEvaluation::get_parametrized_function_node_value(comm(),
1280 : local_node_pf,
1281 : eim_eval.get_interpolation_points_node_id(j),
1282 : eim_eval.get_interpolation_points_comp(j));
1283 : }
1284 0 : eim_solutions[i] = eim_eval.rb_eim_solve(EIM_rhs);
1285 : }
1286 : }
1287 : else
1288 : {
1289 0 : const auto & local_pf = _local_parametrized_functions_for_training[i];
1290 :
1291 0 : if (RB_size > 0)
1292 : {
1293 : // Get the right-hand side vector for the EIM approximation
1294 : // by sampling the parametrized function (stored in solution)
1295 : // at the interpolation points.
1296 0 : DenseVector<Number> EIM_rhs(RB_size);
1297 0 : for (unsigned int j=0; j<RB_size; j++)
1298 : {
1299 0 : EIM_rhs(j) =
1300 0 : RBEIMEvaluation::get_parametrized_function_value(comm(),
1301 : local_pf,
1302 : eim_eval.get_interpolation_points_elem_id(j),
1303 : eim_eval.get_interpolation_points_comp(j),
1304 : eim_eval.get_interpolation_points_qp(j));
1305 : }
1306 0 : eim_solutions[i] = eim_eval.rb_eim_solve(EIM_rhs);
1307 : }
1308 : }
1309 : }
1310 0 : }
1311 :
1312 0 : const RBEIMEvaluation::QpDataMap & RBEIMConstruction::get_parametrized_function_from_training_set(unsigned int training_index) const
1313 : {
1314 0 : libmesh_error_msg_if(training_index >= _local_parametrized_functions_for_training.size(),
1315 : "Invalid index: " << training_index);
1316 0 : return _local_parametrized_functions_for_training[training_index];
1317 : }
1318 :
1319 0 : const RBEIMEvaluation::SideQpDataMap & RBEIMConstruction::get_side_parametrized_function_from_training_set(unsigned int training_index) const
1320 : {
1321 0 : libmesh_error_msg_if(training_index >= _local_side_parametrized_functions_for_training.size(),
1322 : "Invalid index: " << training_index);
1323 0 : return _local_side_parametrized_functions_for_training[training_index];
1324 : }
1325 :
1326 0 : const RBEIMEvaluation::NodeDataMap & RBEIMConstruction::get_node_parametrized_function_from_training_set(unsigned int training_index) const
1327 : {
1328 0 : libmesh_error_msg_if(training_index >= _local_node_parametrized_functions_for_training.size(),
1329 : "Invalid index: " << training_index);
1330 0 : return _local_node_parametrized_functions_for_training[training_index];
1331 : }
1332 :
1333 0 : const std::unordered_map<dof_id_type, std::vector<Real> > & RBEIMConstruction::get_local_quad_point_JxW()
1334 : {
1335 0 : return _local_quad_point_JxW;
1336 : }
1337 :
1338 0 : const std::map<std::pair<dof_id_type,unsigned int>, std::vector<Real> > & RBEIMConstruction::get_local_side_quad_point_JxW()
1339 : {
1340 0 : return _local_side_quad_point_JxW;
1341 : }
1342 :
1343 0 : unsigned int RBEIMConstruction::get_n_parametrized_functions_for_training() const
1344 : {
1345 0 : if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_sides())
1346 0 : return _local_side_parametrized_functions_for_training.size();
1347 0 : else if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_nodes())
1348 0 : return _local_node_parametrized_functions_for_training.size();
1349 : else
1350 0 : return _local_parametrized_functions_for_training.size();
1351 : }
1352 :
1353 0 : void RBEIMConstruction::reinit_eim_projection_matrix()
1354 : {
1355 0 : RBEIMEvaluation & rbe = get_rb_eim_evaluation();
1356 :
1357 : // We need space for one extra interpolation point if we're using the
1358 : // EIM error indicator.
1359 0 : unsigned int max_matrix_size = rbe.use_eim_error_indicator() ? get_Nmax()+1 : get_Nmax();
1360 0 : _eim_projection_matrix.resize(max_matrix_size,max_matrix_size);
1361 0 : }
1362 :
1363 0 : std::pair<Real,unsigned int> RBEIMConstruction::compute_max_eim_error()
1364 : {
1365 0 : LOG_SCOPE("compute_max_eim_error()", "RBEIMConstruction");
1366 :
1367 0 : if (get_n_params() == 0)
1368 : {
1369 : // Just return 0 if we have no parameters.
1370 0 : return std::make_pair(0.,0);
1371 : }
1372 :
1373 : // keep track of the maximum error
1374 0 : unsigned int max_err_index = 0;
1375 0 : Real max_err = 0.;
1376 :
1377 0 : libmesh_error_msg_if(get_n_training_samples() != get_local_n_training_samples(),
1378 : "Error: Training samples should be the same on all procs");
1379 :
1380 0 : const unsigned int RB_size = get_rb_eim_evaluation().get_n_basis_functions();
1381 :
1382 0 : if(best_fit_type_flag == PROJECTION_BEST_FIT)
1383 : {
1384 0 : for (auto training_index : make_range(get_n_training_samples()))
1385 : {
1386 0 : if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_sides())
1387 : {
1388 : // Make a copy of the pre-computed solution for the specified training sample
1389 : // since we will modify it below to compute the best fit error.
1390 0 : SideQpDataMap solution_copy = _local_side_parametrized_functions_for_training[training_index];
1391 :
1392 : // Perform an L2 projection in order to find the best approximation to
1393 : // the parametrized function from the current EIM space.
1394 0 : DenseVector<Number> best_fit_rhs(RB_size);
1395 0 : for (unsigned int i=0; i<RB_size; i++)
1396 : {
1397 0 : best_fit_rhs(i) = side_inner_product(solution_copy,
1398 0 : get_rb_eim_evaluation().get_side_basis_function(i),
1399 : /*apply_comp_scaling*/ false);
1400 : }
1401 :
1402 : // Now compute the best fit by an LU solve
1403 0 : DenseMatrix<Number> RB_inner_product_matrix_N(RB_size);
1404 0 : _eim_projection_matrix.get_principal_submatrix(RB_size, RB_inner_product_matrix_N);
1405 :
1406 0 : DenseVector<Number> best_fit_coeffs;
1407 0 : RB_inner_product_matrix_N.lu_solve(best_fit_rhs, best_fit_coeffs);
1408 :
1409 0 : get_rb_eim_evaluation().side_decrement_vector(solution_copy, best_fit_coeffs);
1410 0 : Real best_fit_error = get_max_abs_value(solution_copy);
1411 :
1412 0 : if (best_fit_error > max_err)
1413 : {
1414 0 : max_err_index = training_index;
1415 0 : max_err = best_fit_error;
1416 : }
1417 0 : }
1418 0 : else if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_nodes())
1419 : {
1420 : // Make a copy of the pre-computed solution for the specified training sample
1421 : // since we will modify it below to compute the best fit error.
1422 0 : NodeDataMap solution_copy = _local_node_parametrized_functions_for_training[training_index];
1423 :
1424 : // Perform an L2 projection in order to find the best approximation to
1425 : // the parametrized function from the current EIM space.
1426 0 : DenseVector<Number> best_fit_rhs(RB_size);
1427 0 : for (unsigned int i=0; i<RB_size; i++)
1428 : {
1429 0 : best_fit_rhs(i) = node_inner_product(solution_copy,
1430 0 : get_rb_eim_evaluation().get_node_basis_function(i),
1431 : /*apply_comp_scaling*/ false);
1432 : }
1433 :
1434 : // Now compute the best fit by an LU solve
1435 0 : DenseMatrix<Number> RB_inner_product_matrix_N(RB_size);
1436 0 : _eim_projection_matrix.get_principal_submatrix(RB_size, RB_inner_product_matrix_N);
1437 :
1438 0 : DenseVector<Number> best_fit_coeffs;
1439 0 : RB_inner_product_matrix_N.lu_solve(best_fit_rhs, best_fit_coeffs);
1440 :
1441 0 : get_rb_eim_evaluation().node_decrement_vector(solution_copy, best_fit_coeffs);
1442 0 : Real best_fit_error = get_node_max_abs_value(solution_copy);
1443 :
1444 0 : if (best_fit_error > max_err)
1445 : {
1446 0 : max_err_index = training_index;
1447 0 : max_err = best_fit_error;
1448 : }
1449 0 : }
1450 : else
1451 : {
1452 : // Make a copy of the pre-computed solution for the specified training sample
1453 : // since we will modify it below to compute the best fit error.
1454 0 : QpDataMap solution_copy = _local_parametrized_functions_for_training[training_index];
1455 :
1456 : // Perform an L2 projection in order to find the best approximation to
1457 : // the parametrized function from the current EIM space.
1458 0 : DenseVector<Number> best_fit_rhs(RB_size);
1459 0 : for (unsigned int i=0; i<RB_size; i++)
1460 : {
1461 0 : best_fit_rhs(i) = inner_product(solution_copy,
1462 0 : get_rb_eim_evaluation().get_basis_function(i),
1463 : /*apply_comp_scaling*/ false);
1464 : }
1465 :
1466 : // Now compute the best fit by an LU solve
1467 0 : DenseMatrix<Number> RB_inner_product_matrix_N(RB_size);
1468 0 : _eim_projection_matrix.get_principal_submatrix(RB_size, RB_inner_product_matrix_N);
1469 :
1470 0 : DenseVector<Number> best_fit_coeffs;
1471 0 : RB_inner_product_matrix_N.lu_solve(best_fit_rhs, best_fit_coeffs);
1472 :
1473 0 : get_rb_eim_evaluation().decrement_vector(solution_copy, best_fit_coeffs);
1474 0 : Real best_fit_error = get_max_abs_value(solution_copy);
1475 :
1476 0 : if (best_fit_error > max_err)
1477 : {
1478 0 : max_err_index = training_index;
1479 0 : max_err = best_fit_error;
1480 : }
1481 0 : }
1482 : }
1483 : }
1484 0 : else if(best_fit_type_flag == EIM_BEST_FIT)
1485 : {
1486 : // Perform EIM solve in order to find the approximation to solution
1487 : // (rb_eim_solve provides the EIM basis function coefficients used below)
1488 :
1489 0 : std::vector<RBParameters> training_parameters_copy(get_n_training_samples());
1490 0 : for (auto training_index : make_range(get_n_training_samples()))
1491 : {
1492 0 : training_parameters_copy[training_index] = get_params_from_training_set(training_index);
1493 : }
1494 :
1495 0 : get_rb_eim_evaluation().rb_eim_solves(training_parameters_copy, RB_size);
1496 0 : const std::vector<DenseVector<Number>> & rb_eim_solutions = get_rb_eim_evaluation().get_rb_eim_solutions();
1497 :
1498 0 : for (auto training_index : make_range(get_n_training_samples()))
1499 : {
1500 0 : const DenseVector<Number> & best_fit_coeffs = rb_eim_solutions[training_index];
1501 :
1502 0 : if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_sides())
1503 : {
1504 0 : SideQpDataMap solution_copy = _local_side_parametrized_functions_for_training[training_index];
1505 0 : get_rb_eim_evaluation().side_decrement_vector(solution_copy, best_fit_coeffs);
1506 0 : Real best_fit_error = get_max_abs_value(solution_copy);
1507 :
1508 0 : if (best_fit_error > max_err)
1509 : {
1510 0 : max_err_index = training_index;
1511 0 : max_err = best_fit_error;
1512 : }
1513 : }
1514 0 : else if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_nodes())
1515 : {
1516 0 : NodeDataMap solution_copy = _local_node_parametrized_functions_for_training[training_index];
1517 0 : get_rb_eim_evaluation().node_decrement_vector(solution_copy, best_fit_coeffs);
1518 0 : Real best_fit_error = get_node_max_abs_value(solution_copy);
1519 :
1520 0 : if (best_fit_error > max_err)
1521 : {
1522 0 : max_err_index = training_index;
1523 0 : max_err = best_fit_error;
1524 : }
1525 : }
1526 : else
1527 : {
1528 0 : QpDataMap solution_copy = _local_parametrized_functions_for_training[training_index];
1529 0 : get_rb_eim_evaluation().decrement_vector(solution_copy, best_fit_coeffs);
1530 0 : Real best_fit_error = get_max_abs_value(solution_copy);
1531 :
1532 0 : if (best_fit_error > max_err)
1533 : {
1534 0 : max_err_index = training_index;
1535 0 : max_err = best_fit_error;
1536 : }
1537 : }
1538 : }
1539 0 : }
1540 : else
1541 : {
1542 0 : libmesh_error_msg("EIM best fit type not recognized");
1543 : }
1544 :
1545 0 : return std::make_pair(max_err,max_err_index);
1546 : }
1547 :
1548 0 : void RBEIMConstruction::initialize_parametrized_functions_in_training_set()
1549 : {
1550 0 : LOG_SCOPE("initialize_parametrized_functions_in_training_set()", "RBEIMConstruction");
1551 :
1552 0 : libmesh_error_msg_if(!serial_training_set,
1553 : "Error: We must have serial_training_set==true in "
1554 : "RBEIMConstruction::initialize_parametrized_functions_in_training_set");
1555 :
1556 0 : libMesh::out << "Initializing parametrized functions in training set..." << std::endl;
1557 :
1558 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
1559 :
1560 0 : if (eim_eval.get_parametrized_function().is_lookup_table)
1561 0 : eim_eval.get_parametrized_function().initialize_lookup_table();
1562 :
1563 : // Store the locations of all quadrature points
1564 0 : initialize_qp_data();
1565 :
1566 : // Keep track of the largest value in our parametrized functions
1567 : // in the training set. We can use this value for normalization
1568 : // purposes, for example.
1569 0 : _max_abs_value_in_training_set = 0.;
1570 :
1571 0 : unsigned int n_comps = eim_eval.get_parametrized_function().get_n_components();
1572 :
1573 : // Keep track of the maximum value per component. This will allow
1574 : // us to scale the components to all have a similar magnitude,
1575 : // which is helpful during the error assessment for the basis
1576 : // enrichment to ensure that components with smaller magnitude
1577 : // are not ignored.
1578 0 : std::vector<Real> max_abs_value_per_component_in_training_set(n_comps);
1579 :
1580 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
1581 : {
1582 0 : _local_side_parametrized_functions_for_training.resize( get_n_training_samples() );
1583 0 : for (auto i : make_range(get_n_training_samples()))
1584 : {
1585 0 : libMesh::out << "Initializing parametrized function for training sample "
1586 0 : << (i+1) << " of " << get_n_training_samples() << std::endl;
1587 :
1588 0 : set_params_from_training_set(i);
1589 :
1590 0 : eim_eval.get_parametrized_function().preevaluate_parametrized_function_on_mesh_sides(get_parameters(),
1591 0 : _local_side_quad_point_locations,
1592 0 : _local_side_quad_point_subdomain_ids,
1593 0 : _local_side_quad_point_boundary_ids,
1594 0 : _local_side_quad_point_side_types,
1595 0 : _local_side_quad_point_locations_perturbations,
1596 0 : *this);
1597 :
1598 0 : for (const auto & [elem_side_pair, xyz_vector] : _local_side_quad_point_locations)
1599 : {
1600 0 : std::vector<std::vector<Number>> comps_and_qps(n_comps);
1601 0 : for (unsigned int comp=0; comp<n_comps; comp++)
1602 : {
1603 0 : comps_and_qps[comp].resize(xyz_vector.size());
1604 0 : for (unsigned int qp : index_range(xyz_vector))
1605 : {
1606 : Number value =
1607 0 : eim_eval.get_parametrized_function().lookup_preevaluated_side_value_on_mesh(comp,
1608 0 : elem_side_pair.first,
1609 0 : elem_side_pair.second,
1610 0 : qp);
1611 0 : comps_and_qps[comp][qp] = value;
1612 :
1613 0 : Real abs_value = std::abs(value);
1614 0 : if (abs_value > _max_abs_value_in_training_set)
1615 : {
1616 0 : _max_abs_value_in_training_set = abs_value;
1617 0 : _max_abs_value_in_training_set_index = i;
1618 : }
1619 :
1620 0 : if (abs_value > max_abs_value_per_component_in_training_set[comp])
1621 0 : max_abs_value_per_component_in_training_set[comp] = abs_value;
1622 : }
1623 : }
1624 :
1625 0 : _local_side_parametrized_functions_for_training[i][elem_side_pair] = comps_and_qps;
1626 0 : }
1627 : }
1628 :
1629 0 : libMesh::out << "Parametrized functions in training set initialized" << std::endl;
1630 :
1631 0 : unsigned int max_id = 0;
1632 0 : comm().maxloc(_max_abs_value_in_training_set, max_id);
1633 0 : comm().broadcast(_max_abs_value_in_training_set_index, max_id);
1634 0 : libMesh::out << "Maximum absolute value in the training set: "
1635 0 : << _max_abs_value_in_training_set << std::endl << std::endl;
1636 :
1637 : // Calculate the maximum value for each component in the training set
1638 : // across all components
1639 0 : comm().max(max_abs_value_per_component_in_training_set);
1640 :
1641 : // We store the maximum value across all components divided by the maximum value for this component
1642 : // so that when we scale using these factors all components should have a magnitude on the same
1643 : // order as the maximum component.
1644 0 : _component_scaling_in_training_set.resize(n_comps);
1645 0 : for (unsigned int i : make_range(n_comps))
1646 : {
1647 0 : if ((eim_eval.scale_components_in_enrichment().count(i) == 0) ||
1648 0 : max_abs_value_per_component_in_training_set[i] == 0.)
1649 0 : _component_scaling_in_training_set[i] = 1.;
1650 : else
1651 0 : _component_scaling_in_training_set[i] = _max_abs_value_in_training_set / max_abs_value_per_component_in_training_set[i];
1652 : }
1653 : }
1654 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
1655 : {
1656 0 : _local_node_parametrized_functions_for_training.resize( get_n_training_samples() );
1657 0 : for (auto i : make_range(get_n_training_samples()))
1658 : {
1659 0 : libMesh::out << "Initializing parametrized function for training sample "
1660 0 : << (i+1) << " of " << get_n_training_samples() << std::endl;
1661 :
1662 0 : set_params_from_training_set(i);
1663 :
1664 0 : eim_eval.get_parametrized_function().preevaluate_parametrized_function_on_mesh_nodes(get_parameters(),
1665 0 : _local_node_locations,
1666 0 : _local_node_boundary_ids,
1667 0 : *this);
1668 :
1669 0 : for (const auto & pr : _local_node_locations)
1670 : {
1671 0 : const auto & node_id = pr.first;
1672 :
1673 0 : std::vector<Number> comps(n_comps);
1674 0 : for (unsigned int comp=0; comp<n_comps; comp++)
1675 : {
1676 : Number value =
1677 0 : eim_eval.get_parametrized_function().lookup_preevaluated_node_value_on_mesh(comp,
1678 0 : node_id);
1679 0 : comps[comp] = value;
1680 :
1681 0 : Real abs_value = std::abs(value);
1682 0 : if (abs_value > _max_abs_value_in_training_set)
1683 : {
1684 0 : _max_abs_value_in_training_set = abs_value;
1685 0 : _max_abs_value_in_training_set_index = i;
1686 : }
1687 :
1688 0 : if (abs_value > max_abs_value_per_component_in_training_set[comp])
1689 0 : max_abs_value_per_component_in_training_set[comp] = abs_value;
1690 : }
1691 :
1692 0 : _local_node_parametrized_functions_for_training[i][node_id] = comps;
1693 : }
1694 : }
1695 :
1696 0 : libMesh::out << "Parametrized functions in training set initialized" << std::endl;
1697 :
1698 0 : unsigned int max_id = 0;
1699 0 : comm().maxloc(_max_abs_value_in_training_set, max_id);
1700 0 : comm().broadcast(_max_abs_value_in_training_set_index, max_id);
1701 0 : libMesh::out << "Maximum absolute value in the training set: "
1702 0 : << _max_abs_value_in_training_set << std::endl << std::endl;
1703 :
1704 : // Calculate the maximum value for each component in the training set
1705 : // across all components
1706 0 : comm().max(max_abs_value_per_component_in_training_set);
1707 :
1708 : // We store the maximum value across all components divided by the maximum value for this component
1709 : // so that when we scale using these factors all components should have a magnitude on the same
1710 : // order as the maximum component.
1711 0 : _component_scaling_in_training_set.resize(n_comps);
1712 0 : for (unsigned int i : make_range(n_comps))
1713 : {
1714 0 : if ((eim_eval.scale_components_in_enrichment().count(i) == 0) ||
1715 0 : max_abs_value_per_component_in_training_set[i] == 0.)
1716 0 : _component_scaling_in_training_set[i] = 1.;
1717 : else
1718 0 : _component_scaling_in_training_set[i] = _max_abs_value_in_training_set / max_abs_value_per_component_in_training_set[i];
1719 : }
1720 : }
1721 : else
1722 : {
1723 0 : _local_parametrized_functions_for_training.resize( get_n_training_samples() );
1724 0 : for (auto i : make_range(get_n_training_samples()))
1725 : {
1726 0 : libMesh::out << "Initializing parametrized function for training sample "
1727 0 : << (i+1) << " of " << get_n_training_samples() << std::endl;
1728 :
1729 0 : set_params_from_training_set(i);
1730 :
1731 0 : eim_eval.get_parametrized_function().preevaluate_parametrized_function_on_mesh(get_parameters(),
1732 0 : _local_quad_point_locations,
1733 0 : _local_quad_point_subdomain_ids,
1734 0 : _local_quad_point_locations_perturbations,
1735 0 : *this);
1736 :
1737 0 : for (const auto & [elem_id, xyz_vector] : _local_quad_point_locations)
1738 : {
1739 0 : std::vector<std::vector<Number>> comps_and_qps(n_comps);
1740 0 : for (unsigned int comp=0; comp<n_comps; comp++)
1741 : {
1742 0 : comps_and_qps[comp].resize(xyz_vector.size());
1743 0 : for (unsigned int qp : index_range(xyz_vector))
1744 : {
1745 : Number value =
1746 0 : eim_eval.get_parametrized_function().lookup_preevaluated_value_on_mesh(comp, elem_id, qp);
1747 0 : comps_and_qps[comp][qp] = value;
1748 :
1749 0 : Real abs_value = std::abs(value);
1750 0 : if (abs_value > _max_abs_value_in_training_set)
1751 : {
1752 0 : _max_abs_value_in_training_set = abs_value;
1753 0 : _max_abs_value_in_training_set_index = i;
1754 : }
1755 :
1756 0 : if (abs_value > max_abs_value_per_component_in_training_set[comp])
1757 0 : max_abs_value_per_component_in_training_set[comp] = abs_value;
1758 : }
1759 : }
1760 :
1761 0 : _local_parametrized_functions_for_training[i][elem_id] = comps_and_qps;
1762 0 : }
1763 : }
1764 :
1765 0 : libMesh::out << "Parametrized functions in training set initialized" << std::endl;
1766 :
1767 0 : unsigned int max_id = 0;
1768 0 : comm().maxloc(_max_abs_value_in_training_set, max_id);
1769 0 : comm().broadcast(_max_abs_value_in_training_set_index, max_id);
1770 0 : libMesh::out << "Maximum absolute value in the training set: "
1771 0 : << _max_abs_value_in_training_set << std::endl << std::endl;
1772 :
1773 : // Calculate the maximum value for each component in the training set
1774 : // across all components
1775 0 : comm().max(max_abs_value_per_component_in_training_set);
1776 :
1777 : // We store the maximum value across all components divided by the maximum value for this component
1778 : // so that when we scale using these factors all components should have a magnitude on the same
1779 : // order as the maximum component.
1780 0 : _component_scaling_in_training_set.resize(n_comps);
1781 0 : for (unsigned int i : make_range(n_comps))
1782 : {
1783 0 : if ((eim_eval.scale_components_in_enrichment().count(i) == 0) ||
1784 0 : max_abs_value_per_component_in_training_set[i] == 0.)
1785 0 : _component_scaling_in_training_set[i] = 1.;
1786 : else
1787 0 : _component_scaling_in_training_set[i] = _max_abs_value_in_training_set / max_abs_value_per_component_in_training_set[i];
1788 : }
1789 : }
1790 : // This function does nothing if rb_property_map from RBParametrizedFunction
1791 : // is empty which would result in an empty rb_property_map in VectorizedEvalInput
1792 : // stored in RBEIMEvaluation.
1793 0 : eim_eval.initialize_rb_property_map();
1794 0 : }
1795 :
1796 0 : void RBEIMConstruction::initialize_qp_data()
1797 : {
1798 0 : LOG_SCOPE("initialize_qp_data()", "RBEIMConstruction");
1799 :
1800 0 : if (!get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
1801 : {
1802 0 : libMesh::out << "Initializing quadrature point locations" << std::endl;
1803 : }
1804 : else
1805 : {
1806 0 : libMesh::out << "Initializing quadrature point and perturbation locations" << std::endl;
1807 : }
1808 :
1809 : // Compute truth representation via L2 projection
1810 0 : const MeshBase & mesh = this->get_mesh();
1811 :
1812 0 : FEMContext context(*this);
1813 0 : init_context(context);
1814 :
1815 0 : if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_sides())
1816 : {
1817 : const std::set<boundary_id_type> & parametrized_function_boundary_ids =
1818 0 : get_rb_eim_evaluation().get_parametrized_function().get_parametrized_function_boundary_ids();
1819 0 : libmesh_error_msg_if (parametrized_function_boundary_ids.empty(),
1820 : "Need to have non-empty boundary IDs to initialize side data");
1821 :
1822 0 : _local_side_quad_point_locations.clear();
1823 0 : _local_side_quad_point_subdomain_ids.clear();
1824 0 : _local_side_quad_point_boundary_ids.clear();
1825 0 : _local_side_quad_point_JxW.clear();
1826 0 : _local_side_quad_point_side_types.clear();
1827 :
1828 0 : _local_side_quad_point_locations_perturbations.clear();
1829 :
1830 : // BoundaryInfo and related data structures
1831 0 : const auto & binfo = mesh.get_boundary_info();
1832 0 : std::vector<boundary_id_type> side_boundary_ids;
1833 :
1834 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
1835 : {
1836 0 : dof_id_type elem_id = elem->id();
1837 :
1838 0 : context.pre_fe_reinit(*this, elem);
1839 :
1840 : // elem_fe is used for shellface data
1841 0 : auto elem_fe = context.get_element_fe(/*var=*/0, elem->dim());
1842 0 : const std::vector<Real> & JxW = elem_fe->get_JxW();
1843 0 : const std::vector<Point> & xyz = elem_fe->get_xyz();
1844 :
1845 : // side_fe is used for element side data
1846 0 : auto side_fe = context.get_side_fe(/*var=*/0, elem->dim());
1847 0 : const std::vector<Real> & JxW_side = side_fe->get_JxW();
1848 0 : const std::vector< Point > & xyz_side = side_fe->get_xyz();
1849 :
1850 0 : for (context.side = 0;
1851 0 : context.side != context.get_elem().n_sides();
1852 0 : ++context.side)
1853 : {
1854 : // skip non-boundary elements
1855 0 : if(!context.get_elem().neighbor_ptr(context.side))
1856 : {
1857 0 : binfo.boundary_ids(elem, context.side, side_boundary_ids);
1858 :
1859 0 : bool has_side_boundary_id = false;
1860 0 : boundary_id_type matching_boundary_id = BoundaryInfo::invalid_id;
1861 0 : for (boundary_id_type side_boundary_id : side_boundary_ids)
1862 0 : if(parametrized_function_boundary_ids.count(side_boundary_id))
1863 : {
1864 0 : has_side_boundary_id = true;
1865 0 : matching_boundary_id = side_boundary_id;
1866 0 : break;
1867 : }
1868 :
1869 0 : if(has_side_boundary_id)
1870 : {
1871 0 : context.get_side_fe(/*var=*/0, elem->dim())->reinit(elem, context.side);
1872 :
1873 0 : auto elem_side_pair = std::make_pair(elem_id, context.side);
1874 :
1875 0 : _local_side_quad_point_locations[elem_side_pair] = xyz_side;
1876 0 : _local_side_quad_point_JxW[elem_side_pair] = JxW_side;
1877 0 : _local_side_quad_point_subdomain_ids[elem_side_pair] = elem->subdomain_id();
1878 0 : _local_side_quad_point_boundary_ids[elem_side_pair] = matching_boundary_id;
1879 :
1880 : // This is a standard side (not a shellface) so set side type to 0
1881 0 : _local_side_quad_point_side_types[elem_side_pair] = 0;
1882 :
1883 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
1884 : {
1885 0 : Real fd_delta = get_rb_eim_evaluation().get_parametrized_function().fd_delta;
1886 :
1887 0 : std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
1888 :
1889 0 : for (const Point & xyz_qp : xyz_side)
1890 : {
1891 0 : std::vector<Point> xyz_perturb_vec;
1892 0 : if (elem->dim() == 3)
1893 : {
1894 : // In this case we have a 3D element, and hence the side is 2D.
1895 : //
1896 : // We use the following approach to perturb xyz:
1897 : // 1) inverse map xyz to the reference element
1898 : // 2) perturb on the reference element in the (xi,eta) "directions"
1899 : // 3) map the perturbed points back to the physical element
1900 : // This approach is necessary to ensure that the perturbed points
1901 : // are still in the element's side.
1902 :
1903 0 : std::unique_ptr<const Elem> elem_side;
1904 0 : elem->build_side_ptr(elem_side, context.side);
1905 :
1906 : Point xi_eta =
1907 0 : FEMap::inverse_map(elem_side->dim(),
1908 : elem_side.get(),
1909 : xyz_qp,
1910 : /*Newton iteration tolerance*/ TOLERANCE,
1911 0 : /*secure*/ true);
1912 :
1913 : // Inverse map should map back to a 2D reference domain
1914 0 : libmesh_assert(std::abs(xi_eta(2)) < TOLERANCE);
1915 :
1916 0 : Point xi_eta_perturb = xi_eta;
1917 :
1918 0 : xi_eta_perturb(0) += fd_delta;
1919 : Point xyz_perturb_0 =
1920 0 : FEMap::map(elem_side->dim(),
1921 : elem_side.get(),
1922 0 : xi_eta_perturb);
1923 0 : xi_eta_perturb(0) -= fd_delta;
1924 :
1925 0 : xi_eta_perturb(1) += fd_delta;
1926 : Point xyz_perturb_1 =
1927 0 : FEMap::map(elem_side->dim(),
1928 : elem_side.get(),
1929 0 : xi_eta_perturb);
1930 0 : xi_eta_perturb(1) -= fd_delta;
1931 :
1932 : // Finally, we rescale xyz_perturb_0 and xyz_perturb_1 so that
1933 : // (xyz_perturb - xyz_qp).norm() == fd_delta, since this is
1934 : // required in order to compute finite differences correctly.
1935 0 : Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
1936 0 : Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
1937 :
1938 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
1939 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
1940 0 : }
1941 : else
1942 : {
1943 : // We current do nothing for sides of dim=2 or dim=1 elements
1944 : // since we have no need for this capability so far.
1945 : // Support for these cases could be added if it is needed.
1946 : }
1947 :
1948 0 : xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
1949 : }
1950 :
1951 0 : _local_side_quad_point_locations_perturbations[elem_side_pair] = xyz_perturb_vec_at_qps;
1952 0 : }
1953 : }
1954 : }
1955 : }
1956 :
1957 : // In the case of 2D elements, we also check the shellfaces
1958 0 : if (elem->dim() == 2)
1959 0 : for (unsigned int shellface_index=0; shellface_index<2; shellface_index++)
1960 : {
1961 0 : binfo.shellface_boundary_ids(elem, shellface_index, side_boundary_ids);
1962 :
1963 0 : bool has_side_boundary_id = false;
1964 0 : boundary_id_type matching_boundary_id = BoundaryInfo::invalid_id;
1965 0 : for (boundary_id_type side_boundary_id : side_boundary_ids)
1966 0 : if(parametrized_function_boundary_ids.count(side_boundary_id))
1967 : {
1968 0 : has_side_boundary_id = true;
1969 0 : matching_boundary_id = side_boundary_id;
1970 0 : break;
1971 : }
1972 :
1973 0 : if(has_side_boundary_id)
1974 : {
1975 0 : context.elem_fe_reinit();
1976 :
1977 : // We use shellface_index as the side_index since shellface boundary conditions
1978 : // are stored separately from side boundary conditions in BoundaryInfo.
1979 0 : auto elem_side_pair = std::make_pair(elem_id, shellface_index);
1980 :
1981 0 : _local_side_quad_point_locations[elem_side_pair] = xyz;
1982 0 : _local_side_quad_point_JxW[elem_side_pair] = JxW;
1983 0 : _local_side_quad_point_subdomain_ids[elem_side_pair] = elem->subdomain_id();
1984 0 : _local_side_quad_point_boundary_ids[elem_side_pair] = matching_boundary_id;
1985 :
1986 : // This is a shellface (not a standard side) so set side type to 1
1987 0 : _local_side_quad_point_side_types[elem_side_pair] = 1;
1988 :
1989 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
1990 : {
1991 0 : Real fd_delta = get_rb_eim_evaluation().get_parametrized_function().fd_delta;
1992 :
1993 0 : std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
1994 :
1995 0 : for (const Point & xyz_qp : xyz)
1996 : {
1997 0 : std::vector<Point> xyz_perturb_vec;
1998 : // Here we follow the same approach as above for getting xyz_perturb_vec,
1999 : // except that we are using the element itself instead of its side.
2000 : {
2001 : Point xi_eta =
2002 0 : FEMap::inverse_map(elem->dim(),
2003 : elem,
2004 : xyz_qp,
2005 : /*Newton iteration tolerance*/ TOLERANCE,
2006 0 : /*secure*/ true);
2007 :
2008 : // Inverse map should map back to a 2D reference domain
2009 0 : libmesh_assert(std::abs(xi_eta(2)) < TOLERANCE);
2010 :
2011 0 : Point xi_eta_perturb = xi_eta;
2012 :
2013 0 : xi_eta_perturb(0) += fd_delta;
2014 : Point xyz_perturb_0 =
2015 0 : FEMap::map(elem->dim(),
2016 : elem,
2017 0 : xi_eta_perturb);
2018 0 : xi_eta_perturb(0) -= fd_delta;
2019 :
2020 0 : xi_eta_perturb(1) += fd_delta;
2021 : Point xyz_perturb_1 =
2022 0 : FEMap::map(elem->dim(),
2023 : elem,
2024 0 : xi_eta_perturb);
2025 0 : xi_eta_perturb(1) -= fd_delta;
2026 :
2027 : // Finally, we rescale xyz_perturb_0 and xyz_perturb_1 so that
2028 : // (xyz_perturb - xyz_qp).norm() == fd_delta, since this is
2029 : // required in order to compute finite differences correctly.
2030 0 : Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
2031 0 : Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
2032 :
2033 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
2034 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
2035 : }
2036 :
2037 0 : xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
2038 : }
2039 :
2040 0 : _local_side_quad_point_locations_perturbations[elem_side_pair] = xyz_perturb_vec_at_qps;
2041 0 : }
2042 : }
2043 : }
2044 0 : }
2045 : }
2046 0 : else if (get_rb_eim_evaluation().get_parametrized_function().on_mesh_nodes())
2047 : {
2048 : const std::set<boundary_id_type> & parametrized_function_boundary_ids =
2049 0 : get_rb_eim_evaluation().get_parametrized_function().get_parametrized_function_boundary_ids();
2050 0 : libmesh_error_msg_if (parametrized_function_boundary_ids.empty(),
2051 : "Need to have non-empty boundary IDs to initialize node data");
2052 :
2053 0 : _local_node_locations.clear();
2054 0 : _local_node_boundary_ids.clear();
2055 :
2056 0 : const auto & binfo = mesh.get_boundary_info();
2057 :
2058 : // Make a set with all the nodes that have nodesets. Use
2059 : // a set so that we don't have any duplicate entries. We
2060 : // deal with duplicate entries below by getting all boundary
2061 : // IDs on each node.
2062 0 : std::set<dof_id_type> nodes_with_nodesets;
2063 0 : for (const auto & t : binfo.build_node_list())
2064 0 : nodes_with_nodesets.insert(std::get<0>(t));
2065 :
2066 : // To be filled in by BoundaryInfo calls in loop below
2067 0 : std::vector<boundary_id_type> node_boundary_ids;
2068 :
2069 0 : for (dof_id_type node_id : nodes_with_nodesets)
2070 : {
2071 0 : const Node * node = mesh.node_ptr(node_id);
2072 :
2073 0 : if (node->processor_id() != mesh.comm().rank())
2074 0 : continue;
2075 :
2076 0 : binfo.boundary_ids(node, node_boundary_ids);
2077 :
2078 0 : bool has_node_boundary_id = false;
2079 0 : boundary_id_type matching_boundary_id = BoundaryInfo::invalid_id;
2080 0 : for (boundary_id_type node_boundary_id : node_boundary_ids)
2081 0 : if(parametrized_function_boundary_ids.count(node_boundary_id))
2082 : {
2083 0 : has_node_boundary_id = true;
2084 0 : matching_boundary_id = node_boundary_id;
2085 0 : break;
2086 : }
2087 :
2088 0 : if(has_node_boundary_id)
2089 : {
2090 0 : _local_node_locations[node_id] = *node;
2091 0 : _local_node_boundary_ids[node_id] = matching_boundary_id;
2092 : }
2093 : }
2094 : }
2095 : else
2096 : {
2097 0 : _local_quad_point_locations.clear();
2098 0 : _local_quad_point_subdomain_ids.clear();
2099 0 : _local_quad_point_JxW.clear();
2100 :
2101 0 : _local_quad_point_locations_perturbations.clear();
2102 :
2103 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
2104 : {
2105 0 : auto elem_fe = context.get_element_fe(/*var=*/0, elem->dim());
2106 0 : const std::vector<Real> & JxW = elem_fe->get_JxW();
2107 0 : const std::vector<Point> & xyz = elem_fe->get_xyz();
2108 :
2109 0 : dof_id_type elem_id = elem->id();
2110 :
2111 0 : context.pre_fe_reinit(*this, elem);
2112 0 : context.elem_fe_reinit();
2113 :
2114 0 : _local_quad_point_locations[elem_id] = xyz;
2115 0 : _local_quad_point_JxW[elem_id] = JxW;
2116 0 : _local_quad_point_subdomain_ids[elem_id] = elem->subdomain_id();
2117 :
2118 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
2119 : {
2120 0 : Real fd_delta = get_rb_eim_evaluation().get_parametrized_function().fd_delta;
2121 :
2122 0 : std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
2123 :
2124 0 : for (const Point & xyz_qp : xyz)
2125 : {
2126 0 : std::vector<Point> xyz_perturb_vec;
2127 0 : if (elem->dim() == 3)
2128 : {
2129 0 : Point xyz_perturb = xyz_qp;
2130 :
2131 0 : xyz_perturb(0) += fd_delta;
2132 0 : xyz_perturb_vec.emplace_back(xyz_perturb);
2133 0 : xyz_perturb(0) -= fd_delta;
2134 :
2135 0 : xyz_perturb(1) += fd_delta;
2136 0 : xyz_perturb_vec.emplace_back(xyz_perturb);
2137 0 : xyz_perturb(1) -= fd_delta;
2138 :
2139 0 : xyz_perturb(2) += fd_delta;
2140 0 : xyz_perturb_vec.emplace_back(xyz_perturb);
2141 0 : xyz_perturb(2) -= fd_delta;
2142 : }
2143 0 : else if (elem->dim() == 2)
2144 : {
2145 : // In this case we assume that we have a 2D element
2146 : // embedded in 3D space. In this case we have to use
2147 : // the following approach to perturb xyz:
2148 : // 1) inverse map xyz to the reference element
2149 : // 2) perturb on the reference element in the (xi,eta) "directions"
2150 : // 3) map the perturbed points back to the physical element
2151 : // This approach is necessary to ensure that the perturbed points
2152 : // are still in the element.
2153 :
2154 : Point xi_eta =
2155 0 : FEMap::inverse_map(elem->dim(),
2156 : elem,
2157 : xyz_qp,
2158 : /*Newton iteration tolerance*/ TOLERANCE,
2159 0 : /*secure*/ true);
2160 :
2161 : // Inverse map should map back to a 2D reference domain
2162 0 : libmesh_assert(std::abs(xi_eta(2)) < TOLERANCE);
2163 :
2164 0 : Point xi_eta_perturb = xi_eta;
2165 :
2166 0 : xi_eta_perturb(0) += fd_delta;
2167 : Point xyz_perturb_0 =
2168 0 : FEMap::map(elem->dim(),
2169 : elem,
2170 0 : xi_eta_perturb);
2171 0 : xi_eta_perturb(0) -= fd_delta;
2172 :
2173 0 : xi_eta_perturb(1) += fd_delta;
2174 : Point xyz_perturb_1 =
2175 0 : FEMap::map(elem->dim(),
2176 : elem,
2177 0 : xi_eta_perturb);
2178 0 : xi_eta_perturb(1) -= fd_delta;
2179 :
2180 : // Finally, we rescale xyz_perturb_0 and xyz_perturb_1 so that
2181 : // (xyz_perturb - xyz_qp).norm() == fd_delta, since this is
2182 : // required in order to compute finite differences correctly.
2183 0 : Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
2184 0 : Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
2185 :
2186 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
2187 0 : xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
2188 : }
2189 : else
2190 : {
2191 : // We current do nothing in the dim=1 case since
2192 : // we have no need for this capability so far.
2193 : // Support for this case could be added if it is
2194 : // needed.
2195 : }
2196 :
2197 0 : xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
2198 : }
2199 :
2200 0 : _local_quad_point_locations_perturbations[elem_id] = xyz_perturb_vec_at_qps;
2201 0 : }
2202 0 : }
2203 : }
2204 0 : }
2205 :
2206 : Number
2207 0 : RBEIMConstruction::inner_product(const QpDataMap & v, const QpDataMap & w, bool apply_comp_scaling)
2208 : {
2209 0 : LOG_SCOPE("inner_product()", "RBEIMConstruction");
2210 :
2211 0 : Number val = 0.;
2212 :
2213 0 : for (const auto & [elem_id, v_comp_and_qp] : v)
2214 : {
2215 0 : const auto & w_comp_and_qp = libmesh_map_find(w, elem_id);
2216 0 : const auto & JxW = libmesh_map_find(_local_quad_point_JxW, elem_id);
2217 :
2218 0 : for (const auto & comp : index_range(v_comp_and_qp))
2219 : {
2220 0 : const std::vector<Number> & v_qp = v_comp_and_qp[comp];
2221 0 : const std::vector<Number> & w_qp = w_comp_and_qp[comp];
2222 :
2223 0 : Real comp_scaling = 1.;
2224 0 : if (apply_comp_scaling)
2225 : {
2226 : // We square the component scaling here because it occurs twice in
2227 : // the inner product calculation below.
2228 0 : comp_scaling = std::pow(_component_scaling_in_training_set[comp], 2.);
2229 : }
2230 :
2231 0 : for (unsigned int qp : index_range(JxW))
2232 0 : val += JxW[qp] * comp_scaling * v_qp[qp] * libmesh_conj(w_qp[qp]);
2233 : }
2234 : }
2235 :
2236 0 : comm().sum(val);
2237 0 : return val;
2238 : }
2239 :
2240 : Number
2241 0 : RBEIMConstruction::side_inner_product(const SideQpDataMap & v, const SideQpDataMap & w, bool apply_comp_scaling)
2242 : {
2243 0 : LOG_SCOPE("side_inner_product()", "RBEIMConstruction");
2244 :
2245 0 : Number val = 0.;
2246 :
2247 0 : for (const auto & [elem_and_side, v_comp_and_qp] : v)
2248 : {
2249 0 : const auto & w_comp_and_qp = libmesh_map_find(w, elem_and_side);
2250 0 : const auto & JxW = libmesh_map_find(_local_side_quad_point_JxW, elem_and_side);
2251 :
2252 0 : for (const auto & comp : index_range(v_comp_and_qp))
2253 : {
2254 0 : const std::vector<Number> & v_qp = v_comp_and_qp[comp];
2255 0 : const std::vector<Number> & w_qp = w_comp_and_qp[comp];
2256 :
2257 0 : Real comp_scaling = 1.;
2258 0 : if (apply_comp_scaling)
2259 : {
2260 : // We square the component scaling here because it occurs twice in
2261 : // the inner product calculation below.
2262 0 : comp_scaling = std::pow(_component_scaling_in_training_set[comp], 2.);
2263 : }
2264 :
2265 0 : for (unsigned int qp : index_range(JxW))
2266 0 : val += JxW[qp] * comp_scaling * v_qp[qp] * libmesh_conj(w_qp[qp]);
2267 : }
2268 : }
2269 :
2270 0 : comm().sum(val);
2271 0 : return val;
2272 : }
2273 :
2274 : Number
2275 0 : RBEIMConstruction::node_inner_product(const NodeDataMap & v, const NodeDataMap & w, bool apply_comp_scaling)
2276 : {
2277 0 : LOG_SCOPE("node_inner_product()", "RBEIMConstruction");
2278 :
2279 0 : Number val = 0.;
2280 :
2281 0 : for (const auto & [node_id, v_comps] : v)
2282 : {
2283 0 : const auto & w_comps = libmesh_map_find(w, node_id);
2284 :
2285 0 : for (const auto & comp : index_range(v_comps))
2286 : {
2287 : // There is no quadrature rule on nodes, so we just multiply the values directly.
2288 : // Hence we effectively work with the Euclidean inner product in this case.
2289 :
2290 0 : Real comp_scaling = 1.;
2291 0 : if (apply_comp_scaling)
2292 : {
2293 : // We square the component scaling here because it occurs twice in
2294 : // the inner product calculation below.
2295 0 : comp_scaling = std::pow(_component_scaling_in_training_set[comp], 2.);
2296 : }
2297 :
2298 0 : val += comp_scaling * v_comps[comp] * libmesh_conj(w_comps[comp]);
2299 : }
2300 : }
2301 :
2302 0 : comm().sum(val);
2303 0 : return val;
2304 : }
2305 :
2306 0 : Real RBEIMConstruction::get_node_max_abs_value(const NodeDataMap & v) const
2307 : {
2308 0 : Real max_value = 0.;
2309 :
2310 0 : for (const auto & pr : v)
2311 : {
2312 0 : const auto & values = pr.second;
2313 0 : for (const auto & comp : index_range(values))
2314 : {
2315 0 : const auto & value = values[comp];
2316 :
2317 0 : Real comp_scaling = 1.;
2318 0 : if (get_rb_eim_evaluation().scale_components_in_enrichment().count(comp))
2319 : {
2320 : // Make sure that _component_scaling_in_training_set is initialized
2321 0 : libmesh_error_msg_if(comp >= _component_scaling_in_training_set.size(),
2322 : "Invalid vector index");
2323 0 : comp_scaling = _component_scaling_in_training_set[comp];
2324 : }
2325 :
2326 0 : max_value = std::max(max_value, std::abs(value * comp_scaling));
2327 : }
2328 : }
2329 :
2330 0 : comm().max(max_value);
2331 0 : return max_value;
2332 : }
2333 :
2334 0 : void RBEIMConstruction::enrich_eim_approximation(unsigned int training_index,
2335 : bool add_basis_function,
2336 : EimPointData * eim_point_data)
2337 : {
2338 0 : LOG_SCOPE("enrich_eim_approximation()", "RBEIMConstruction");
2339 :
2340 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
2341 :
2342 0 : set_params_from_training_set(training_index);
2343 :
2344 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
2345 0 : enrich_eim_approximation_on_sides(_local_side_parametrized_functions_for_training[training_index],
2346 : add_basis_function,
2347 : eim_point_data);
2348 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
2349 0 : enrich_eim_approximation_on_nodes(_local_node_parametrized_functions_for_training[training_index],
2350 : add_basis_function,
2351 : eim_point_data);
2352 : else
2353 : {
2354 0 : enrich_eim_approximation_on_interiors(_local_parametrized_functions_for_training[training_index],
2355 : add_basis_function,
2356 : eim_point_data);
2357 : }
2358 0 : }
2359 :
2360 0 : bool RBEIMConstruction::enrich_eim_approximation_on_sides(const SideQpDataMap & side_pf,
2361 : bool add_basis_function,
2362 : EimPointData * eim_point_data)
2363 : {
2364 : // Make a copy of the input parametrized function, since we will modify this below
2365 : // to give us a new basis function.
2366 0 : SideQpDataMap local_pf = side_pf;
2367 :
2368 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
2369 :
2370 : // If we have at least one basis function, then we need to use
2371 : // rb_eim_solve() to find the EIM interpolation error. Otherwise,
2372 : // just use solution as is.
2373 0 : if (!eim_point_data && (eim_eval.get_n_basis_functions() > 0))
2374 : {
2375 : // Get the right-hand side vector for the EIM approximation
2376 : // by sampling the parametrized function (stored in solution)
2377 : // at the interpolation points.
2378 0 : unsigned int RB_size = eim_eval.get_n_basis_functions();
2379 0 : DenseVector<Number> EIM_rhs(RB_size);
2380 0 : for (unsigned int i=0; i<RB_size; i++)
2381 : {
2382 0 : EIM_rhs(i) =
2383 0 : RBEIMEvaluation::get_parametrized_function_side_value(comm(),
2384 : local_pf,
2385 : eim_eval.get_interpolation_points_elem_id(i),
2386 : eim_eval.get_interpolation_points_side_index(i),
2387 : eim_eval.get_interpolation_points_comp(i),
2388 : eim_eval.get_interpolation_points_qp(i));
2389 : }
2390 :
2391 0 : eim_eval.set_parameters( get_parameters() );
2392 0 : DenseVector<Number> rb_eim_solution = eim_eval.rb_eim_solve(EIM_rhs);
2393 :
2394 : // Load the "EIM residual" into solution by subtracting
2395 : // the EIM approximation
2396 0 : eim_eval.side_decrement_vector(local_pf, rb_eim_solution);
2397 : }
2398 :
2399 : // Find the quadrature point at which local_pf (which now stores
2400 : // the "EIM residual") has maximum absolute value
2401 0 : Number optimal_value = 0.;
2402 0 : Point optimal_point;
2403 0 : unsigned int optimal_comp = 0;
2404 0 : dof_id_type optimal_elem_id = DofObject::invalid_id;
2405 0 : unsigned int optimal_side_index = 0;
2406 0 : subdomain_id_type optimal_subdomain_id = 0;
2407 0 : boundary_id_type optimal_boundary_id = 0;
2408 0 : unsigned int optimal_qp = 0;
2409 0 : std::vector<Point> optimal_point_perturbs;
2410 0 : std::vector<Real> optimal_point_phi_i_qp;
2411 :
2412 : // Initialize largest_abs_value to be negative so that it definitely gets updated.
2413 0 : Real largest_abs_value = -1.;
2414 :
2415 : // In order to compute phi_i_qp, we initialize a FEMContext
2416 0 : FEMContext con(*this);
2417 0 : init_context(con);
2418 :
2419 0 : for (const auto & [elem_and_side, comp_and_qp] : local_pf)
2420 : {
2421 0 : dof_id_type elem_id = elem_and_side.first;
2422 0 : unsigned int side_index = elem_and_side.second;
2423 :
2424 0 : const Elem & elem_ref = get_mesh().elem_ref(elem_id);
2425 0 : con.pre_fe_reinit(*this, &elem_ref);
2426 :
2427 0 : unsigned int side_type = libmesh_map_find(_local_side_quad_point_side_types, elem_and_side);
2428 :
2429 0 : std::vector<std::vector<Real>> phi;
2430 : // side_type == 0 --> standard side
2431 : // side_type == 1 --> shellface
2432 0 : if (side_type == 0)
2433 : {
2434 : // TODO: We only want the "dofs on side" entries
2435 : // from phi_side. Could do this by initing an FE object
2436 : // on the side itself, rather than using get_side_fe().
2437 0 : auto side_fe = con.get_side_fe(/*var=*/ 0);
2438 0 : side_fe->reinit(&elem_ref, side_index);
2439 :
2440 0 : phi = side_fe->get_phi();
2441 : }
2442 0 : else if (side_type == 1)
2443 : {
2444 0 : con.elem_fe_reinit();
2445 :
2446 0 : auto elem_fe = con.get_element_fe(/*var=*/0, elem_ref.dim());
2447 0 : phi = elem_fe->get_phi();
2448 : }
2449 : else
2450 0 : libmesh_error_msg ("Unrecognized side_type: " << side_type);
2451 :
2452 0 : for (const auto & comp : index_range(comp_and_qp))
2453 : {
2454 0 : const std::vector<Number> & qp_values = comp_and_qp[comp];
2455 :
2456 0 : for (auto qp : index_range(qp_values))
2457 : {
2458 0 : Number value = qp_values[qp];
2459 0 : Real abs_value = std::abs(value);
2460 :
2461 0 : bool update_optimal_point = false;
2462 0 : if (!eim_point_data)
2463 0 : update_optimal_point = (abs_value > largest_abs_value);
2464 : else
2465 0 : update_optimal_point = (elem_id == eim_point_data->elem_id) &&
2466 0 : (side_index == eim_point_data->side_index) &&
2467 0 : (comp == eim_point_data->comp_index) &&
2468 0 : (qp == eim_point_data->qp_index);
2469 :
2470 0 : if (update_optimal_point)
2471 : {
2472 0 : largest_abs_value = abs_value;
2473 0 : optimal_value = value;
2474 0 : optimal_comp = comp;
2475 0 : optimal_elem_id = elem_id;
2476 0 : optimal_side_index = side_index;
2477 0 : optimal_qp = qp;
2478 :
2479 0 : optimal_point_phi_i_qp.resize(phi.size());
2480 0 : for (auto i : index_range(phi))
2481 0 : optimal_point_phi_i_qp[i] = phi[i][qp];
2482 :
2483 : const auto & point_list =
2484 0 : libmesh_map_find(_local_side_quad_point_locations, elem_and_side);
2485 :
2486 0 : libmesh_error_msg_if(qp >= point_list.size(), "Error: Invalid qp");
2487 :
2488 0 : optimal_point = point_list[qp];
2489 :
2490 0 : optimal_subdomain_id = libmesh_map_find(_local_side_quad_point_subdomain_ids, elem_and_side);
2491 0 : optimal_boundary_id = libmesh_map_find(_local_side_quad_point_boundary_ids, elem_and_side);
2492 :
2493 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
2494 : {
2495 : const auto & perturb_list =
2496 0 : libmesh_map_find(_local_side_quad_point_locations_perturbations, elem_and_side);
2497 :
2498 0 : libmesh_error_msg_if(qp >= perturb_list.size(), "Error: Invalid qp");
2499 :
2500 0 : optimal_point_perturbs = perturb_list[qp];
2501 : }
2502 : }
2503 : }
2504 : }
2505 0 : }
2506 :
2507 : // Find out which processor has the largest of the abs values
2508 : // and broadcast from that processor.
2509 : unsigned int proc_ID_index;
2510 0 : this->comm().maxloc(largest_abs_value, proc_ID_index);
2511 :
2512 0 : this->comm().broadcast(optimal_value, proc_ID_index);
2513 0 : this->comm().broadcast(optimal_point, proc_ID_index);
2514 0 : this->comm().broadcast(optimal_comp, proc_ID_index);
2515 0 : this->comm().broadcast(optimal_elem_id, proc_ID_index);
2516 0 : this->comm().broadcast(optimal_side_index, proc_ID_index);
2517 0 : this->comm().broadcast(optimal_subdomain_id, proc_ID_index);
2518 0 : this->comm().broadcast(optimal_boundary_id, proc_ID_index);
2519 0 : this->comm().broadcast(optimal_qp, proc_ID_index);
2520 0 : this->comm().broadcast(optimal_point_perturbs, proc_ID_index);
2521 0 : this->comm().broadcast(optimal_point_phi_i_qp, proc_ID_index);
2522 :
2523 0 : libmesh_error_msg_if(optimal_elem_id == DofObject::invalid_id, "Error: Invalid element ID");
2524 :
2525 0 : if (add_basis_function)
2526 : {
2527 0 : if (optimal_value == 0.)
2528 : {
2529 0 : libMesh::out << "Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2530 0 : return true;
2531 : }
2532 :
2533 : // Scale local_pf so that its largest value is 1.0
2534 0 : scale_parametrized_function(local_pf, 1./optimal_value);
2535 :
2536 : // Add local_pf as the new basis function and store data
2537 : // associated with the interpolation point.
2538 0 : eim_eval.add_side_basis_function(local_pf);
2539 : }
2540 :
2541 0 : eim_eval.add_side_interpolation_data(optimal_point,
2542 : optimal_comp,
2543 : optimal_elem_id,
2544 : optimal_side_index,
2545 : optimal_subdomain_id,
2546 : optimal_boundary_id,
2547 : optimal_qp,
2548 : optimal_point_perturbs,
2549 : optimal_point_phi_i_qp);
2550 :
2551 : // In this case we did not encounter a linearly dependent basis function, so return false
2552 0 : return false;
2553 0 : }
2554 :
2555 0 : bool RBEIMConstruction::enrich_eim_approximation_on_nodes(const NodeDataMap & node_pf,
2556 : bool add_basis_function,
2557 : EimPointData * eim_point_data)
2558 : {
2559 : // Make a copy of the input parametrized function, since we will modify this below
2560 : // to give us a new basis function.
2561 0 : NodeDataMap local_pf = node_pf;
2562 :
2563 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
2564 :
2565 : // If we have at least one basis function, then we need to use
2566 : // rb_eim_solve() to find the EIM interpolation error. Otherwise,
2567 : // just use solution as is.
2568 0 : if (!eim_point_data && (eim_eval.get_n_basis_functions() > 0))
2569 : {
2570 : // Get the right-hand side vector for the EIM approximation
2571 : // by sampling the parametrized function (stored in solution)
2572 : // at the interpolation points.
2573 0 : unsigned int RB_size = eim_eval.get_n_basis_functions();
2574 0 : DenseVector<Number> EIM_rhs(RB_size);
2575 0 : for (unsigned int i=0; i<RB_size; i++)
2576 : {
2577 0 : EIM_rhs(i) =
2578 0 : RBEIMEvaluation::get_parametrized_function_node_value(comm(),
2579 : local_pf,
2580 : eim_eval.get_interpolation_points_node_id(i),
2581 : eim_eval.get_interpolation_points_comp(i));
2582 : }
2583 :
2584 0 : eim_eval.set_parameters( get_parameters() );
2585 0 : DenseVector<Number> rb_eim_solution = eim_eval.rb_eim_solve(EIM_rhs);
2586 :
2587 : // Load the "EIM residual" into solution by subtracting
2588 : // the EIM approximation
2589 0 : eim_eval.node_decrement_vector(local_pf, rb_eim_solution);
2590 : }
2591 :
2592 : // Find the quadrature point at which local_pf (which now stores
2593 : // the "EIM residual") has maximum absolute value
2594 0 : Number optimal_value = 0.;
2595 0 : Point optimal_point;
2596 0 : unsigned int optimal_comp = 0;
2597 0 : dof_id_type optimal_node_id = DofObject::invalid_id;
2598 0 : boundary_id_type optimal_boundary_id = 0;
2599 :
2600 : // Initialize largest_abs_value to be negative so that it definitely gets updated.
2601 0 : Real largest_abs_value = -1.;
2602 :
2603 0 : for (const auto & [node_id, values] : local_pf)
2604 : {
2605 0 : for (unsigned int comp : index_range(values))
2606 : {
2607 0 : Number value = values[comp];
2608 0 : Real abs_value = std::abs(value);
2609 :
2610 0 : bool update_optimal_point = false;
2611 0 : if (!eim_point_data)
2612 0 : update_optimal_point = (abs_value > largest_abs_value);
2613 : else
2614 0 : update_optimal_point = (node_id == eim_point_data->node_id) &&
2615 0 : (comp == eim_point_data->comp_index);
2616 :
2617 0 : if (update_optimal_point)
2618 : {
2619 0 : largest_abs_value = abs_value;
2620 0 : optimal_value = value;
2621 0 : optimal_comp = comp;
2622 0 : optimal_node_id = node_id;
2623 :
2624 0 : optimal_point = libmesh_map_find(_local_node_locations, node_id);
2625 :
2626 0 : optimal_boundary_id = libmesh_map_find(_local_node_boundary_ids, node_id);
2627 : }
2628 : }
2629 : }
2630 :
2631 : // Find out which processor has the largest of the abs values
2632 : // and broadcast from that processor.
2633 : unsigned int proc_ID_index;
2634 0 : this->comm().maxloc(largest_abs_value, proc_ID_index);
2635 :
2636 0 : this->comm().broadcast(optimal_value, proc_ID_index);
2637 0 : this->comm().broadcast(optimal_point, proc_ID_index);
2638 0 : this->comm().broadcast(optimal_comp, proc_ID_index);
2639 0 : this->comm().broadcast(optimal_node_id, proc_ID_index);
2640 0 : this->comm().broadcast(optimal_boundary_id, proc_ID_index);
2641 :
2642 0 : libmesh_error_msg_if(optimal_node_id == DofObject::invalid_id, "Error: Invalid node ID");
2643 :
2644 0 : if (add_basis_function)
2645 : {
2646 0 : if (optimal_value == 0.)
2647 : {
2648 0 : libMesh::out << "Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2649 0 : return true;
2650 : }
2651 :
2652 : // Scale local_pf so that its largest value is 1.0
2653 0 : scale_node_parametrized_function(local_pf, 1./optimal_value);
2654 :
2655 : // Add local_pf as the new basis function and store data
2656 : // associated with the interpolation point.
2657 0 : eim_eval.add_node_basis_function(local_pf);
2658 : }
2659 :
2660 0 : eim_eval.add_node_interpolation_data(optimal_point,
2661 : optimal_comp,
2662 : optimal_node_id,
2663 : optimal_boundary_id);
2664 :
2665 : // In this case we did not encounter a linearly dependent basis function, so return false
2666 0 : return false;
2667 : }
2668 :
2669 0 : bool RBEIMConstruction::enrich_eim_approximation_on_interiors(const QpDataMap & interior_pf,
2670 : bool add_basis_function,
2671 : EimPointData * eim_point_data)
2672 : {
2673 : // Make a copy of the input parametrized function, since we will modify this below
2674 : // to give us a new basis function.
2675 0 : QpDataMap local_pf = interior_pf;
2676 :
2677 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
2678 :
2679 : // If we have at least one basis function, then we need to use
2680 : // rb_eim_solve() to find the EIM interpolation error. Otherwise,
2681 : // just use solution as is.
2682 0 : if (!eim_point_data && (eim_eval.get_n_basis_functions() > 0))
2683 : {
2684 : // Get the right-hand side vector for the EIM approximation
2685 : // by sampling the parametrized function (stored in solution)
2686 : // at the interpolation points.
2687 0 : unsigned int RB_size = eim_eval.get_n_basis_functions();
2688 0 : DenseVector<Number> EIM_rhs(RB_size);
2689 0 : for (unsigned int i=0; i<RB_size; i++)
2690 : {
2691 0 : EIM_rhs(i) =
2692 0 : RBEIMEvaluation::get_parametrized_function_value(comm(),
2693 : local_pf,
2694 : eim_eval.get_interpolation_points_elem_id(i),
2695 : eim_eval.get_interpolation_points_comp(i),
2696 : eim_eval.get_interpolation_points_qp(i));
2697 : }
2698 :
2699 0 : eim_eval.set_parameters( get_parameters() );
2700 0 : DenseVector<Number> rb_eim_solution = eim_eval.rb_eim_solve(EIM_rhs);
2701 :
2702 : // Load the "EIM residual" into solution by subtracting
2703 : // the EIM approximation
2704 0 : eim_eval.decrement_vector(local_pf, rb_eim_solution);
2705 : }
2706 :
2707 : // Find the quadrature point at which local_pf (which now stores
2708 : // the "EIM residual") has maximum absolute value
2709 0 : Number optimal_value = 0.;
2710 0 : Point optimal_point;
2711 0 : unsigned int optimal_comp = 0;
2712 0 : dof_id_type optimal_elem_id = DofObject::invalid_id;
2713 0 : subdomain_id_type optimal_subdomain_id = 0;
2714 0 : unsigned int optimal_qp = 0;
2715 0 : std::vector<Point> optimal_point_perturbs;
2716 0 : std::vector<Real> optimal_point_phi_i_qp;
2717 0 : ElemType optimal_elem_type = INVALID_ELEM;
2718 0 : std::vector<Real> optimal_JxW_all_qp;
2719 0 : std::vector<std::vector<Real>> optimal_phi_i_all_qp;
2720 0 : Order optimal_qrule_order = INVALID_ORDER;
2721 0 : Point optimal_dxyzdxi_elem_center;
2722 0 : Point optimal_dxyzdeta_elem_center;
2723 :
2724 : // Initialize largest_abs_value to be negative so that it definitely gets updated.
2725 0 : Real largest_abs_value = -1.;
2726 :
2727 : // In order to compute phi_i_qp, we initialize a FEMContext
2728 0 : FEMContext con(*this);
2729 0 : for (auto dim : con.elem_dimensions())
2730 : {
2731 0 : auto fe = con.get_element_fe(/*var=*/0, dim);
2732 0 : fe->get_phi();
2733 0 : fe->get_JxW();
2734 0 : fe->get_dxyzdxi();
2735 0 : fe->get_dxyzdeta();
2736 : }
2737 :
2738 0 : for (const auto & [elem_id, comp_and_qp] : local_pf)
2739 : {
2740 : // Also initialize phi in order to compute phi_i_qp
2741 0 : const Elem & elem_ref = get_mesh().elem_ref(elem_id);
2742 0 : con.pre_fe_reinit(*this, &elem_ref);
2743 :
2744 0 : auto elem_fe = con.get_element_fe(/*var=*/0, elem_ref.dim());
2745 0 : const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
2746 0 : const auto & JxW = elem_fe->get_JxW();
2747 0 : const auto & dxyzdxi = elem_fe->get_dxyzdxi();
2748 0 : const auto & dxyzdeta = elem_fe->get_dxyzdeta();
2749 :
2750 0 : elem_fe->reinit(&elem_ref);
2751 :
2752 0 : for (const auto & comp : index_range(comp_and_qp))
2753 : {
2754 0 : const std::vector<Number> & qp_values = comp_and_qp[comp];
2755 :
2756 0 : for (auto qp : index_range(qp_values))
2757 : {
2758 0 : Number value = qp_values[qp];
2759 0 : Real abs_value = std::abs(value);
2760 :
2761 0 : if (get_rb_eim_evaluation().scale_components_in_enrichment().count(comp))
2762 0 : abs_value *= _component_scaling_in_training_set[comp];
2763 :
2764 0 : bool update_optimal_point = false;
2765 0 : if (!eim_point_data)
2766 0 : update_optimal_point = (abs_value > largest_abs_value);
2767 : else
2768 0 : update_optimal_point = (elem_id == eim_point_data->elem_id) &&
2769 0 : (comp == eim_point_data->comp_index) &&
2770 0 : (qp == eim_point_data->qp_index);
2771 :
2772 0 : if (update_optimal_point)
2773 : {
2774 0 : largest_abs_value = abs_value;
2775 0 : optimal_value = value;
2776 0 : optimal_comp = comp;
2777 0 : optimal_elem_id = elem_id;
2778 0 : optimal_qp = qp;
2779 0 : optimal_elem_type = elem_ref.type();
2780 :
2781 0 : optimal_point_phi_i_qp.resize(phi.size());
2782 0 : for (auto i : index_range(phi))
2783 0 : optimal_point_phi_i_qp[i] = phi[i][qp];
2784 :
2785 : const auto & point_list =
2786 0 : libmesh_map_find(_local_quad_point_locations, elem_id);
2787 :
2788 0 : libmesh_error_msg_if(qp >= point_list.size(), "Error: Invalid qp");
2789 :
2790 0 : optimal_point = point_list[qp];
2791 :
2792 0 : optimal_subdomain_id = libmesh_map_find(_local_quad_point_subdomain_ids, elem_id);
2793 :
2794 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_xyz_perturbations)
2795 : {
2796 : const auto & perturb_list =
2797 0 : libmesh_map_find(_local_quad_point_locations_perturbations, elem_id);
2798 :
2799 0 : libmesh_error_msg_if(qp >= perturb_list.size(), "Error: Invalid qp");
2800 :
2801 0 : optimal_point_perturbs = perturb_list[qp];
2802 : }
2803 :
2804 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_all_elem_qp_data)
2805 : {
2806 0 : optimal_JxW_all_qp = JxW;
2807 0 : optimal_phi_i_all_qp = phi;
2808 : }
2809 :
2810 0 : if (get_rb_eim_evaluation().get_parametrized_function().requires_all_elem_center_data)
2811 : {
2812 0 : optimal_qrule_order = con.get_element_qrule().get_order();
2813 : // Get data derivatives at vertex average
2814 0 : std::vector<Point> nodes = { elem_ref.reference_elem()->vertex_average() };
2815 0 : elem_fe->reinit (&elem_ref, &nodes);
2816 :
2817 0 : Point dxyzdxi_pt, dxyzdeta_pt;
2818 0 : if (con.get_elem_dim()>0)
2819 0 : dxyzdxi_pt = dxyzdxi[0];
2820 0 : if (con.get_elem_dim()>1)
2821 0 : dxyzdeta_pt = dxyzdeta[0];
2822 :
2823 0 : optimal_dxyzdxi_elem_center = dxyzdxi_pt;
2824 0 : optimal_dxyzdeta_elem_center = dxyzdeta_pt;
2825 :
2826 0 : elem_fe->reinit(&elem_ref);
2827 : }
2828 : }
2829 : }
2830 : }
2831 : }
2832 :
2833 : // Find out which processor has the largest of the abs values
2834 : // and broadcast from that processor.
2835 : unsigned int proc_ID_index;
2836 0 : this->comm().maxloc(largest_abs_value, proc_ID_index);
2837 :
2838 0 : this->comm().broadcast(optimal_value, proc_ID_index);
2839 0 : this->comm().broadcast(optimal_point, proc_ID_index);
2840 0 : this->comm().broadcast(optimal_comp, proc_ID_index);
2841 0 : this->comm().broadcast(optimal_elem_id, proc_ID_index);
2842 0 : this->comm().broadcast(optimal_subdomain_id, proc_ID_index);
2843 0 : this->comm().broadcast(optimal_qp, proc_ID_index);
2844 0 : this->comm().broadcast(optimal_point_perturbs, proc_ID_index);
2845 0 : this->comm().broadcast(optimal_point_phi_i_qp, proc_ID_index);
2846 0 : this->comm().broadcast(optimal_JxW_all_qp, proc_ID_index);
2847 0 : this->comm().broadcast(optimal_phi_i_all_qp, proc_ID_index);
2848 0 : this->comm().broadcast(optimal_dxyzdxi_elem_center, proc_ID_index);
2849 0 : this->comm().broadcast(optimal_dxyzdeta_elem_center, proc_ID_index);
2850 :
2851 : // Cast optimal_elem_type to an int in order to broadcast it
2852 : {
2853 0 : int optimal_elem_type_int = static_cast<int>(optimal_elem_type);
2854 0 : this->comm().broadcast(optimal_elem_type_int, proc_ID_index);
2855 0 : optimal_elem_type = static_cast<ElemType>(optimal_elem_type_int);
2856 : }
2857 :
2858 : // Cast optimal_qrule_order to an int in order to broadcast it
2859 : {
2860 0 : int optimal_qrule_order_int = static_cast<int>(optimal_qrule_order);
2861 0 : this->comm().broadcast(optimal_qrule_order_int, proc_ID_index);
2862 0 : optimal_qrule_order = static_cast<Order>(optimal_qrule_order_int);
2863 : }
2864 :
2865 0 : libmesh_error_msg_if(optimal_elem_id == DofObject::invalid_id, "Error: Invalid element ID");
2866 :
2867 0 : if (add_basis_function)
2868 : {
2869 0 : if (optimal_value == 0.)
2870 : {
2871 0 : libMesh::out << "Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2872 0 : return true;
2873 : }
2874 :
2875 : // Scale local_pf so that its largest value is 1.0
2876 0 : scale_parametrized_function(local_pf, 1./optimal_value);
2877 :
2878 : // Add local_pf as the new basis function and store data
2879 : // associated with the interpolation point.
2880 0 : eim_eval.add_basis_function(local_pf);
2881 : }
2882 :
2883 0 : eim_eval.add_interpolation_data(optimal_point,
2884 : optimal_comp,
2885 : optimal_elem_id,
2886 : optimal_subdomain_id,
2887 : optimal_qp,
2888 : optimal_point_perturbs,
2889 : optimal_point_phi_i_qp,
2890 : optimal_elem_type,
2891 : optimal_JxW_all_qp,
2892 : optimal_phi_i_all_qp,
2893 : optimal_qrule_order,
2894 : optimal_dxyzdxi_elem_center,
2895 : optimal_dxyzdeta_elem_center);
2896 :
2897 : // In this case we did not encounter a linearly dependent basis function, so return false
2898 0 : return false;
2899 0 : }
2900 :
2901 0 : void RBEIMConstruction::update_eim_matrices(bool set_eim_error_indicator)
2902 : {
2903 0 : LOG_SCOPE("update_eim_matrices()", "RBEIMConstruction");
2904 :
2905 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
2906 0 : unsigned int RB_size = eim_eval.get_n_basis_functions();
2907 :
2908 0 : libmesh_assert_msg(RB_size >= 1, "Must have at least 1 basis function.");
2909 :
2910 0 : if (set_eim_error_indicator)
2911 : {
2912 : // Here we have RB_size EIM basis functions, and RB_size+1 interpolation points,
2913 : // since we should have added one extra interpolation point for the EIM error
2914 : // indicator. As a result, we use RB_size as the index to access the (RB_size+1)^th
2915 : // interpolation point in the calls to eim_eval.get_interpolation_points_*.
2916 0 : DenseVector<Number> extra_point_row(RB_size);
2917 :
2918 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
2919 : {
2920 : // update the EIM interpolation matrix
2921 0 : for (unsigned int j=0; j<RB_size; j++)
2922 : {
2923 : // Evaluate the basis functions at the new interpolation point in order
2924 : // to update the interpolation matrix
2925 : Number value =
2926 0 : eim_eval.get_eim_basis_function_side_value(j,
2927 : eim_eval.get_interpolation_points_elem_id(RB_size),
2928 : eim_eval.get_interpolation_points_side_index(RB_size),
2929 : eim_eval.get_interpolation_points_comp(RB_size),
2930 : eim_eval.get_interpolation_points_qp(RB_size));
2931 0 : extra_point_row(j) = value;
2932 : }
2933 : }
2934 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
2935 : {
2936 : // update the EIM interpolation matrix
2937 0 : for (unsigned int j=0; j<RB_size; j++)
2938 : {
2939 : // Evaluate the basis functions at the new interpolation point in order
2940 : // to update the interpolation matrix
2941 : Number value =
2942 0 : eim_eval.get_eim_basis_function_node_value(j,
2943 : eim_eval.get_interpolation_points_node_id(RB_size),
2944 : eim_eval.get_interpolation_points_comp(RB_size));
2945 0 : extra_point_row(j) = value;
2946 : }
2947 : }
2948 : else
2949 : {
2950 : // update the EIM interpolation matrix
2951 0 : for (unsigned int j=0; j<RB_size; j++)
2952 : {
2953 : // Evaluate the basis functions at the new interpolation point in order
2954 : // to update the interpolation matrix
2955 : Number value =
2956 0 : eim_eval.get_eim_basis_function_value(j,
2957 : eim_eval.get_interpolation_points_elem_id(RB_size),
2958 : eim_eval.get_interpolation_points_comp(RB_size),
2959 : eim_eval.get_interpolation_points_qp(RB_size));
2960 0 : extra_point_row(j) = value;
2961 : }
2962 : }
2963 :
2964 0 : eim_eval.set_error_indicator_interpolation_row(extra_point_row);
2965 0 : return;
2966 : }
2967 :
2968 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
2969 : {
2970 : // update the matrix that is used to evaluate L2 projections
2971 : // into the EIM approximation space
2972 0 : for (unsigned int i=(RB_size-1); i<RB_size; i++)
2973 : {
2974 0 : for (unsigned int j=0; j<RB_size; j++)
2975 : {
2976 0 : Number value = side_inner_product(eim_eval.get_side_basis_function(j),
2977 0 : eim_eval.get_side_basis_function(i),
2978 : /*apply_comp_scaling*/ false);
2979 :
2980 0 : _eim_projection_matrix(i,j) = value;
2981 0 : if (i!=j)
2982 : {
2983 : // The inner product matrix is assumed to be hermitian
2984 0 : _eim_projection_matrix(j,i) = libmesh_conj(value);
2985 : }
2986 : }
2987 : }
2988 :
2989 : // update the EIM interpolation matrix
2990 0 : for (unsigned int j=0; j<RB_size; j++)
2991 : {
2992 : // Evaluate the basis functions at the new interpolation point in order
2993 : // to update the interpolation matrix
2994 : Number value =
2995 0 : eim_eval.get_eim_basis_function_side_value(j,
2996 : eim_eval.get_interpolation_points_elem_id(RB_size-1),
2997 : eim_eval.get_interpolation_points_side_index(RB_size-1),
2998 : eim_eval.get_interpolation_points_comp(RB_size-1),
2999 : eim_eval.get_interpolation_points_qp(RB_size-1));
3000 0 : eim_eval.set_interpolation_matrix_entry(RB_size-1, j, value);
3001 : }
3002 : }
3003 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
3004 : {
3005 : // update the matrix that is used to evaluate L2 projections
3006 : // into the EIM approximation space
3007 0 : for (unsigned int i=(RB_size-1); i<RB_size; i++)
3008 : {
3009 0 : for (unsigned int j=0; j<RB_size; j++)
3010 : {
3011 0 : Number value = node_inner_product(eim_eval.get_node_basis_function(j),
3012 0 : eim_eval.get_node_basis_function(i),
3013 : /*apply_comp_scaling*/ false);
3014 :
3015 0 : _eim_projection_matrix(i,j) = value;
3016 0 : if (i!=j)
3017 : {
3018 : // The inner product matrix is assumed to be hermitian
3019 0 : _eim_projection_matrix(j,i) = libmesh_conj(value);
3020 : }
3021 : }
3022 : }
3023 :
3024 : // update the EIM interpolation matrix
3025 0 : for (unsigned int j=0; j<RB_size; j++)
3026 : {
3027 : // Evaluate the basis functions at the new interpolation point in order
3028 : // to update the interpolation matrix
3029 : Number value =
3030 0 : eim_eval.get_eim_basis_function_node_value(j,
3031 : eim_eval.get_interpolation_points_node_id(RB_size-1),
3032 : eim_eval.get_interpolation_points_comp(RB_size-1));
3033 0 : eim_eval.set_interpolation_matrix_entry(RB_size-1, j, value);
3034 : }
3035 : }
3036 : else
3037 : {
3038 : // update the matrix that is used to evaluate L2 projections
3039 : // into the EIM approximation space
3040 0 : for (unsigned int i=(RB_size-1); i<RB_size; i++)
3041 : {
3042 0 : for (unsigned int j=0; j<RB_size; j++)
3043 : {
3044 0 : Number value = inner_product(eim_eval.get_basis_function(j),
3045 0 : eim_eval.get_basis_function(i),
3046 : /*apply_comp_scaling*/ false);
3047 :
3048 0 : _eim_projection_matrix(i,j) = value;
3049 0 : if (i!=j)
3050 : {
3051 : // The inner product matrix is assumed to be hermitian
3052 0 : _eim_projection_matrix(j,i) = libmesh_conj(value);
3053 : }
3054 : }
3055 : }
3056 :
3057 : // update the EIM interpolation matrix
3058 0 : for (unsigned int j=0; j<RB_size; j++)
3059 : {
3060 : // Evaluate the basis functions at the new interpolation point in order
3061 : // to update the interpolation matrix
3062 : Number value =
3063 0 : eim_eval.get_eim_basis_function_value(j,
3064 : eim_eval.get_interpolation_points_elem_id(RB_size-1),
3065 : eim_eval.get_interpolation_points_comp(RB_size-1),
3066 : eim_eval.get_interpolation_points_qp(RB_size-1));
3067 0 : eim_eval.set_interpolation_matrix_entry(RB_size-1, j, value);
3068 : }
3069 : }
3070 : }
3071 :
3072 0 : void RBEIMConstruction::scale_node_parametrized_function(NodeDataMap & local_pf,
3073 : Number scaling_factor)
3074 : {
3075 0 : for (auto & pr : local_pf)
3076 : {
3077 0 : auto & values = pr.second;
3078 0 : for ( auto & value : values)
3079 0 : value *= scaling_factor;
3080 : }
3081 0 : }
3082 :
3083 0 : unsigned int RBEIMConstruction::get_random_int_0_to_n(unsigned int n)
3084 : {
3085 : // std::random_device seed;
3086 : // std::mt19937 gen{seed()};
3087 : // We do not use a random seed here, since we generally prefer our results
3088 : // to reproducible, rather than fully random. If desired we could provide an
3089 : // option to use the random seed approach (commented out above).
3090 0 : std::default_random_engine gen;
3091 0 : std::uniform_int_distribution<> dist{0, static_cast<int>(n)};
3092 0 : return dist(gen);
3093 : }
3094 :
3095 0 : EimPointData RBEIMConstruction::get_random_point(const QpDataMap & v)
3096 : {
3097 : EimPointData eim_point_data;
3098 :
3099 : // If we have more than one process, then we need to do a parallel union
3100 : // of v to make sure that we have data from all processors. Our approach
3101 : // here is to set v_ptr to either v or global_v, depending on whether we
3102 : // are in parallel or serial. The purpose of this approach is to avoid
3103 : // making a copy of v in the case that this is a serial job.
3104 0 : QpDataMap const * v_ptr = nullptr;
3105 0 : QpDataMap global_v;
3106 0 : if (comm().size() > 1)
3107 : {
3108 0 : global_v = v;
3109 :
3110 : // We only use global_v on proc 0, so we set the second argument of
3111 : // set_union() to zero here to indicate that we only need the result
3112 : // on proc 0.
3113 0 : comm().set_union(global_v, 0);
3114 0 : v_ptr = &global_v;
3115 : }
3116 : else
3117 : {
3118 0 : v_ptr = &v;
3119 : }
3120 :
3121 0 : bool error_finding_new_element = false;
3122 0 : if (comm().rank() == 0)
3123 : {
3124 0 : const VectorizedEvalInput & vec_eval_input = get_rb_eim_evaluation().get_vec_eval_input();
3125 :
3126 : {
3127 0 : std::set<dof_id_type> previous_elem_ids(vec_eval_input.elem_ids.begin(), vec_eval_input.elem_ids.end());
3128 :
3129 : // We ensure that we select a point that has not been selected previously
3130 : // by setting up new_elem_ids to contain only elements that are not in
3131 : // previous_elem_ids, and then selecting the elem_id at random from new_elem_ids.
3132 : // We give an error if there are no elements in new_elem_ids. This is potentially
3133 : // an overzealous assertion since we could pick an element that has already
3134 : // been selected as long as we pick a (comp_index, qp_index) that has not already
3135 : // been selected for that element.
3136 : //
3137 : // However, in general we do not expect all elements to be selected in the EIM
3138 : // training, so it is reasonable to use the simple assertion below. Moreover, by
3139 : // ensuring that we choose a new element we should typically ensure that the
3140 : // randomly selected point has some separation from the previous EIM points, which
3141 : // is typically desirable if we want EIM evaluations that are independent from
3142 : // the EIM points (e.g. for EIM error indicator purposes).
3143 0 : std::set<dof_id_type> new_elem_ids;
3144 0 : for (const auto & v_pair : *v_ptr)
3145 0 : if (previous_elem_ids.count(v_pair.first) == 0)
3146 0 : new_elem_ids.insert(v_pair.first);
3147 :
3148 : // If new_elem_ids is empty then we set error_finding_new_element to true.
3149 : // We then broadcast the value of error_finding_new_element to all processors
3150 : // below in order to ensure that all processors agree on whether or not
3151 : // there was an error.
3152 0 : error_finding_new_element = (new_elem_ids.empty());
3153 :
3154 0 : if (!error_finding_new_element)
3155 : {
3156 0 : unsigned int random_elem_idx = get_random_int_0_to_n(new_elem_ids.size()-1);
3157 :
3158 0 : auto item = new_elem_ids.begin();
3159 0 : std::advance(item, random_elem_idx);
3160 0 : eim_point_data.elem_id = *item;
3161 : }
3162 : }
3163 :
3164 0 : if (!error_finding_new_element)
3165 : {
3166 : {
3167 0 : const auto & vars_and_qps = libmesh_map_find(*v_ptr,eim_point_data.elem_id);
3168 0 : eim_point_data.comp_index = get_random_int_0_to_n(vars_and_qps.size()-1);
3169 : }
3170 :
3171 : {
3172 0 : const auto & qps = libmesh_map_find(*v_ptr,eim_point_data.elem_id)[eim_point_data.comp_index];
3173 0 : eim_point_data.qp_index = get_random_int_0_to_n(qps.size()-1);
3174 : }
3175 : }
3176 : }
3177 :
3178 0 : comm().broadcast(error_finding_new_element);
3179 0 : libmesh_error_msg_if(error_finding_new_element, "Could not find new element in get_random_point()");
3180 :
3181 : // Broadcast the values computed above from rank 0
3182 0 : comm().broadcast(eim_point_data.elem_id);
3183 0 : comm().broadcast(eim_point_data.comp_index);
3184 0 : comm().broadcast(eim_point_data.qp_index);
3185 :
3186 0 : return eim_point_data;
3187 : }
3188 :
3189 0 : EimPointData RBEIMConstruction::get_random_point(const SideQpDataMap & v)
3190 : {
3191 : EimPointData eim_point_data;
3192 :
3193 : // If we have more than one process, then we need to do a parallel union
3194 : // of v to make sure that we have data from all processors. Our approach
3195 : // here is to set v_ptr to either v or global_v, depending on whether we
3196 : // are in parallel or serial. The purpose of this approach is to avoid
3197 : // making a copy of v in the case that this is a serial job.
3198 0 : SideQpDataMap const * v_ptr = nullptr;
3199 0 : SideQpDataMap global_v;
3200 0 : if (comm().size() > 1)
3201 : {
3202 0 : global_v = v;
3203 :
3204 : // We only use global_v on proc 0, so we set the second argument of
3205 : // set_union() to zero here to indicate that we only need the result
3206 : // on proc 0.
3207 0 : comm().set_union(global_v, 0);
3208 0 : v_ptr = &global_v;
3209 : }
3210 : else
3211 : {
3212 0 : v_ptr = &v;
3213 : }
3214 :
3215 0 : bool error_finding_new_element_and_side = false;
3216 0 : if (comm().rank() == 0)
3217 : {
3218 0 : const VectorizedEvalInput & vec_eval_input = get_rb_eim_evaluation().get_vec_eval_input();
3219 :
3220 0 : std::pair<dof_id_type,unsigned int> elem_and_side;
3221 : {
3222 0 : std::set<std::pair<dof_id_type,unsigned int>> previous_elem_and_side_ids;
3223 0 : for (const auto idx : index_range(vec_eval_input.elem_ids))
3224 : {
3225 0 : previous_elem_and_side_ids.insert(
3226 0 : std::make_pair(vec_eval_input.elem_ids[idx],
3227 0 : vec_eval_input.side_indices[idx]));
3228 : }
3229 :
3230 : // See discussion above in the QpDataMap case for the justification
3231 : // of how we set up new_elem_and_side_ids below.
3232 0 : std::set<std::pair<dof_id_type,unsigned int>> new_elem_and_side_ids;
3233 0 : for (const auto & v_pair : *v_ptr)
3234 0 : if (previous_elem_and_side_ids.count(v_pair.first) == 0)
3235 0 : new_elem_and_side_ids.insert(v_pair.first);
3236 :
3237 : // If new_elem_and_side_ids is empty then we set error_finding_new_element_and_side
3238 : // to true. We then broadcast the value of error_finding_new_element_and_side to all
3239 : // processors below in order to ensure that all processors agree on whether
3240 : // or not there was an error.
3241 0 : error_finding_new_element_and_side = (new_elem_and_side_ids.empty());
3242 :
3243 0 : if (!error_finding_new_element_and_side)
3244 : {
3245 0 : unsigned int random_elem_and_side_idx = get_random_int_0_to_n(new_elem_and_side_ids.size()-1);
3246 :
3247 0 : auto item = new_elem_and_side_ids.begin();
3248 0 : std::advance(item, random_elem_and_side_idx);
3249 0 : elem_and_side = *item;
3250 0 : eim_point_data.elem_id = elem_and_side.first;
3251 0 : eim_point_data.side_index = elem_and_side.second;
3252 : }
3253 : }
3254 :
3255 0 : if (!error_finding_new_element_and_side)
3256 : {
3257 : {
3258 0 : const auto & vars_and_qps = libmesh_map_find(*v_ptr,elem_and_side);
3259 0 : eim_point_data.comp_index = get_random_int_0_to_n(vars_and_qps.size()-1);
3260 : }
3261 :
3262 : {
3263 0 : const auto & qps = libmesh_map_find(*v_ptr,elem_and_side)[eim_point_data.comp_index];
3264 0 : eim_point_data.qp_index = get_random_int_0_to_n(qps.size()-1);
3265 : }
3266 : }
3267 : }
3268 :
3269 0 : comm().broadcast(error_finding_new_element_and_side);
3270 0 : libmesh_error_msg_if(error_finding_new_element_and_side, "Could not find new (element,side) in get_random_point()");
3271 :
3272 : // Broadcast the values computed above from rank 0
3273 0 : comm().broadcast(eim_point_data.elem_id);
3274 0 : comm().broadcast(eim_point_data.side_index);
3275 0 : comm().broadcast(eim_point_data.comp_index);
3276 0 : comm().broadcast(eim_point_data.qp_index);
3277 :
3278 0 : return eim_point_data;
3279 : }
3280 :
3281 0 : EimPointData RBEIMConstruction::get_random_point(const NodeDataMap & v)
3282 : {
3283 : EimPointData eim_point_data;
3284 :
3285 : // If we have more than one process, then we need to do a parallel union
3286 : // of v to make sure that we have data from all processors. Our approach
3287 : // here is to set v_ptr to either v or global_v, depending on whether we
3288 : // are in parallel or serial. The purpose of this approach is to avoid
3289 : // making a copy of v in the case that this is a serial job.
3290 0 : NodeDataMap const * v_ptr = nullptr;
3291 0 : NodeDataMap global_v;
3292 0 : if (comm().size() > 1)
3293 : {
3294 0 : global_v = v;
3295 :
3296 : // We only use global_v on proc 0, so we set the second argument of
3297 : // set_union() to zero here to indicate that we only need the result
3298 : // on proc 0.
3299 0 : comm().set_union(global_v, 0);
3300 0 : v_ptr = &global_v;
3301 : }
3302 : else
3303 : {
3304 0 : v_ptr = &v;
3305 : }
3306 :
3307 0 : bool error_finding_new_node = false;
3308 0 : if (comm().rank() == 0)
3309 : {
3310 0 : const VectorizedEvalInput & vec_eval_input = get_rb_eim_evaluation().get_vec_eval_input();
3311 :
3312 : {
3313 0 : std::set<dof_id_type> previous_node_ids(vec_eval_input.node_ids.begin(), vec_eval_input.node_ids.end());
3314 :
3315 : // See discussion above in the QpDataMap case for the justification
3316 : // of how we set up new_node_ids below.
3317 0 : std::set<dof_id_type> new_node_ids;
3318 0 : for (const auto & v_pair : *v_ptr)
3319 0 : if (previous_node_ids.count(v_pair.first) == 0)
3320 0 : new_node_ids.insert(v_pair.first);
3321 :
3322 : // If new_node_ids is empty then we set error_finding_new_node
3323 : // to true. We then broadcast the value of error_finding_new_node to all
3324 : // processors below in order to ensure that all processors agree on whether
3325 : // or not there was an error.
3326 0 : error_finding_new_node = (new_node_ids.empty());
3327 :
3328 0 : if (!error_finding_new_node)
3329 : {
3330 0 : unsigned int random_node_idx = get_random_int_0_to_n(new_node_ids.size()-1);
3331 :
3332 0 : auto item = new_node_ids.begin();
3333 0 : std::advance(item, random_node_idx);
3334 0 : eim_point_data.node_id = *item;
3335 : }
3336 : }
3337 :
3338 0 : if (!error_finding_new_node)
3339 : {
3340 0 : const auto & vars = libmesh_map_find(*v_ptr,eim_point_data.node_id);
3341 0 : eim_point_data.comp_index = get_random_int_0_to_n(vars.size()-1);
3342 : }
3343 : }
3344 :
3345 0 : comm().broadcast(error_finding_new_node);
3346 0 : libmesh_error_msg_if(error_finding_new_node, "Could not find new node in get_random_point()");
3347 :
3348 : // Broadcast the values computed above from rank 0
3349 0 : comm().broadcast(eim_point_data.node_id);
3350 0 : comm().broadcast(eim_point_data.comp_index);
3351 :
3352 0 : return eim_point_data;
3353 : }
3354 :
3355 0 : EimPointData RBEIMConstruction::get_random_point_from_training_sample()
3356 : {
3357 0 : RBEIMEvaluation & eim_eval = get_rb_eim_evaluation();
3358 :
3359 0 : if (eim_eval.get_parametrized_function().on_mesh_sides())
3360 0 : return get_random_point(_local_side_parametrized_functions_for_training[0]);
3361 0 : else if (eim_eval.get_parametrized_function().on_mesh_nodes())
3362 0 : return get_random_point(_local_node_parametrized_functions_for_training[0]);
3363 : else
3364 0 : return get_random_point(_local_parametrized_functions_for_training[0]);
3365 : }
3366 :
3367 : } // namespace libMesh
|
