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