libMesh
transient_rb_evaluation.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_evaluation.h"
22 #include "libmesh/transient_rb_theta_expansion.h"
23 
24 // libMesh includes
25 #include "libmesh/numeric_vector.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/xdr_cxx.h"
28 #include "libmesh/parallel.h"
29 #include "libmesh/getpot.h"
30 
31 // C++ includes
32 #include <fstream>
33 #include <sstream>
34 #include <iomanip>
35 
36 namespace libMesh
37 {
38 
40  RBEvaluation(comm_in),
41  _rb_solve_data_cached(false)
42 {
43  // Indicate that we need to compute the RB
44  // inner product matrix in this case
46 }
47 
49 {
50  clear();
51 }
52 
54 {
55  Parent::clear();
56 
58 }
59 
61 {
63 
64  // Delete the M_q representors
65  for (std::size_t q_m=0; q_m<M_q_representor.size(); q_m++)
66  for (std::size_t i=0; i<M_q_representor[q_m].size(); i++)
67  {
68  if (M_q_representor[q_m][i])
69  {
70  M_q_representor[q_m][i]->clear();
71  delete M_q_representor[q_m][i];
73  }
74  }
75 }
76 
77 void TransientRBEvaluation::resize_data_structures(const unsigned int Nmax,
78  bool resize_error_bound_data)
79 {
80  LOG_SCOPE("resize_data_structures()", "TransientRBEvaluation");
81 
82  Parent::resize_data_structures(Nmax, resize_error_bound_data);
83 
84  RB_L2_matrix.resize(Nmax,Nmax);
85  RB_LHS_matrix.resize(Nmax,Nmax);
86  RB_RHS_matrix.resize(Nmax,Nmax);
87  RB_RHS_save.resize(Nmax);
88 
89  TransientRBThetaExpansion & trans_theta_expansion =
90  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
91  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
92  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
93  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
94 
95  // Allocate dense matrices for RB solves
96  RB_M_q_vector.resize(Q_m);
97  for (unsigned int q=0; q<Q_m; q++)
98  {
99  // Initialize the memory for the RB matrices
100  RB_M_q_vector[q].resize(Nmax,Nmax);
101  }
102 
103  // Initialize the initial condition storage
104  RB_initial_condition_all_N.resize(Nmax);
105  for (std::size_t i=0; i<RB_initial_condition_all_N.size(); i++)
106  {
107  // The i^th row holds a vector of lenght i+1
108  RB_initial_condition_all_N[i].resize(i+1);
109  }
110 
111  initial_L2_error_all_N.resize(Nmax, 0.);
112 
113 
114  if (resize_error_bound_data)
115  {
116  // Initialize vectors for the norms of the representors
117  Fq_Mq_representor_innerprods.resize(Q_f);
118  for (unsigned int i=0; i<Q_f; i++)
119  {
120  Fq_Mq_representor_innerprods[i].resize(Q_m);
121  for (unsigned int j=0; j<Q_m; j++)
122  {
123  Fq_Mq_representor_innerprods[i][j].resize(Nmax, 0.);
124  }
125  }
126 
127  unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
128  Mq_Mq_representor_innerprods.resize(Q_m_hat);
129  for (unsigned int i=0; i<Q_m_hat; i++)
130  {
131  Mq_Mq_representor_innerprods[i].resize(Nmax);
132  for (unsigned int j=0; j<Nmax; j++)
133  {
134  Mq_Mq_representor_innerprods[i][j].resize(Nmax, 0.);
135  }
136  }
137 
138  Aq_Mq_representor_innerprods.resize(Q_a);
139  for (unsigned int i=0; i<Q_a; i++)
140  {
141  Aq_Mq_representor_innerprods[i].resize(Q_m);
142  for (unsigned int j=0; j<Q_m; j++)
143  {
144  Aq_Mq_representor_innerprods[i][j].resize(Nmax);
145  for (unsigned int k=0; k<Nmax; k++)
146  {
147  Aq_Mq_representor_innerprods[i][j][k].resize(Nmax, 0.);
148  }
149  }
150  }
151 
152  // Resize M_q_representor
153  // This is cleared in the call to clear_riesz_representors
154  // in Parent::resize_RB_data, so just resize here
155  M_q_representor.resize(Q_m);
156  for (unsigned int q_m=0; q_m<Q_m; q_m++)
157  {
158  M_q_representor[q_m].resize(Nmax);
159  }
160  }
161 }
162 
164 {
165  LOG_SCOPE("rb_solve()", "TransientRBEvaluation");
166 
167  if (N > get_n_basis_functions())
168  libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");
169 
170  const RBParameters & mu = get_parameters();
171 
172  TransientRBThetaExpansion & trans_theta_expansion =
173  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
174  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
175  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
176  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
177 
178  const unsigned int n_time_steps = get_n_time_steps();
179  const Real dt = get_delta_t();
180  const Real euler_theta = get_euler_theta();
181 
182  // Resize the RB and error bound vectors
183  error_bound_all_k.resize(n_time_steps+1);
184  RB_outputs_all_k.resize(trans_theta_expansion.get_n_outputs());
185  RB_output_error_bounds_all_k.resize(trans_theta_expansion.get_n_outputs());
186  for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
187  {
188  RB_outputs_all_k[n].resize(n_time_steps+1, 0.);
189  RB_output_error_bounds_all_k[n].resize(n_time_steps+1, 0.);
190  }
191 
192  // First assemble the mass matrix
193  DenseMatrix<Number> RB_mass_matrix_N(N,N);
194  RB_mass_matrix_N.zero();
195  DenseMatrix<Number> RB_M_q_m;
196  for (unsigned int q_m=0; q_m<Q_m; q_m++)
197  {
198  RB_M_q_vector[q_m].get_principal_submatrix(N, RB_M_q_m);
199  RB_mass_matrix_N.add(trans_theta_expansion.eval_M_theta(q_m, mu), RB_M_q_m);
200  }
201 
202  RB_LHS_matrix.resize(N,N);
204 
205  RB_RHS_matrix.resize(N,N);
207 
208  RB_LHS_matrix.add(1./dt, RB_mass_matrix_N);
209  RB_RHS_matrix.add(1./dt, RB_mass_matrix_N);
210 
211  DenseMatrix<Number> RB_Aq_a;
212  for (unsigned int q_a=0; q_a<Q_a; q_a++)
213  {
214  RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);
215 
216  RB_LHS_matrix.add( euler_theta*trans_theta_expansion.eval_A_theta(q_a,mu), RB_Aq_a);
217  RB_RHS_matrix.add( -(1.-euler_theta)*trans_theta_expansion.eval_A_theta(q_a,mu), RB_Aq_a);
218  }
219 
220  // Add forcing terms
221  DenseVector<Number> RB_Fq_f;
222  RB_RHS_save.resize(N);
223  RB_RHS_save.zero();
224  for (unsigned int q_f=0; q_f<Q_f; q_f++)
225  {
226  RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);
227  RB_RHS_save.add(trans_theta_expansion.eval_F_theta(q_f,mu), RB_Fq_f);
228  }
229 
230  // Set system time level to 0
231  set_time_step(0);
232 
233  // Resize/clear the solution vector
234  RB_solution.resize(N);
235 
236  // Load the initial condition into RB_solution
237  if (N > 0)
238  {
240  }
241 
242  // Resize/clear the old solution vector
244 
245  // Initialize the RB rhs
246  DenseVector<Number> RB_rhs(N);
247  RB_rhs.zero();
248 
249  // Initialize the vectors storing solution data
250  RB_temporal_solution_data.resize(n_time_steps+1);
251  for (unsigned int time_level=0; time_level<=n_time_steps; time_level++)
252  {
253  RB_temporal_solution_data[time_level].resize(N);
254  }
255  // and load the initial data
257 
258  // Set outputs at initial time
259  {
260  DenseVector<Number> RB_output_vector_N;
261  for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
262  {
263  RB_outputs_all_k[n][0] = 0.;
264  for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
265  {
266  RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
267  RB_outputs_all_k[n][0] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
268  }
269  }
270  }
271 
272  // Initialize error bounds, if necessary
273  Real error_bound_sum = 0.;
274  Real alpha_LB = 0.;
276  {
277  if (N > 0)
278  {
279  error_bound_sum += pow( initial_L2_error_all_N[N-1], 2.);
280  }
281 
282  // Set error bound at the initial time
283  error_bound_all_k[get_time_step()] = std::sqrt(error_bound_sum);
284 
285  // Compute the outputs and associated error bounds at the initial time
286  DenseVector<Number> RB_output_vector_N;
287  for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
288  {
289  RB_outputs_all_k[n][0] = 0.;
290  for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
291  {
292  RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
293  RB_outputs_all_k[n][0] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
294  }
295 
297  }
298 
299  alpha_LB = get_stability_lower_bound();
300 
301  // Precompute time-invariant parts of the dual norm of the residual.
303  }
304 
305  for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
306  {
307  set_time_step(time_level);
309 
310  // Compute RB_rhs, as RB_LHS_matrix x old_RB_solution
312 
313  // Add forcing terms
314  RB_rhs.add(get_control(time_level), RB_RHS_save);
315 
316  if (N > 0)
317  {
319  }
320 
321  // Save RB_solution for current time level
323 
324  // Evaluate outputs
325  DenseVector<Number> RB_output_vector_N;
326  for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
327  {
328  RB_outputs_all_k[n][time_level] = 0.;
329  for (unsigned int q_l=0; q_l<trans_theta_expansion.get_n_output_terms(n); q_l++)
330  {
331  RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
332  RB_outputs_all_k[n][time_level] += trans_theta_expansion.eval_output_theta(n,q_l,mu)*
333  RB_output_vector_N.dot(RB_solution);
334  }
335  }
336 
337  // Calculate RB error bounds
339  {
340  // Evaluate the dual norm of the residual for RB_solution_vector
341  // Real epsilon_N = uncached_compute_residual_dual_norm(N);
342  Real epsilon_N = compute_residual_dual_norm(N);
343 
344  error_bound_sum += residual_scaling_numer(alpha_LB) * pow(epsilon_N, 2.);
345 
346  // store error bound at time-level _k
347  error_bound_all_k[time_level] = std::sqrt(error_bound_sum/residual_scaling_denom(alpha_LB));
348 
349  // Now evaluated output error bounds
350  for (unsigned int n=0; n<trans_theta_expansion.get_n_outputs(); n++)
351  {
352  RB_output_error_bounds_all_k[n][time_level] = error_bound_all_k[time_level] *
353  eval_output_dual_norm(n,mu);
354  }
355  }
356  }
357 
358  _rb_solve_data_cached = true ;
359 
360  if (evaluate_RB_error_bound) // Calculate the error bounds
361  {
362  return error_bound_all_k[n_time_steps];
363  }
364  else // Don't calculate the error bounds
365  {
366  // Just return -1. if we did not compute the error bound
367  return -1.;
368  }
369 }
370 
372 {
374 
375  const unsigned int n_time_steps = get_n_time_steps();
376  // Set system time level to 0
377  set_time_step(0);
378 
379  // Resize/clear the solution vector
380  const unsigned int N = RB_RHS_save.size();
381  RB_solution.resize(N);
382 
383  // Load the initial condition into RB_solution
384  if (N > 0)
386 
387  // Resize/clear the old solution vector
389 
390  // Initialize the RB rhs
391  DenseVector<Number> RB_rhs(N);
392  RB_rhs.zero();
393 
394  for (unsigned int time_level=1; time_level<=n_time_steps; time_level++)
395  {
396  set_time_step(time_level);
398 
399  // Compute RB_rhs, as *RB_lhs_matrix x old_RB_solution
401 
402  // Add forcing terms
403  RB_rhs.add(get_control(time_level), RB_RHS_save);
404 
405  if (N > 0)
407  }
408 
409  {
410  // Just return -1. We did not compute the error bound
411  return -1.;
412  }
413 }
414 
416 {
417  // Just set the normalization factor to 1 in this case.
418  // Users can override this method if specific behavior
419  // is required.
420 
421  return 1.;
422 }
423 
425 {
426  return get_delta_t();
427 }
428 
430 {
431  LOG_SCOPE("cache_online_residual_terms()", "TransientRBEvaluation");
432 
433  const RBParameters mu = get_parameters();
434 
435  TransientRBThetaExpansion & trans_theta_expansion =
436  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
437  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
438  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
439  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
440 
441  cached_Fq_term = 0.;
442  unsigned int q=0;
443  for (unsigned int q_f1=0; q_f1<Q_f; q_f1++)
444  {
445  Number cached_theta_q_f1 = trans_theta_expansion.eval_F_theta(q_f1,mu);
446  for (unsigned int q_f2=q_f1; q_f2<Q_f; q_f2++)
447  {
448  Real delta = (q_f1==q_f2) ? 1. : 2.;
449  cached_Fq_term += delta*cached_theta_q_f1*trans_theta_expansion.eval_F_theta(q_f2,mu) *
451 
452  q++;
453  }
454  }
455 
457  for (unsigned int q_f=0; q_f<Q_f; q_f++)
458  {
459  Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
460  for (unsigned int q_a=0; q_a<Q_a; q_a++)
461  {
462  Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
463  for (unsigned int i=0; i<N; i++)
464  {
465  cached_Fq_Aq_vector(i) += 2.*cached_theta_q_f*cached_theta_q_a*
466  Fq_Aq_representor_innerprods[q_f][q_a][i];
467  }
468  }
469  }
470 
472  q=0;
473  for (unsigned int q_a1=0; q_a1<Q_a; q_a1++)
474  {
475  Number cached_theta_q_a1 = trans_theta_expansion.eval_A_theta(q_a1,mu);
476  for (unsigned int q_a2=q_a1; q_a2<Q_a; q_a2++)
477  {
478  Number cached_theta_q_a2 = trans_theta_expansion.eval_A_theta(q_a2,mu);
479  Real delta = (q_a1==q_a2) ? 1. : 2.;
480 
481  for (unsigned int i=0; i<N; i++)
482  {
483  for (unsigned int j=0; j<N; j++)
484  {
485  cached_Aq_Aq_matrix(i,j) += delta*
486  cached_theta_q_a1*cached_theta_q_a2*
488  }
489  }
490  q++;
491  }
492  }
493 
495  for (unsigned int q_f=0; q_f<Q_f; q_f++)
496  {
497  Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
498  for (unsigned int q_m=0; q_m<Q_m; q_m++)
499  {
500  Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
501  for (unsigned int i=0; i<N; i++)
502  {
503  cached_Fq_Mq_vector(i) += 2.*cached_theta_q_f * cached_theta_q_m * Fq_Mq_representor_innerprods[q_f][q_m][i];
504  }
505  }
506  }
507 
509  for (unsigned int q_a=0; q_a<Q_a; q_a++)
510  {
511  Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
512 
513  for (unsigned int q_m=0; q_m<Q_m; q_m++)
514  {
515  Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
516 
517  for (unsigned int i=0; i<N; i++)
518  {
519  for (unsigned int j=0; j<N; j++)
520  {
521  cached_Aq_Mq_matrix(i,j) += 2.*cached_theta_q_a*cached_theta_q_m*Aq_Mq_representor_innerprods[q_a][q_m][i][j];
522  }
523  }
524  }
525  }
526 
528  q=0;
529  for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
530  {
531  Number cached_theta_q_m1 = trans_theta_expansion.eval_M_theta(q_m1,mu);
532  for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
533  {
534  Number cached_theta_q_m2 = trans_theta_expansion.eval_M_theta(q_m2,mu);
535  Real delta = (q_m1==q_m2) ? 1. : 2.;
536 
537  for (unsigned int i=0; i<N; i++)
538  {
539  for (unsigned int j=0; j<N; j++)
540  {
541  cached_Mq_Mq_matrix(i,j) += delta*
542  cached_theta_q_m1*cached_theta_q_m2*
544  }
545  }
546  q++;
547  }
548  }
549 }
550 
552 {
553  LOG_SCOPE("compute_residual_dual_norm()", "TransientRBEvaluation");
554 
555  // This assembly assumes we have already called cache_online_residual_terms
556  // and that the rb_solve parameter is constant in time
557 
558  const Real dt = get_delta_t();
559  const Real euler_theta = get_euler_theta();
560  const Real current_control = get_control(get_time_step());
561 
562  DenseVector<Number> RB_u_euler_theta(N);
563  DenseVector<Number> mass_coeffs(N);
564  for (unsigned int i=0; i<N; i++)
565  {
566  RB_u_euler_theta(i) = euler_theta*RB_solution(i) +
567  (1.-euler_theta)*old_RB_solution(i);
568  mass_coeffs(i) = -(RB_solution(i) - old_RB_solution(i))/dt;
569  }
570 
571  Number residual_norm_sq = current_control*current_control*cached_Fq_term;
572 
573  residual_norm_sq += current_control*RB_u_euler_theta.dot(cached_Fq_Aq_vector);
574  residual_norm_sq += current_control*mass_coeffs.dot(cached_Fq_Mq_vector);
575 
576  for (unsigned int i=0; i<N; i++)
577  for (unsigned int j=0; j<N; j++)
578  {
579  residual_norm_sq += RB_u_euler_theta(i)*RB_u_euler_theta(j)*cached_Aq_Aq_matrix(i,j);
580  residual_norm_sq += mass_coeffs(i)*mass_coeffs(j)*cached_Mq_Mq_matrix(i,j);
581  residual_norm_sq += RB_u_euler_theta(i)*mass_coeffs(j)*cached_Aq_Mq_matrix(i,j);
582  }
583 
584 
585  if (libmesh_real(residual_norm_sq) < 0)
586  {
587  libMesh::out << "Warning: Square of residual norm is negative "
588  << "in TransientRBEvaluation::compute_residual_dual_norm()" << std::endl;
589 
590  // Sometimes this is negative due to rounding error,
591  // but error is on the order of 1.e-10, so shouldn't
592  // affect result
593  residual_norm_sq = std::abs(residual_norm_sq);
594  }
595 
596  return libmesh_real(std::sqrt( residual_norm_sq ));
597 }
598 
600 {
601  LOG_SCOPE("uncached_compute_residual_dual_norm()", "TransientRBEvaluation");
602 
603  // Use the stored representor inner product values
604  // to evaluate the residual norm
605 
606  const RBParameters & mu = get_parameters();
607 
608  TransientRBThetaExpansion & trans_theta_expansion =
609  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
610  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
611  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
612  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
613 
614  const Real dt = get_delta_t();
615  const Real euler_theta = get_euler_theta();
616 
617  std::vector<Number> RB_u_euler_theta(N);
618  std::vector<Number> mass_coeffs(N);
619  for (unsigned int i=0; i<N; i++)
620  {
621  RB_u_euler_theta[i] = euler_theta*RB_solution(i) +
622  (1.-euler_theta)*old_RB_solution(i);
623  mass_coeffs[i] = -(RB_solution(i) - old_RB_solution(i))/dt;
624  }
625 
626  Number residual_norm_sq = 0.;
627 
628  unsigned int q=0;
629  for (unsigned int q_f1=0; q_f1<Q_f; q_f1++)
630  {
631  Number cached_theta_q_f1 = trans_theta_expansion.eval_F_theta(q_f1,mu);
632  for (unsigned int q_f2=q_f1; q_f2<Q_f; q_f2++)
633  {
634  Real delta = (q_f1==q_f2) ? 1. : 2.;
635  residual_norm_sq += delta*cached_theta_q_f1*trans_theta_expansion.eval_F_theta(q_f2,mu) * Fq_representor_innerprods[q];
636 
637  q++;
638  }
639  }
640 
641  for (unsigned int q_f=0; q_f<Q_f; q_f++)
642  {
643  Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
644  for (unsigned int q_a=0; q_a<Q_a; q_a++)
645  {
646  Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
647  for (unsigned int i=0; i<N; i++)
648  {
649  residual_norm_sq += 2.*RB_u_euler_theta[i]*cached_theta_q_f*cached_theta_q_a*
650  Fq_Aq_representor_innerprods[q_f][q_a][i];
651  }
652  }
653  }
654 
655  q=0;
656  for (unsigned int q_a1=0; q_a1<Q_a; q_a1++)
657  {
658  Number cached_theta_q_a1 = trans_theta_expansion.eval_A_theta(q_a1,mu);
659  for (unsigned int q_a2=q_a1; q_a2<Q_a; q_a2++)
660  {
661  Number cached_theta_q_a2 = trans_theta_expansion.eval_A_theta(q_a2,mu);
662  Real delta = (q_a1==q_a2) ? 1. : 2.;
663 
664  for (unsigned int i=0; i<N; i++)
665  {
666  for (unsigned int j=0; j<N; j++)
667  {
668  residual_norm_sq += delta*RB_u_euler_theta[i]*RB_u_euler_theta[j]*
669  cached_theta_q_a1*cached_theta_q_a2*
671  }
672  }
673  q++;
674  }
675  }
676 
677  // Now add the terms due to the time-derivative
678  q=0;
679  for (unsigned int q_m1=0; q_m1<Q_m; q_m1++)
680  {
681  Number cached_theta_q_m1 = trans_theta_expansion.eval_M_theta(q_m1,mu);
682  for (unsigned int q_m2=q_m1; q_m2<Q_m; q_m2++)
683  {
684  Number cached_theta_q_m2 = trans_theta_expansion.eval_M_theta(q_m2,mu);
685  Real delta = (q_m1==q_m2) ? 1. : 2.;
686 
687  for (unsigned int i=0; i<N; i++)
688  {
689  for (unsigned int j=0; j<N; j++)
690  {
691  residual_norm_sq += delta*mass_coeffs[i]*mass_coeffs[j]*
692  cached_theta_q_m1*cached_theta_q_m2*
694  }
695  }
696  q++;
697  }
698  }
699 
700  for (unsigned int q_f=0; q_f<Q_f; q_f++)
701  {
702  Number cached_theta_q_f = trans_theta_expansion.eval_F_theta(q_f,mu);
703  for (unsigned int q_m=0; q_m<Q_m; q_m++)
704  {
705  Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
706  for (unsigned int i=0; i<N; i++)
707  {
708  residual_norm_sq += 2.*mass_coeffs[i]*cached_theta_q_f * cached_theta_q_m * Fq_Mq_representor_innerprods[q_f][q_m][i];
709  }
710  }
711  }
712 
713  for (unsigned int q_a=0; q_a<Q_a; q_a++)
714  {
715  Number cached_theta_q_a = trans_theta_expansion.eval_A_theta(q_a,mu);
716 
717  for (unsigned int q_m=0; q_m<Q_m; q_m++)
718  {
719  Number cached_theta_q_m = trans_theta_expansion.eval_M_theta(q_m,mu);
720 
721  for (unsigned int i=0; i<N; i++)
722  {
723  for (unsigned int j=0; j<N; j++)
724  {
725  residual_norm_sq += 2.*RB_u_euler_theta[i]*mass_coeffs[j]*
726  cached_theta_q_a*cached_theta_q_m*
727  Aq_Mq_representor_innerprods[q_a][q_m][i][j];
728  }
729  }
730  }
731  }
732 
733  if (libmesh_real(residual_norm_sq) < 0)
734  {
735  libMesh::out << "Warning: Square of residual norm is negative "
736  << "in TransientRBEvaluation::compute_residual_dual_norm()" << std::endl;
737 
738  // Sometimes this is negative due to rounding error,
739  // but error is on the order of 1.e-10, so shouldn't
740  // affect result
741  residual_norm_sq = std::abs(residual_norm_sq);
742  }
743 
744  // libMesh::out << "slow residual_sq = " << slow_residual_norm_sq
745  // << ", fast residual_sq = " << residual_norm_sq << std::endl;
746 
747  return libmesh_real(std::sqrt( residual_norm_sq ));
748 }
749 
750 void TransientRBEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
751  const bool write_binary_data)
752 {
753  LOG_SCOPE("legacy_write_offline_data_to_files()", "TransientRBEvaluation");
754 
756 
757  TransientRBThetaExpansion & trans_theta_expansion =
758  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
759  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
760  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
761  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
762 
763  const unsigned int n_bfs = get_n_basis_functions();
764 
765  // The writing mode: ENCODE for binary, WRITE for ASCII
766  XdrMODE mode = write_binary_data ? ENCODE : WRITE;
767 
768  // The suffix to use for all the files that are written out
769  const std::string suffix = write_binary_data ? ".xdr" : ".dat";
770 
771  if (this->processor_id() == 0)
772  {
773  std::ostringstream file_name;
774 
775  // Write out the temporal discretization data
776  file_name.str("");
777  file_name << directory_name << "/temporal_discretization_data" << suffix;
778  Xdr temporal_discretization_data_out(file_name.str(), mode);
779 
780  Real real_value; unsigned int int_value;
781  real_value = get_delta_t(); temporal_discretization_data_out << real_value;
782  real_value = get_euler_theta(); temporal_discretization_data_out << real_value;
783  int_value = get_n_time_steps(); temporal_discretization_data_out << int_value;
784  int_value = get_time_step(); temporal_discretization_data_out << int_value;
785  temporal_discretization_data_out.close();
786 
787 
788  // Write out the L2 matrix
789  file_name.str("");
790  file_name << directory_name << "/RB_L2_matrix" << suffix;
791  Xdr RB_L2_matrix_out(file_name.str(), mode);
792 
793  for (unsigned int i=0; i<n_bfs; i++)
794  {
795  for (unsigned int j=0; j<n_bfs; j++)
796  {
797  RB_L2_matrix_out << RB_L2_matrix(i,j);
798  }
799  }
800  RB_L2_matrix_out.close();
801 
802  // Write out the M_q matrices
803  for (unsigned int q_m=0; q_m<Q_m; q_m++)
804  {
805  file_name.str("");
806  file_name << directory_name << "/RB_M_";
807  file_name << std::setw(3)
808  << std::setprecision(0)
809  << std::setfill('0')
810  << std::right
811  << q_m;
812  file_name << suffix;
813  Xdr RB_M_q_m_out(file_name.str(), mode);
814 
815  for (unsigned int i=0; i<n_bfs; i++)
816  {
817  for (unsigned int j=0; j<n_bfs; j++)
818  {
819  RB_M_q_m_out << RB_M_q_vector[q_m](i,j);
820  }
821  }
822  RB_M_q_m_out.close();
823  }
824 
825  // Write out the initial condition data
826  // and the initial L2 error for all N
827  file_name.str("");
828  file_name << directory_name << "/initial_conditions" << suffix;
829  Xdr initial_conditions_out(file_name.str(), mode);
830  file_name.str("");
831  file_name << directory_name << "/initial_L2_error" << suffix;
832  Xdr initial_L2_error_out(file_name.str(), mode);
833 
834  for (unsigned int i=0; i<n_bfs; i++)
835  {
836  initial_L2_error_out << initial_L2_error_all_N[i];
837  for (unsigned int j=0; j<=i; j++)
838  {
839  initial_conditions_out << RB_initial_condition_all_N[i](j);
840  }
841  }
842  initial_conditions_out.close();
843  initial_L2_error_out.close();
844 
845  // Next write out the Fq_Mq representor norm data
846  file_name.str("");
847  file_name << directory_name << "/Fq_Mq_terms" << suffix;
848  Xdr RB_Fq_Mq_terms_out(file_name.str(), mode);
849 
850  for (unsigned int q_f=0; q_f<Q_f; q_f++)
851  {
852  for (unsigned int q_m=0; q_m<Q_m; q_m++)
853  {
854  for (unsigned int i=0; i<n_bfs; i++)
855  {
856  RB_Fq_Mq_terms_out << Fq_Mq_representor_innerprods[q_f][q_m][i];
857  }
858  }
859  }
860  RB_Fq_Mq_terms_out.close();
861 
862  // Next write out the Mq_Mq representor norm data
863  file_name.str("");
864  file_name << directory_name << "/Mq_Mq_terms" << suffix;
865  Xdr RB_Mq_Mq_terms_out(file_name.str(), mode);
866 
867  unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
868  for (unsigned int q=0; q<Q_m_hat; q++)
869  {
870  for (unsigned int i=0; i<n_bfs; i++)
871  {
872  for (unsigned int j=0; j<n_bfs; j++)
873  {
874  RB_Mq_Mq_terms_out << Mq_Mq_representor_innerprods[q][i][j];
875  }
876  }
877  }
878  RB_Mq_Mq_terms_out.close();
879 
880  // Next write out the Aq_Mq representor norm data
881  file_name.str("");
882  file_name << directory_name << "/Aq_Mq_terms" << suffix;
883  Xdr RB_Aq_Mq_terms_out(file_name.str(), mode);
884 
885  for (unsigned int q_a=0; q_a<Q_a; q_a++)
886  {
887  for (unsigned int q_m=0; q_m<Q_m; q_m++)
888  {
889  for (unsigned int i=0; i<n_bfs; i++)
890  {
891  for (unsigned int j=0; j<n_bfs; j++)
892  {
893  RB_Aq_Mq_terms_out << Aq_Mq_representor_innerprods[q_a][q_m][i][j];
894  }
895  }
896  }
897  }
898  RB_Aq_Mq_terms_out.close();
899  }
900 }
901 
902 void TransientRBEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
903  bool read_error_bound_data,
904  const bool read_binary_data)
905 {
906  LOG_SCOPE("legacy_read_offline_data_from_files()", "TransientRBEvaluation");
907 
909 
910  TransientRBThetaExpansion & trans_theta_expansion =
911  cast_ref<TransientRBThetaExpansion &>(get_rb_theta_expansion());
912  const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
913  const unsigned int Q_a = trans_theta_expansion.get_n_A_terms();
914  const unsigned int Q_f = trans_theta_expansion.get_n_F_terms();
915 
916  // First, find out how many basis functions we had when Greedy terminated
917  // This was set in RBSystem::read_offline_data_from_files
918  unsigned int n_bfs = this->get_n_basis_functions();
919 
920  // The reading mode: DECODE for binary, READ for ASCII
921  XdrMODE mode = read_binary_data ? DECODE : READ;
922 
923  // The suffix to use for all the files that are written out
924  const std::string suffix = read_binary_data ? ".xdr" : ".dat";
925 
926  // The string stream we'll use to make the file names
927  std::ostringstream file_name;
928 
929  // Write out the temporal discretization data
930  file_name.str("");
931  file_name << directory_name << "/temporal_discretization_data" << suffix;
932  assert_file_exists(file_name.str());
933 
934  Xdr temporal_discretization_data_in(file_name.str(), mode);
935 
936  Real real_value; unsigned int int_value;
937  temporal_discretization_data_in >> real_value; set_delta_t(real_value);
938  temporal_discretization_data_in >> real_value; set_euler_theta(real_value);
939  temporal_discretization_data_in >> int_value; set_n_time_steps(int_value);
940  temporal_discretization_data_in >> int_value; set_time_step(int_value);
941  temporal_discretization_data_in.close();
942 
943  file_name.str("");
944  file_name << directory_name << "/RB_L2_matrix" << suffix;
945  assert_file_exists(file_name.str());
946 
947  Xdr RB_L2_matrix_in(file_name.str(), mode);
948 
949  for (unsigned int i=0; i<n_bfs; i++)
950  {
951  for (unsigned int j=0; j<n_bfs; j++)
952  {
953  Number value;
954  RB_L2_matrix_in >> value;
955  RB_L2_matrix(i,j) = value;
956  }
957  }
958  RB_L2_matrix_in.close();
959 
960  // Read in the M_q matrices
961  for (unsigned int q_m=0; q_m<Q_m; q_m++)
962  {
963  file_name.str("");
964  file_name << directory_name << "/RB_M_";
965  file_name << std::setw(3)
966  << std::setprecision(0)
967  << std::setfill('0')
968  << std::right
969  << q_m;
970 
971  file_name << suffix;
972  assert_file_exists(file_name.str());
973 
974  Xdr RB_M_q_m_in(file_name.str(), mode);
975 
976  for (unsigned int i=0; i<n_bfs; i++)
977  {
978  for (unsigned int j=0; j<n_bfs; j++)
979  {
980  Number value;
981  RB_M_q_m_in >> value;
982  RB_M_q_vector[q_m](i,j) = value;
983  }
984  }
985  RB_M_q_m_in.close();
986  }
987 
988 
989  // Read in the initial condition data
990  // and the initial L2 error for all N
991  file_name.str("");
992  file_name << directory_name << "/initial_conditions" << suffix;
993  assert_file_exists(file_name.str());
994 
995  Xdr initial_conditions_in(file_name.str(), mode);
996 
997  file_name.str("");
998  file_name << directory_name << "/initial_L2_error" << suffix;
999  assert_file_exists(file_name.str());
1000 
1001  Xdr initial_L2_error_in(file_name.str(), mode);
1002 
1003  for (unsigned int i=0; i<n_bfs; i++)
1004  {
1005  initial_L2_error_in >> initial_L2_error_all_N[i];
1006  for (unsigned int j=0; j<=i; j++)
1007  {
1008  initial_conditions_in >> RB_initial_condition_all_N[i](j);
1009  }
1010  }
1011  initial_conditions_in.close();
1012  initial_L2_error_in.close();
1013 
1014 
1015  if (read_error_bound_data)
1016  {
1017  // Next read in the Fq_Mq representor norm data
1018  file_name.str("");
1019  file_name << directory_name << "/Fq_Mq_terms" << suffix;
1020  assert_file_exists(file_name.str());
1021 
1022  Xdr RB_Fq_Mq_terms_in(file_name.str(), mode);
1023 
1024  for (unsigned int q_f=0; q_f<Q_f; q_f++)
1025  {
1026  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1027  {
1028  for (unsigned int i=0; i<n_bfs; i++)
1029  {
1030  RB_Fq_Mq_terms_in >> Fq_Mq_representor_innerprods[q_f][q_m][i];
1031  }
1032  }
1033  }
1034  RB_Fq_Mq_terms_in.close();
1035 
1036  // Next read in the Mq_Mq representor norm data
1037  file_name.str("");
1038  file_name << directory_name << "/Mq_Mq_terms" << suffix;
1039  assert_file_exists(file_name.str());
1040 
1041  Xdr RB_Mq_Mq_terms_in(file_name.str(), mode);
1042 
1043  unsigned int Q_m_hat = Q_m*(Q_m+1)/2;
1044  for (unsigned int q=0; q<Q_m_hat; q++)
1045  {
1046  for (unsigned int i=0; i<n_bfs; i++)
1047  {
1048  for (unsigned int j=0; j<n_bfs; j++)
1049  {
1050  RB_Mq_Mq_terms_in >> Mq_Mq_representor_innerprods[q][i][j];
1051  }
1052  }
1053  }
1054  RB_Mq_Mq_terms_in.close();
1055 
1056  // Next read in the Aq_Mq representor norm data
1057  file_name.str("");
1058  file_name << directory_name << "/Aq_Mq_terms" << suffix;
1059  assert_file_exists(file_name.str());
1060 
1061  Xdr RB_Aq_Mq_terms_in(file_name.str(), mode);
1062 
1063  for (unsigned int q_a=0; q_a<Q_a; q_a++)
1064  {
1065  for (unsigned int q_m=0; q_m<Q_m; q_m++)
1066  {
1067  for (unsigned int i=0; i<n_bfs; i++)
1068  {
1069  for (unsigned int j=0; j<n_bfs; j++)
1070  {
1071  RB_Aq_Mq_terms_in >> Aq_Mq_representor_innerprods[q_a][q_m][i][j];
1072  }
1073  }
1074  }
1075  }
1076  RB_Aq_Mq_terms_in.close();
1077  }
1078 }
1079 
1080 }
T libmesh_real(T a)
TransientRBEvaluation(const Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
std::vector< std::vector< std::vector< Number > > > Mq_Mq_representor_innerprods
virtual void clear_riesz_representors() libmesh_override
Clear all the Riesz representors that are used to compute the RB residual (and hence error bound)...
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool evaluate_RB_error_bound
Boolean to indicate whether we evaluate a posteriori error bounds when rb_solve is called...
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
Adds factor times vec to this vector.
Definition: dense_vector.h:431
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
virtual Number eval_A_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_a at the current parameter.
virtual Real uncached_compute_residual_dual_norm(const unsigned int N)
Compute the dual norm of the residual for the solution saved in RB_solution.
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
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.
DenseVector< Number > RB_solution
The RB solution vector.
Number cached_Fq_term
Cached residual terms.
DenseMatrix< Number > RB_LHS_matrix
Cached data for subsequent solves.
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
std::vector< DenseMatrix< Number > > RB_M_q_vector
Dense matrices for the RB mass matrices.
virtual Number eval_M_theta(unsigned int q, const RBParameters &mu)
Evaluate theta at the current parameter.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
void cache_online_residual_terms(const unsigned int N)
Helper function for caching the terms in the online residual assembly that do not change in time...
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
DenseVector< Number > old_RB_solution
The RB solution at the previous time-level.
virtual Real compute_residual_dual_norm(const unsigned int N) libmesh_override
Compute the dual norm of the residual for the solution saved in RB_solution.
const class libmesh_nullptr_t libmesh_nullptr
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", bool read_error_bound_data=true, const bool read_binary_data=true) libmesh_override
Read in the saved Offline reduced basis data to initialize the system for Online solves.
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
libmesh_assert(j)
virtual Number eval_output_theta(unsigned int output_index, unsigned int q_l, const RBParameters &mu)
Evaluate theta_q_l at the current parameter.
std::vector< std::vector< NumericVector< Number > * > > M_q_representor
Vector storing the mass matrix representors.
Real get_euler_theta() const
Get/set euler_theta, parameter that determines the temporal discretization.
std::vector< DenseVector< Number > > RB_Fq_vector
Dense vector for the RHS.
bool _rb_solve_data_cached
Check that the data has been cached in case of using rb_solve_again.
std::vector< Real > initial_L2_error_all_N
Vector storing initial L2 error for all 1 <= N <= RB_size.
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:32
virtual Real get_stability_lower_bound()
Get a lower bound for the stability constant (e.g.
virtual void legacy_write_offline_data_to_files(const std::string &directory_name="offline_data", const bool write_binary_data=true)
Write out all the data to text files in order to segregate the Offline stage from the Online stage...
virtual void legacy_read_offline_data_from_files(const std::string &directory_name="offline_data", bool read_error_bound_data=true, const bool read_binary_data=true)
Read in the saved Offline reduced basis data to initialize the system for Online solves.
virtual Real residual_scaling_numer(Real alpha_LB)
Specifies the residual scaling on the numerator to be used in the a posteriori error bound...
Real get_control(const unsigned int k) const
Get/set the RHS control.
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
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...
std::vector< std::vector< Real > > RB_output_error_bounds_all_k
The error bounds for each RB output for all time-levels from the most recent rb_solve.
std::vector< Real > error_bound_all_k
The error bound data for all time-levels from the most recent rb_solve.
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
virtual unsigned int get_n_M_terms()
Get Q_m, the number of terms in the affine expansion for the mass operator.
virtual Real residual_scaling_denom(Real alpha_LB)
Specifies the residual scaling on the denominator to be used in the a posteriori error bound...
double pow(double a, int b)
This class implements a C++ interface to the XDR (eXternal Data Representation) format.
Definition: xdr_cxx.h:68
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_again()
If a solve has already been performed, then we cached some data and we can perform a new solve much m...
std::vector< DenseVector< Number > > RB_initial_condition_all_N
The RB initial conditions (i.e.
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.
This class is part of the rbOOmit framework.
Definition: rb_evaluation.h:50
virtual Number eval_F_theta(unsigned int q, const RBParameters &mu)
Evaluate theta_q_f at the current parameter.
virtual void clear_riesz_representors()
Clear all the Riesz representors that are used to compute the RB residual (and hence error bound)...
std::vector< DenseMatrix< Number > > RB_Aq_vector
Dense matrices for the RB computations.
virtual void clear() libmesh_override
Clear this RBEvaluation object.
Definition: rb_evaluation.C:60
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int get_time_step() const
Get/set the current time-step.
Real eval_output_dual_norm(unsigned int n, const RBParameters &mu)
Evaluate the dual norm of output n for the current parameters.
void assert_file_exists(const std::string &file_name)
Helper function that checks if file_name exists.
virtual Real get_error_bound_normalization() libmesh_override
OStreamProxy out
std::vector< DenseVector< Number > > RB_temporal_solution_data
Array storing the solution data at each time level from the most recent solve.
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
virtual void legacy_write_offline_data_to_files(const std::string &directory_name="offline_data", const bool write_binary_data=true) libmesh_override
Write out all the data to text files in order to segregate the Offline stage from the Online stage...
void set_euler_theta(const Real euler_theta_in)
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true) libmesh_override
Resize and clear the data vectors corresponding to the value of Nmax.
std::vector< std::vector< Number > > RB_outputs_all_k
The RB outputs for all time-levels from the most recent rb_solve.
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
virtual Real rb_solve(unsigned int N) libmesh_override
Perform online solve for current_params with the N basis functions.
void set_n_time_steps(const unsigned int K)
const RBParameters & get_parameters() const
Get the current parameters.
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.
Definition: dense_matrix.C:382
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the rb_theta_expansion.
Definition: rb_evaluation.C:89
DenseMatrix< Number > RB_L2_matrix
Dense RB L2 matrix.
std::vector< std::vector< std::vector< Number > > > Aq_Aq_representor_innerprods
processor_id_type processor_id() const
virtual void clear() libmesh_override
Clear this TransientRBEvaluation object.
std::vector< std::vector< DenseVector< Number > > > RB_output_vectors
The vectors storing the RB output vectors.