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 
44 #include "libmesh/rb_eim_construction.h"
45 #include "libmesh/rb_eim_evaluation.h"
46 
47 namespace libMesh
48 {
49 
51  const std::string & name_in,
52  const unsigned int number_in)
53  : Parent(es, name_in, number_in),
54  best_fit_type_flag(PROJECTION_BEST_FIT),
55  _parametrized_functions_in_training_set_initialized(false),
56  _point_locator_tol(TOLERANCE)
57 {
58  _explicit_system_name = name_in + "_explicit_sys";
59 
60  // We cannot do rb_solve with an empty
61  // "rb space" with EIM
63 
64  // Indicate that we need to compute the RB
65  // inner product matrix in this case
67 
68  // Indicate that we need the training set
69  // for the Greedy to be the same on all
70  // processors
71  serial_training_set = true;
72 
73  // attach empty RBAssemblyExpansion object
75 
76  // We set implicit_neighbor_dofs = false. This is important when we use
77  // DISCONTINUOUS basis functions, since by default libMesh sets
78  // implicit_neighbor_dofs = true for "all discontinuous" systems, which
79  // results in a larger sparsity: dofs on neighboring elements are added
80  // to the sparsity pattern since this is typically required for DG or FV
81  // discretizations. Since we're only doing L2 projects here, we do not
82  // need extra dofs in the sparsity pattern, so we set implicit neighbor
83  // dofs to false.
85 }
86 
88 {
89  this->clear();
90 }
91 
93 {
94  Parent::clear();
95 
96  // clear the mesh function
97  _mesh_function.reset();
98 
99  // clear the eim assembly vector
100  for (std::size_t i=0; i<_rb_eim_assembly_objects.size(); i++)
101  delete _rb_eim_assembly_objects[i];
102  _rb_eim_assembly_objects.clear();
103 
104  // clear the parametrized functions from the training set
105  for (std::size_t i=0; i<_parametrized_functions_in_training_set.size(); i++)
106  {
108  {
112  }
113  }
115 
116  for (std::size_t i=0; i<_matrix_times_bfs.size(); i++)
117  {
118  if (_matrix_times_bfs[i])
119  {
120  _matrix_times_bfs[i]->clear();
121  delete _matrix_times_bfs[i];
123  }
124  }
125 }
126 
127 void RBEIMConstruction::process_parameters_file (const std::string & parameters_filename)
128 {
129  Parent::process_parameters_file(parameters_filename);
130 
131  GetPot infile(parameters_filename);
132 
133  std::string best_fit_type_string = infile("best_fit_type","projection");
134  set_best_fit_type_flag(best_fit_type_string);
135 }
136 
137 void RBEIMConstruction::set_best_fit_type_flag (const std::string & best_fit_type_string)
138 {
139  if (best_fit_type_string == "projection")
140  {
142  }
143  else
144  if (best_fit_type_string == "eim")
145  {
147  }
148  else
149  libmesh_error_msg("Error: invalid best_fit_type in input file");
150 }
151 
153 {
155 
156  // Print out setup info
157  libMesh::out << std::endl << "RBEIMConstruction parameters:" << std::endl;
159  {
160  libMesh::out << "best fit type: projection" << std::endl;
161  }
162  else
164  {
165  libMesh::out << "best fit type: eim" << std::endl;
166  }
167  libMesh::out << std::endl;
168 }
169 
171 {
172  // Add the ExplicitSystem that we use to store the EIM basis functions
174 
177 
179 }
180 
181 void RBEIMConstruction::initialize_rb_construction(bool skip_matrix_assembly,
182  bool skip_vector_assembly)
183 {
184  Parent::initialize_rb_construction(skip_matrix_assembly, skip_vector_assembly);
185 
186  // initialize a serial vector that we will use for MeshFunction evaluations
189  get_explicit_system().get_dof_map().get_send_list(), false,
190  GHOSTED);
191 
192  // Initialize the MeshFunction for interpolating the
193  // solution vector at quadrature points
194  std::vector<unsigned int> vars;
199  vars));
200  _mesh_function->init();
201 
202  // inner_product_solver performs solves with the same matrix every time
203  // hence we can set reuse_preconditioner(true).
204  inner_product_solver->reuse_preconditioner(true);
205 
207 }
208 
209 Real RBEIMConstruction::train_reduced_basis(const bool resize_rb_eval_data)
210 {
211  // precompute all the parametrized functions that we'll use in the greedy
213 
214  return Parent::train_reduced_basis(resize_rb_eval_data);
215 }
216 
218  Point p)
219 {
220  // Set default values to be an empty vector so that we can use it
221  // below
222  // to check if this processor returned valid data.
223  DenseVector<Number> default_values;
224  _mesh_function->enable_out_of_mesh_mode(default_values);
225 
226  _mesh_function->set_point_locator_tolerance( get_point_locator_tol() );
227 
228  DenseVector<Number> values;
229  (*_mesh_function)(p,
230  /*time*/ 0.,
231  values);
232 
233  // We evaluated the mesh function, but it will only return a valid set of values on one processor
234  // (values will be empty on all other processors) so we need to broadcast those valid values
235  // to all processors.
236  Number value = 0;
237  unsigned int root_id=0;
238  unsigned int check_for_valid_value = 0;
239  if (values.size() != 0)
240  {
241  root_id = this->processor_id();
242  value = values(var_number);
243  check_for_valid_value = 1;
244  }
245 
246  // If this sum is zero, then we didn't enter the if block above on any processor. In that
247  // case we should throw an error.
248  this->comm().sum(check_for_valid_value);
249  if (check_for_valid_value == 0)
250  {
251  libmesh_error_msg("MeshFunction evaluation failed on all processors");
252  }
253 
254  // root_id may be non-zero on more than one processor due to ghost elements
255  // so use this->comm().max to get just one proc id
256  this->comm().max(root_id);
257 
258  // Then broadcast the result
259  this->comm().broadcast(value, root_id);
260 
261  return value;
262 }
263 
265 {
266  _point_locator_tol = point_locator_tol;
267 }
268 
270 {
271  return _point_locator_tol;
272 }
273 
275 {
276  _rb_eim_assembly_objects.clear();
277  for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
278  {
279  _rb_eim_assembly_objects.push_back( build_eim_assembly(i).release() );
280  }
281 }
282 
284 {
286 }
287 
289 {
290  LOG_SCOPE("load_basis_function()", "RBEIMConstruction");
291 
292  libmesh_assert_less (i, get_rb_evaluation().get_n_basis_functions());
293 
295 
297 }
298 
300 {
301  LOG_SCOPE("load_rb_solution()", "RBEIMConstruction");
302 
303  solution->zero();
304 
305  if (get_rb_evaluation().RB_solution.size() > get_rb_evaluation().get_n_basis_functions())
306  libmesh_error_msg("ERROR: System contains " << get_rb_evaluation().get_n_basis_functions() << " basis functions." \
307  << " RB_solution vector constains " << get_rb_evaluation().RB_solution.size() << " entries." \
308  << " RB_solution in RBConstruction::load_rb_solution is too long!");
309 
310  for (unsigned int i=0; i<get_rb_evaluation().RB_solution.size(); i++)
311  get_explicit_system().solution->add(get_rb_evaluation().RB_solution(i),
312  get_rb_evaluation().get_basis_function(i));
313 
315 }
316 
317 std::vector<ElemAssembly *> RBEIMConstruction::get_eim_assembly_objects()
318 {
320 }
321 
323 {
324  LOG_SCOPE("enrich_RB_space()", "RBEIMConstruction");
325 
326  // put solution in _ghosted_meshfunction_vector so we can access it from the mesh function
327  // this allows us to compute EIM_rhs appropriately
329  get_explicit_system().get_dof_map().get_send_list());
330 
331  RBEIMEvaluation & eim_eval = cast_ref<RBEIMEvaluation &>(get_rb_evaluation());
332 
333  // If we have at least one basis function we need to use
334  // rb_solve to find the EIM interpolation error, otherwise just use solution as is
336  {
337  // get the right-hand side vector for the EIM approximation
338  // by sampling the parametrized function (stored in solution)
339  // at the interpolation points
340  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
341  DenseVector<Number> EIM_rhs(RB_size);
342  for (unsigned int i=0; i<RB_size; i++)
343  {
344  EIM_rhs(i) = evaluate_mesh_function( eim_eval.interpolation_points_var[i],
345  eim_eval.interpolation_points[i] );
346  }
347 
348  eim_eval.set_parameters( get_parameters() );
349  eim_eval.rb_solve(EIM_rhs);
350 
351  // Load the "EIM residual" into solution by subtracting
352  // the EIM approximation
353  for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
354  {
356  }
357  }
358 
359  // need to update since context uses current_local_solution
361 
362  // Find the quadrature point at which solution (which now stores
363  // the "EIM residual") has maximum absolute value
364  // by looping over the mesh
365  Point optimal_point;
366  Number optimal_value = 0.;
367  unsigned int optimal_var = 0;
368  dof_id_type optimal_elem_id = DofObject::invalid_id;
369 
370  // Initialize largest_abs_value to be negative so that it definitely gets updated.
371  Real largest_abs_value = -1.;
372 
373  // Compute truth representation via projection
374  MeshBase & mesh = this->get_mesh();
375 
377  DGFEMContext & explicit_context = cast_ref<DGFEMContext &>(*explicit_c);
378  init_context_with_sys(explicit_context, get_explicit_system());
379 
380  for (const auto & elem : mesh.active_local_element_ptr_range())
381  {
382  explicit_context.pre_fe_reinit(get_explicit_system(), elem);
383  explicit_context.elem_fe_reinit();
384 
385  for (unsigned int var=0; var<get_explicit_system().n_vars(); var++)
386  {
387  unsigned int n_qpoints = explicit_context.get_element_qrule().n_points();
388 
389  for (unsigned int qp=0; qp<n_qpoints; qp++)
390  {
391  Number value = explicit_context.interior_value(var, qp);
392  Real abs_value = std::abs(value);
393 
394  if (abs_value > largest_abs_value)
395  {
396  optimal_value = value;
397  largest_abs_value = abs_value;
398  optimal_var = var;
399  optimal_elem_id = elem->id();
400 
401  FEBase * elem_fe = libmesh_nullptr;
402  explicit_context.get_element_fe( var, elem_fe );
403  optimal_point = elem_fe->get_xyz()[qp];
404  }
405 
406  }
407  }
408  }
409 
410  // Find out which processor has the largest of the abs values
411  unsigned int proc_ID_index;
412  this->comm().maxloc(largest_abs_value, proc_ID_index);
413 
414  // Broadcast the optimal point from proc_ID_index
415  this->comm().broadcast(optimal_point, proc_ID_index);
416 
417  // Also broadcast the corresponding optimal_var, optimal_value, and optimal_elem_id
418  this->comm().broadcast(optimal_var, proc_ID_index);
419  this->comm().broadcast(optimal_value, proc_ID_index);
420  this->comm().broadcast(optimal_elem_id, proc_ID_index);
421 
422  // In debug mode, assert that we found an optimal_elem_id
423  libmesh_assert_not_equal_to(optimal_elem_id, DofObject::invalid_id);
424 
425  // Scale the solution
426  get_explicit_system().solution->scale(1./optimal_value);
427 
428  // Store optimal point in interpolation_points
429  eim_eval.interpolation_points.push_back(optimal_point);
430  eim_eval.interpolation_points_var.push_back(optimal_var);
431  Elem * elem_ptr = mesh.elem_ptr(optimal_elem_id);
432  eim_eval.interpolation_points_elem.push_back( elem_ptr );
433 
434  NumericVector<Number> * new_bf = NumericVector<Number>::build(this->comm()).release();
436  *new_bf = *get_explicit_system().solution;
437  get_rb_evaluation().basis_functions.push_back( new_bf );
438 
440  {
441  // In order to speed up dot products, we store the product
442  // of the basis function and the inner product matrix
443 
444  UniquePtr<NumericVector<Number>> implicit_sys_temp1 = this->solution->zero_clone();
445  UniquePtr<NumericVector<Number>> implicit_sys_temp2 = this->solution->zero_clone();
446  NumericVector<Number>* matrix_times_new_bf =
447  get_explicit_system().solution->zero_clone().release();
448 
449  // We must localize new_bf before calling get_explicit_sys_subvector
450  UniquePtr<NumericVector<Number>> localized_new_bf =
452  localized_new_bf->init(get_explicit_system().n_dofs(), false, SERIAL);
453  new_bf->localize(*localized_new_bf);
454 
455  for (unsigned int var=0; var<get_explicit_system().n_vars(); var++)
456  {
457  get_explicit_sys_subvector(*implicit_sys_temp1,
458  var,
459  *localized_new_bf);
460 
461  inner_product_matrix->vector_mult(*implicit_sys_temp2, *implicit_sys_temp1);
462 
463  set_explicit_sys_subvector(*matrix_times_new_bf,
464  var,
465  *implicit_sys_temp2);
466  }
467 
468  _matrix_times_bfs.push_back(matrix_times_new_bf);
469  }
470 }
471 
473 {
474  if (!serial_training_set)
475  libmesh_error_msg("Error: We must have serial_training_set==true in " \
476  << "RBEIMConstruction::initialize_parametrized_functions_in_training_set");
477 
478  libMesh::out << "Initializing parametrized functions in training set..." << std::endl;
479  // initialize rb_eval's parameters
481 
483  for (unsigned int i=0; i<get_n_training_samples(); i++)
484  {
486  truth_solve(-1);
487 
489 
490  libMesh::out << "Completed solve for training sample " << (i+1) << " of " << get_n_training_samples() << std::endl;
491  }
492 
494 
495  libMesh::out << "Parametrized functions in training set initialized" << std::endl << std::endl;
496 }
497 
499 {
500  // Ignore unused parameter warnings when Exodus is not available.
501  libmesh_ignore(pathname);
502 
504 
505  for (std::size_t i=0; i<_parametrized_functions_in_training_set.size(); i++)
506  {
507 #ifdef LIBMESH_HAVE_EXODUS_API
509 
510  std::stringstream pathname_i;
511  pathname_i << pathname << "_" << i << ".exo";
512 
513  std::set<std::string> system_names;
514  system_names.insert(get_explicit_system().name());
515  ExodusII_IO(get_mesh()).write_equation_systems (pathname_i.str(),
516  this->get_equation_systems(),
517  &system_names);
518  libMesh::out << "Plotted parameterized function " << i << std::endl;
519 #endif
520  }
521 }
522 
523 
524 
526 {
527  LOG_SCOPE("compute_best_fit_error()", "RBEIMConstruction");
528 
529  const unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
530 
531  // load up the parametrized function for the current parameters
532  truth_solve(-1);
533 
534  switch(best_fit_type_flag)
535  {
536  // Perform an L2 projection in order to find an approximation to solution (from truth_solve above)
537  case(PROJECTION_BEST_FIT):
538  {
539  // We have pre-stored inner_product_matrix * basis_function[i] for each i
540  // so we can just evaluate the dot product here.
541  DenseVector<Number> best_fit_rhs(RB_size);
542  for (unsigned int i=0; i<RB_size; i++)
543  {
544  best_fit_rhs(i) = get_explicit_system().solution->dot(*_matrix_times_bfs[i]);
545  }
546 
547  // Now compute the best fit by an LU solve
549  DenseMatrix<Number> RB_inner_product_matrix_N(RB_size);
550  get_rb_evaluation().RB_inner_product_matrix.get_principal_submatrix(RB_size, RB_inner_product_matrix_N);
551 
552  RB_inner_product_matrix_N.lu_solve(best_fit_rhs, get_rb_evaluation().RB_solution);
553  break;
554  }
555  // Perform EIM solve in order to find the approximation to solution
556  // (rb_solve provides the EIM basis function coefficients used below)
557  case(EIM_BEST_FIT):
558  {
559  // Turn off error estimation for this rb_solve, we use the linfty norm instead
562  get_rb_evaluation().rb_solve(RB_size);
564  break;
565  }
566  default:
567  libmesh_error_msg("Should not reach here");
568  }
569 
570  // load the error into solution
571  for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
572  get_explicit_system().solution->add(-get_rb_evaluation().RB_solution(i),
573  get_rb_evaluation().get_basis_function(i));
574 
575  Real best_fit_error = get_explicit_system().solution->linfty_norm();
576 
577  return best_fit_error;
578 }
579 
581 {
582  LOG_SCOPE("truth_solve()", "RBEIMConstruction");
583 
584  int training_parameters_found_index = -1;
586  {
587  // Check if parameters are in the training set. If so, we can just load the
588  // solution from _parametrized_functions_in_training_set
589 
590  for (unsigned int i=0; i<get_n_training_samples(); i++)
591  {
593  {
594  training_parameters_found_index = i;
595  break;
596  }
597  }
598  }
599 
600  // If the parameters are in the training set, just copy the solution vector
601  if (training_parameters_found_index >= 0)
602  {
604  *_parametrized_functions_in_training_set[training_parameters_found_index];
605  get_explicit_system().update(); // put the solution into current_local_solution as well
606  }
607  // Otherwise, we have to compute the projection
608  else
609  {
610  if (this->n_vars() != 1)
611  {
612  libmesh_error_msg("The system that we use to perform EIM L2 solves should have one variable");
613  }
614 
615  RBEIMEvaluation & eim_eval = cast_ref<RBEIMEvaluation &>(get_rb_evaluation());
616  eim_eval.set_parameters( get_parameters() );
617 
618  // Compute truth representation via L2 projection
619  const MeshBase & mesh = this->get_mesh();
620 
621  UniquePtr<DGFEMContext> c(new DGFEMContext( *this ));
622  DGFEMContext & context = cast_ref<DGFEMContext &>(*c);
623  init_context_with_sys(context, *this);
624 
625  // First cache all the element data
626  std::vector<std::vector<std::vector<Number>>> parametrized_fn_vals(mesh.n_elem());
627  std::vector<std::vector<Real>> JxW_values(mesh.n_elem());
628  std::vector<std::vector<std::vector<Real>>> phi_values(mesh.n_elem());
629 
630  for (const auto & elem : mesh.active_local_element_ptr_range())
631  {
632  dof_id_type elem_id = elem->id();
633 
634  context.pre_fe_reinit(*this, elem);
635  context.elem_fe_reinit();
636 
637  FEBase * elem_fe = libmesh_nullptr;
638  context.get_element_fe( 0, elem_fe );
639  unsigned int n_qpoints = context.get_element_qrule().n_points();
640  const std::vector<Real> & JxW = elem_fe->get_JxW();
641  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
642  const std::vector<Point> & xyz = elem_fe->get_xyz();
643 
644  // Loop over qp before var because parametrized functions often use
645  // some caching based on qp.
646  parametrized_fn_vals[elem_id].resize(n_qpoints);
647  JxW_values[elem_id].resize(n_qpoints);
648  phi_values[elem_id].resize(n_qpoints);
649  for (unsigned int qp=0; qp<n_qpoints; qp++)
650  {
651  JxW_values[elem_id][qp] = JxW[qp];
652 
653  unsigned int n_var_dofs = cast_int<unsigned int>(context.get_dof_indices().size());
654  phi_values[elem_id][qp].resize(n_var_dofs);
655  for (unsigned int i=0; i != n_var_dofs; i++)
656  {
657  phi_values[elem_id][qp][i] = phi[i][qp];
658  }
659 
660  parametrized_fn_vals[elem_id][qp].resize(get_explicit_system().n_vars());
661  for (unsigned int var=0; var<get_explicit_system().n_vars(); var++)
662  {
663  Number eval_result = eim_eval.evaluate_parametrized_function(var, xyz[qp], *elem);
664  parametrized_fn_vals[elem_id][qp][var] = eval_result;
665  }
666  }
667  }
668 
669  // We do a distinct solve for each variable in the ExplicitSystem
670  for (unsigned int var=0; var<get_explicit_system().n_vars(); var++)
671  {
672  rhs->zero();
673 
674  for (const auto & elem : mesh.active_local_element_ptr_range())
675  {
676  dof_id_type elem_id = elem->id();
677 
678  context.pre_fe_reinit(*this, elem);
679  //context.elem_fe_reinit(); <--- skip this because we cached all the FE data
680 
681  // Loop over qp before var because parametrized functions often use
682  // some caching based on qp.
683  for (std::size_t qp=0; qp<JxW_values[elem_id].size(); qp++)
684  {
685  unsigned int n_var_dofs = phi_values[elem_id][qp].size();
686 
687  Number eval_result = parametrized_fn_vals[elem_id][qp][var];
688  for (unsigned int i=0; i != n_var_dofs; i++)
689  {
690  context.get_elem_residual()(i) +=
691  JxW_values[elem_id][qp] * eval_result * phi_values[elem_id][qp][i];
692  }
693  }
694 
695  // Apply constraints, e.g. periodic constraints
697 
698  // Add element vector to global vector
699  rhs->add_vector(context.get_elem_residual(), context.get_dof_indices() );
700  }
701 
702  // Solve to find the best fit, then solution stores the truth representation
703  // of the function to be approximated
705 
706  if (assert_convergence)
708 
709  // Now copy the solution to the explicit system's solution.
711  }
713  }
714 
715  if (plot_solution > 0)
716  {
717 #ifdef LIBMESH_HAVE_EXODUS_API
719  this->get_equation_systems());
720 #endif
721  }
722 
723  return 0.;
724 }
725 
727 {
728  // default implementation of init_context
729  // for compute_best_fit
730  for (unsigned int var=0; var<sys.n_vars(); var++)
731  {
732  FEBase * elem_fe = libmesh_nullptr;
733  c.get_element_fe( var, elem_fe );
734  elem_fe->get_JxW();
735  elem_fe->get_phi();
736  elem_fe->get_xyz();
737  }
738 }
739 
741 {
742  LOG_SCOPE("update_RB_system_matrices()", "RBEIMConstruction");
743 
744  // First, update the inner product matrix
745  {
746  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
747 
748  UniquePtr<NumericVector<Number>> explicit_sys_temp =
749  get_explicit_system().solution->zero_clone();
750 
751  UniquePtr<NumericVector<Number>> temp1 = this->solution->zero_clone();
752  UniquePtr<NumericVector<Number>> temp2 = this->solution->zero_clone();
753 
754  for (unsigned int i=(RB_size-1); i<RB_size; i++)
755  {
756  for (unsigned int j=0; j<RB_size; j++)
757  {
758  // We must localize get_rb_evaluation().get_basis_function(j) before calling
759  // get_explicit_sys_subvector
760  UniquePtr<NumericVector<Number>> localized_basis_function =
762  localized_basis_function->init(get_explicit_system().n_dofs(), false, SERIAL);
763  get_rb_evaluation().get_basis_function(j).localize(*localized_basis_function);
764 
765  // Compute reduced inner_product_matrix via a series of matvecs
766  for (unsigned int var=0; var<get_explicit_system().n_vars(); var++)
767  {
768  get_explicit_sys_subvector(*temp1, var, *localized_basis_function);
769  inner_product_matrix->vector_mult(*temp2, *temp1);
770  set_explicit_sys_subvector(*explicit_sys_temp, var, *temp2);
771  }
772 
773  Number value = explicit_sys_temp->dot( get_rb_evaluation().get_basis_function(i) );
775  if (i!=j)
776  {
777  // The inner product matrix is assumed
778  // to be hermitian
780  }
781  }
782  }
783  }
784 
785  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
786 
787  RBEIMEvaluation & eim_eval = cast_ref<RBEIMEvaluation &>(get_rb_evaluation());
788 
789  // update the EIM interpolation matrix
790  for (unsigned int j=0; j<RB_size; j++)
791  {
792  // Sample the basis functions at the
793  // new interpolation point
795  get_explicit_system().get_dof_map().get_send_list());
796 
797  eim_eval.interpolation_matrix(RB_size-1,j) =
799  eim_eval.interpolation_points[RB_size-1] );
800  }
801 }
802 
804 {
805  Real best_fit_error = compute_best_fit_error();
806  return best_fit_error;
807 }
808 
810 {
811  libMesh::out << "Updating RB matrices" << std::endl;
813 }
814 
816  unsigned int var,
817  NumericVector<Number> & source)
818 {
819  LOG_SCOPE("set_explicit_sys_subvector()", "RBEIMConstruction");
820 
821  // For convenience we localize the source vector first to make it easier to
822  // copy over (no need to do distinct send/receives).
823  UniquePtr<NumericVector<Number>> localized_source =
825  localized_source->init(this->n_dofs(), false, SERIAL);
826  source.localize(*localized_source);
827 
828  for (std::size_t i=0; i<_dof_map_between_systems[var].size(); i++)
829  {
830  dof_id_type implicit_sys_dof_index = i;
831  dof_id_type explicit_sys_dof_index = _dof_map_between_systems[var][i];
832 
833  if ((dest.first_local_index() <= explicit_sys_dof_index) &&
834  (explicit_sys_dof_index < dest.last_local_index()))
835  dest.set(explicit_sys_dof_index,
836  (*localized_source)(implicit_sys_dof_index));
837  }
838 
839  dest.close();
840 }
841 
843  unsigned int var,
844  NumericVector<Number> & localized_source)
845 {
846  LOG_SCOPE("get_explicit_sys_subvector()", "RBEIMConstruction");
847 
848  for (std::size_t i=0; i<_dof_map_between_systems[var].size(); i++)
849  {
850  dof_id_type implicit_sys_dof_index = i;
851  dof_id_type explicit_sys_dof_index = _dof_map_between_systems[var][i];
852 
853  if ((dest.first_local_index() <= implicit_sys_dof_index) &&
854  (implicit_sys_dof_index < dest.last_local_index()))
855  dest.set(implicit_sys_dof_index,
856  localized_source(explicit_sys_dof_index));
857  }
858 
859  dest.close();
860 }
861 
863 {
864  LOG_SCOPE("init_dof_map_between_systems()", "RBEIMConstruction");
865 
866  unsigned int n_vars = get_explicit_system().n_vars();
867  unsigned int n_sys_dofs = this->n_dofs();
868 
869  _dof_map_between_systems.resize(n_vars);
870  for (unsigned int var=0; var<n_vars; var++)
871  {
872  _dof_map_between_systems[var].resize(n_sys_dofs);
873  }
874 
875  std::vector<dof_id_type> implicit_sys_dof_indices;
876  std::vector<dof_id_type> explicit_sys_dof_indices;
877 
878  for (const auto & elem : get_mesh().active_element_ptr_range())
879  {
880  this->get_dof_map().dof_indices (elem, implicit_sys_dof_indices);
881 
882  const unsigned int n_dofs = implicit_sys_dof_indices.size();
883 
884  for (unsigned int var=0; var<n_vars; var++)
885  {
886  get_explicit_system().get_dof_map().dof_indices (elem, explicit_sys_dof_indices, var);
887 
888  libmesh_assert(explicit_sys_dof_indices.size() == n_dofs);
889 
890  for (unsigned int i=0; i<n_dofs; i++)
891  {
892  dof_id_type implicit_sys_dof_index = implicit_sys_dof_indices[i];
893  dof_id_type explicit_sys_dof_index = explicit_sys_dof_indices[i];
894 
895  _dof_map_between_systems[var][implicit_sys_dof_index] =
896  explicit_sys_dof_index;
897  }
898  }
899  }
900 }
901 
902 } // namespace libMesh
DenseMatrix< Number > interpolation_matrix
Dense matrix that stores the lower triangular interpolation matrix that can be used.
std::vector< NumericVector< Number > * > basis_functions
The libMesh vectors storing the finite element coefficients of the RB basis functions.
void get_explicit_sys_subvector(NumericVector< Number > &dest, unsigned int var, NumericVector< Number > &localized_source)
Load the subvector of localized_source corresponding to variable var into dest.
virtual void initialize_rb_construction(bool skip_matrix_assembly=false, bool skip_vector_assembly=false)
Allocate all the data structures necessary for the construction stage of the RB method.
double abs(double a)
This is the EquationSystems class.
UniquePtr< NumericVector< Number > > _ghosted_meshfunction_vector
We also need an extra vector in which we can store a ghosted copy of the vector that we wish to use M...
bool evaluate_RB_error_bound
Boolean to indicate whether we evaluate a posteriori error bounds when rb_solve is called...
virtual Real get_RB_error_bound() libmesh_override
Override to return the best fit error.
void set_point_locator_tol(Real point_locator_tol)
Set a point locator tolerance to be used in this class&#39;s MeshFunction, and other operations that requ...
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1582
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true) libmesh_override
Override train_reduced_basis to first initialize _parametrized_functions_in_training_set.
virtual numeric_index_type last_local_index() const =0
T libmesh_conj(T a)
bool assert_convergence
A boolean flag to indicate whether to check for proper convergence after each solve.
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true)
Train the reduced basis.
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
DenseVector< Number > RB_solution
The RB solution vector.
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
virtual void print_info() libmesh_override
Print out info that describes the current setup of this RBConstruction.
void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs)
Allow the implicit_neighbor_dofs flag to be set programmatically.
Definition: dof_map.C:1671
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
ImplicitSystem & sys
void maxloc(T &r, unsigned int &max_id) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors, returning the minimum rank of a processor which originally held the maximum value.
std::string _explicit_system_name
We use an ExplicitSystem to store the EIM basis functions.
DenseMatrix< Number > RB_inner_product_matrix
The inner product matrix.
virtual void enrich_RB_space() libmesh_override
Add a new basis function to the RB space.
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
std::vector< NumericVector< Number > * > _parametrized_functions_in_training_set
The libMesh vectors storing the finite element coefficients of the RB basis functions.
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:1802
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
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:262
static const Real TOLERANCE
UniquePtr< MeshFunction > _mesh_function
A mesh function to interpolate on the mesh.
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual UniquePtr< ElemAssembly > build_eim_assembly(unsigned int bf_index)=0
Build an element assembly object that will access basis function bf_index.
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
UniquePtr< LinearSolver< Number > > inner_product_solver
We store an extra linear solver object which we can optionally use for solving all systems in which t...
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
const std::string & name() const
Definition: system.h:1998
virtual void clear() libmesh_override
Clear this object.
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
virtual void init_explicit_system()=0
Add variables to the ExplicitSystem that is used to store the basis functions.
libmesh_assert(j)
numeric_index_type get_n_training_samples() const
Get the total number of training samples.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual void update_system() libmesh_override
Update the system after enriching the RB space; this calls a series of functions to update the system...
std::vector< ElemAssembly * > get_eim_assembly_objects()
virtual Real compute_best_fit_error()
We compute the best fit of parametrized_function into the EIM space and then evaluate the error in th...
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
virtual void solve_for_matrix_and_rhs(LinearSolver< Number > &input_solver, SparseMatrix< Number > &input_matrix, NumericVector< Number > &input_rhs)
Assembles & solves the linear system A*x=b for the specified matrix input_matrix and right-hand side ...
bool use_empty_rb_solve_in_greedy
A boolean flag to indicate whether or not we initialize the Greedy algorithm by performing rb_solves ...
const MeshBase & get_mesh() const
Definition: system.h:2014
virtual Real truth_solve(int plot_solution) libmesh_override
Load the truth representation of the parametrized function at the current parameters into the solutio...
const DofMap & get_dof_map() const
Definition: system.h:2030
void set_rb_assembly_expansion(RBAssemblyExpansion &rb_assembly_expansion_in)
Set the rb_assembly_expansion object.
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
std::vector< NumericVector< Number > * > _matrix_times_bfs
This vector is used to store inner_product_matrix * basis_function[i] for each i, since we frequently...
virtual Real rb_solve(unsigned int N) libmesh_override
Calculate the EIM approximation to parametrized_function using the first N EIM basis functions...
void set_parameters(const RBParameters &params)
Set the current parameters to params.
virtual void init_implicit_system()=0
Add one variable to the ImplicitSystem (i.e.
void set_params_from_training_set(unsigned int index)
Set parameters to the RBParameters stored in index index of the training set.
ExplicitSystem & get_explicit_system()
Get the ExplicitSystem associated with this system.
Real _point_locator_tol
The point locator tolerance.
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.
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
NumericVector< Number > & get_basis_function(unsigned int i)
Get a reference to the i^th basis function.
virtual void init_data() libmesh_override
Override to initialize the coupling matrix to decouple variables in this system.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
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.
Definition: dense_matrix.C:553
virtual void update_RB_system_matrices() libmesh_override
Compute the reduced basis matrices for the current basis.
virtual void load_basis_function(unsigned int i) libmesh_override
Load the i^th RB function into the RBConstruction solution vector.
This class extends FEMContext in order to provide extra data required to perform local element residu...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:324
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
virtual void print_info()
Print out info that describes the current setup of this RBConstruction.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
std::vector< ElemAssembly * > _rb_eim_assembly_objects
The vector of assembly objects that are created to point to this RBEIMConstruction.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void broadcast(T &data, const unsigned int root_id=0) const
Take a local value and broadcast it to all processors.
std::vector< unsigned int > interpolation_points_var
The corresponding list of variables indices at which the interpolation points were identified...
std::vector< Elem * > interpolation_points_elem
The corresponding list of elements at which the interpolation points were identified.
virtual void elem_fe_reinit(const std::vector< Point > *const pts=libmesh_nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1382
BEST_FIT_TYPE best_fit_type_flag
Enum that indicates which type of "best fit" algorithm we should use.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void initialize_rb_construction(bool skip_matrix_assembly=false, bool skip_vector_assembly=false) libmesh_override
Initialize this system so that we can perform the Construction stage of the RB method.
void initialize_parametrized_functions_in_training_set()
Loop over the training set and compute the parametrized function for each training index...
virtual Real rb_solve(unsigned int N)
Perform online solve with the N RB basis functions, for the set of parameters in current_params, where 0 <= N <= RB_size.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
dof_id_type n_local_dofs() const
Definition: system.C:185
void libmesh_ignore(const T &)
Number evaluate_parametrized_function(unsigned int var_index, const Point &p, const Elem &elem)
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
virtual void process_parameters_file(const std::string &parameters_filename) libmesh_override
Read parameters in from file and set up this system accordingly.
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number evaluate_mesh_function(unsigned int var_number, Point p)
Evaluate the mesh function at the specified point and for the specified variable. ...
void check_convergence(LinearSolver< Number > &input_solver)
Check if the linear solver reports convergence.
RBParameters get_params_from_training_set(unsigned int index)
Return the RBParameters in index index of training set.
const Parallel::Communicator & comm() const
OStreamProxy out
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
static const bool value
Definition: xdr_io.C:108
const EquationSystems & get_equation_systems() const
Definition: system.h:712
virtual numeric_index_type first_local_index() const =0
virtual void init_context_with_sys(FEMContext &c, System &sys)
Initialize c based on sys.
bool serial_training_set
This boolean flag indicates whether or not the training set should be the same on all processors...
void plot_parametrized_functions_in_training_set(const std::string &pathname)
Plot all the parameterized functions that we are storing in _parametrized_functions_in_training_set.
virtual void process_parameters_file(const std::string &parameters_filename)
Read in from the file specified by parameters_filename and set the this system&#39;s member variables acc...
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
dof_id_type n_dofs() const
Definition: system.C:148
std::vector< std::vector< dof_id_type > > _dof_map_between_systems
The index map between the explicit system and the implicit system.
const T_sys & get_system(const std::string &name) const
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
unsigned int n_vars() const
Definition: system.h:2086
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
RBAssemblyExpansion _empty_rb_assembly_expansion
We initialize RBEIMConstruction so that it has an "empty" RBAssemblyExpansion, because this isn&#39;t use...
virtual ~RBEIMConstruction()
Destructor.
This class provides function-like objects for data distributed over a mesh.
Definition: mesh_function.h:53
const RBParameters & get_parameters() const
Get the current parameters.
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
void init_dof_map_between_systems()
Set up the index map between the implicit and explicit systems.
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.
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual dof_id_type n_elem() const =0
std::vector< Point > interpolation_points
The list of interpolation points, i.e.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
This class is part of the rbOOmit framework.
void get_all_variable_numbers(std::vector< unsigned int > &all_variable_numbers) const
Fills all_variable_numbers with all the variable numbers for the variables that have been added to th...
Definition: system.C:1278
bool _parametrized_functions_in_training_set_initialized
Boolean flag to indicate whether or not we have called compute_parametrized_functions_in_training_set...
RBEIMConstruction(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
unsigned int n_points() const
Definition: quadrature.h:113
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
virtual void initialize_eim_assembly_objects()
Build a vector of ElemAssembly objects that accesses the basis functions stored in this RBEIMConstruc...
virtual void load_rb_solution() libmesh_override
Load the RB solution from the most recent solve with rb_eval into this system&#39;s solution vector...
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
void set_explicit_sys_subvector(NumericVector< Number > &dest, unsigned int var, NumericVector< Number > &source)
Load source into the subvector of dest corresponding to var var.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917