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