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