libMesh
rb_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 // rbOOmit includes
21 #include "libmesh/rb_construction.h"
22 #include "libmesh/rb_assembly_expansion.h"
23 #include "libmesh/rb_evaluation.h"
24 #include "libmesh/elem_assembly.h"
25 
26 // LibMesh includes
27 #include "libmesh/numeric_vector.h"
28 #include "libmesh/sparse_matrix.h"
29 #include "libmesh/dof_map.h"
30 #include "libmesh/libmesh_logging.h"
31 #include "libmesh/equation_systems.h"
32 #include "libmesh/exodusII_io.h"
33 #include "libmesh/gmv_io.h"
34 #include "libmesh/linear_solver.h"
35 #include "libmesh/getpot.h"
36 #include "libmesh/mesh_base.h"
37 #include "libmesh/parallel.h"
38 #include "libmesh/xdr_cxx.h"
39 #include "libmesh/timestamp.h"
40 #include "libmesh/petsc_linear_solver.h"
41 #include "libmesh/dg_fem_context.h"
42 #include "libmesh/dirichlet_boundaries.h"
43 #include "libmesh/zero_function.h"
44 #include "libmesh/coupling_matrix.h"
45 #include "libmesh/face_tri3_subdivision.h"
46 #include "libmesh/quadrature.h"
47 #include "libmesh/utility.h"
48 
49 // C++ includes
50 #include <sys/types.h>
51 #include <sys/stat.h>
52 #include <errno.h>
53 #include <fstream>
54 #include <sstream>
55 #include <limits>
56 
57 namespace libMesh
58 {
59 
61  const std::string & name_in,
62  const unsigned int number_in)
63  : Parent(es, name_in, number_in),
64  inner_product_solver(LinearSolver<Number>::build(es.comm())),
65  extra_linear_solver(libmesh_nullptr),
66  inner_product_matrix(SparseMatrix<Number>::build(es.comm())),
67  exit_on_repeated_greedy_parameters(true),
68  impose_internal_fluxes(false),
69  compute_RB_inner_product(false),
70  store_non_dirichlet_operators(false),
71  use_empty_rb_solve_in_greedy(true),
72  Fq_representor_innerprods_computed(false),
73  Nmax(0),
74  delta_N(1),
75  quiet_mode(true),
76  output_dual_innerprods_computed(false),
77  assert_convergence(true),
78  rb_eval(libmesh_nullptr),
79  inner_product_assembly(libmesh_nullptr),
80  rel_training_tolerance(1.e-4),
81  abs_training_tolerance(1.e-12),
82  normalize_rb_bound_in_greedy(false)
83 {
84  // set assemble_before_solve flag to false
85  // so that we control matrix assembly.
86  assemble_before_solve = false;
87 }
88 
89 
91 {
92  this->clear();
93 }
94 
95 
97 {
98  LOG_SCOPE("clear()", "RBConstruction");
99 
100  Parent::clear();
101 
102  for (std::size_t q=0; q<Aq_vector.size(); q++)
103  {
104  delete Aq_vector[q];
106  }
107 
108  for (std::size_t q=0; q<Fq_vector.size(); q++)
109  {
110  delete Fq_vector[q];
112  }
113 
114  for (std::size_t i=0; i<outputs_vector.size(); i++)
115  for (std::size_t q_l=0; q_l<outputs_vector[i].size(); q_l++)
116  {
117  delete outputs_vector[i][q_l];
119  }
120 
122  {
123  for (std::size_t q=0; q<non_dirichlet_Aq_vector.size(); q++)
124  {
125  delete non_dirichlet_Aq_vector[q];
127  }
128 
129  for (std::size_t q=0; q<non_dirichlet_Fq_vector.size(); q++)
130  {
131  delete non_dirichlet_Fq_vector[q];
133  }
134 
135  for (std::size_t i=0; i<non_dirichlet_outputs_vector.size(); i++)
136  for (std::size_t q_l=0; q_l<non_dirichlet_outputs_vector[i].size(); q_l++)
137  {
138  delete non_dirichlet_outputs_vector[i][q_l];
140  }
141  }
142 
143  // Also delete the Fq representors
144  for (std::size_t q_f=0; q_f<Fq_representor.size(); q_f++)
145  {
146  delete Fq_representor[q_f];
148  }
149  // Set Fq_representor_innerprods_computed flag to false now
150  // that we've cleared the Fq representors
152 }
153 
154 std::string RBConstruction::system_type () const
155 {
156  return "RBConstruction";
157 }
158 
160  SparseMatrix<Number> & input_matrix,
161  NumericVector<Number> & input_rhs)
162 {
163  // This is similar to LinearImplicitSysmte::solve()
164 
165  // Get a reference to the EquationSystems
166  const EquationSystems & es =
167  this->get_equation_systems();
168 
169  // If the linear solver hasn't been initialized, we do so here.
170  input_solver.init();
171 
172  // Get the user-specifiied linear solver tolerance
173  const Real tol =
174  es.parameters.get<Real>("linear solver tolerance");
175 
176  // Get the user-specified maximum # of linear solver iterations
177  const unsigned int maxits =
178  es.parameters.get<unsigned int>("linear solver maximum iterations");
179 
180  // Solve the linear system. Several cases:
181  std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
182 
183  // It's good practice to clear the solution vector first since it can
184  // affect convergence of iterative solvers
185  solution->zero();
186  rval = input_solver.solve (input_matrix, *solution, input_rhs, tol, maxits);
187 
188  // Store the number of linear iterations required to
189  // solve and the final residual.
190  _n_linear_iterations = rval.first;
191  _final_linear_residual = rval.second;
192 
193  // Update the system after the solve
194  this->update();
195 }
196 
198 {
199  rb_eval = &rb_eval_in;
200 }
201 
203 {
204  if (!rb_eval)
205  libmesh_error_msg("Error: RBEvaluation object hasn't been initialized yet");
206 
207  return *rb_eval;
208 }
209 
211 {
212  return (rb_eval != libmesh_nullptr);
213 }
214 
216 {
218 }
219 
220 void RBConstruction::process_parameters_file (const std::string & parameters_filename)
221 {
222  // First read in data from input_filename
223  GetPot infile(parameters_filename);
224 
225  const unsigned int n_training_samples = infile("n_training_samples",0);
226  const bool deterministic_training = infile("deterministic_training",false);
227  unsigned int training_parameters_random_seed_in =
228  static_cast<unsigned int>(-1);
229  training_parameters_random_seed_in = infile("training_parameters_random_seed",
230  training_parameters_random_seed_in);
231  const bool quiet_mode_in = infile("quiet_mode", quiet_mode);
232  const unsigned int Nmax_in = infile("Nmax", Nmax);
233  const Real rel_training_tolerance_in = infile("rel_training_tolerance",
235  const Real abs_training_tolerance_in = infile("abs_training_tolerance",
237 
238  // Initialize value to false, let the input file value override.
239  bool normalize_rb_bound_in_greedy = false;
240  const bool normalize_rb_bound_in_greedy_in = infile("normalize_rb_bound_in_greedy",
241  normalize_rb_bound_in_greedy);
242 
243  // Read in the parameters from the input file too
244  unsigned int n_continuous_parameters = infile.vector_variable_size("parameter_names");
245  RBParameters mu_min_in;
246  RBParameters mu_max_in;
247  for (unsigned int i=0; i<n_continuous_parameters; i++)
248  {
249  // Read in the parameter names
250  std::string param_name = infile("parameter_names", "NONE", i);
251 
252  {
253  Real min_val = infile(param_name, 0., 0);
254  mu_min_in.set_value(param_name, min_val);
255  }
256 
257  {
258  Real max_val = infile(param_name, 0., 1);
259  mu_max_in.set_value(param_name, max_val);
260  }
261  }
262 
263  std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
264 
265  unsigned int n_discrete_parameters = infile.vector_variable_size("discrete_parameter_names");
266  for (unsigned int i=0; i<n_discrete_parameters; i++)
267  {
268  std::string param_name = infile("discrete_parameter_names", "NONE", i);
269 
270  unsigned int n_vals_for_param = infile.vector_variable_size(param_name);
271  std::vector<Real> vals_for_param(n_vals_for_param);
272  for (std::size_t j=0; j<vals_for_param.size(); j++)
273  vals_for_param[j] = infile(param_name, 0., j);
274 
275  discrete_parameter_values_in[param_name] = vals_for_param;
276  }
277 
278  std::map<std::string,bool> log_scaling_in;
279  RBParameters::const_iterator it = mu_min_in.begin();
280  RBParameters::const_iterator it_end = mu_min_in.end();
281  for ( ; it != it_end; ++it)
282  {
283  std::string param_name = it->first;
284 
285  // For now, just set all entries to false.
286  // TODO: Implement a decent way to specify log-scaling true/false
287  // in the input text file
288  log_scaling_in[param_name] = false;
289  }
290 
291  // Set the parameters that have been read in
292  set_rb_construction_parameters(n_training_samples,
293  deterministic_training,
294  training_parameters_random_seed_in,
295  quiet_mode_in,
296  Nmax_in,
297  rel_training_tolerance_in,
298  abs_training_tolerance_in,
299  normalize_rb_bound_in_greedy_in,
300  mu_min_in,
301  mu_max_in,
302  discrete_parameter_values_in,
303  log_scaling_in);
304 }
305 
307  unsigned int n_training_samples_in,
308  bool deterministic_training_in,
309  unsigned int training_parameters_random_seed_in,
310  bool quiet_mode_in,
311  unsigned int Nmax_in,
312  Real rel_training_tolerance_in,
313  Real abs_training_tolerance_in,
314  bool normalize_rb_bound_in_greedy_in,
315  RBParameters mu_min_in,
316  RBParameters mu_max_in,
317  std::map<std::string, std::vector<Real>> discrete_parameter_values_in,
318  std::map<std::string,bool> log_scaling_in)
319 {
320  // Read in training_parameters_random_seed value. This is used to
321  // seed the RNG when picking the training parameters. By default the
322  // value is -1, which means use std::time to seed the RNG.
323  set_training_random_seed(training_parameters_random_seed_in);
324 
325  // Set quiet mode
326  set_quiet_mode(quiet_mode_in);
327 
328  // Initialize RB parameters
329  set_Nmax(Nmax_in);
330 
331  set_rel_training_tolerance(rel_training_tolerance_in);
332  set_abs_training_tolerance(abs_training_tolerance_in);
333 
334  set_normalize_rb_bound_in_greedy(normalize_rb_bound_in_greedy_in);
335 
336  // Initialize the parameter ranges and the parameters themselves
337  initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_in);
338 
340  this->get_parameters_max(),
341  n_training_samples_in,
342  log_scaling_in,
343  deterministic_training_in); // use deterministic parameters
344 }
345 
347 {
348  // Print out info that describes the current setup
349  libMesh::out << std::endl << "RBConstruction parameters:" << std::endl;
350  libMesh::out << "system name: " << this->name() << std::endl;
351  libMesh::out << "Nmax: " << Nmax << std::endl;
352  libMesh::out << "Greedy relative error tolerance: " << get_rel_training_tolerance() << std::endl;
353  libMesh::out << "Greedy absolute error tolerance: " << get_abs_training_tolerance() << std::endl;
354  libMesh::out << "Do we normalize RB error bound in greedy? " << get_normalize_rb_bound_in_greedy() << std::endl;
356  {
357  libMesh::out << "Aq operators attached: " << get_rb_theta_expansion().get_n_A_terms() << std::endl;
358  libMesh::out << "Fq functions attached: " << get_rb_theta_expansion().get_n_F_terms() << std::endl;
359  libMesh::out << "n_outputs: " << get_rb_theta_expansion().get_n_outputs() << std::endl;
360  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
361  libMesh::out << "output " << n << ", Q_l = " << get_rb_theta_expansion().get_n_output_terms(n) << std::endl;
362  }
363  else
364  {
365  libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
366  }
367  libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
370  for ( ; it != it_end; ++it)
371  {
372  std::string param_name = it->first;
373  if (!is_discrete_parameter(param_name))
374  {
375  libMesh::out << "Parameter " << param_name
376  << ": Min = " << get_parameter_min(param_name)
377  << ", Max = " << get_parameter_max(param_name) << std::endl;
378  }
379  }
381  libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
382  libMesh::out << "quiet mode? " << is_quiet() << std::endl;
383  libMesh::out << std::endl;
384 }
385 
387 {
388  UniquePtr<NumericVector<Number>> temp = solution->clone();
389 
390  for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
391  {
392  for (unsigned int j=0; j<get_rb_evaluation().get_n_basis_functions(); j++)
393  {
394  inner_product_matrix->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
395  Number value = temp->dot( get_rb_evaluation().get_basis_function(i) );
396 
397  libMesh::out << value << " ";
398  }
399  libMesh::out << std::endl;
400  }
401  libMesh::out << std::endl;
402 }
403 
405 {
406  rb_assembly_expansion = &rb_assembly_expansion_in;
407 }
408 
410 {
412  libmesh_error_msg("Error: RBAssemblyExpansion object hasn't been initialized yet");
413 
414  return *rb_assembly_expansion;
415 }
416 
418 {
419  inner_product_assembly = &inner_product_assembly_in;
420 }
421 
423 {
425  libmesh_error_msg("Error: inner_product_assembly hasn't been initialized yet");
426 
427  return *inner_product_assembly;
428 }
429 
431 {
432  const DofMap & dof_map = get_dof_map();
433 
434  for (dof_id_type i=dof_map.first_dof(); i<dof_map.end_dof(); i++)
435  {
437  {
438  vector.set(i, 0.);
439  }
440  }
441  vector.close();
442 }
443 
444 void RBConstruction::initialize_rb_construction(bool skip_matrix_assembly,
445  bool skip_vector_assembly)
446 {
447  if (!skip_matrix_assembly && !skip_vector_assembly)
448  {
449  // Check that the theta and assembly objects are consistently sized
450  libmesh_assert_equal_to (get_rb_theta_expansion().get_n_A_terms(), get_rb_assembly_expansion().get_n_A_terms());
451  libmesh_assert_equal_to (get_rb_theta_expansion().get_n_F_terms(), get_rb_assembly_expansion().get_n_F_terms());
452  libmesh_assert_equal_to (get_rb_theta_expansion().get_n_outputs(), get_rb_assembly_expansion().get_n_outputs());
453  for (unsigned int i=0; i<get_rb_theta_expansion().get_n_outputs(); i++)
454  {
455  libmesh_assert_equal_to (get_rb_theta_expansion().get_n_output_terms(i),
456  get_rb_assembly_expansion().get_n_output_terms(i));
457  }
458  }
459 
460  // Perform the initialization
462  assemble_affine_expansion(skip_matrix_assembly, skip_vector_assembly);
463 
464  // inner_product_solver performs solves with the same matrix every time
465  // hence we can set reuse_preconditioner(true).
466  inner_product_solver->reuse_preconditioner(true);
467 
468  // The primary solver is used for truth solves and other solves that
469  // require different matrices, so set reuse_preconditioner(false).
471 
472 }
473 
474 void RBConstruction::assemble_affine_expansion(bool skip_matrix_assembly,
475  bool skip_vector_assembly)
476 {
477  if (!skip_matrix_assembly)
478  {
479  // Assemble and store all of the matrices
480  this->assemble_misc_matrices();
482  }
483 
484  if (!skip_vector_assembly)
485  {
486  // Assemble and store all of the vectors
489  }
490 }
491 
493 {
494  // Resize vectors for storing mesh-dependent data
495  Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
496  Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
497 
498  // Resize the Fq_representors and initialize each to NULL
499  // These are basis independent and hence stored here, whereas
500  // the Aq_representors are stored in RBEvaluation
501  Fq_representor.resize(get_rb_theta_expansion().get_n_F_terms());
502 
503  // Initialize vectors for the inner products of the Fq representors
504  // These are basis independent and therefore stored here.
505  unsigned int Q_f_hat = get_rb_theta_expansion().get_n_F_terms()*(get_rb_theta_expansion().get_n_F_terms()+1)/2;
506  Fq_representor_innerprods.resize(Q_f_hat);
507 
508  // Resize the output vectors
509  outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
510  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
512 
513  // Resize the output dual norm vectors
514  output_dual_innerprods.resize(get_rb_theta_expansion().get_n_outputs());
515  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
516  {
518  output_dual_innerprods[n].resize(Q_l_hat);
519  }
520 
521  {
522  DofMap & dof_map = this->get_dof_map();
523 
525  inner_product_matrix->init();
526  inner_product_matrix->zero();
527 
528  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
529  {
530  // Initialize the memory for the matrices
531  Aq_vector[q] = SparseMatrix<Number>::build(this->comm()).release();
532  dof_map.attach_matrix(*Aq_vector[q]);
533  Aq_vector[q]->init();
534  Aq_vector[q]->zero();
535  }
536 
537  // We also need to initialize a second set of non-Dirichlet operators
539  {
540  non_dirichlet_Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
541  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
542  {
543  // Initialize the memory for the matrices
546  non_dirichlet_Aq_vector[q]->init();
547  non_dirichlet_Aq_vector[q]->zero();
548  }
549  }
550  }
551 
552  // Initialize the vectors
553  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
554  {
555  // Initialize the memory for the vectors
556  Fq_vector[q] = NumericVector<Number>::build(this->comm()).release();
557  Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
558  }
559 
560  // We also need to initialize a second set of non-Dirichlet operators
562  {
563  non_dirichlet_Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
564  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
565  {
566  // Initialize the memory for the vectors
568  non_dirichlet_Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
569  }
570  }
571 
572  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
573  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
574  {
575  // Initialize the memory for the truth output vectors
576  outputs_vector[n][q_l] = (NumericVector<Number>::build(this->comm()).release());
577  outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
578  }
579 
581  {
582  non_dirichlet_outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
583  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
584  {
585  non_dirichlet_outputs_vector[n].resize( get_rb_theta_expansion().get_n_output_terms(n) );
586  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
587  {
588  // Initialize the memory for the truth output vectors
589  non_dirichlet_outputs_vector[n][q_l] = (NumericVector<Number>::build(this->comm()).release());
590  non_dirichlet_outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
591  }
592  }
593  }
594 
595  // Resize truth_outputs vector
596  truth_outputs.resize(this->get_rb_theta_expansion().get_n_outputs());
597 }
598 
600 {
601  return UniquePtr<DGFEMContext>(new DGFEMContext(*this));
602 }
603 
605  ElemAssembly * elem_assembly,
606  SparseMatrix<Number> * input_matrix,
607  NumericVector<Number> * input_vector,
608  bool symmetrize,
609  bool apply_dof_constraints)
610 {
611  LOG_SCOPE("add_scaled_matrix_and_vector()", "RBConstruction");
612 
613  bool assemble_matrix = (input_matrix != libmesh_nullptr);
614  bool assemble_vector = (input_vector != libmesh_nullptr);
615 
616  if (!assemble_matrix && !assemble_vector)
617  return;
618 
619  const MeshBase & mesh = this->get_mesh();
620 
621  // First add any node-based terms (e.g. point loads)
622  // We only enter this loop if we have at least one
623  // nodeset, since we use nodesets to indicate
624  // where to impose the node-based terms.
625  if (mesh.get_boundary_info().n_nodeset_conds() > 0)
626  {
627  std::vector<numeric_index_type> node_id_list;
628  std::vector<boundary_id_type> bc_id_list;
629 
630  // Get the list of nodes with boundary IDs
631  mesh.get_boundary_info().build_node_list(node_id_list, bc_id_list);
632 
633  for (std::size_t i=0; i<node_id_list.size(); i++)
634  {
635  const Node & node = mesh.node_ref(node_id_list[i]);
636 
637  // If node is on this processor, then all dofs on node are too
638  // so we can do the add below safely
639  if (node.processor_id() == this->comm().rank())
640  {
641  // Get the values to add to the rhs vector
642  std::vector<dof_id_type> nodal_dof_indices;
643  DenseMatrix<Number> nodal_matrix;
644  DenseVector<Number> nodal_rhs;
645  elem_assembly->get_nodal_values(
646  nodal_dof_indices, nodal_matrix, nodal_rhs, *this, node);
647 
648  if(!nodal_dof_indices.empty())
649  {
650  if(assemble_vector)
651  {
652  input_vector->add_vector(nodal_rhs, nodal_dof_indices);
653  }
654 
655  if(assemble_matrix)
656  {
657  input_matrix->add_matrix(nodal_matrix, nodal_dof_indices);
658  }
659  }
660 #ifdef LIBMESH_ENABLE_DEPRECATED
661  else if(elem_assembly->is_nodal_rhs_values_overriden)
662  {
663  // This is the "old" implementation, to be deprecated soon.
664  // Only enter this if ElemAssembly::get_nodal_rhs_values has
665  // a non-default implementation.
666  libmesh_deprecated();
667 
668  std::map<numeric_index_type, Number> rhs_values;
669  elem_assembly->get_nodal_rhs_values(rhs_values, *this, node);
670 
671  std::map<numeric_index_type, Number>::const_iterator it =
672  rhs_values.begin();
673  const std::map<numeric_index_type, Number>::const_iterator it_end =
674  rhs_values.end();
675  for ( ; it != it_end; ++it)
676  {
677  numeric_index_type dof_index = it->first;
678  Number value = it->second;
679  input_vector->add( dof_index, value);
680  }
681  }
682 #endif
683  }
684  }
685  }
686 
688  DGFEMContext & context = cast_ref<DGFEMContext &>(*c);
689 
690  this->init_context(context);
691 
692  for (const auto & elem : mesh.active_local_element_ptr_range())
693  {
694  // Subdivision elements need special care:
695  // - skip ghost elements
696  // - init special quadrature rule
697  UniquePtr<QBase> qrule;
698  if (elem->type() == TRI3SUBDIVISION)
699  {
700  const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
701  if (gh_elem->is_ghost())
702  continue ;
703  // A Gauss quadrature rule for numerical integration.
704  // For subdivision shell elements, a single Gauss point per
705  // element is sufficient, hence we use extraorder = 0.
706  const int extraorder = 0;
707  FEBase * elem_fe = libmesh_nullptr;
708  context.get_element_fe( 0, elem_fe );
709 
710  qrule = elem_fe->get_fe_type().default_quadrature_rule (2, extraorder);
711 
712  // Tell the finite element object to use our quadrature rule.
713  elem_fe->attach_quadrature_rule (qrule.get());
714  }
715 
716  context.pre_fe_reinit(*this, elem);
717  context.elem_fe_reinit();
718  elem_assembly->interior_assembly(context);
719 
720  for (context.side = 0;
721  context.side != context.get_elem().n_sides();
722  ++context.side )
723  {
724  // May not need to apply fluxes on non-boundary elements
725  if ((context.get_elem().neighbor_ptr(context.get_side()) != libmesh_nullptr) && !impose_internal_fluxes)
726  continue;
727 
728  bool reinit_succeeded = false;
729 
730  // Exceptions can be thrown in side_fe_reinit, e.g. due to a zero or negative
731  // Jacobian on an element side. This is sometimes intentional, e.g. some people
732  // use "degenerate HEX" elements when generating meshes. In order to support this
733  // case we just print a message and skip the side assembly on this element
734  // if an exception occurs during side_fe_reinit.
735  try
736  {
737  context.side_fe_reinit();
738  reinit_succeeded = true;
739  }
740  catch(std::exception& e)
741  {
742  libMesh::err << std::endl << "Error detected when computing element data on side "
743  << static_cast<int>(context.side) << " of element "
744  << context.get_elem().id() << std::endl;
745  libMesh::err << "Skipping assembly on this element side" << std::endl << std::endl;
746  }
747 
748  if(reinit_succeeded)
749  {
750  elem_assembly->boundary_assembly(context);
751  }
752 
753  if (context.dg_terms_are_active())
754  {
755  input_matrix->add_matrix (context.get_elem_elem_jacobian(),
756  context.get_dof_indices(),
757  context.get_dof_indices());
758 
759  input_matrix->add_matrix (context.get_elem_neighbor_jacobian(),
760  context.get_dof_indices(),
761  context.get_neighbor_dof_indices());
762 
763  input_matrix->add_matrix (context.get_neighbor_elem_jacobian(),
764  context.get_neighbor_dof_indices(),
765  context.get_dof_indices());
766 
767  input_matrix->add_matrix (context.get_neighbor_neighbor_jacobian(),
768  context.get_neighbor_dof_indices(),
769  context.get_neighbor_dof_indices());
770  }
771  }
772 
773  // Need to symmetrize before imposing
774  // periodic constraints
775  if (assemble_matrix && symmetrize)
776  {
777  DenseMatrix<Number> Ke_transpose;
778  context.get_elem_jacobian().get_transpose(Ke_transpose);
779  context.get_elem_jacobian() += Ke_transpose;
780  context.get_elem_jacobian() *= 0.5;
781  }
782 
783  if (apply_dof_constraints)
784  {
785  // Apply constraints, e.g. Dirichlet and periodic constraints
787  (context.get_elem_jacobian(), context.get_elem_residual(), context.get_dof_indices() );
788  }
789 
790  // Scale and add to global matrix and/or vector
791  context.get_elem_jacobian() *= scalar;
792  context.get_elem_residual() *= scalar;
793 
794  if (assemble_matrix)
795  {
796 
797  CouplingMatrix * coupling_matrix = get_dof_map()._dof_coupling;
798  if (!coupling_matrix)
799  {
800  // If we haven't defined a _dof_coupling matrix then just add
801  // the whole matrix
802  input_matrix->add_matrix (context.get_elem_jacobian(),
803  context.get_dof_indices() );
804  }
805  else
806  {
807  // Otherwise we should only add the relevant submatrices
808  for (unsigned int var1=0; var1<n_vars(); var1++)
809  {
810  ConstCouplingRow ccr(var1, *coupling_matrix);
813  ccr.begin(); it != end; ++it)
814  {
815  unsigned int var2 = *it;
816 
817  unsigned int sub_m = context.get_elem_jacobian( var1, var2 ).m();
818  unsigned int sub_n = context.get_elem_jacobian( var1, var2 ).n();
819  DenseMatrix<Number> sub_jac(sub_m, sub_n);
820  for (unsigned int row=0; row<sub_m; row++)
821  for (unsigned int col=0; col<sub_n; col++)
822  {
823  sub_jac(row,col) = context.get_elem_jacobian( var1, var2 ).el(row,col);
824  }
825  input_matrix->add_matrix (sub_jac,
826  context.get_dof_indices(var1),
827  context.get_dof_indices(var2) );
828  }
829  }
830  }
831 
832  }
833 
834  if (assemble_vector)
835  input_vector->add_vector (context.get_elem_residual(),
836  context.get_dof_indices() );
837  }
838 
839  if (assemble_matrix)
840  input_matrix->close();
841  if (assemble_vector)
842  input_vector->close();
843 }
844 
846 {
847  // Set current_local_solution = vec so that we can access
848  // vec from DGFEMContext during assembly
849  vec.localize
850  (*current_local_solution, this->get_dof_map().get_send_list());
851 }
852 
854 {
855  LOG_SCOPE("truth_assembly()", "RBConstruction");
856 
857  const RBParameters & mu = get_parameters();
858 
859  this->matrix->zero();
860  this->rhs->zero();
861 
862  this->matrix->close();
863  this->rhs->close();
864 
865  {
866  // We should have already assembled the matrices
867  // and vectors in the affine expansion, so
868  // just use them
869 
870  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
871  {
872  matrix->add(get_rb_theta_expansion().eval_A_theta(q_a, mu), *get_Aq(q_a));
873  }
874 
876  temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
877  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
878  {
879  *temp_vec = *get_Fq(q_f);
880  temp_vec->scale( get_rb_theta_expansion().eval_F_theta(q_f, mu) );
881  rhs->add(*temp_vec);
882  }
883  }
884 
885  this->matrix->close();
886  this->rhs->close();
887 }
888 
890  bool apply_dof_constraints)
891 {
892  input_matrix->zero();
895  input_matrix,
897  false, /* symmetrize */
898  apply_dof_constraints);
899 }
900 
902  SparseMatrix<Number> * input_matrix,
903  bool apply_dof_constraints)
904 {
905  if (q >= get_rb_theta_expansion().get_n_A_terms())
906  libmesh_error_msg("Error: We must have q < Q_a in assemble_Aq_matrix.");
907 
908  input_matrix->zero();
909 
912  input_matrix,
914  false, /* symmetrize */
915  apply_dof_constraints);
916 }
917 
919  unsigned int q_a,
920  SparseMatrix<Number> * input_matrix,
921  bool symmetrize)
922 {
923  LOG_SCOPE("add_scaled_Aq()", "RBConstruction");
924 
925  if (q_a >= get_rb_theta_expansion().get_n_A_terms())
926  libmesh_error_msg("Error: We must have q < Q_a in add_scaled_Aq.");
927 
928  if (!symmetrize)
929  {
930  input_matrix->add(scalar, *get_Aq(q_a));
931  input_matrix->close();
932  }
933  else
934  {
937  input_matrix,
939  symmetrize);
940  }
941 }
942 
944 {
945  libMesh::out << "Assembling inner product matrix" << std::endl;
947 }
948 
950 {
951  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
952  {
953  libMesh::out << "Assembling affine operator " << (q_a+1) << " of "
954  << get_rb_theta_expansion().get_n_A_terms() << std::endl;
955  assemble_Aq_matrix(q_a, get_Aq(q_a));
956  }
957 
959  {
960  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
961  {
962  libMesh::out << "Assembling non-Dirichlet affine operator " << (q_a+1) << " of "
963  << get_rb_theta_expansion().get_n_A_terms() << std::endl;
964  assemble_Aq_matrix(q_a, get_non_dirichlet_Aq(q_a), false);
965  }
966  }
967 }
968 
970 {
971  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
972  {
973  libMesh::out << "Assembling affine vector " << (q_f+1) << " of "
974  << get_rb_theta_expansion().get_n_F_terms() << std::endl;
975  assemble_Fq_vector(q_f, get_Fq(q_f));
976  }
977 
979  {
980  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
981  {
982  libMesh::out << "Assembling non-Dirichlet affine vector " << (q_f+1) << " of "
983  << get_rb_theta_expansion().get_n_F_terms() << std::endl;
984  assemble_Fq_vector(q_f, get_non_dirichlet_Fq(q_f), false);
985  }
986  }
987 
988 }
989 
991  NumericVector<Number> * input_vector,
992  bool apply_dof_constraints)
993 {
994  if (q >= get_rb_theta_expansion().get_n_F_terms())
995  libmesh_error_msg("Error: We must have q < Q_f in assemble_Fq_vector.");
996 
997  input_vector->zero();
998 
1002  input_vector,
1003  false, /* symmetrize */
1004  apply_dof_constraints /* apply_dof_constraints */);
1005 }
1006 
1008 {
1009  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1010  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1011  {
1012  libMesh::out << "Assembling output vector, (" << (n+1) << "," << (q_l+1)
1013  << ") of (" << get_rb_theta_expansion().get_n_outputs()
1014  << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1015  << std::endl;
1016  get_output_vector(n, q_l)->zero();
1019  get_output_vector(n,q_l),
1020  false, /* symmetrize */
1021  true /* apply_dof_constraints */);
1022  }
1023 
1025  {
1026  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1027  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1028  {
1029  libMesh::out << "Assembling non-Dirichlet output vector, (" << (n+1) << "," << (q_l+1)
1030  << ") of (" << get_rb_theta_expansion().get_n_outputs()
1031  << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1032  << std::endl;
1037  false, /* symmetrize */
1038  false /* apply_dof_constraints */);
1039  }
1040  }
1041 }
1042 
1043 Real RBConstruction::train_reduced_basis(const bool resize_rb_eval_data)
1044 {
1045  LOG_SCOPE("train_reduced_basis()", "RBConstruction");
1046 
1047  int count = 0;
1048 
1049  // initialize rb_eval's parameters
1051 
1052  // possibly resize data structures according to Nmax
1053  if (resize_rb_eval_data)
1054  {
1056  }
1057 
1058  // Clear the Greedy param list
1059  for (std::size_t i=0; i<get_rb_evaluation().greedy_param_list.size(); i++)
1060  get_rb_evaluation().greedy_param_list[i].clear();
1061 
1063 
1064  Real training_greedy_error;
1065 
1066 
1067  // If we are continuing from a previous training run,
1068  // we might already be at the max number of basis functions.
1069  // If so, we can just return.
1070  if (get_rb_evaluation().get_n_basis_functions() >= get_Nmax())
1071  {
1072  libMesh::out << "Maximum number of basis functions reached: Nmax = "
1073  << get_Nmax() << std::endl;
1074  return 0.;
1075  }
1076 
1077 
1078  // Compute the dual norms of the outputs if we haven't already done so
1080 
1081  // Compute the Fq Riesz representor dual norms if we haven't already done so
1083 
1084  libMesh::out << std::endl << "---- Performing Greedy basis enrichment ----" << std::endl;
1085  Real initial_greedy_error = 0.;
1086  bool initial_greedy_error_initialized = false;
1087  while (true)
1088  {
1089  libMesh::out << std::endl << "---- Basis dimension: "
1090  << get_rb_evaluation().get_n_basis_functions() << " ----" << std::endl;
1091 
1092  if (count > 0 || (count==0 && use_empty_rb_solve_in_greedy))
1093  {
1094  libMesh::out << "Performing RB solves on training set" << std::endl;
1095  training_greedy_error = compute_max_error_bound();
1096 
1097  libMesh::out << "Maximum error bound is " << training_greedy_error << std::endl << std::endl;
1098 
1099  // record the initial error
1100  if (!initial_greedy_error_initialized)
1101  {
1102  initial_greedy_error = training_greedy_error;
1103  initial_greedy_error_initialized = true;
1104  }
1105 
1106  // Break out of training phase if we have reached Nmax
1107  // or if the training_tolerance is satisfied.
1108  if (greedy_termination_test(training_greedy_error, initial_greedy_error, count))
1109  break;
1110  }
1111 
1112  libMesh::out << "Performing truth solve at parameter:" << std::endl;
1113  print_parameters();
1114 
1115  // Update the list of Greedily selected parameters
1116  this->update_greedy_param_list();
1117 
1118  // Perform an Offline truth solve for the current parameter
1119  truth_solve(-1);
1120 
1121  // Add orthogonal part of the snapshot to the RB space
1122  libMesh::out << "Enriching the RB space" << std::endl;
1123  enrich_RB_space();
1124 
1125  update_system();
1126 
1127  // Increment counter
1128  count++;
1129  }
1130  this->update_greedy_param_list();
1131 
1132  return training_greedy_error;
1133 }
1134 
1136  Real initial_error,
1137  int)
1138 {
1139  if (abs_greedy_error < this->abs_training_tolerance)
1140  {
1141  libMesh::out << "Absolute error tolerance reached." << std::endl;
1142  return true;
1143  }
1144 
1145  Real rel_greedy_error = abs_greedy_error/initial_error;
1146  if (rel_greedy_error < this->rel_training_tolerance)
1147  {
1148  libMesh::out << "Relative error tolerance reached." << std::endl;
1149  return true;
1150  }
1151 
1152  if (get_rb_evaluation().get_n_basis_functions() >= this->get_Nmax())
1153  {
1154  libMesh::out << "Maximum number of basis functions reached: Nmax = "
1155  << get_Nmax() << std::endl;
1156  return true;
1157  }
1158 
1160  {
1161  for (std::size_t i=0; i<get_rb_evaluation().greedy_param_list.size(); i++)
1162  {
1163  RBParameters & previous_parameters = get_rb_evaluation().greedy_param_list[i];
1164  if (previous_parameters == get_parameters())
1165  {
1166  libMesh::out << "Exiting greedy because the same parameters were selected twice"
1167  << std::endl;
1168  return true;
1169  }
1170  }
1171  }
1172 
1173  return false;
1174 }
1175 
1177 {
1179 }
1180 
1182 {
1183  if (i >= get_rb_evaluation().greedy_param_list.size())
1184  libmesh_error_msg("Error: Argument in RBConstruction::get_greedy_parameter is too large.");
1185 
1187 }
1188 
1190 {
1191  LOG_SCOPE("truth_solve()", "RBConstruction");
1192 
1193  truth_assembly();
1194 
1195  // truth_assembly assembles into matrix and rhs, so use those for the solve
1196  if (extra_linear_solver)
1197  {
1198  // If extra_linear_solver has been initialized, then we use it for the
1199  // truth solves.
1201 
1202  if (assert_convergence)
1204  }
1205  else
1206  {
1208 
1209  if (assert_convergence)
1211  }
1212 
1213 
1214 
1215  const RBParameters & mu = get_parameters();
1216 
1217  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1218  {
1219  truth_outputs[n] = 0.;
1220  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1222  get_output_vector(n,q_l)->dot(*solution);
1223  }
1224 
1225  if (plot_solution > 0)
1226  {
1227 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
1228  GMVIO(get_mesh()).write_equation_systems ("truth.gmv",
1229  this->get_equation_systems());
1230 #else
1231 #ifdef LIBMESH_HAVE_EXODUS_API
1233  this->get_equation_systems());
1234 #endif
1235 #endif
1236  }
1237 
1238  // Get the X norm of the truth solution
1239  // Useful for normalizing our true error data
1241  Number truth_X_norm = std::sqrt(inner_product_storage_vector->dot(*solution));
1242 
1243  return libmesh_real(truth_X_norm);
1244 }
1245 
1246 void RBConstruction::set_Nmax(unsigned int Nmax_in)
1247 {
1248  this->Nmax = Nmax_in;
1249 }
1250 
1252 {
1253  LOG_SCOPE("load_basis_function()", "RBConstruction");
1254 
1255  libmesh_assert_less (i, get_rb_evaluation().get_n_basis_functions());
1256 
1258 
1259  this->update();
1260 }
1261 
1263 {
1264  LOG_SCOPE("enrich_RB_space()", "RBConstruction");
1265 
1266  NumericVector<Number> * new_bf = NumericVector<Number>::build(this->comm()).release();
1267  new_bf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1268  *new_bf = *solution;
1269 
1270  for (unsigned int index=0; index<get_rb_evaluation().get_n_basis_functions(); index++)
1271  {
1272  inner_product_matrix->vector_mult(*inner_product_storage_vector, *new_bf);
1273 
1274  Number scalar =
1275  inner_product_storage_vector->dot(get_rb_evaluation().get_basis_function(index));
1276  new_bf->add(-scalar, get_rb_evaluation().get_basis_function(index));
1277  }
1278 
1279  // Normalize new_bf
1280  inner_product_matrix->vector_mult(*inner_product_storage_vector, *new_bf);
1281  Number new_bf_norm = std::sqrt( inner_product_storage_vector->dot(*new_bf) );
1282 
1283  if (new_bf_norm == 0.)
1284  {
1285  new_bf->zero(); // avoid potential nan's
1286  }
1287  else
1288  {
1289  new_bf->scale(1./new_bf_norm);
1290  }
1291 
1292  // load the new basis function into the basis_functions vector.
1293  get_rb_evaluation().basis_functions.push_back( new_bf );
1294 }
1295 
1297 {
1298  libMesh::out << "Updating RB matrices" << std::endl;
1300 
1301  libMesh::out << "Updating RB residual terms" << std::endl;
1302 
1304 }
1305 
1307 {
1309 
1310  Real error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions());
1311 
1313  {
1314  Real error_bound_normalization = get_rb_evaluation().get_error_bound_normalization();
1315 
1316  if ((error_bound < abs_training_tolerance) ||
1317  (error_bound_normalization < abs_training_tolerance))
1318  {
1319  // We don't want to normalize this error bound if the bound or the
1320  // normalization value are below the absolute tolerance. Hence do nothing
1321  // in this case.
1322  }
1323  else
1324  error_bound /= error_bound_normalization;
1325  }
1326 
1327  return error_bound;
1328 }
1329 
1330 void RBConstruction::recompute_all_residual_terms(bool compute_inner_products)
1331 {
1332  // Compute the basis independent terms
1334  compute_Fq_representor_innerprods(compute_inner_products);
1335 
1336  // and all the basis dependent terms
1337  unsigned int saved_delta_N = delta_N;
1339 
1340  update_residual_terms(compute_inner_products);
1341 
1342  delta_N = saved_delta_N;
1343 }
1344 
1346 {
1347  LOG_SCOPE("compute_max_error_bound()", "RBConstruction");
1348 
1349  // Treat the case with no parameters in a special way
1350  if (get_n_params() == 0)
1351  {
1352  Real max_val;
1353  if (std::numeric_limits<Real>::has_infinity)
1354  {
1355  max_val = std::numeric_limits<Real>::infinity();
1356  }
1357  else
1358  {
1359  max_val = std::numeric_limits<Real>::max();
1360  }
1361 
1362  // Make sure we do at least one solve, but otherwise return a zero error bound
1363  // when we have no parameters
1364  return (get_rb_evaluation().get_n_basis_functions() == 0) ? max_val : 0.;
1365  }
1366 
1368 
1369  // keep track of the maximum error
1370  unsigned int max_err_index = 0;
1371  Real max_err = 0.;
1372 
1374  for (unsigned int i=0; i<get_local_n_training_samples(); i++)
1375  {
1376  // Load training parameter i, this is only loaded
1377  // locally since the RB solves are local.
1378  set_params_from_training_set( first_index+i );
1379 
1381 
1382  if (training_error_bounds[i] > max_err)
1383  {
1384  max_err_index = i;
1385  max_err = training_error_bounds[i];
1386  }
1387  }
1388 
1389  std::pair<numeric_index_type, Real> error_pair(first_index+max_err_index, max_err);
1390  get_global_max_error_pair(this->comm(),error_pair);
1391 
1392  // If we have a serial training set (i.e. a training set that is the same on all processors)
1393  // just set the parameters on all processors
1394  if (serial_training_set)
1395  {
1396  set_params_from_training_set( error_pair.first );
1397  }
1398  // otherwise, broadcast the parameter that produced the maximum error
1399  else
1400  {
1401  unsigned int root_id=0;
1402  if ((get_first_local_training_index() <= error_pair.first) &&
1403  (error_pair.first < get_last_local_training_index()))
1404  {
1405  set_params_from_training_set( error_pair.first );
1406  root_id = this->processor_id();
1407  }
1408 
1409  this->comm().sum(root_id); // root_id is only non-zero on one processor
1410  broadcast_parameters(root_id);
1411  }
1412 
1413  return error_pair.second;
1414 }
1415 
1417 {
1418  LOG_SCOPE("update_RB_system_matrices()", "RBConstruction");
1419 
1420  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1421 
1423  temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1424 
1425  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1426  {
1427  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1428  {
1429  get_rb_evaluation().RB_Fq_vector[q_f](i) = get_Fq(q_f)->dot(get_rb_evaluation().get_basis_function(i));
1430  }
1431  }
1432 
1433  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1434  {
1435  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1436  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1437  {
1438  get_rb_evaluation().RB_output_vectors[n][q_l](i) =
1439  get_output_vector(n,q_l)->dot(get_rb_evaluation().get_basis_function(i));
1440  }
1441 
1442  for (unsigned int j=0; j<RB_size; j++)
1443  {
1444  Number value = 0.;
1445 
1447  {
1448  // Compute reduced inner_product_matrix
1449  temp->zero();
1450  inner_product_matrix->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1451 
1452  value = temp->dot( get_rb_evaluation().get_basis_function(i) );
1454  if (i!=j)
1455  {
1456  // The inner product matrix is assumed
1457  // to be hermitian
1459  }
1460  }
1461 
1462  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1463  {
1464  // Compute reduced Aq matrix
1465  temp->zero();
1466  get_Aq(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1467 
1468  value = (*temp).dot(get_rb_evaluation().get_basis_function(i));
1469  get_rb_evaluation().RB_Aq_vector[q_a](i,j) = value;
1470 
1471  if (i!=j)
1472  {
1473  temp->zero();
1474  get_Aq(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(i));
1475 
1476  value = (*temp).dot(get_rb_evaluation().get_basis_function(j));
1477  get_rb_evaluation().RB_Aq_vector[q_a](j,i) = value;
1478  }
1479  }
1480  }
1481  }
1482 }
1483 
1484 
1485 void RBConstruction::update_residual_terms(bool compute_inner_products)
1486 {
1487  LOG_SCOPE("update_residual_terms()", "RBConstruction");
1488 
1489  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1490 
1491  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1492  {
1493  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1494  {
1495  // Initialize the vector in which we'll store the representor
1496  if (!get_rb_evaluation().Aq_representor[q_a][i])
1497  {
1498  get_rb_evaluation().Aq_representor[q_a][i] = (NumericVector<Number>::build(this->comm()).release());
1499  get_rb_evaluation().Aq_representor[q_a][i]->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1500  }
1501 
1502  libmesh_assert(get_rb_evaluation().Aq_representor[q_a][i]->size() == this->n_dofs() &&
1503  get_rb_evaluation().Aq_representor[q_a][i]->local_size() == this->n_local_dofs() );
1504 
1505  rhs->zero();
1506  get_Aq(q_a)->vector_mult(*rhs, get_rb_evaluation().get_basis_function(i));
1507  rhs->scale(-1.);
1508 
1509  if (!is_quiet())
1510  {
1511  libMesh::out << "Starting solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
1512  << Utility::get_timestamp() << std::endl;
1513  }
1514 
1516 
1517  if (assert_convergence)
1519 
1520  if (!is_quiet())
1521  {
1522  libMesh::out << "Finished solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
1523  << Utility::get_timestamp() << std::endl;
1524  libMesh::out << this->n_linear_iterations() << " iterations, final residual "
1525  << this->final_linear_residual() << std::endl;
1526  }
1527 
1528  // Store the representor
1530  }
1531  }
1532 
1533  // Now compute and store the inner products (if requested)
1534  if (compute_inner_products)
1535  {
1536 
1537  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1538  {
1540 
1541  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1542  {
1543  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1544  {
1546  inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a][i]);
1547  }
1548  }
1549  }
1550 
1551  unsigned int q=0;
1552  for (unsigned int q_a1=0; q_a1<get_rb_theta_expansion().get_n_A_terms(); q_a1++)
1553  {
1554  for (unsigned int q_a2=q_a1; q_a2<get_rb_theta_expansion().get_n_A_terms(); q_a2++)
1555  {
1556  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1557  {
1558  for (unsigned int j=0; j<RB_size; j++)
1559  {
1560  inner_product_matrix->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][j]);
1562  inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][i]);
1563 
1564  if (i != j)
1565  {
1566  inner_product_matrix->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][i]);
1568  inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][j]);
1569  }
1570  }
1571  }
1572  q++;
1573  }
1574  }
1575  } // end if (compute_inner_products)
1576 }
1577 
1579 {
1580  return *inner_product_matrix;
1581 }
1582 
1584 {
1585  // Skip calculations if we've already computed the output dual norms
1587  {
1588  // Short circuit if we don't have any outputs
1589  if (get_rb_theta_expansion().get_n_outputs() == 0)
1590  {
1592  return;
1593  }
1594 
1595  // Only log if we get to here
1596  LOG_SCOPE("compute_output_dual_innerprods()", "RBConstruction");
1597 
1598  libMesh::out << "Compute output dual inner products" << std::endl;
1599 
1600  // Find out the largest value of Q_l
1601  unsigned int max_Q_l = 0;
1602  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1603  max_Q_l = (get_rb_theta_expansion().get_n_output_terms(n) > max_Q_l) ? get_rb_theta_expansion().get_n_output_terms(n) : max_Q_l;
1604 
1605  std::vector<NumericVector<Number> *> L_q_representor(max_Q_l);
1606  for (unsigned int q=0; q<max_Q_l; q++)
1607  {
1608  L_q_representor[q] = (NumericVector<Number>::build(this->comm()).release());
1609  L_q_representor[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1610  }
1611 
1612  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1613  {
1614  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1615  {
1616  rhs->zero();
1617  rhs->add(1., *get_output_vector(n,q_l));
1618 
1619  if (!is_quiet())
1620  libMesh::out << "Starting solve n=" << n << ", q_l=" << q_l
1621  << " in RBConstruction::compute_output_dual_innerprods() at "
1622  << Utility::get_timestamp() << std::endl;
1623 
1624  // Use the main linear solver here instead of the inner_product solver, since
1625  // get_matrix_for_output_dual_solves() may not return the inner product matrix.
1627 
1628  // We possibly perform multiple solves here with the same matrix, hence
1629  // set reuse_preconditioner(true) (and set it back to false again below
1630  // at the end of this function).
1631  linear_solver->reuse_preconditioner(true);
1632 
1633  if (assert_convergence)
1635 
1636  if (!is_quiet())
1637  {
1638  libMesh::out << "Finished solve n=" << n << ", q_l=" << q_l
1639  << " in RBConstruction::compute_output_dual_innerprods() at "
1640  << Utility::get_timestamp() << std::endl;
1641 
1643  << " iterations, final residual "
1644  << this->final_linear_residual() << std::endl;
1645  }
1646 
1647  *L_q_representor[q_l] = *solution;
1648  }
1649 
1650  unsigned int q=0;
1651  for (unsigned int q_l1=0; q_l1<get_rb_theta_expansion().get_n_output_terms(n); q_l1++)
1652  {
1654 
1655  for (unsigned int q_l2=q_l1; q_l2<get_rb_theta_expansion().get_n_output_terms(n); q_l2++)
1656  {
1657  output_dual_innerprods[n][q] = L_q_representor[q_l2]->dot(*inner_product_storage_vector);
1658  libMesh::out << "output_dual_innerprods[" << n << "][" << q << "] = " << output_dual_innerprods[n][q] << std::endl;
1659 
1660  q++;
1661  }
1662  }
1663  }
1664 
1665  // Finally clear the L_q_representor vectors
1666  for (unsigned int q=0; q<max_Q_l; q++)
1667  {
1668  if (L_q_representor[q])
1669  {
1670  delete L_q_representor[q];
1671  L_q_representor[q] = libmesh_nullptr;
1672  }
1673  }
1674 
1675  // We may not need to use linear_solver again (e.g. this would happen if we use
1676  // extra_linear_solver for the truth_solves). As a result, let's clear linear_solver
1677  // to release any memory it may be taking up. If we do need it again, it will
1678  // be initialized when necessary.
1679  linear_solver->clear();
1680  linear_solver->reuse_preconditioner(false);
1681 
1683  }
1684 
1686 }
1687 
1688 void RBConstruction::compute_Fq_representor_innerprods(bool compute_inner_products)
1689 {
1690 
1691  // Skip calculations if we've already computed the Fq_representors
1693  {
1694  // Only log if we get to here
1695  LOG_SCOPE("compute_Fq_representor_innerprods()", "RBConstruction");
1696 
1697  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1698  {
1699  if (!Fq_representor[q_f])
1700  {
1701  Fq_representor[q_f] = (NumericVector<Number>::build(this->comm()).release());
1702  Fq_representor[q_f]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1703  }
1704 
1705  libmesh_assert(Fq_representor[q_f]->size() == this->n_dofs() &&
1706  Fq_representor[q_f]->local_size() == this->n_local_dofs() );
1707 
1708  rhs->zero();
1709  rhs->add(1., *get_Fq(q_f));
1710 
1711  if (!is_quiet())
1712  libMesh::out << "Starting solve q_f=" << q_f
1713  << " in RBConstruction::update_residual_terms() at "
1714  << Utility::get_timestamp() << std::endl;
1715 
1717 
1718  if (assert_convergence)
1720 
1721  if (!is_quiet())
1722  {
1723  libMesh::out << "Finished solve q_f=" << q_f
1724  << " in RBConstruction::update_residual_terms() at "
1725  << Utility::get_timestamp() << std::endl;
1726 
1728  << " iterations, final residual "
1729  << this->final_linear_residual() << std::endl;
1730  }
1731 
1732  *Fq_representor[q_f] = *solution;
1733  }
1734 
1735  if (compute_inner_products)
1736  {
1737  unsigned int q=0;
1738 
1739  for (unsigned int q_f1=0; q_f1<get_rb_theta_expansion().get_n_F_terms(); q_f1++)
1740  {
1742 
1743  for (unsigned int q_f2=q_f1; q_f2<get_rb_theta_expansion().get_n_F_terms(); q_f2++)
1744  {
1746 
1747  q++;
1748  }
1749  }
1750  } // end if (compute_inner_products)
1751 
1753  }
1754 
1756 }
1757 
1759 {
1760  LOG_SCOPE("load_rb_solution()", "RBConstruction");
1761 
1762  solution->zero();
1763 
1764  if (get_rb_evaluation().RB_solution.size() > get_rb_evaluation().get_n_basis_functions())
1765  libmesh_error_msg("ERROR: System contains " << get_rb_evaluation().get_n_basis_functions() << " basis functions." \
1766  << " RB_solution vector constains " << get_rb_evaluation().RB_solution.size() << " entries." \
1767  << " RB_solution in RBConstruction::load_rb_solution is too long!");
1768 
1769  for (std::size_t i=0; i<get_rb_evaluation().RB_solution.size(); i++)
1770  solution->add(get_rb_evaluation().RB_solution(i), get_rb_evaluation().get_basis_function(i));
1771 
1772  update();
1773 }
1774 
1775 // The slow (but simple, non-error prone) way to compute the residual dual norm
1776 // Useful for error checking
1777 //Real RBConstruction::compute_residual_dual_norm(const unsigned int N)
1778 //{
1779 // LOG_SCOPE("compute_residual_dual_norm()", "RBConstruction");
1780 //
1781 // // Put the residual in rhs in order to compute the norm of the Riesz representor
1782 // // Note that this only works in serial since otherwise each processor will
1783 // // have a different parameter value during the Greedy training.
1784 //
1785 // UniquePtr<NumericVector<Number>> RB_sol = NumericVector<Number>::build();
1786 // RB_sol->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1787 //
1788 // UniquePtr<NumericVector<Number>> temp = NumericVector<Number>::build();
1789 // temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1790 //
1791 // for (unsigned int i=0; i<N; i++)
1792 // {
1793 // RB_sol->add(RB_solution(i), get_rb_evaluation().get_basis_function(i));
1794 // }
1795 //
1796 // this->truth_assembly();
1797 // matrix->vector_mult(*temp, *RB_sol);
1798 // rhs->add(-1., *temp);
1799 //
1800 // // Then solve to get the Reisz representor
1801 // matrix->zero();
1802 // matrix->add(1., *inner_product_matrix);
1803 // if (constrained_problem)
1804 // matrix->add(1., *constraint_matrix);
1805 //
1806 // solution->zero();
1807 // solve();
1808 // // Make sure we didn't max out the number of iterations
1809 // if ((this->n_linear_iterations() >=
1810 // this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations")) &&
1811 // (this->final_linear_residual() >
1812 // this->get_equation_systems().parameters.get<Real>("linear solver tolerance")))
1813 // {
1814 // libmesh_error_msg("Warning: Linear solver may not have converged! Final linear residual = "
1815 // << this->final_linear_residual() << ", number of iterations = "
1816 // << this->n_linear_iterations());
1817 // }
1818 //
1819 // inner_product_matrix->vector_mult(*inner_product_storage_vector, *solution);
1820 //
1821 // Real slow_residual_norm_sq = solution->dot(*inner_product_storage_vector);
1822 //
1823 // return std::sqrt( libmesh_real(slow_residual_norm_sq) );
1824 //}
1825 
1827 {
1828  return inner_product_matrix.get();
1829 }
1830 
1832 {
1833  if (q >= get_rb_theta_expansion().get_n_A_terms())
1834  libmesh_error_msg("Error: We must have q < Q_a in get_Aq.");
1835 
1836  return Aq_vector[q];
1837 }
1838 
1840 {
1842  libmesh_error_msg("Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Aq.");
1843 
1844  if (q >= get_rb_theta_expansion().get_n_A_terms())
1845  libmesh_error_msg("Error: We must have q < Q_a in get_Aq.");
1846 
1847  return non_dirichlet_Aq_vector[q];
1848 }
1849 
1851 {
1852  if (q >= get_rb_theta_expansion().get_n_F_terms())
1853  libmesh_error_msg("Error: We must have q < Q_f in get_Fq.");
1854 
1855  return Fq_vector[q];
1856 }
1857 
1859 {
1861  libmesh_error_msg("Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Fq.");
1862 
1863  if (q >= get_rb_theta_expansion().get_n_F_terms())
1864  libmesh_error_msg("Error: We must have q < Q_f in get_Fq.");
1865 
1866  return non_dirichlet_Fq_vector[q];
1867 }
1868 
1870 {
1871  if ((n >= get_rb_theta_expansion().get_n_outputs()) || (q_l >= get_rb_theta_expansion().get_n_output_terms(n)))
1872  libmesh_error_msg("Error: We must have n < n_outputs and " \
1873  << "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_output_vector.");
1874 
1875  return outputs_vector[n][q_l];
1876 }
1877 
1879 {
1880  if ((n >= get_rb_theta_expansion().get_n_outputs()) || (q_l >= get_rb_theta_expansion().get_n_output_terms(n)))
1881  libmesh_error_msg("Error: We must have n < n_outputs and " \
1882  << "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_non_dirichlet_output_vector.");
1883 
1884  return non_dirichlet_outputs_vector[n][q_l];
1885 }
1886 
1887 void RBConstruction::get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices)
1888 {
1889  all_matrices.clear();
1890 
1891  all_matrices["inner_product"] = get_inner_product_matrix();
1892 
1893  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1894  {
1895  std::stringstream matrix_name;
1896  matrix_name << "A" << q_a;
1897  all_matrices[matrix_name.str()] = get_Aq(q_a);
1898 
1900  {
1901  matrix_name << "_non_dirichlet";
1902  all_matrices[matrix_name.str()] = get_non_dirichlet_Aq(q_a);
1903  }
1904  }
1905 }
1906 
1907 void RBConstruction::get_all_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors)
1908 {
1909  all_vectors.clear();
1910 
1911  get_output_vectors(all_vectors);
1912 
1913  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1914  {
1915  std::stringstream F_vector_name;
1916  F_vector_name << "F" << q_f;
1917  all_vectors[F_vector_name.str()] = get_Fq(q_f);
1918 
1920  {
1921  F_vector_name << "_non_dirichlet";
1922  all_vectors[F_vector_name.str()] = get_non_dirichlet_Fq(q_f);
1923  }
1924  }
1925 }
1926 
1927 void RBConstruction::get_output_vectors(std::map<std::string, NumericVector<Number> *> & output_vectors)
1928 {
1929  output_vectors.clear();
1930 
1931  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1932  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1933  {
1934  std::stringstream output_name;
1935  output_name << "output_" << n << "_"<< q_l;
1936  output_vectors[output_name.str()] = get_output_vector(n,q_l);
1937 
1939  {
1940  output_name << "_non_dirichlet";
1941  output_vectors[output_name.str()] = get_non_dirichlet_output_vector(n,q_l);
1942  }
1943  }
1944 }
1945 
1947 {
1948  ZeroFunction<> zf;
1949 
1950  std::set<boundary_id_type> dirichlet_ids;
1951  std::vector<unsigned int> variables;
1952 
1953  // The DirichletBoundary constructor clones zf, so it's OK that zf is only in local scope
1954  return UniquePtr<DirichletBoundary> (new DirichletBoundary(dirichlet_ids, variables, &zf));
1955 }
1956 
1957 void RBConstruction::write_riesz_representors_to_files(const std::string & riesz_representors_dir,
1958  const bool write_binary_residual_representors)
1959 {
1960  LOG_SCOPE("write_riesz_representors_to_files()", "RBConstruction");
1961 
1962  // Write out Riesz representors. These are useful to have when restarting,
1963  // so you don't have to recompute them all over again.
1964 
1965  // First we write out the Fq representors, these are independent of an RBEvaluation object.
1966  libMesh::out << "Writing out the Fq_representors..." << std::endl;
1967 
1968  std::ostringstream file_name;
1969  const std::string riesz_representor_suffix = (write_binary_residual_representors ? ".xdr" : ".dat");
1970  struct stat stat_info;
1971 
1972  // Residual representors written out to their own separate directory
1973  if ( this->processor_id() == 0)
1974  if ( Utility::mkdir(riesz_representors_dir.c_str()) != 0)
1975  libMesh::out << "Skipping creating residual_representors directory: " << strerror(errno) << std::endl;
1976 
1977  for (std::size_t i=0; i<Fq_representor.size(); ++i)
1978  {
1979  if (Fq_representor[i] != libmesh_nullptr)
1980  {
1981  file_name.str(""); // reset filename
1982  file_name << riesz_representors_dir << "/Fq_representor" << i << riesz_representor_suffix;
1983 
1984  // Check to see if file exists, if so, don't overwrite it, we assume it was
1985  // there from a previous call to this function. Note: if stat returns zero
1986  // it means it successfully got the file attributes (and therefore the file
1987  // exists). Because of the following factors:
1988  // 1.) write_serialized_data takes longer for proc 0 than it does for others,
1989  // so processors can get out of sync.
1990  // 2.) The constructor for Xdr opens a (0 length) file on *all* processors,
1991  // there are typically hundreds of 0-length files created during this loop,
1992  // and that screws up checking for the existence of files. One way to stay
1993  // in sync is to but a barrier at each iteration of the loop -- not sure how
1994  // bad this will affect performance, but it can't be much worse than serialized
1995  // I/O already is :)
1996  int stat_result = stat(file_name.str().c_str(), &stat_info);
1997 
1998  if ( (stat_result != 0) || // file definitely doesn't already exist
1999  (stat_info.st_size == 0)) // file exists, but has zero length (can happen if another proc already opened it!)
2000  {
2001  // No need to copy!
2002  // *solution = *(Fq_representor[i]);
2003  // std::swap doesn't work on pointers
2004  //std::swap(solution.get(), Fq_representor[i]);
2005  Fq_representor[i]->swap(*solution);
2006 
2007  Xdr fqr_data(file_name.str(),
2008  write_binary_residual_representors ? ENCODE : WRITE);
2009 
2010  write_serialized_data(fqr_data, false);
2011 
2012  // Synchronize before moving on
2013  this->comm().barrier();
2014 
2015  // Swap back.
2016  Fq_representor[i]->swap(*solution);
2017 
2018  // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2019  // for the system call, be sure to do it only on one processor, etc.
2020  }
2021  }
2022  }
2023 
2024 
2025  // Next, write out the Aq representors associated with rb_eval.
2026  libMesh::out << "Writing out the Aq_representors..." << std::endl;
2027 
2028  const unsigned int jstop = get_rb_evaluation().get_n_basis_functions();
2029  const unsigned int jstart = jstop-get_delta_N();
2030  for (std::size_t i=0; i<get_rb_evaluation().Aq_representor.size(); ++i)
2031  for (unsigned int j=jstart; j<jstop; ++j)
2032  {
2033  libMesh::out << "Writing out Aq_representor[" << i << "][" << j << "]..." << std::endl;
2034  libmesh_assert(get_rb_evaluation().Aq_representor[i][j]);
2035 
2036  file_name.str(""); // reset filename
2037  file_name << riesz_representors_dir
2038  << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2039 
2040  {
2041  // No need to copy! Use swap instead.
2042  // *solution = *(Aq_representor[i][j]);
2043  get_rb_evaluation().Aq_representor[i][j]->swap(*solution);
2044 
2045  Xdr aqr_data(file_name.str(),
2046  write_binary_residual_representors ? ENCODE : WRITE);
2047 
2048  write_serialized_data(aqr_data, false);
2049 
2050  // Synchronize before moving on
2051  this->comm().barrier();
2052 
2053  // Swap back.
2054  get_rb_evaluation().Aq_representor[i][j]->swap(*solution);
2055 
2056  // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2057  // for the system call, be sure to do it only on one processor, etc.
2058  }
2059  }
2060 }
2061 
2062 
2063 
2064 void RBConstruction::read_riesz_representors_from_files(const std::string & riesz_representors_dir,
2065  const bool read_binary_residual_representors)
2066 {
2067  LOG_SCOPE("read_riesz_representors_from_files()", "RBConstruction");
2068 
2069  libMesh::out << "Reading in the Fq_representors..." << std::endl;
2070 
2071  const std::string riesz_representor_suffix = (read_binary_residual_representors ? ".xdr" : ".dat");
2072  std::ostringstream file_name;
2073  struct stat stat_info;
2074 
2075  // Read in the Fq_representors. There should be Q_f of these. FIXME:
2076  // should we be worried about leaks here?
2077  for (std::size_t i=0; i<Fq_representor.size(); ++i)
2078  {
2079  if (Fq_representor[i] != libmesh_nullptr)
2080  libmesh_error_msg("Error, must delete existing Fq_representor before reading in from file.");
2081  }
2082 
2083  for (std::size_t i=0; i<Fq_representor.size(); i++)
2084  {
2085  file_name.str(""); // reset filename
2086  file_name << riesz_representors_dir
2087  << "/Fq_representor" << i << riesz_representor_suffix;
2088 
2089  // On processor zero check to be sure the file exists
2090  if (this->processor_id() == 0)
2091  {
2092  int stat_result = stat(file_name.str().c_str(), &stat_info);
2093 
2094  if (stat_result != 0)
2095  libmesh_error_msg("File does not exist: " << file_name.str());
2096  }
2097 
2098  Xdr fqr_data(file_name.str(),
2099  read_binary_residual_representors ? DECODE : READ);
2100 
2101  read_serialized_data(fqr_data, false);
2102 
2103  Fq_representor[i] = NumericVector<Number>::build(this->comm()).release();
2104  Fq_representor[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2105 
2106  // No need to copy, just swap
2107  // *Fq_representor[i] = *solution;
2108  Fq_representor[i]->swap(*solution);
2109  }
2110 
2111  // Alert the update_residual_terms() function that we don't need to recompute
2112  // the Fq_representors as we have already read them in from file!
2114 
2115 
2116  libMesh::out << "Reading in the Aq_representors..." << std::endl;
2117 
2118  // Read in the Aq representors. The class makes room for [Q_a][Nmax] of these. We are going to
2119  // read in [Q_a][get_rb_evaluation().get_n_basis_functions()]. FIXME:
2120  // should we be worried about leaks in the locations where we're about to fill entries?
2121  for (std::size_t i=0; i<get_rb_evaluation().Aq_representor.size(); ++i)
2122  for (std::size_t j=0; j<get_rb_evaluation().Aq_representor[i].size(); ++j)
2123  {
2125  libmesh_error_msg("Error, must delete existing Aq_representor before reading in from file.");
2126  }
2127 
2128  // Now ready to read them in from file!
2129  for (std::size_t i=0; i<get_rb_evaluation().Aq_representor.size(); ++i)
2130  for (unsigned int j=0; j<get_rb_evaluation().get_n_basis_functions(); ++j)
2131  {
2132  file_name.str(""); // reset filename
2133  file_name << riesz_representors_dir
2134  << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2135 
2136  // On processor zero check to be sure the file exists
2137  if (this->processor_id() == 0)
2138  {
2139  int stat_result = stat(file_name.str().c_str(), &stat_info);
2140 
2141  if (stat_result != 0)
2142  libmesh_error_msg("File does not exist: " << file_name.str());
2143  }
2144 
2145  Xdr aqr_data(file_name.str(), read_binary_residual_representors ? DECODE : READ);
2146 
2147  read_serialized_data(aqr_data, false);
2148 
2151  false, PARALLEL);
2152 
2153  // No need to copy, just swap
2154  //*Aq_representor[i][j] = *solution;
2155  get_rb_evaluation().Aq_representor[i][j]->swap(*solution);
2156  }
2157 }
2158 
2160 {
2162 
2163  conv_flag = input_solver.get_converged_reason();
2164 
2165  if (conv_flag < 0)
2166  {
2167  std::stringstream err_msg;
2168  err_msg << "Convergence error. Error id: " << conv_flag;
2169  libmesh_error_msg(err_msg.str());
2170  }
2171 }
2172 
2174 {
2175  return assert_convergence;
2176 }
2177 
2179 {
2180  assert_convergence = flag;
2181 }
2182 
2183 } // namespace libMesh
virtual void read_riesz_representors_from_files(const std::string &riesz_representors_dir, const bool write_binary_residual_representors)
Read in all the Riesz representor data from files.
T libmesh_real(T a)
OStreamProxy err
unsigned int n_linear_iterations() const
std::vector< NumericVector< Number > * > basis_functions
The libMesh vectors storing the finite element coefficients of the RB basis functions.
numeric_index_type get_local_n_training_samples() const
Get the total number of training samples local to this processor.
bool normalize_rb_bound_in_greedy
This boolean indicates if we normalize the RB error in the greedy using RBEvaluation::get_error_bound...
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.
const RBParameters & get_parameters_max() const
Get an RBParameters object that specifies the maximum allowable value for each parameter.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
std::vector< SparseMatrix< Number > * > non_dirichlet_Aq_vector
We may also need a second set of matrices/vectors that do not have the Dirichlet boundary conditions ...
unsigned char get_side() const
Accessor for current side of Elem object.
Definition: fem_context.h:885
This is the EquationSystems class.
std::vector< Real > training_error_bounds
Vector storing the values of the error bound for each parameter in the training set — the parameter ...
virtual void get_all_vectors(std::map< std::string, NumericVector< Number > * > &all_vectors)
Get a map that stores pointers to all of the vectors.
A Node is like a Point, but with more information.
Definition: node.h:52
virtual void truth_assembly()
Assemble the truth matrix and right-hand side for current_parameters.
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
unsigned int n() const
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
UniquePtr< NumericVector< Number > > inner_product_storage_vector
We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary me...
T libmesh_conj(T a)
ElemAssembly & get_F_assembly(unsigned int q)
Return a reference to the specified F_assembly object.
bool assert_convergence
A boolean flag to indicate whether to check for proper convergence after each solve.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true)
Train the reduced basis.
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true)
Resize and clear the data vectors corresponding to the value of Nmax.
virtual ~RBConstruction()
Destructor.
void write_serialized_data(Xdr &io, const bool write_additional_data=true) const
Writes additional data, namely vectors, for this System.
Definition: system_io.C:1728
void broadcast_parameters(unsigned int proc_id)
Broadcasts parameters on processor proc_id to all processors.
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.
virtual Real truth_solve(int plot_solution)
Perform a "truth" solve, i.e.
RBConstruction(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
unsigned int _n_linear_iterations
The number of linear iterations required to solve the linear system Ax=b.
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
std::string get_timestamp()
Definition: timestamp.C:37
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
std::vector< std::vector< NumericVector< Number > * > > Aq_representor
Vector storing the residual representors associated with the left-hand side.
DenseMatrix< Number > RB_inner_product_matrix
The inner product matrix.
bool get_convergence_assertion_flag() const
Getter for the flag determining if convergence should be checked after each solve.
void print_discrete_parameter_values() const
Print out all the discrete parameter values.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:420
SparseMatrix< Number > * get_Aq(unsigned int q)
Get a pointer to Aq.
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
virtual void set_context_solution_vec(NumericVector< Number > &vec)
Set current_local_solution = vec so that we can access vec from FEMContext during assembly...
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]...
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
virtual void interior_assembly(FEMContext &)
Perform the element interior assembly.
Definition: elem_assembly.h:59
unsigned int get_n_params() const
Get the number of parameters.
void attach_matrix(SparseMatrix< Number > &matrix)
Additional matrices may be handled with this DofMap.
Definition: dof_map.C:246
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
MeshBase & mesh
ElemAssembly & get_inner_product_assembly()
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
This proxy class acts like a container of indices from a single coupling row.
bool is_discrete_parameter(const std::string &mu_name) const
Is parameter mu_name discrete?
const T & get(const std::string &) const
Definition: parameters.h:430
virtual bool greedy_termination_test(Real abs_greedy_error, Real initial_greedy_error, int count)
Function that indicates when to terminate the Greedy basis training.
void get_transpose(DenseMatrix< T > &dest) const
Put the tranposed matrix into dest.
Definition: dense_matrix.C:576
ElemAssembly & get_output_assembly(unsigned int output_index, unsigned int q_l)
Return a reference to the specified output assembly object.
void add_scaled_Aq(Number scalar, unsigned int q_a, SparseMatrix< Number > *input_matrix, bool symmetrize)
Add the scaled q^th affine matrix to input_matrix.
const RBParameters & get_greedy_parameter(unsigned int i)
Return the parameters chosen during the i^th step of the Greedy algorithm.
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
virtual LinearConvergenceReason get_converged_reason() const =0
int mkdir(const char *pathname)
Create a directory.
Definition: utility.C:140
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
static void get_global_max_error_pair(const Parallel::Communicator &communicator, std::pair< numeric_index_type, Real > &error_pair)
Static function to return the error pair (index,error) that is corresponds to the largest error on al...
static UniquePtr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
virtual Real get_error_bound_normalization()
const_iterator end() const
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
NumericVector< Number > * get_Fq(unsigned int q)
Get a pointer to Fq.
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
std::vector< RBParameters > greedy_param_list
The list of parameters selected by the Greedy algorithm in generating the Reduced Basis associated wi...
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
FEType get_fe_type() const
Definition: fe_abstract.h:421
virtual std::string system_type() const libmesh_override
virtual void get_nodal_values(std::vector< dof_id_type > &, DenseMatrix< Number > &, DenseVector< Number > &, const System &, const Node &)
Get values to add to the matrix or rhs vector based on node.
Definition: elem_assembly.h:71
long double max(long double a, double b)
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
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...
const std::string & name() const
Definition: system.h:1998
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
virtual void init_context(FEMContext &)
Initialize the FEMContext prior to performing an element loop.
virtual Real compute_max_error_bound()
(i) Compute the a posteriori error bound for each set of parameters in the training set...
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
numeric_index_type get_n_training_samples() const
Get the total number of training samples.
unsigned int get_Nmax() const
Get/set Nmax, the maximum number of RB functions we are willing to compute.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void set_quiet_mode(bool quiet_mode_in)
Set the quiet_mode flag.
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...
void assemble_inner_product_matrix(SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
Assemble the inner product matrix and store it in input_matrix.
virtual void get_output_vectors(std::map< std::string, NumericVector< Number > * > &all_vectors)
Get a map that stores pointers to all of the vectors.
LinearSolver< Number > * extra_linear_solver
Also, we store a pointer to an extra linear solver.
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
virtual void boundary_assembly(FEMContext &)
Perform the element boundary assembly.
Definition: elem_assembly.h:64
virtual void get_nodal_rhs_values(std::map< numeric_index_type, Number > &values, const System &, const Node &)
Get values to add to the RHS vector based on node.
Definition: elem_assembly.h:88
virtual Number eval_output_theta(unsigned int output_index, unsigned int q_l, const RBParameters &mu)
Evaluate theta_q_l at the current parameter.
Generic sparse matrix.
Definition: dof_map.h:65
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual void scale(const T factor)=0
Scale each element of the vector by the given factor.
virtual void recompute_all_residual_terms(const bool compute_inner_products=true)
This function computes all of the residual representors, can be useful when restarting a basis traini...
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
virtual void load_basis_function(unsigned int i)
Load the i^th RB function into the RBConstruction solution vector.
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:58
std::vector< DenseVector< Number > > RB_Fq_vector
Dense vector for the RHS.
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
bool is_rb_eval_initialized() const
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 impose_internal_fluxes
Boolean flag to indicate whether we impose "fluxes" (i.e.
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
bool Fq_representor_innerprods_computed
A boolean flag to indicate whether or not the Fq representor norms have already been computed — used...
bool use_empty_rb_solve_in_greedy
A boolean flag to indicate whether or not we initialize the Greedy algorithm by performing rb_solves ...
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1734
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
std::vector< std::vector< NumericVector< Number > * > > non_dirichlet_outputs_vector
bool dg_terms_are_active() const
Are the DG terms active, i.e.
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.
Real get_parameter_max(const std::string &param_name) const
Get maximum allowable value of parameter param_name.
virtual void get_all_matrices(std::map< std::string, SparseMatrix< Number > * > &all_matrices)
Get a map that stores pointers to all of the matrices.
virtual void zero()=0
Set all entries to 0.
dof_id_type numeric_index_type
Definition: id_types.h:92
void update_greedy_param_list()
Update the list of Greedily chosen parameters with current_parameters.
void set_parameters(const RBParameters &params)
Set the current parameters to params.
Real rel_training_tolerance
Relative and absolute tolerances for training reduced basis using the Greedy scheme.
virtual void init(const char *name=libmesh_nullptr)=0
Initialize data structures if not done so already.
unsigned int m() const
void print_basis_function_orthogonality()
Print out a matrix that shows the orthogonality of the RB basis functions.
void set_params_from_training_set(unsigned int index)
Set parameters to the RBParameters stored in index index of the training set.
RBAssemblyExpansion & get_rb_assembly_expansion()
void set_rb_construction_parameters(unsigned int n_training_samples_in, bool deterministic_training_in, unsigned int training_parameters_random_seed_in, bool quiet_mode_in, unsigned int Nmax_in, Real rel_training_tolerance_in, Real abs_training_tolerance_in, bool normalize_rb_error_bound_in_greedy_in, RBParameters mu_min_in, RBParameters mu_max_in, std::map< std::string, std::vector< Real >> discrete_parameter_values_in, std::map< std::string, bool > log_scaling)
Set the state of this RBConstruction object based on the arguments to this function.
SparseMatrix< Number > * get_inner_product_matrix()
Get a pointer to inner_product_matrix.
virtual void update_RB_system_matrices()
Compute the reduced basis matrices for the current basis.
NumericVector< Number > & get_basis_function(unsigned int i)
Get a reference to the i^th basis function.
virtual void update_residual_terms(bool compute_inner_products=true)
Compute the terms that are combined `online&#39; to determine the dual norm of the residual.
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
void assemble_Fq_vector(unsigned int q, NumericVector< Number > *input_vector, bool apply_dof_constraints=true)
Assemble the q^th affine vector and store it in input_matrix.
virtual void update_system()
Update the system after enriching the RB space; this calls a series of functions to update the system...
This class extends FEMContext in order to provide extra data required to perform local element residu...
const_iterator end() const
Get a constant iterator to the end of this RBParameters object.
Definition: rb_parameters.C:85
void barrier() const
Pause execution until all processors reach a certain point.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
std::vector< std::vector< NumericVector< Number > * > > outputs_vector
The libMesh vectors that define the output functionals.
virtual SparseMatrix< Number > & get_matrix_for_output_dual_solves()
Return the matrix for the output residual dual norm solves.
virtual void print_info()
Print out info that describes the current setup of this RBConstruction.
NumericVector< Number > * get_non_dirichlet_Fq(unsigned int q)
Get a pointer to non-Dirichlet Fq.
virtual T el(const unsigned int i, const unsigned int j) const libmesh_override
Definition: dense_matrix.h:84
virtual Real get_RB_error_bound()
virtual void initialize_training_parameters(const RBParameters &mu_min, const RBParameters &mu_max, unsigned int n_training_parameters, std::map< std::string, bool > log_param_scale, bool deterministic=true)
Initialize the parameter ranges and indicate whether deterministic or random training parameters shou...
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0
This function calls the solver _solver_type preconditioned with the _preconditioner_type precondition...
std::vector< NumericVector< Number > * > Fq_representor
Vector storing the residual representors associated with the right-hand side.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
CouplingMatrix * _dof_coupling
Degree of freedom coupling.
Definition: dof_map.h:1240
SparseMatrix< Number > * get_non_dirichlet_Aq(unsigned int q)
Get a pointer to non_dirichlet_Aq.
void set_normalize_rb_bound_in_greedy(bool normalize_rb_bound_in_greedy)
Get/set the boolean to indicate if we normalize the RB error in the greedy.
void set_inner_product_assembly(ElemAssembly &inner_product_assembly_in)
Set the rb_assembly_expansion object.
virtual unsigned int n_sides() const =0
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
void set_convergence_assertion_flag(bool flag)
Setter for the flag determining if convergence should be checked after each solve.
const DenseMatrix< Number > & get_elem_neighbor_jacobian() const
Const accessor for element-neighbor Jacobian.
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
unsigned int delta_N
The number of basis functions that we add at each greedy step.
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
void build_node_list(std::vector< dof_id_type > &node_id_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of nodes and ids for those nodes.
RBEvaluation * rb_eval
The current RBEvaluation object we are using to perform the Evaluation stage of the reduced basis met...
const_iterator begin() const
const DenseMatrix< Number > & get_neighbor_elem_jacobian() const
Const accessor for element-neighbor Jacobian.
const DenseMatrix< Number > & get_elem_elem_jacobian() const
Const accessor for element-element Jacobian.
const std::vector< dof_id_type > & get_neighbor_dof_indices() const
Accessor for neighbor dof indices.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
const RBParameters & get_parameters_min() const
Get an RBParameters object that specifies the minimum allowable value for each parameter.
Real _final_linear_residual
The final residual for the linear system Ax=b.
std::vector< std::vector< std::vector< Number > > > Fq_Aq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
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
unsigned int get_n_outputs() const
Get n_outputs, the number output functionals.
This class is part of the rbOOmit framework.
Definition: rb_evaluation.h:50
UniquePtr< LinearSolver< Number > > linear_solver
The LinearSolver defines the interface used to solve the linear_implicit system.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
void set_abs_training_tolerance(Real new_training_tolerance)
Get/set the absolute tolerance for the basis training.
std::vector< SparseMatrix< Number > * > Aq_vector
Vector storing the Q_a matrices from the affine expansion.
std::vector< DenseMatrix< Number > > RB_Aq_vector
Dense matrices for the RB computations.
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool is_quiet() const
Is the system in quiet mode?
numeric_index_type get_first_local_training_index() const
Get the first local index of the training parameters.
void set_training_random_seed(unsigned int seed)
Set the seed that is used to randomly generate training parameters.
ElemAssembly * inner_product_assembly
Pointer to inner product assembly.
virtual void reuse_preconditioner(bool)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
bool is_nodal_rhs_values_overriden
Temporary flag to help us figure out if we should call the deprecated get_nodal_rhs_values method or ...
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
void check_convergence(LinearSolver< Number > &input_solver)
Check if the linear solver reports convergence.
static UniquePtr< DirichletBoundary > build_zero_dirichlet_boundary_object()
It&#39;s helpful to be able to generate a DirichletBoundary that stores a ZeroFunction in order to impose...
unsigned int get_delta_N() const
Get delta_N, the number of basis functions we add to the RB space per iteration of the greedy algorit...
void set_value(const std::string &param_name, Real value)
Set the value of the specified parameter.
Definition: rb_parameters.C:57
void print_parameters() const
Print the current parameters.
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
This class stores the set of ElemAssembly functor objects that define the "parameter-independent expa...
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
ElemAssembly & get_A_assembly(unsigned int q)
Return a reference to the specified A_assembly object.
ElemAssembly provides a per-element (interior and boundary) assembly functionality.
Definition: elem_assembly.h:38
std::vector< Number > truth_outputs
Vector storing the truth output values from the most recent truth solve.
static const bool value
Definition: xdr_io.C:108
virtual void compute_Fq_representor_innerprods(bool compute_inner_products=true)
Compute the terms that are combined `online&#39; to determine the dual norm of the residual.
virtual void set_Nmax(unsigned int Nmax)
void set_rel_training_tolerance(Real new_training_tolerance)
Get/set the relative tolerance for the basis training.
virtual void write_riesz_representors_to_files(const std::string &riesz_representors_dir, const bool write_binary_residual_representors)
Write out all the Riesz representor data to files.
numeric_index_type get_last_local_training_index() const
Get the last local index of the training parameters.
virtual void assemble_all_output_vectors()
Assemble and store the output vectors.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
bool serial_training_set
This boolean flag indicates whether or not the training set should be the same on all processors...
void assemble_Aq_matrix(unsigned int q, SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
Assemble the q^th affine matrix and store it in input_matrix.
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags)
virtual T dot(const NumericVector< T > &v) const =0
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...
RBAssemblyExpansion * rb_assembly_expansion
This member holds the (parameter independent) assembly functors that define the "affine expansion" of...
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
std::size_t n_nodeset_conds() const
NumericVector< Number > * get_output_vector(unsigned int n, unsigned int q_l)
Get a pointer to the n^th output.
dof_id_type n_dofs() const
Definition: system.C:148
unsigned int rank() const
Definition: parallel.h:724
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
unsigned int n_vars() const
Definition: system.h:2086
bool quiet_mode
Flag to indicate whether we print out extra information during the Offline stage. ...
virtual void assemble_misc_matrices()
Assemble and store all the inner-product matrix, the constraint matrix (for constrained problems) and...
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
const_iterator begin() const
Get a constant iterator to beginning of this RBParameters object.
Definition: rb_parameters.C:80
const RBParameters & get_parameters() const
Get the current parameters.
virtual void load_rb_solution()
Load the RB solution from the most recent solve with rb_eval into this system&#39;s solution vector...
void set_rb_evaluation(RBEvaluation &rb_eval_in)
Set the RBEvaluation object.
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
virtual void compute_output_dual_innerprods()
Compute and store the dual norm of each output functional.
dof_id_type id() const
Definition: dof_object.h:632
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
virtual void side_fe_reinit() libmesh_override
Override side_fe_reinit to set a boolean flag so that by default DG terms are assumed to be inactive...
std::map< std::string, Real >::const_iterator const_iterator
Definition: rb_parameters.h:57
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.
bool exit_on_repeated_greedy_parameters
Boolean flag to indicate whether we exit the greedy if we select the same parameters twice in a row...
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
std::vector< NumericVector< Number > * > non_dirichlet_Fq_vector
virtual UniquePtr< DGFEMContext > build_context()
Builds a DGFEMContext object with enough information to do evaluations on each element.
bool output_dual_innerprods_computed
A boolean flag to indicate whether or not the output dual norms have already been computed — used to...
std::vector< NumericVector< Number > * > Fq_vector
Vector storing the Q_f vectors in the affine decomposition of the right-hand side.
virtual void enrich_RB_space()
Add a new basis function to the RB space.
virtual void assemble_all_affine_vectors()
Assemble and store the affine RHS vectors.
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
void read_serialized_data(Xdr &io, const bool read_additional_data=true)
Reads additional data, namely vectors, for this System.
Definition: system_io.C:729
virtual void attach_quadrature_rule(QBase *q)=0
Provides the class with the quadrature rule.
NumericVector< Number > * get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l)
Get a pointer to non-Dirichlet output vector.
virtual void allocate_data_structures()
Helper function that actually allocates all the data structures required by this class.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the rb_theta_expansion.
Definition: rb_evaluation.C:89
const DenseMatrix< Number > & get_neighbor_neighbor_jacobian() const
Const accessor for element-neighbor Jacobian.
This class forms the foundation from which generic finite elements may be derived.
virtual void clear()
Clear all the data structures associated with the system.
unsigned int Nmax
Maximum number of reduced basis functions we are willing to use.
virtual void assemble_affine_expansion(bool skip_matrix_assembly, bool skip_vector_assembly)
Assemble the matrices and vectors for this system.
std::vector< std::vector< std::vector< Number > > > Aq_Aq_representor_innerprods
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694
This class defines a coupling matrix.
Real get_parameter_min(const std::string &param_name) const
Get minimum allowable value of parameter param_name.
void zero_constrained_dofs_on_vector(NumericVector< Number > &vector)
It is sometimes useful to be able to zero vector entries that correspond to constrained dofs...
virtual void assemble_all_affine_operators()
Assemble and store all Q_a affine operators as well as the inner-product matrix.
std::vector< std::vector< DenseVector< Number > > > RB_output_vectors
The vectors storing the RB output vectors.