libMesh
transient_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/transient_rb_construction.h"
22 #include "libmesh/transient_rb_evaluation.h"
23 #include "libmesh/transient_rb_theta_expansion.h"
24 #include "libmesh/transient_rb_assembly_expansion.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/linear_solver.h"
32 #include "libmesh/equation_systems.h"
33 #include "libmesh/exodusII_io.h"
34 #include "libmesh/getpot.h"
35 #include "libmesh/dense_matrix.h"
36 #include "libmesh/dense_vector.h"
37 #include "libmesh/xdr_cxx.h"
38 #include "libmesh/timestamp.h"
39 #include "libmesh/parallel.h"
40 
41 // For checking for the existence of files
42 #include <sys/stat.h>
43 
44 #include <fstream>
45 #include <sstream>
46 
47 // Need SLEPc to get the POD eigenvalues
48 #if defined(LIBMESH_HAVE_SLEPC)
49 // LAPACK include (via SLEPc)
50 #include <petscsys.h>
51 #include <slepcblaslapack.h>
52 #endif // LIBMESH_HAVE_SLEPC
53 
54 namespace libMesh
55 {
56 
58  const std::string & name_in,
59  const unsigned int number_in)
60  : Parent(es, name_in, number_in),
61  L2_matrix(SparseMatrix<Number>::build(es.comm())),
62  non_dirichlet_L2_matrix(SparseMatrix<Number>::build(es.comm())),
63  nonzero_initialization(false),
64  compute_truth_projection_error(false),
65  init_filename(""),
66  POD_tol(-1.),
67  max_truth_solves(-1),
68  L2_assembly(libmesh_nullptr)
69 {
70  // Indicate that we need to compute the RB
71  // inner product matrix in this case
73 
74  temporal_data.resize(0);
75 
76  // We should not necessarily exit the greedy due to repeated parameters in
77  // the transient case
79 }
80 
81 
82 
84 {
85  this->clear();
86 }
87 
88 
90 {
91  Parent::clear();
92 
93  // clear the mass matrices
94  for (std::size_t q=0; q<M_q_vector.size(); q++)
95  {
96  delete M_q_vector[q];
98  }
99 
101  {
102  for (std::size_t q=0; q<non_dirichlet_M_q_vector.size(); q++)
103  {
104  delete non_dirichlet_M_q_vector[q];
106  }
107  }
108 
109  // clear the temporal_data
110  for (std::size_t i=0; i<temporal_data.size(); i++)
111  {
112  if (temporal_data[i])
113  {
114  temporal_data[i]->clear();
115  delete temporal_data[i];
117  }
118  }
119  temporal_data.resize(0);
120 }
121 
123  bool skip_vector_assembly)
124 {
125  // Check that the theta and assembly objects are consistently sized
126 #ifndef NDEBUG
127  TransientRBThetaExpansion & trans_theta_expansion =
128  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
129 
130  TransientRBAssemblyExpansion & trans_assembly_expansion =
131  cast_ref<TransientRBAssemblyExpansion &>(get_rb_assembly_expansion());
132 #endif
133  // This assert only gets called if DEBUG is on
134  libmesh_assert_equal_to (trans_theta_expansion.get_n_M_terms(), trans_assembly_expansion.get_n_M_terms());
135 
136  Parent::initialize_rb_construction(skip_matrix_assembly, skip_vector_assembly);
137 }
138 
139 void TransientRBConstruction::process_parameters_file (const std::string & parameters_filename)
140 {
141  Parent::process_parameters_file(parameters_filename);
142 
143  // Read in data from parameters_filename
144  GetPot infile(parameters_filename);
145 
146  // Read in the generic temporal discretization data
147  process_temporal_parameters_file(parameters_filename);
148 
149  // Read in the data specific to Construction
150  nonzero_initialization = infile("nonzero_initialization",nonzero_initialization);
151  init_filename = infile("init_filename",init_filename);
152 
153  const Real POD_tol_in = infile("POD_tol", POD_tol);
154  const int max_truth_solves_in = infile("max_truth_solves", max_truth_solves);
155  const unsigned int delta_N_in = infile("delta_N", delta_N);
156 
157  set_POD_tol(POD_tol_in);
158  set_max_truth_solves(max_truth_solves_in);
159  set_delta_N(delta_N_in);
160 
161  // Pass the temporal discretization data to the RBEvaluation
162  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
163  trans_rb_eval.pull_temporal_discretization_data( *this );
164 }
165 
167 {
169 
170  libMesh::out << std::endl << "TransientRBConstruction parameters:" << std::endl;
171 
173  {
174  // Print out info that describes the current setup
175  TransientRBThetaExpansion & trans_theta_expansion =
176  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
177  libMesh::out << "Q_m: " << trans_theta_expansion.get_n_M_terms() << std::endl;
178  }
179  else
180  {
181  libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
182  }
183  libMesh::out << "Number of time-steps: " << get_n_time_steps() << std::endl;
184  libMesh::out << "dt: " << get_delta_t() << std::endl;
185  libMesh::out << "euler_theta (time discretization parameter): " << get_euler_theta() << std::endl;
186  if (get_POD_tol() > 0.)
187  libMesh::out << "POD_tol: " << get_POD_tol() << std::endl;
188  if (max_truth_solves > 0)
189  libMesh::out << "Maximum number of truth solves: " << max_truth_solves << std::endl;
190  libMesh::out << "delta_N (number of basis functions to add each POD-Greedy step): " << get_delta_N() << std::endl;
192  {
193  libMesh::out << "Reading initial condition from " << init_filename << std::endl;
194  }
195  else
196  {
197  libMesh::out << "Using zero initial condition" << std::endl;
198  }
199  libMesh::out << std::endl;
200 }
201 
203 {
205 
206  TransientRBThetaExpansion & trans_theta_expansion =
207  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
208  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
209  const unsigned int n_outputs = trans_theta_expansion.get_n_outputs();
210 
211  // Resize and allocate vectors for storing mesh-dependent data
212  const unsigned int n_time_levels = get_n_time_steps()+1;
213  temporal_data.resize(n_time_levels);
214 
215  // Resize vectors for storing mesh-dependent data but only
216  // initialize if initialize_mesh_dependent_data == true
217  M_q_vector.resize(Q_m);
218 
219  // Only initialize the mass matrices if we
220  // are not in single-matrix mode
221  {
222  DofMap & dof_map = this->get_dof_map();
223 
224  dof_map.attach_matrix(*L2_matrix);
225  L2_matrix->init();
226  L2_matrix->zero();
227 
228  for (unsigned int q=0; q<Q_m; q++)
229  {
230  // Initialize the memory for the matrices
231  M_q_vector[q] = SparseMatrix<Number>::build(this->comm()).release();
232  dof_map.attach_matrix(*M_q_vector[q]);
233  M_q_vector[q]->init();
234  M_q_vector[q]->zero();
235  }
236 
237  // We also need to initialize a second set of non-Dirichlet operators
239  {
241  non_dirichlet_L2_matrix->init();
242  non_dirichlet_L2_matrix->zero();
243 
244  non_dirichlet_M_q_vector.resize(Q_m);
245  for (unsigned int q=0; q<Q_m; q++)
246  {
247  // Initialize the memory for the matrices
250  non_dirichlet_M_q_vector[q]->init();
251  non_dirichlet_M_q_vector[q]->zero();
252  }
253  }
254  }
255 
256  for (unsigned int i=0; i<n_time_levels; i++)
257  {
258  temporal_data[i] = (NumericVector<Number>::build(this->comm()).release());
259  temporal_data[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
260  }
261 
262  // and the truth output vectors
263  truth_outputs_all_k.resize(n_outputs);
264  for (unsigned int n=0; n<n_outputs; n++)
265  {
266  truth_outputs_all_k[n].resize(n_time_levels);
267  }
268 
269  // This vector is for storing rhs entries for
270  // computing the projection of the initial condition
271  // into the RB space
273 }
274 
276  bool skip_vector_assembly)
277 {
278  // Call parent's assembly functions
279  Parent::assemble_affine_expansion(skip_matrix_assembly, skip_vector_assembly);
280 
281  // Now update RB_ic_proj_rhs_all_N if necessary.
282  // This allows us to compute the L2 projection
283  // of the initial condition into the RB space
284  // so that we can continue to enrich a given RB
285  // space.
286  if (get_rb_evaluation().get_n_basis_functions() > 0)
287  {
288  // Load the initial condition into the solution vector
290 
292  temp1->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
293 
294  // First compute the right-hand side vector for the L2 projection
295  L2_matrix->vector_mult(*temp1, *solution);
296 
297  for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
298  {
299  RB_ic_proj_rhs_all_N(i) = temp1->dot(get_rb_evaluation().get_basis_function(i));
300  }
301  }
302 }
303 
304 Real TransientRBConstruction::train_reduced_basis(const bool resize_rb_eval_data)
305 {
307  Real value = Parent::train_reduced_basis(resize_rb_eval_data);
309 
310  return value;
311 }
312 
314 {
315  TransientRBThetaExpansion & trans_theta_expansion =
316  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
317 
318  if (q >= trans_theta_expansion.get_n_M_terms())
319  libmesh_error_msg("Error: We must have q < Q_m in get_M_q.");
320 
321  return M_q_vector[q];
322 }
323 
325 {
327  libmesh_error_msg("Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_M_q.");
328 
329  TransientRBThetaExpansion & trans_theta_expansion =
330  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
331 
332  if (q >= trans_theta_expansion.get_n_M_terms())
333  libmesh_error_msg("Error: We must have q < Q_m in get_M_q.");
334 
335  return non_dirichlet_M_q_vector[q];
336 }
337 
338 void TransientRBConstruction::get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices)
339 {
340  Parent::get_all_matrices(all_matrices);
341 
342  all_matrices["L2_matrix"] = L2_matrix.get();
343 
345  all_matrices["L2_matrix_non_dirichlet"] = non_dirichlet_L2_matrix.get();
346 
347  TransientRBThetaExpansion & trans_theta_expansion =
348  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
349  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
350 
351  for (unsigned int q_m=0; q_m<Q_m; q_m++)
352  {
353  std::stringstream matrix_name;
354  matrix_name << "M" << q_m;
355  all_matrices[matrix_name.str()] = get_M_q(q_m);
356 
358  {
359  matrix_name << "_non_dirichlet";
360  all_matrices[matrix_name.str()] = get_non_dirichlet_M_q(q_m);
361  }
362  }
363 }
364 
365 void TransientRBConstruction::assemble_L2_matrix(SparseMatrix<Number> * input_matrix, bool apply_dirichlet_bc)
366 {
367  input_matrix->zero();
369  L2_assembly,
370  input_matrix,
372  false, /* symmetrize */
373  apply_dirichlet_bc);
374 }
375 
377 {
378  input_matrix->zero();
379  add_scaled_mass_matrix(1., input_matrix);
380 }
381 
383 {
384  const RBParameters & mu = get_parameters();
385 
386  TransientRBThetaExpansion & trans_theta_expansion =
387  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
388 
389  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
390 
391  for (unsigned int q=0; q<Q_m; q++)
392  input_matrix->add(scalar * trans_theta_expansion.eval_M_theta(q,mu), *get_M_q(q));
393 }
394 
396  NumericVector<Number> & dest,
397  NumericVector<Number> & arg)
398 {
399  LOG_SCOPE("mass_matrix_scaled_matvec()", "TransientRBConstruction");
400 
401  dest.zero();
402 
403  const RBParameters & mu = get_parameters();
404 
405  TransientRBThetaExpansion & trans_theta_expansion =
406  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
407 
408  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
409 
411  temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
412 
413  for (unsigned int q=0; q<Q_m; q++)
414  {
415  get_M_q(q)->vector_mult(*temp_vec, arg);
416  dest.add(scalar * trans_theta_expansion.eval_M_theta(q,mu), *temp_vec);
417  }
418 }
419 
421 {
422  LOG_SCOPE("truth_assembly()", "TransientRBConstruction");
423 
424  this->matrix->close();
425 
426  this->matrix->zero();
427  this->rhs->zero();
428 
429  const RBParameters & mu = get_parameters();
430 
431  TransientRBThetaExpansion & trans_theta_expansion =
432  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
433 
434  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
435  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
436 
437  const Real dt = get_delta_t();
438  const Real euler_theta = get_euler_theta();
439 
440  {
441  // We should have already assembled the matrices
442  // and vectors in the affine expansion, so
443  // just use them
444 
447 
449  temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
450 
451  for (unsigned int q_a=0; q_a<Q_a; q_a++)
452  {
453  matrix->add(euler_theta*trans_theta_expansion.eval_A_theta(q_a,mu), *get_Aq(q_a));
454 
455  get_Aq(q_a)->vector_mult(*temp_vec, *current_local_solution);
456  temp_vec->scale( -(1.-euler_theta)*trans_theta_expansion.eval_A_theta(q_a,mu) );
457  rhs->add(*temp_vec);
458  }
459 
460  for (unsigned int q_f=0; q_f<Q_f; q_f++)
461  {
462  *temp_vec = *get_Fq(q_f);
463  temp_vec->scale( get_control(get_time_step())*trans_theta_expansion.eval_F_theta(q_f,mu) );
464  rhs->add(*temp_vec);
465  }
466 
467  }
468 
469  this->matrix->close();
470  this->rhs->close();
471 }
472 
474 {
475  L2_assembly = &L2_assembly_in;
476 }
477 
479 {
480  if (!L2_assembly)
481  libmesh_error_msg("Error: L2_assembly hasn't been initialized yet");
482 
483  return *L2_assembly;
484 }
485 
486 void TransientRBConstruction::assemble_Mq_matrix(unsigned int q, SparseMatrix<Number> * input_matrix, bool apply_dirichlet_bc)
487 {
488  TransientRBThetaExpansion & trans_theta_expansion =
489  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
490 
491  TransientRBAssemblyExpansion & trans_assembly_expansion =
492  cast_ref<TransientRBAssemblyExpansion &>(get_rb_assembly_expansion());
493 
494  if (q >= trans_theta_expansion.get_n_M_terms())
495  libmesh_error_msg("Error: We must have q < Q_m in assemble_Mq_matrix.");
496 
497  input_matrix->zero();
499  &trans_assembly_expansion.get_M_assembly(q),
500  input_matrix,
502  false, /* symmetrize */
503  apply_dirichlet_bc);
504 }
505 
507 {
509 
510  TransientRBThetaExpansion & trans_theta_expansion =
511  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
512 
513  for (unsigned int q=0; q<trans_theta_expansion.get_n_M_terms(); q++)
515 
517  {
518  for (unsigned int q=0; q<trans_theta_expansion.get_n_M_terms(); q++)
520  }
521 }
522 
524 {
525  libMesh::out << "Assembling L2 matrix" << std::endl;
527 
529  {
530  libMesh::out << "Assembling non-Dirichlet L2 matrix" << std::endl;
531  assemble_L2_matrix(non_dirichlet_L2_matrix.get(), /* apply_dirichlet_bc = */ false);
532  }
533 
535 }
536 
538 {
539  LOG_SCOPE("truth_solve()", "TransientRBConstruction");
540 
541  const RBParameters & mu = get_parameters();
542  const unsigned int n_time_steps = get_n_time_steps();
543 
544  // // NumericVector for computing true L2 error
545  // UniquePtr<NumericVector<Number>> temp = NumericVector<Number>::build();
546  // temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
547 
548  // Apply initial condition again.
550  set_time_step(0);
551 
552  // Now compute the truth outputs
553  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
554  {
555  truth_outputs_all_k[n][0] = 0.;
556  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
557  {
559  get_output_vector(n,q_l)->dot(*solution);
560  }
561  }
562 
563  // Load initial projection error into temporal_data dense matrix
566 
567  for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
568  {
569  set_time_step(time_level);
570 
572 
573  // We assume that the truth assembly has been attached to the system
574  truth_assembly();
575 
576  // truth_assembly assembles into matrix and rhs, so use those for the solve
578 
579  // The matrix doesn't change at each timestep, so we
580  // can set reuse_preconditioner == true
581  linear_solver->reuse_preconditioner(true);
582 
583  if (assert_convergence)
584  {
586  }
587 
588  // Now compute the truth outputs
589  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
590  {
591  truth_outputs_all_k[n][time_level] = 0.;
592  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
593  {
594  truth_outputs_all_k[n][time_level] +=
596  }
597  }
598 
599  // load projection error into column _k of temporal_data matrix
602 
603  if ((write_interval > 0) && (time_level%write_interval == 0))
604  {
605  libMesh::out << std::endl << "Truth solve, plotting time step " << time_level << std::endl;
606 
607  std::ostringstream file_name;
608 
609  file_name << "truth.e.";
610  file_name << std::setw(3)
611  << std::setprecision(0)
612  << std::setfill('0')
613  << std::right
614  << time_level;
615 
616 #ifdef LIBMESH_HAVE_EXODUS_API
617  ExodusII_IO(get_mesh()).write_equation_systems (file_name.str(),
618  this->get_equation_systems());
619 #endif
620  }
621  }
622 
623  // Set reuse_preconditioner back to false for subsequent solves.
624  linear_solver->reuse_preconditioner(false);
625 
626  // Get the L2 norm of the truth solution at time-level _K
627  // Useful for normalizing our true error data
629  Real final_truth_L2_norm = libmesh_real(std::sqrt(inner_product_storage_vector->dot(*solution)));
630 
631 
632  return final_truth_L2_norm;
633 }
634 
635 bool
637  Real initial_greedy_error,
638  int count)
639 {
640  if ((get_max_truth_solves()>0) && (count >= get_max_truth_solves()))
641  {
642  libMesh::out << "Maximum number of truth solves reached: max = "
643  << count << std::endl;
644  return true;
645  }
646 
647  return Parent::greedy_termination_test(abs_greedy_error, initial_greedy_error, count);
648 }
649 
651 {
652  LOG_SCOPE("set_error_temporal_data()", "TransientRBConstruction");
653 
654  // first compute the projection of solution onto the current
655  // RB space
656 
657  const unsigned int time_step = get_time_step();
658 
659  if (get_rb_evaluation().get_n_basis_functions() == 0)
660  {
661  // If the basis is empty, then the error is the solution itself
662  temporal_data[time_step]->zero();
663  temporal_data[time_step]->add(1., *solution);
664  }
665  else
666  {
667  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
668 
670  temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
671 
672  // First compute the right-hand side vector for the projection
673  inner_product_matrix->vector_mult(*temp, *solution);
674 
675  // zero_dirichlet_dofs_on_vector(*temp);
676 
677  // Do not assume that RB_stiffness matrix is diagonal,
678  // diagonality degrades as N increases
679 
680  // Get an appropriately sized copy of RB_inner_product_matrix
681  DenseMatrix<Number> RB_inner_product_matrix_N(RB_size,RB_size);
682  for (unsigned int i=0; i<RB_size; i++)
683  for (unsigned int j=0; j<RB_size; j++)
684  {
685  RB_inner_product_matrix_N(i,j) = get_rb_evaluation().RB_inner_product_matrix(i,j);
686  }
687 
688  // Compute the projection RHS
689  DenseVector<Number> RB_proj_rhs(RB_size);
690  for (unsigned int i=0; i<RB_size; i++)
691  {
692  RB_proj_rhs(i) = temp->dot(get_rb_evaluation().get_basis_function(i));
693  }
694 
695  DenseVector<Number> RB_proj(RB_size);
696 
697  // Now solve the linear system
698  RB_inner_product_matrix_N.lu_solve(RB_proj_rhs, RB_proj);
699 
700  // Load the RB projection into temp
701  temp->zero();
702  for (unsigned int i=0; i<RB_size; i++)
703  {
704  temp->add(RB_proj(i), get_rb_evaluation().get_basis_function(i));
705  }
706 
707  temp->add(-1., *solution);
708 
709  // Now temp holds the projection error, store in temporal_data
710  *(temporal_data[time_step]) = *temp;
711  }
712 
713  // return the square of the X norm of the truth solution
715 
717 }
718 
720 {
721  LOG_SCOPE("get_error_temporal_data()", "TransientRBConstruction");
722 
723  const unsigned int time_step = get_time_step();
724 
725  return *temporal_data[time_step];
726 }
727 
729 {
731  {
732  // Use System::read_serialized_data to read the initial condition
733  // into this->solution
734  Xdr IC_data(init_filename, READ);
735  read_serialized_data(IC_data, false);
736  }
737  else
738  {
739  // Otherwise zero out the solution as a default
740  this->solution->zero();
741  }
742  this->solution->close();
743  this->update();
744 }
745 
747 {
748  LOG_SCOPE("add_IC_to_RB_space()", "TransientRBConstruction");
749 
750  if (get_rb_evaluation().get_n_basis_functions() > 0)
751  libmesh_error_msg("Error: Should not call TransientRBConstruction::add_IC_to_RB_space() " \
752  << "on a system that already contains basis functions.");
753 
755  libmesh_error_msg("Error: Should not call TransientRBConstruction::add_IC_to_RB_space() " \
756  << "when nonzero_initialization==false.");
757 
759 
760  // load the new basis function into the basis_functions vector.
761  get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build(this->comm()).release() );
762  NumericVector<Number> & current_bf = get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
763  current_bf.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
764  current_bf = *solution;
765 
766  // We can just set the norm to 1.
768 
769  Real current_bf_norm = libmesh_real(std::sqrt( current_bf.dot(*inner_product_storage_vector) ));
770  current_bf.scale(1./current_bf_norm);
771 
772  unsigned int saved_delta_N = get_delta_N();
773  set_delta_N(1);
774  update_system();
775  set_delta_N(saved_delta_N);
776 }
777 
779 {
780  // Need SLEPc to get the POD eigenvalues
781 #if defined(LIBMESH_HAVE_SLEPC)
782  LOG_SCOPE("enrich_RB_space()", "TransientRBConstruction");
783 
784  // With the "method of snapshots", the size of
785  // the eigenproblem is determined by the number
786  // of time-steps (rather than the number of spatial dofs).
787  PetscBLASInt eigen_size = temporal_data.size();
788  PetscBLASInt LDA = eigen_size; // The leading order of correlation_matrix
789  std::vector<Number> correlation_matrix(LDA*eigen_size);
790 
791  for (int i=0; i<eigen_size; i++)
792  {
794 
795  for (int j=i; j<eigen_size; j++)
796  {
797  // Scale the inner products by the number of time-steps to normalize the
798  // POD energy norm appropriately
799  Number inner_prod = (temporal_data[j]->dot(*inner_product_storage_vector)) /
800  (Real)(get_n_time_steps()+1);
801 
802  // Fill upper triangular part of correlation_matrix
803  correlation_matrix[j*eigen_size+i] = inner_prod;
804  }
805  }
806 
807  // Call the LAPACK eigensolver
808  char JOBZ = 'V'; // Compute eigenvectors and eigenvalues
809  char RANGE = 'A'; // Compute all eigenvalues
810  char UPLO = 'U'; // Upper triangular symmetric matrix
811 
812  Real VL = 0.; // Not used when RANGE = A
813  Real VU = 0.; // Not used when RANGE = A
814 
815  PetscBLASInt IL = 0; // Not used when RANGE = A
816  PetscBLASInt IU = 0; // Not used when RANGE = A
817 
818  Real ABSTOL = 1.e-14; // Absolute tolerance for eigensolver
819 
820  PetscBLASInt M = 0; // (output) The total number of evals found
821 
822  std::vector<Real> W(eigen_size); // (output) the eigenvalues
823 
824  PetscBLASInt LDZ = eigen_size; // The leading order of Z
825  std::vector<Number> Z(LDZ*eigen_size); // (output) the eigenvectors
826 
827  std::vector<PetscBLASInt> ISUPPZ(2*eigen_size); // Indicates which evecs in Z are nonzero
828 
829  // Work array, sized according to lapack documentation
830  PetscBLASInt LWORK = 26*eigen_size;
831  std::vector<Number> WORK(LWORK);
832 
833  // Work array, sized according to lapack documentation
834  PetscBLASInt LIWORK = 10*eigen_size;
835  std::vector<PetscBLASInt> IWORK(LIWORK);
836 
837  PetscBLASInt INFO = 0;
838 
839 #ifdef LIBMESH_USE_REAL_NUMBERS // Use real numbers
840  // Call the eigensolver for symmetric eigenvalue problems.
841  // NOTE: evals in W are in ascending order
842  LAPACKsyevr_(&JOBZ, &RANGE, &UPLO, &eigen_size, &correlation_matrix[0],
843  &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, &W[0], &Z[0], &LDZ,
844  &ISUPPZ[0], &WORK[0], &LWORK, &IWORK[0], &LIWORK, &INFO );
845 #elif LIBMESH_USE_COMPLEX_NUMBERS // Use complex numbers
846  // Need some extra data in the complex case
847 
848  // Work array, sized according to lapack documentation
849  PetscBLASInt LRWORK = 24*eigen_size;
850  std::vector<Real> RWORK(LRWORK);
851 
852  // Call the eigensolver for symmetric eigenvalue problems.
853  // NOTE: evals in W are in ascending order
854  LAPACKsyevr_(&JOBZ, &RANGE, &UPLO, &eigen_size, &correlation_matrix[0],
855  &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, &W[0], &Z[0], &LDZ,
856  &ISUPPZ[0], &WORK[0], &LWORK, &RWORK[0], &LRWORK, &IWORK[0],
857  &LIWORK, &INFO );
858 #else
859 #error libMesh does not yet support quaternions!
860 #endif
861 
862  if (INFO != 0)
863  libmesh_error_msg("Error in LAPACK syev eigensolver routine, INFO = " << INFO);
864 
865  // eval and evec now hold the sorted eigenvalues/eigenvectors
866  libMesh::out << std::endl << "POD Eigenvalues:" << std::endl;
867  for (unsigned int i=0; i<=1; i++)
868  {
869  libMesh::out << "eigenvalue " << i << " = " << W[eigen_size-1-i] << std::endl;
870  }
871  libMesh::out << "..." << std::endl;// << "." << std::endl << "." << std::endl;
872  libMesh::out << "last eigenvalue = " << W[0] << std::endl;
873  libMesh::out << std::endl;
874 
875  // Now load the new basis functions
876  unsigned int count = 0;
877  while (true)
878  {
879  // load the new basis function into the basis_functions vector.
880  get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build(this->comm()).release() );
881  NumericVector<Number> & current_bf = get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
882  current_bf.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
883  current_bf.zero();
884 
885  // Perform the matrix multiplication of temporal data with
886  // the next POD eigenvector
887  for (int j=0; j<eigen_size; j++)
888  {
889  int Z_row = j;
890  int Z_col = (eigen_size-1-count);
891  int Z_index = Z_col*eigen_size + Z_row;
892  current_bf.add(Z[Z_index], *temporal_data[j]);
893  }
894 
895  // We just set the norm to 1.
896  inner_product_matrix->vector_mult(*inner_product_storage_vector,current_bf);
897 
898  Real current_bf_norm = std::abs( std::sqrt( current_bf.dot(*inner_product_storage_vector) ) );
899  current_bf.scale(1./current_bf_norm);
900 
901  // Increment count here since we use the incremented counter
902  // in the if clauses below
903  count++;
904 
905  // If positive POD_tol, we use it to determine the number of basis functions
906  // to add, and then break the loop when POD_tol is satisfied, or after Nmax
907  // basis functions have been added. Else we break the loop after delta_N
908  // (or Nmax) new basis functions.
909  if (POD_tol > 0.)
910  {
911  set_delta_N(1);
912 
913  // We need to define the updated RB system matrices before the RB solve
914  update_system();
915  Real error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions());
916 
917  if ((error_bound <= POD_tol) || (get_rb_evaluation().get_n_basis_functions()==get_Nmax()))
918  {
919  set_delta_N(0);
920  break;
921  }
922  }
923  else
924  {
925  if (count == get_delta_N())
926  {
927  break;
928  }
929  else
930  if (get_rb_evaluation().get_n_basis_functions()==get_Nmax())
931  {
932  set_delta_N(count);
933  break;
934  }
935  }
936  }
937 
938 #else
939  libmesh_not_implemented();
940 #endif
941 }
942 
943 
945 {
946  // If delta_N is set to zero, there is nothing to update
947  if (get_delta_N() == 0)
948  return;
949 
951 
952  libMesh::out << "Updating RB initial conditions" << std::endl;
954 }
955 
957 {
958  return *L2_matrix;
959 }
960 
962 {
963  LOG_SCOPE("load_rb_solution()", "TransientRBConstruction");
964 
965  solution->zero();
966 
967  const unsigned int time_step = get_time_step();
968 
969  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
970  DenseVector<Number> RB_solution_vector_k = trans_rb_eval.RB_temporal_solution_data[time_step];
971 
972  if (RB_solution_vector_k.size() > get_rb_evaluation().get_n_basis_functions())
973  libmesh_error_msg("ERROR: rb_eval object contains " \
975  << " basis functions. RB_solution vector constains " \
976  << RB_solution_vector_k.size() \
977  << " entries. RB_solution in TransientRBConstruction::load_rb_solution is too long!");
978 
979  for (unsigned int i=0; i<RB_solution_vector_k.size(); i++)
980  solution->add(RB_solution_vector_k(i), get_rb_evaluation().get_basis_function(i));
981 
982  update();
983 }
984 
986 {
987  LOG_SCOPE("update_RB_system_matrices()", "TransientRBConstruction");
988 
990 
991  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
992 
993  TransientRBThetaExpansion & trans_theta_expansion =
994  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
995  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
996 
997  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
998 
1000  temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1001 
1002  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1003  {
1004  for (unsigned int j=0; j<RB_size; j++)
1005  {
1006  // Compute reduced L2 matrix
1007  temp->zero();
1008  L2_matrix->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1009 
1011  trans_rb_eval.RB_L2_matrix(i,j) = value;
1012  if (i!=j)
1013  {
1014  // The L2 matrix is assumed
1015  // to be symmetric
1016  trans_rb_eval.RB_L2_matrix(j,i) = value;
1017  }
1018 
1019  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1020  {
1021  // Compute reduced M_q matrix
1022  temp->zero();
1023  get_M_q(q_m)->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1024 
1025  value = (get_rb_evaluation().get_basis_function(i)).dot(*temp);
1026  trans_rb_eval.RB_M_q_vector[q_m](i,j) = value;
1027 
1028  if (i!=j)
1029  {
1030  // Each mass matrix term is assumed
1031  // to be symmetric
1032  trans_rb_eval.RB_M_q_vector[q_m](j,i) = value;
1033  }
1034  }
1035 
1036  }
1037  }
1038 }
1039 
1040 
1041 
1042 
1043 void TransientRBConstruction::update_residual_terms(bool compute_inner_products)
1044 {
1045  LOG_SCOPE("update_residual_terms()", "TransientRBConstruction");
1046 
1047  Parent::update_residual_terms(compute_inner_products);
1048 
1049  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1050 
1051  TransientRBThetaExpansion & trans_theta_expansion =
1052  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
1053 
1054  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
1055  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
1056  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
1057 
1058  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1059 
1060  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1061  {
1062  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1063  {
1064  // Initialize the vectors when we need them
1065  if (!trans_rb_eval.M_q_representor[q_m][i])
1066  {
1067  trans_rb_eval.M_q_representor[q_m][i] = (NumericVector<Number>::build(this->comm()).release());
1068  trans_rb_eval.M_q_representor[q_m][i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1069  }
1070 
1071  libmesh_assert(trans_rb_eval.M_q_representor[q_m][i]->size() == this->n_dofs() &&
1072  trans_rb_eval.M_q_representor[q_m][i]->local_size() == this->n_local_dofs() );
1073 
1074  rhs->zero();
1075  M_q_vector[q_m]->vector_mult(*rhs, get_rb_evaluation().get_basis_function(i));
1076 
1077  if (!is_quiet())
1078  libMesh::out << "Starting solve i="
1079  << i << " in TransientRBConstruction::update_residual_terms() at "
1080  << Utility::get_timestamp() << std::endl;
1081 
1083 
1084  if (assert_convergence)
1086 
1087  if (!is_quiet())
1088  {
1089  libMesh::out << "Finished solve i="
1090  << i << " in TransientRBConstruction::update_residual_terms() at "
1091  << Utility::get_timestamp() << std::endl;
1092 
1094  << " iterations, final residual "
1095  << this->final_linear_residual() << std::endl;
1096  }
1097 
1098  *trans_rb_eval.M_q_representor[q_m][i] = *solution;
1099  }
1100  }
1101 
1102  // Now compute and store the inner products if requested
1103  if (compute_inner_products)
1104  {
1105  for (unsigned int q_f=0; q_f<Q_f; q_f++)
1106  {
1108 
1109  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1110  {
1111  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1112  {
1113  trans_rb_eval.Fq_Mq_representor_innerprods[q_f][q_m][i] =
1114  trans_rb_eval.M_q_representor[q_m][i]->dot(*inner_product_storage_vector);
1115  } // end for q_m
1116  } // end for i
1117  } // end for q_f
1118 
1119  unsigned int q=0;
1120  for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
1121  {
1122  for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
1123  {
1124  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1125  {
1126  for (unsigned int j=0; j<RB_size; j++)
1127  {
1128  inner_product_matrix->vector_mult(*inner_product_storage_vector, *trans_rb_eval.M_q_representor[q_m2][j]);
1129 
1130  trans_rb_eval.Mq_Mq_representor_innerprods[q][i][j] =
1131  trans_rb_eval.M_q_representor[q_m1][i]->dot(*inner_product_storage_vector);
1132 
1133  if (i != j)
1134  {
1136  *trans_rb_eval.M_q_representor[q_m2][i]);
1137 
1138  trans_rb_eval.Mq_Mq_representor_innerprods[q][j][i] =
1139  trans_rb_eval.M_q_representor[q_m1][j]->dot(*inner_product_storage_vector);
1140  }
1141  } // end for j
1142  } // end for i
1143  q++;
1144  } // end for q_m2
1145  } // end for q_m1
1146 
1147 
1148  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1149  {
1150  for (unsigned int j=0; j<RB_size; j++)
1151  {
1152  for (unsigned int q_a=0; q_a<Q_a; q_a++)
1153  {
1154  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1155  {
1157  *trans_rb_eval.M_q_representor[q_m][j]);
1158 
1159  trans_rb_eval.Aq_Mq_representor_innerprods[q_a][q_m][i][j] =
1160  trans_rb_eval.Aq_representor[q_a][i]->dot(*inner_product_storage_vector);
1161 
1162  if (i != j)
1163  {
1165  *trans_rb_eval.M_q_representor[q_m][i]);
1166 
1167  trans_rb_eval.Aq_Mq_representor_innerprods[q_a][q_m][j][i] =
1168  trans_rb_eval.Aq_representor[q_a][j]->dot(*inner_product_storage_vector);
1169  }
1170  } // end for q_m
1171  } // end for q_a
1172  } // end for j
1173  } // end for i
1174  } // end if (compute_inner_products)
1175 }
1176 
1177 
1179 {
1180  LOG_SCOPE("update_RB_initial_condition_all_N()", "TransientRBConstruction");
1181 
1182  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1183 
1184  // Load the initial condition into the solution vector
1185  initialize_truth();
1186 
1188  temp1->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1189 
1191  temp2->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1192 
1193 
1194  unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1195 
1196  // First compute the right-hand side vector for the L2 projection
1197  L2_matrix->vector_mult(*temp1, *solution);
1198 
1199  for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1200  {
1201  RB_ic_proj_rhs_all_N(i) = temp1->dot(get_rb_evaluation().get_basis_function(i));
1202  }
1203 
1204 
1205  // Now compute the projection for each N
1206  DenseMatrix<Number> RB_L2_matrix_N;
1207  DenseVector<Number> RB_rhs_N;
1208  for (unsigned int N=(RB_size-delta_N); N<RB_size; N++)
1209  {
1210  // We have to index here by N+1 since the loop index is zero-based.
1211  trans_rb_eval.RB_L2_matrix.get_principal_submatrix(N+1, RB_L2_matrix_N);
1212 
1214 
1215  DenseVector<Number> RB_ic_N(N+1);
1216 
1217  // Now solve the linear system
1218  RB_L2_matrix_N.lu_solve(RB_rhs_N, RB_ic_N);
1219 
1220  // Load RB_ic_N into RB_initial_condition_all_N
1221  trans_rb_eval.RB_initial_condition_all_N[N] = RB_ic_N;
1222 
1223  // Compute the L2 error for the RB initial condition
1224  // This part is dependent on the truth space.
1225 
1226  // load the RB solution into temp1
1227  temp1->zero();
1228  for (unsigned int i=0; i<N+1; i++)
1229  {
1230  temp1->add(RB_ic_N(i), get_rb_evaluation().get_basis_function(i));
1231  }
1232 
1233  // subtract truth initial condition from RB_ic_N
1234  temp1->add(-1., *solution);
1235 
1236  // Compute L2 norm error, i.e. sqrt(M(solution,solution))
1237  temp2->zero();
1238  L2_matrix->vector_mult(*temp2, *temp1);
1239 
1240  trans_rb_eval.initial_L2_error_all_N[N] = libmesh_real(std::sqrt(temp2->dot(*temp1)));
1241  }
1242 }
1243 
1244 //Real TransientRBConstruction::uncached_compute_residual_dual_norm(const unsigned int N)
1245 //{
1246 // LOG_SCOPE("uncached_compute_residual_dual_norm()", "TransientRBConstruction");
1247 //
1248 // // This is the "slow" way of computing the residual, but it is useful
1249 // // for validating the "fast" way.
1250 // // Note that this only works in serial since otherwise each processor will
1251 // // have a different parameter value during the Greedy training.
1252 //
1253 // // Assemble the right-hand side to find the Reisz representor
1254 // // of the residual in the X norm
1255 // UniquePtr<NumericVector<Number>> RB_sol = NumericVector<Number>::build();
1256 // RB_sol->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1257 // RB_sol->zero();
1258 //
1259 // UniquePtr<NumericVector<Number>> ghosted_temp = NumericVector<Number>::build();
1260 // ghosted_temp->init (this->n_dofs(), this->n_local_dofs(),
1261 // this->get_dof_map().get_send_list(), false,
1262 // GHOSTED);
1263 //
1264 // UniquePtr<NumericVector<Number>> parallel_temp = NumericVector<Number>::build();
1265 // parallel_temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1266 //
1267 // // Store current_local_solution, since we don't want to corrupt it
1268 // *ghosted_temp = *current_local_solution;
1269 //
1270 // for (unsigned int i=0; i<N; i++)
1271 // {
1272 // RB_sol->add(RB_solution(i), rb_eval->get_basis_function(i));
1273 // parallel_temp->add(old_RB_solution(i), rb_eval->get_basis_function(i));
1274 // }
1275 //
1276 // // Load parallel_temp into current_local_solution in order to do assembly
1277 // const std::vector<unsigned int> & send_list = this->get_dof_map().get_send_list();
1278 // parallel_temp->localize (*current_local_solution, send_list);
1279 //
1280 // // Load the system_matrix
1281 // this->truth_assembly();
1282 //
1283 // // Restore current_local_solution
1284 // *current_local_solution = *ghosted_temp;
1285 //
1286 // matrix->vector_mult(*parallel_temp, *RB_sol);
1287 // rhs->add(-1., *parallel_temp);
1288 // rhs->close();
1289 //
1290 // zero_dirichlet_dofs_on_rhs();
1291 //
1292 // // Then solve the system to get the Reisz representor
1293 // matrix->zero();
1294 // {
1295 // matrix->add(1., *inner_product_matrix);
1296 // if (constrained_problem)
1297 // matrix->add(1., *constraint_matrix);
1298 // }
1299 //
1300 //
1301 // solution->zero();
1302 // solve();
1303 // // Make sure we didn't max out the number of iterations
1304 // if ((this->n_linear_iterations() >=
1305 // this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations")) &&
1306 // (this->final_linear_residual() >
1307 // this->get_equation_systems().parameters.get<Real>("linear solver tolerance")))
1308 // {
1309 // libmesh_error_msg("Warning: Linear solver may not have converged! Final linear residual = "
1310 // << this->final_linear_residual() << ", number of iterations = "
1311 // << this->n_linear_iterations());
1312 // }
1313 //
1314 // inner_product_matrix->vector_mult(*inner_product_storage_vector, *solution);
1315 //
1316 // Real slow_residual_norm_sq = solution->dot(*inner_product_storage_vector);
1317 //
1318 // return libmesh_real(std::sqrt( slow_residual_norm_sq ));
1319 //}
1320 
1321 void TransientRBConstruction::write_riesz_representors_to_files(const std::string & riesz_representors_dir,
1322  const bool write_binary_residual_representors)
1323 {
1324  LOG_SCOPE("write_riesz_representors_to_files()", "TransientRBConstruction");
1325 
1326  // Write out the M_q_representors. These are useful to have when restarting,
1327  // so you don't have to recompute them all over again. There should be
1328  // this->rb_eval->get_n_basis_functions() of these.
1329  libMesh::out << "Writing out the M_q_representors..." << std::endl;
1330 
1331  std::ostringstream file_name;
1332  const std::string riesz_representor_suffix = (write_binary_residual_representors ? ".xdr" : ".dat");
1333 
1334  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1335 
1336  const unsigned int istop = trans_rb_eval.get_n_basis_functions();
1337  const unsigned int istart = istop-get_delta_N();
1338 
1339  for (std::size_t q=0; q<trans_rb_eval.M_q_representor.size(); ++q)
1340  for (unsigned int i=istart; i<istop; ++i)
1341  {
1342  libMesh::out << "Writing out M_q_representor[" << q << "][" << i << "]..." << std::endl;
1343  libmesh_assert(trans_rb_eval.M_q_representor[q][i]);
1344 
1345  file_name.str(""); // reset filename
1346  file_name << riesz_representors_dir << "/M_q_representor" << i << riesz_representor_suffix;
1347 
1348  {
1349  // No need to copy!
1350  //*solution = *(M_q_representor[q][i]);
1351  trans_rb_eval.M_q_representor[q][i]->swap(*solution);
1352 
1353  Xdr mr_data(file_name.str(),
1354  write_binary_residual_representors ? ENCODE : WRITE);
1355 
1356  write_serialized_data(mr_data, false);
1357 
1358  // Synchronize before moving on
1359  this->comm().barrier();
1360 
1361  // Swap back.
1362  trans_rb_eval.M_q_representor[q][i]->swap(*solution);
1363 
1364  // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
1365  // for the system call, be sure to do it only on one processor, etc.
1366  }
1367  }
1368 }
1369 
1370 void TransientRBConstruction::read_riesz_representors_from_files(const std::string & riesz_representors_dir,
1371  const bool read_binary_residual_representors)
1372 {
1373  LOG_SCOPE("read_riesz_representors_from_files()", "TransientRBConstruction");
1374 
1375  const std::string riesz_representor_suffix =
1376  (read_binary_residual_representors ? ".xdr" : ".dat");
1377 
1378  std::ostringstream file_name;
1379  struct stat stat_info;
1380 
1381  TransientRBEvaluation & trans_rb_eval = cast_ref<TransientRBEvaluation &>(get_rb_evaluation());
1382 
1383  libMesh::out << "Reading in the M_q_representors..." << std::endl;
1384 
1385  // Read in the Aq representors. The class makes room for [Q_m][Nmax] of these. We are going to
1386  // read in [Q_m][this->rb_eval->get_n_basis_functions()]. FIXME:
1387  // should we be worried about leaks in the locations where we're about to fill entries?
1388  for (std::size_t i=0; i<trans_rb_eval.M_q_representor.size(); ++i)
1389  for (std::size_t j=0; j<trans_rb_eval.M_q_representor[i].size(); ++j)
1390  {
1391  if (trans_rb_eval.M_q_representor[i][j] != libmesh_nullptr)
1392  libmesh_error_msg("Error, must delete existing M_q_representor before reading in from file.");
1393  }
1394 
1395  // Now ready to read them in from file!
1396  for (std::size_t i=0; i<trans_rb_eval.M_q_representor.size(); ++i)
1397  for (std::size_t j=0; j<trans_rb_eval.get_n_basis_functions(); ++j)
1398  {
1399  file_name.str(""); // reset filename
1400  file_name << riesz_representors_dir
1401  << "/M_q_representor" << i << "_" << j << riesz_representor_suffix;
1402 
1403  // On processor zero check to be sure the file exists
1404  if (this->processor_id() == 0)
1405  {
1406  int stat_result = stat(file_name.str().c_str(), &stat_info);
1407 
1408  if (stat_result != 0)
1409  libmesh_error_msg("File does not exist: " << file_name.str());
1410  }
1411 
1412  Xdr aqr_data(file_name.str(),
1413  read_binary_residual_representors ? DECODE : READ);
1414 
1415  read_serialized_data(aqr_data, false);
1416 
1417  trans_rb_eval.M_q_representor[i][j] = NumericVector<Number>::build(this->comm()).release();
1418  trans_rb_eval.M_q_representor[i][j]->init (n_dofs(), n_local_dofs(),
1419  false, PARALLEL);
1420 
1421  // No need to copy, just swap
1422  //*M_q_representor[i][j] = *solution;
1423  trans_rb_eval.M_q_representor[i][j]->swap(*solution);
1424  }
1425 }
1426 
1427 } // namespace libMesh
T libmesh_real(T a)
std::vector< std::vector< std::vector< Number > > > Mq_Mq_representor_innerprods
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.
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.
virtual void write_riesz_representors_to_files(const std::string &riesz_representors_dir, const bool write_binary_residual_representors) libmesh_override
Write out all the Riesz representor data to files.
virtual void update_system() libmesh_override
Update the system after enriching the RB space.
virtual void read_riesz_representors_from_files(const std::string &riesz_representors_dir, const bool write_binary_residual_representors) libmesh_override
Write out all the Riesz representor data to files.
double abs(double a)
This is the EquationSystems class.
int max_truth_solves
Maximum number of truth solves in the POD-Greedy.
UniquePtr< NumericVector< Number > > old_local_solution
All the values I need to compute my contribution to the simulation at hand.
void set_max_truth_solves(int max_truth_solves_in)
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
virtual Number eval_A_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_a at the current parameter.
virtual void load_rb_solution() libmesh_override
Load the RB solution from the current time-level into the libMesh solution vector.
UniquePtr< NumericVector< Number > > inner_product_storage_vector
We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary me...
std::vector< SparseMatrix< Number > * > non_dirichlet_M_q_vector
We sometimes also need a second set of M_q matrices that do not have the Dirichlet boundary condition...
virtual bool greedy_termination_test(Real abs_greedy_error, Real initial_greedy_error, int count) libmesh_override
Function that indicates when to terminate the Greedy basis training.
void assemble_mass_matrix(SparseMatrix< Number > *input_matrix)
Assemble the mass matrix at the current parameter and store it in input_matrix.
bool assert_convergence
A boolean flag to indicate whether to check for proper convergence after each solve.
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true)
Train the reduced basis.
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
UniquePtr< SparseMatrix< Number > > L2_matrix
The L2 matrix.
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.
void process_temporal_parameters_file(const std::string &parameters_filename)
Read in and initialize parameters from parameters_filename.
SparseMatrix< Number > * get_M_q(unsigned int q)
Get a pointer to M_q.
virtual void enrich_RB_space() libmesh_override
Add a new basis functions to the RB space.
This extends RBAssemblyExpansion to provide an assembly expansion for the case of time-dependent PDEs...
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
virtual void print_info() libmesh_override
Print out info that describes the current setup of this RBConstruction.
virtual SparseMatrix< Number > & get_matrix_for_output_dual_solves() libmesh_override
Override to return the L2 product matrix for output dual norm solves for transient state problems...
std::vector< DenseMatrix< Number > > RB_M_q_vector
Dense matrices for the RB mass matrices.
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.
virtual Number eval_M_theta(unsigned int q, const RBParameters &mu)
Evaluate theta at the current parameter.
DenseMatrix< Number > RB_inner_product_matrix
The inner product matrix.
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.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
CompareTypes< T, T2 >::supertype dot(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:446
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
virtual void initialize_truth()
This function imposes a truth initial condition, defaults to zero initial condition if the flag nonze...
virtual void assemble_all_affine_operators() libmesh_override
Assemble and store all the affine operators.
void attach_matrix(SparseMatrix< Number > &matrix)
Additional matrices may be handled with this DofMap.
Definition: dof_map.C:246
virtual void initialize_rb_construction(bool skip_matrix_assembly=false, bool skip_vector_assembly=false) libmesh_override
Allocate all the data structures necessary for the construction stage of the RB method.
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.
Real POD_tol
If positive, this tolerance determines the number of POD modes we add to the space on a call to enric...
This class is part of the rbOOmit framework.
std::vector< std::vector< Number > > truth_outputs_all_k
The truth outputs for all time-levels from the most recent truth_solve.
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.
virtual void allocate_data_structures() libmesh_override
Helper function that actually allocates all the data structures required by this class.
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
void update_RB_initial_condition_all_N()
Compute the L2 projection of the initial condition onto the RB space for 1 <= N <= RB_size and store ...
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.
The libMesh namespace provides an interface to certain functionality in the library.
NumericVector< Number > * get_Fq(unsigned int q)
Get a pointer to Fq.
std::string init_filename
The filename of the file containing the initial condition projected onto the truth mesh...
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
UniquePtr< LinearSolver< Number > > inner_product_solver
We store an extra linear solver object which we can optionally use for solving all systems in which t...
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
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
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual void assemble_misc_matrices() libmesh_override
Override to assemble the L2 matrix as well.
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
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
std::vector< std::vector< NumericVector< Number > * > > M_q_representor
Vector storing the mass matrix representors.
std::vector< SparseMatrix< Number > * > M_q_vector
Vector storing the Q_m matrices from the mass operator.
virtual void scale(const T factor)=0
Scale each element of the vector by the given factor.
Real get_euler_theta() const
Get/set euler_theta, parameter that determines the temporal discretization.
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 ...
void assemble_L2_matrix(SparseMatrix< Number > *input_matrix, bool apply_dirichlet_bc=true)
Assemble the L2 matrix.
void pull_temporal_discretization_data(RBTemporalDiscretization &other)
Pull the temporal discretization data from other.
const MeshBase & get_mesh() const
Definition: system.h:2014
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
const DofMap & get_dof_map() const
Definition: system.h:2030
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.
std::vector< Real > initial_L2_error_all_N
Vector storing initial L2 error for all 1 <= N <= RB_size.
virtual void process_parameters_file(const std::string &parameters_filename) libmesh_override
Read in the parameters from file and set up the system accordingly.
void set_L2_assembly(ElemAssembly &L2_assembly_in)
Set the L2 object.
RBAssemblyExpansion & get_rb_assembly_expansion()
Real get_control(const unsigned int k) const
Get/set the RHS control.
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.
Number set_error_temporal_data()
Set column k (i.e.
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
void get_principal_submatrix(unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
Put the sub_m x sub_n principal submatrix into dest.
Definition: dense_matrix.C:553
virtual void update_system()
Update the system after enriching the RB space; this calls a series of functions to update the system...
void barrier() const
Pause execution until all processors reach a certain point.
virtual void truth_assembly() libmesh_override
Assemble the truth system in the transient linear case.
std::vector< std::vector< std::vector< Number > > > Fq_Mq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
virtual void print_info()
Print out info that describes the current setup of this RBConstruction.
bool nonzero_initialization
Boolean flag to indicate whether we are using a non-zero initialization.
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
ElemAssembly * L2_assembly
Function pointer for assembling the L2 matrix.
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...
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 unsigned int get_n_M_terms()
Get Q_m, the number of terms in the affine expansion for the mass operator.
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
std::vector< DenseVector< Number > > RB_initial_condition_all_N
The RB initial conditions (i.e.
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.
Real get_delta_t() const
Get/set delta_t, the time-step size.
virtual Number eval_F_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_f at the current parameter.
TransientRBConstruction(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
UniquePtr< LinearSolver< Number > > linear_solver
The LinearSolver defines the interface used to solve the linear_implicit system.
void set_POD_tol(const Real POD_tol_in)
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
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?
std::vector< NumericVector< Number > * > temporal_data
Dense matrix to store the data that we use for the temporal POD.
virtual void update_residual_terms(bool compute_inner_products) libmesh_override
Compute the terms that are combined `online&#39; to determine the dual norm of the residual.
void check_convergence(LinearSolver< Number > &input_solver)
Check if the linear solver reports convergence.
unsigned int get_time_step() const
Get/set the current time-step.
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...
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
std::vector< DenseVector< Number > > RB_temporal_solution_data
Array storing the solution data at each time level from the most recent solve.
ElemAssembly provides a per-element (interior and boundary) assembly functionality.
Definition: elem_assembly.h:38
UniquePtr< SparseMatrix< Number > > non_dirichlet_L2_matrix
The L2 matrix without Dirichlet conditions enforced.
static const bool value
Definition: xdr_io.C:108
void assemble_Mq_matrix(unsigned int q, SparseMatrix< Number > *input_matrix, bool apply_dirichlet_bc=true)
Assemble the q^th affine term of the mass matrix and store it in input_matrix.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
void add_scaled_mass_matrix(Number scalar, SparseMatrix< Number > *input_matrix)
Add the scaled mass matrix (assembled for the current parameter) to input_matrix. ...
const EquationSystems & get_equation_systems() const
Definition: system.h:712
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...
int get_max_truth_solves() const
Get/set max_truth_solves, the maximum number of RB truth solves we are willing to compute in the tran...
const NumericVector< Number > & get_error_temporal_data()
Get the column of temporal_data corresponding to the current time level.
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
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
virtual void assemble_misc_matrices()
Assemble and store all the inner-product matrix, the constraint matrix (for constrained problems) and...
unsigned int get_n_time_steps() const
Get/set the total number of time-steps.
std::vector< std::vector< std::vector< std::vector< Number > > > > Aq_Mq_representor_innerprods
bool compute_truth_projection_error
Boolean flag that indicates whether we will compute the projection error for the truth solution into ...
void mass_matrix_scaled_matvec(Number scalar, NumericVector< Number > &dest, NumericVector< Number > &arg)
Perform a matrix-vector multiplication with the current mass matrix and store the result in dest...
virtual Real truth_solve(int write_interval) libmesh_override
Perform a truth solve at the current parameter.
virtual Real train_reduced_basis(const bool resize_rb_eval_data=true) libmesh_override
Train the reduced basis.
const RBParameters & get_parameters() const
Get the current parameters.
void set_delta_N(const unsigned int new_delta_N)
Set delta_N, the number of basis functions we add to the RB space from each POD.
void add_IC_to_RB_space()
Initialize RB space by adding the truth initial condition as the first RB basis function.
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.
SparseMatrix< Number > * get_non_dirichlet_M_q(unsigned int q)
Get a pointer to non_dirichlet_M_q.
void get_principal_subvector(unsigned int sub_n, DenseVector< T > &dest) const
Puts the principal subvector of size sub_n (i.e.
Definition: dense_vector.h:667
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...
virtual void assemble_affine_expansion(bool skip_matrix_assembly, bool skip_vector_assembly) libmesh_override
Override assemble_affine_expansion to also initialize RB_ic_proj_rhs_all_N, if necessary.
virtual void update_RB_system_matrices() libmesh_override
Compute the reduced basis matrices for the current basis.
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 get_all_matrices(std::map< std::string, SparseMatrix< Number > * > &all_matrices) libmesh_override
Get a map that stores pointers to all of the matrices.
virtual void allocate_data_structures()
Helper function that actually allocates all the data structures required by this class.
Real get_POD_tol() const
Get/set POD_tol.
DenseMatrix< Number > RB_L2_matrix
Dense RB L2 matrix.
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.
processor_id_type processor_id() const
DenseVector< Number > RB_ic_proj_rhs_all_N
The vector that stores the right-hand side for the initial condition projections. ...
virtual void assemble_all_affine_operators()
Assemble and store all Q_a affine operators as well as the inner-product matrix.