libMesh
rb_data_deserialization.C
Go to the documentation of this file.
1 // rbOOmit: An implementation of the Certified Reduced Basis method.
2 // Copyright (C) 2009, 2010, 2015 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 #include "libmesh/libmesh_config.h"
21 #if defined(LIBMESH_HAVE_CAPNPROTO)
22 
23 //libMesh includes
24 #include "libmesh/rb_eim_evaluation.h"
25 #include "libmesh/string_to_enum.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/mesh.h"
28 #include "libmesh/rb_data_deserialization.h"
29 #include "libmesh/transient_rb_theta_expansion.h"
30 #include "libmesh/transient_rb_evaluation.h"
31 #include "libmesh/rb_eim_evaluation.h"
32 #include "libmesh/rb_scm_evaluation.h"
33 
34 // Cap'n'Proto includes
35 #include "capnp/serialize.h"
36 
37 // C++ includes
38 #include <unistd.h>
39 #include <iostream>
40 #include <fstream>
41 #include <fcntl.h>
42 
43 namespace libMesh
44 {
45 
46 namespace
47 {
48 
53 template <typename T>
54 Number load_scalar_value(const T & value)
55 {
56 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
57  return Number(value.getReal(), value.getImag());
58 #else
59  return value;
60 #endif
61 }
62 
63 }
64 
65 namespace RBDataDeserialization
66 {
67 
68 // ---- RBEvaluationDeserialization (BEGIN) ----
69 
71  :
72  _rb_eval(rb_eval)
73 {}
74 
76 {}
77 
78 void RBEvaluationDeserialization::read_from_file(const std::string & path,
79  bool read_error_bound_data)
80 {
81  LOG_SCOPE("read_from_file()", "RBEvaluationDeserialization");
82 
83  int fd = open(path.c_str(), O_RDONLY);
84  if (!fd)
85  libmesh_error_msg("Couldn't open the buffer file: " + path);
86 
87  // Turn off the limit to the amount of data we can read in
88  capnp::ReaderOptions reader_options;
89  reader_options.traversalLimitInWords = std::numeric_limits<uint64_t>::max();
90 
91  capnp::StreamFdMessageReader message(fd, reader_options);
92 
93 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
94  RBData::RBEvaluationReal::Reader rb_eval_reader =
95  message.getRoot<RBData::RBEvaluationReal>();
96 #else
97  RBData::RBEvaluationComplex::Reader rb_eval_reader =
98  message.getRoot<RBData::RBEvaluationComplex>();
99 #endif
100 
101  load_rb_evaluation_data(_rb_eval, rb_eval_reader, read_error_bound_data);
102 }
103 
104 // ---- RBEvaluationDeserialization (END) ----
105 
106 
107 // ---- TransientRBEvaluationDeserialization (BEGIN) ----
108 
111  _trans_rb_eval(trans_rb_eval)
112 {}
113 
115 {}
116 
118  bool read_error_bound_data)
119 {
120  LOG_SCOPE("read_from_file()", "TransientRBEvaluationDeserialization");
121 
122  int fd = open(path.c_str(), O_RDONLY);
123  if (!fd)
124  libmesh_error_msg("Couldn't open the buffer file: " + path);
125 
126  // Turn off the limit to the amount of data we can read in
127  capnp::ReaderOptions reader_options;
128  reader_options.traversalLimitInWords = std::numeric_limits<uint64_t>::max();
129 
130  capnp::StreamFdMessageReader message(fd, reader_options);
131 
132 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
133  RBData::TransientRBEvaluationReal::Reader trans_rb_eval_reader =
134  message.getRoot<RBData::TransientRBEvaluationReal>();
135  RBData::RBEvaluationReal::Reader rb_eval_reader =
136  trans_rb_eval_reader.getRbEvaluation();
137 #else
138  RBData::TransientRBEvaluationComplex::Reader trans_rb_eval_reader =
139  message.getRoot<RBData::TransientRBEvaluationComplex>();
140  RBData::RBEvaluationComplex::Reader rb_eval_reader =
141  trans_rb_eval_reader.getRbEvaluation();
142 #endif
143 
145  rb_eval_reader,
146  trans_rb_eval_reader,
147  read_error_bound_data);
148 }
149 
150 // ---- TransientRBEvaluationDeserialization (END) ----
151 
152 
153 // ---- RBEIMEvaluationDeserialization (BEGIN) ----
154 
157  _rb_eim_eval(rb_eim_eval)
158 {}
159 
161 {}
162 
164 {
165  LOG_SCOPE("read_from_file()", "RBEIMEvaluationDeserialization");
166 
167  int fd = open(path.c_str(), O_RDONLY);
168  if (!fd)
169  libmesh_error_msg("Couldn't open the buffer file: " + path);
170 
171  // Turn off the limit to the amount of data we can read in
172  capnp::ReaderOptions reader_options;
173  reader_options.traversalLimitInWords = std::numeric_limits<uint64_t>::max();
174 
175  capnp::StreamFdMessageReader message(fd, reader_options);
176 
177 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
178  RBData::RBEIMEvaluationReal::Reader rb_eim_eval_reader =
179  message.getRoot<RBData::RBEIMEvaluationReal>();
180  RBData::RBEvaluationReal::Reader rb_eval_reader =
181  rb_eim_eval_reader.getRbEvaluation();
182 #else
183  RBData::RBEIMEvaluationComplex::Reader rb_eim_eval_reader =
184  message.getRoot<RBData::RBEIMEvaluationComplex>();
185  RBData::RBEvaluationComplex::Reader rb_eval_reader =
186  rb_eim_eval_reader.getRbEvaluation();
187 #endif
188 
190  rb_eval_reader,
191  rb_eim_eval_reader);
192 }
193 
194 // ---- RBEIMEvaluationDeserialization (END) ----
195 
196 
197 // ---- RBSCMEvaluationDeserialization (BEGIN) ----
198 // RBSCMEvaluationDeserialization is only available if both SLEPC and
199 // GLPK are available.
200 
201 #if defined(LIBMESH_HAVE_SLEPC) && (LIBMESH_HAVE_GLPK)
202 
205  _rb_scm_eval(rb_scm_eval)
206 {}
207 
209 {}
210 
212 {
213  LOG_SCOPE("read_from_file()", "RBSCMEvaluationDeserialization");
214 
215  int fd = open(path.c_str(), O_RDONLY);
216  if (!fd)
217  libmesh_error_msg("Couldn't open the buffer file: " + path);
218 
219  // Turn off the limit to the amount of data we can read in
220  capnp::ReaderOptions reader_options;
221  reader_options.traversalLimitInWords = std::numeric_limits<uint64_t>::max();
222 
223  capnp::StreamFdMessageReader message(fd, reader_options);
224 
225  RBData::RBSCMEvaluation::Reader rb_scm_eval_reader =
226  message.getRoot<RBData::RBSCMEvaluation>();
227 
229  rb_scm_eval_reader);
230 }
231 
232 #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
233 
234 // ---- RBSCMEvaluationDeserialization (END) ----
235 
236 
237 // ---- Helper functions for loading data from buffers (BEGIN) ----
238 
240  RBData::ParameterRanges::Reader & parameter_ranges,
241  RBData::DiscreteParameterList::Reader & discrete_parameters_list)
242 {
243  // Continuous parameters
244  RBParameters parameters_min;
245  RBParameters parameters_max;
246  {
247  unsigned int n_parameter_ranges = parameter_ranges.getNames().size();
248 
249  for (unsigned int i=0; i<n_parameter_ranges; ++i)
250  {
251  std::string parameter_name = parameter_ranges.getNames()[i];
252  Real min_value = parameter_ranges.getMinValues()[i];
253  Real max_value = parameter_ranges.getMaxValues()[i];
254 
255  parameters_min.set_value(parameter_name, min_value);
256  parameters_max.set_value(parameter_name, max_value);
257  }
258  }
259 
260  // Discrete parameters
261  std::map<std::string, std::vector<Real>> discrete_parameter_values;
262  {
263  unsigned int n_discrete_parameters = discrete_parameters_list.getNames().size();
264 
265  for (unsigned int i=0; i<n_discrete_parameters; ++i)
266  {
267  std::string parameter_name = discrete_parameters_list.getNames()[i];
268 
269  auto value_list = discrete_parameters_list.getValues()[i];
270  unsigned int n_values = value_list.size();
271  std::vector<Real> values(n_values);
272  for (unsigned int j=0; j<n_values; ++j)
273  values[j] = value_list[j];
274 
275  discrete_parameter_values[parameter_name] = values;
276  }
277  }
278 
279  rb_evaluation.initialize_parameters(parameters_min,
280  parameters_max,
281  discrete_parameter_values);
282 }
283 
284 template <typename RBEvaluationReaderNumber>
286  RBEvaluationReaderNumber & rb_evaluation_reader,
287  bool read_error_bound_data)
288 {
289  // Set number of basis functions
290  unsigned int n_bfs = rb_evaluation_reader.getNBfs();
291  rb_evaluation.set_n_basis_functions(n_bfs);
292 
293  rb_evaluation.resize_data_structures(n_bfs, read_error_bound_data);
294 
295  auto parameter_ranges =
296  rb_evaluation_reader.getParameterRanges();
297  auto discrete_parameters_list =
298  rb_evaluation_reader.getDiscreteParameters();
299 
300  load_parameter_ranges(rb_evaluation,
301  parameter_ranges,
302  discrete_parameters_list);
303 
304  const RBThetaExpansion & rb_theta_expansion = rb_evaluation.get_rb_theta_expansion();
305 
306  unsigned int n_F_terms = rb_theta_expansion.get_n_F_terms();
307  unsigned int n_A_terms = rb_theta_expansion.get_n_A_terms();
308 
309  if (read_error_bound_data)
310  {
311 
312  // Fq representor inner-product data
313  {
314  unsigned int Q_f_hat = n_F_terms*(n_F_terms+1)/2;
315 
316  auto fq_innerprods_list = rb_evaluation_reader.getFqInnerprods();
317  if (fq_innerprods_list.size() != Q_f_hat)
318  libmesh_error_msg("Size error while reading Fq representor norm data from buffer.");
319 
320  for (unsigned int i=0; i < Q_f_hat; ++i)
321  rb_evaluation.Fq_representor_innerprods[i] = load_scalar_value(fq_innerprods_list[i]);
322  }
323 
324  // Fq_Aq representor inner-product data
325  {
326  auto fq_aq_innerprods_list = rb_evaluation_reader.getFqAqInnerprods();
327  if (fq_aq_innerprods_list.size() != n_F_terms*n_A_terms*n_bfs)
328  libmesh_error_msg("Size error while reading Fq-Aq representor norm data from buffer.");
329 
330  for (unsigned int q_f=0; q_f<n_F_terms; ++q_f)
331  for (unsigned int q_a=0; q_a<n_A_terms; ++q_a)
332  for (unsigned int i=0; i<n_bfs; ++i)
333  {
334  unsigned int offset = q_f*n_A_terms*n_bfs + q_a*n_bfs + i;
335  rb_evaluation.Fq_Aq_representor_innerprods[q_f][q_a][i] =
336  load_scalar_value(fq_aq_innerprods_list[offset]);
337  }
338  }
339 
340  // Aq_Aq representor inner-product data
341  {
342  unsigned int Q_a_hat = n_A_terms*(n_A_terms+1)/2;
343  auto aq_aq_innerprods_list = rb_evaluation_reader.getAqAqInnerprods();
344  if (aq_aq_innerprods_list.size() != Q_a_hat*n_bfs*n_bfs)
345  libmesh_error_msg("Size error while reading Aq-Aq representor norm data from buffer.");
346 
347  for (unsigned int i=0; i<Q_a_hat; ++i)
348  for (unsigned int j=0; j<n_bfs; ++j)
349  for (unsigned int l=0; l<n_bfs; ++l)
350  {
351  unsigned int offset = i*n_bfs*n_bfs + j*n_bfs + l;
352  rb_evaluation.Aq_Aq_representor_innerprods[i][j][l] =
353  load_scalar_value(aq_aq_innerprods_list[offset]);
354  }
355  }
356 
357  // Output dual inner-product data
358  {
359  unsigned int n_outputs = rb_theta_expansion.get_n_outputs();
360  auto output_innerprod_outer = rb_evaluation_reader.getOutputDualInnerprods();
361 
362  if (output_innerprod_outer.size() != n_outputs)
363  libmesh_error_msg("Incorrect number of outputs detected in the buffer");
364 
365  for (unsigned int output_id=0; output_id<n_outputs; ++output_id)
366  {
367  unsigned int n_output_terms = rb_theta_expansion.get_n_output_terms(output_id);
368 
369  unsigned int Q_l_hat = n_output_terms*(n_output_terms+1)/2;
370  auto output_innerprod_inner = output_innerprod_outer[output_id];
371 
372  if (output_innerprod_inner.size() != Q_l_hat)
373  libmesh_error_msg("Incorrect number of output terms detected in the buffer");
374 
375  for (unsigned int q=0; q<Q_l_hat; ++q)
376  {
377  rb_evaluation.output_dual_innerprods[output_id][q] =
378  load_scalar_value(output_innerprod_inner[q]);
379  }
380  }
381  }
382  }
383 
384  // Output vectors
385  {
386  unsigned int n_outputs = rb_theta_expansion.get_n_outputs();
387  auto output_vector_outer = rb_evaluation_reader.getOutputVectors();
388 
389  if (output_vector_outer.size() != n_outputs)
390  libmesh_error_msg("Incorrect number of outputs detected in the buffer");
391 
392  for (unsigned int output_id=0; output_id<n_outputs; ++output_id)
393  {
394  unsigned int n_output_terms = rb_theta_expansion.get_n_output_terms(output_id);
395 
396  auto output_vector_middle = output_vector_outer[output_id];
397  if (output_vector_middle.size() != n_output_terms)
398  libmesh_error_msg("Incorrect number of output terms detected in the buffer");
399 
400  for (unsigned int q_l=0; q_l<n_output_terms; ++q_l)
401  {
402  auto output_vectors_inner_list = output_vector_middle[q_l];
403 
404  if (output_vectors_inner_list.size() != n_bfs)
405  libmesh_error_msg("Incorrect number of output terms detected in the buffer");
406 
407  for (unsigned int j=0; j<n_bfs; ++j)
408  {
409  rb_evaluation.RB_output_vectors[output_id][q_l](j) =
410  load_scalar_value(output_vectors_inner_list[j]);
411  }
412  }
413  }
414  }
415 
416  // Fq vectors and Aq matrices
417  {
418  auto rb_fq_vectors_outer_list = rb_evaluation_reader.getRbFqVectors();
419  if (rb_fq_vectors_outer_list.size() != n_F_terms)
420  libmesh_error_msg("Incorrect number of Fq vectors detected in the buffer");
421 
422  for (unsigned int q_f=0; q_f<n_F_terms; ++q_f)
423  {
424  auto rb_fq_vectors_inner_list = rb_fq_vectors_outer_list[q_f];
425  if (rb_fq_vectors_inner_list.size() != n_bfs)
426  libmesh_error_msg("Incorrect Fq vector size detected in the buffer");
427 
428  for (unsigned int i=0; i < n_bfs; ++i)
429  {
430  rb_evaluation.RB_Fq_vector[q_f](i) =
431  load_scalar_value(rb_fq_vectors_inner_list[i]);
432  }
433  }
434 
435  auto rb_Aq_matrices_outer_list = rb_evaluation_reader.getRbAqMatrices();
436  if (rb_Aq_matrices_outer_list.size() != n_A_terms)
437  libmesh_error_msg("Incorrect number of Aq matrices detected in the buffer");
438 
439  for (unsigned int q_a=0; q_a<n_A_terms; ++q_a)
440  {
441  auto rb_Aq_matrices_inner_list = rb_Aq_matrices_outer_list[q_a];
442  if (rb_Aq_matrices_inner_list.size() != n_bfs*n_bfs)
443  libmesh_error_msg("Incorrect Aq matrix size detected in the buffer");
444 
445  for (unsigned int i=0; i<n_bfs; ++i)
446  for (unsigned int j=0; j<n_bfs; ++j)
447  {
448  unsigned int offset = i*n_bfs+j;
449  rb_evaluation.RB_Aq_vector[q_a](i,j) =
450  load_scalar_value(rb_Aq_matrices_inner_list[offset]);
451  }
452  }
453  }
454 
455  // Inner-product matrix
456  if (rb_evaluation.compute_RB_inner_product)
457  {
458  auto rb_inner_product_matrix_list =
459  rb_evaluation_reader.getRbInnerProductMatrix();
460 
461  if (rb_inner_product_matrix_list.size() != n_bfs*n_bfs)
462  libmesh_error_msg("Size error while reading the inner product matrix.");
463 
464  for (unsigned int i=0; i<n_bfs; ++i)
465  for (unsigned int j=0; j<n_bfs; ++j)
466  {
467  unsigned int offset = i*n_bfs + j;
468  rb_evaluation.RB_inner_product_matrix(i,j) =
469  load_scalar_value(rb_inner_product_matrix_list[offset]);
470  }
471  }
472 }
473 
474 template <typename RBEvaluationReaderNumber, typename TransRBEvaluationReaderNumber>
476  RBEvaluationReaderNumber & rb_eval_reader,
477  TransRBEvaluationReaderNumber & trans_rb_eval_reader,
478  bool read_error_bound_data)
479 {
480  load_rb_evaluation_data(trans_rb_eval,
481  rb_eval_reader,
482  read_error_bound_data);
483 
484  trans_rb_eval.set_delta_t( trans_rb_eval_reader.getDeltaT() );
485  trans_rb_eval.set_euler_theta( trans_rb_eval_reader.getEulerTheta() );
486  trans_rb_eval.set_n_time_steps( trans_rb_eval_reader.getNTimeSteps() );
487  trans_rb_eval.set_time_step( trans_rb_eval_reader.getTimeStep() );
488 
489  unsigned int n_bfs = rb_eval_reader.getNBfs();
490  unsigned int n_F_terms = trans_rb_eval.get_rb_theta_expansion().get_n_F_terms();
491  unsigned int n_A_terms = trans_rb_eval.get_rb_theta_expansion().get_n_A_terms();
492 
493  TransientRBThetaExpansion & trans_theta_expansion =
494  cast_ref<TransientRBThetaExpansion &>(trans_rb_eval.get_rb_theta_expansion());
495  unsigned int n_M_terms = trans_theta_expansion.get_n_M_terms();
496 
497  // L2 matrix
498  {
499  auto rb_L2_matrix_list =
500  trans_rb_eval_reader.getRbL2Matrix();
501 
502  if (rb_L2_matrix_list.size() != n_bfs*n_bfs)
503  libmesh_error_msg("Size error while reading the L2 matrix.");
504 
505  for (unsigned int i=0; i<n_bfs; ++i)
506  for (unsigned int j=0; j<n_bfs; ++j)
507  {
508  unsigned int offset = i*n_bfs + j;
509  trans_rb_eval.RB_L2_matrix(i,j) =
510  load_scalar_value(rb_L2_matrix_list[offset]);
511  }
512  }
513 
514  // Mq matrices
515  {
516  auto rb_Mq_matrices_outer_list = trans_rb_eval_reader.getRbMqMatrices();
517 
518  if (rb_Mq_matrices_outer_list.size() != n_M_terms)
519  libmesh_error_msg("Incorrect number of Mq matrices detected in the buffer");
520 
521  for (unsigned int q_m=0; q_m < n_M_terms; ++q_m)
522  {
523  auto rb_Mq_matrices_inner_list = rb_Mq_matrices_outer_list[q_m];
524  if (rb_Mq_matrices_inner_list.size() != n_bfs*n_bfs)
525  libmesh_error_msg("Incorrect Mq matrix size detected in the buffer");
526 
527  for (unsigned int i=0; i<n_bfs; ++i)
528  for (unsigned int j=0; j<n_bfs; ++j)
529  {
530  unsigned int offset = i*n_bfs+j;
531  trans_rb_eval.RB_M_q_vector[q_m](i,j) =
532  load_scalar_value(rb_Mq_matrices_inner_list[offset]);
533  }
534  }
535  }
536 
537  // The initial condition and L2 error at t=0.
538  {
539  auto initial_l2_errors_reader =
540  trans_rb_eval_reader.getInitialL2Errors();
541  if (initial_l2_errors_reader.size() != n_bfs)
542  libmesh_error_msg("Incorrect number of initial L2 error terms detected in the buffer");
543 
544  auto initial_conditions_outer_list =
545  trans_rb_eval_reader.getInitialConditions();
546  if (initial_conditions_outer_list.size() != n_bfs)
547  libmesh_error_msg("Incorrect number of outer initial conditions detected in the buffer");
548 
549  for (unsigned int i=0; i<n_bfs; i++)
550  {
551  trans_rb_eval.initial_L2_error_all_N[i] =
552  initial_l2_errors_reader[i];
553 
554  auto initial_conditions_inner_list = initial_conditions_outer_list[i];
555  if (initial_conditions_inner_list.size() != (i+1))
556  libmesh_error_msg("Incorrect number of inner initial conditions detected in the buffer");
557 
558  for (unsigned int j=0; j<=i; j++)
559  {
560  trans_rb_eval.RB_initial_condition_all_N[i](j) =
561  load_scalar_value(initial_conditions_inner_list[j]);
562  }
563  }
564  }
565 
566 
567  if (read_error_bound_data)
568  {
569  // Fq_Mq data
570  {
571  auto fq_mq_innerprods_list = trans_rb_eval_reader.getFqMqInnerprods();
572  if (fq_mq_innerprods_list.size() != n_F_terms*n_M_terms*n_bfs)
573  libmesh_error_msg("Size error while reading Fq-Mq representor data from buffer.");
574 
575  for (unsigned int q_f=0; q_f<n_F_terms; ++q_f)
576  for (unsigned int q_m=0; q_m<n_M_terms; ++q_m)
577  for (unsigned int i=0; i<n_bfs; ++i)
578  {
579  unsigned int offset = q_f*n_M_terms*n_bfs + q_m*n_bfs + i;
580  trans_rb_eval.Fq_Mq_representor_innerprods[q_f][q_m][i] =
581  load_scalar_value(fq_mq_innerprods_list[offset]);
582  }
583  }
584 
585 
586  // Mq_Mq representor inner-product data
587  {
588  unsigned int Q_m_hat = n_M_terms*(n_M_terms+1)/2;
589  auto mq_mq_innerprods_list = trans_rb_eval_reader.getMqMqInnerprods();
590  if (mq_mq_innerprods_list.size() != Q_m_hat*n_bfs*n_bfs)
591  libmesh_error_msg("Size error while reading Mq-Mq representor data from buffer.");
592 
593  for (unsigned int i=0; i<Q_m_hat; ++i)
594  for (unsigned int j=0; j<n_bfs; ++j)
595  for (unsigned int l=0; l<n_bfs; ++l)
596  {
597  unsigned int offset = i*n_bfs*n_bfs + j*n_bfs + l;
598  trans_rb_eval.Mq_Mq_representor_innerprods[i][j][l] =
599  load_scalar_value(mq_mq_innerprods_list[offset]);
600  }
601  }
602 
603  // Aq_Mq representor inner-product data
604  {
605  auto aq_mq_innerprods_list =
606  trans_rb_eval_reader.getAqMqInnerprods();
607  if (aq_mq_innerprods_list.size() != n_A_terms*n_M_terms*n_bfs*n_bfs)
608  libmesh_error_msg("Size error while reading Aq-Mq representor data from buffer.");
609 
610  for (unsigned int q_a=0; q_a<n_A_terms; q_a++)
611  for (unsigned int q_m=0; q_m<n_M_terms; q_m++)
612  for (unsigned int i=0; i<n_bfs; i++)
613  for (unsigned int j=0; j<n_bfs; j++)
614  {
615  unsigned int offset =
616  q_a*(n_M_terms*n_bfs*n_bfs) + q_m*(n_bfs*n_bfs) + i*n_bfs + j;
617 
618  trans_rb_eval.Aq_Mq_representor_innerprods[q_a][q_m][i][j] =
619  load_scalar_value(aq_mq_innerprods_list[offset]);
620  }
621  }
622 
623  }
624 }
625 
626 template <typename RBEvaluationReaderNumber, typename RBEIMEvaluationReaderNumber>
628  RBEvaluationReaderNumber & rb_evaluation_reader,
629  RBEIMEvaluationReaderNumber & rb_eim_evaluation_reader)
630 {
631  // We use read_error_bound_data=false here, since the RBEvaluation error bound data
632  // is not relevant to EIM.
633  load_rb_evaluation_data(rb_eim_evaluation,
634  rb_evaluation_reader,
635  /*read_error_bound_data*/ false);
636 
637  unsigned int n_bfs = rb_eim_evaluation.get_n_basis_functions();
638 
639  // EIM interpolation matrix
640  {
641  auto interpolation_matrix_list = rb_eim_evaluation_reader.getInterpolationMatrix();
642 
643  if (interpolation_matrix_list.size() != n_bfs*(n_bfs+1)/2)
644  libmesh_error_msg("Size error while reading the eim inner product matrix.");
645 
646  for (unsigned int i=0; i<n_bfs; ++i)
647  for (unsigned int j=0; j<=i; ++j)
648  {
649  unsigned int offset = i*(i+1)/2 + j;
650  rb_eim_evaluation.interpolation_matrix(i,j) =
651  load_scalar_value(interpolation_matrix_list[offset]);
652  }
653  }
654 
655  // Interpolation points
656  {
657  auto interpolation_points_list =
658  rb_eim_evaluation_reader.getInterpolationPoints();
659 
660  if (interpolation_points_list.size() != n_bfs)
661  libmesh_error_msg("Size error while reading the eim interpolation points.");
662 
663  rb_eim_evaluation.interpolation_points.resize(n_bfs);
664  for (unsigned int i=0; i<n_bfs; ++i)
665  {
666  load_point(interpolation_points_list[i],
667  rb_eim_evaluation.interpolation_points[i]);
668  }
669  }
670 
671  // Interpolation points variables
672  {
673  auto interpolation_points_var_list =
674  rb_eim_evaluation_reader.getInterpolationPointsVar();
675  rb_eim_evaluation.interpolation_points_var.resize(n_bfs);
676 
677  if (interpolation_points_var_list.size() != n_bfs)
678  libmesh_error_msg("Size error while reading the eim interpolation variables.");
679 
680  for (unsigned int i=0; i<n_bfs; ++i)
681  {
682  rb_eim_evaluation.interpolation_points_var[i] =
683  interpolation_points_var_list[i];
684  }
685  }
686 
687  // Interpolation elements
688  {
689  libMesh::dof_id_type elem_id = 0;
690  libMesh::ReplicatedMesh & interpolation_points_mesh =
691  rb_eim_evaluation.get_interpolation_points_mesh();
692  interpolation_points_mesh.clear();
693 
694  auto interpolation_points_elem_list =
695  rb_eim_evaluation_reader.getInterpolationPointsElems();
696  unsigned int n_interpolation_elems = interpolation_points_elem_list.size();
697 
698  if (n_interpolation_elems != n_bfs)
699  libmesh_error_msg("The number of elements should match the number of basis functions");
700 
701  rb_eim_evaluation.interpolation_points_elem.resize(n_interpolation_elems);
702 
703  for (unsigned int i=0; i<n_interpolation_elems; ++i)
704  {
705  auto mesh_elem_reader = interpolation_points_elem_list[i];
706  std::string elem_type_name = mesh_elem_reader.getType().cStr();
707  libMesh::ElemType elem_type =
708  libMesh::Utility::string_to_enum<libMesh::ElemType>(elem_type_name);
709 
710  libMesh::Elem * elem = libMesh::Elem::build(elem_type).release();
711  elem->set_id(elem_id++);
712  load_elem_into_mesh(mesh_elem_reader, elem, interpolation_points_mesh);
713 
714  rb_eim_evaluation.interpolation_points_elem[i] = elem;
715  }
716  }
717 }
718 
719 #if defined(LIBMESH_HAVE_SLEPC) && (LIBMESH_HAVE_GLPK)
721  RBData::RBSCMEvaluation::Reader & rb_scm_evaluation_reader)
722 {
723  auto parameter_ranges =
724  rb_scm_evaluation_reader.getParameterRanges();
725  auto discrete_parameters_list =
726  rb_scm_evaluation_reader.getDiscreteParameters();
727  load_parameter_ranges(rb_scm_eval,
728  parameter_ranges,
729  discrete_parameters_list);
730 
731  unsigned int n_A_terms = rb_scm_eval.get_rb_theta_expansion().get_n_A_terms();
732 
733  {
734  auto b_min_list = rb_scm_evaluation_reader.getBMin();
735 
736  if (b_min_list.size() != n_A_terms)
737  libmesh_error_msg("Size error while reading B_min");
738 
739  rb_scm_eval.B_min.clear();
740  for (unsigned int i=0; i<n_A_terms; i++)
741  rb_scm_eval.B_min.push_back(b_min_list[i]);
742  }
743 
744  {
745  auto b_max_list = rb_scm_evaluation_reader.getBMax();
746 
747  if (b_max_list.size() != n_A_terms)
748  libmesh_error_msg("Size error while reading B_max");
749 
750  rb_scm_eval.B_max.clear();
751  for (unsigned int i=0; i<n_A_terms; i++)
752  rb_scm_eval.B_max.push_back(b_max_list[i]);
753  }
754 
755  {
756  auto cJ_stability_vector =
757  rb_scm_evaluation_reader.getCJStabilityVector();
758 
759  rb_scm_eval.C_J_stability_vector.clear();
760  for (unsigned int i=0; i<cJ_stability_vector.size(); i++)
761  rb_scm_eval.C_J_stability_vector.push_back( cJ_stability_vector[i] );
762  }
763 
764  {
765  auto cJ_parameters_outer =
766  rb_scm_evaluation_reader.getCJ();
767 
768  rb_scm_eval.C_J.resize(cJ_parameters_outer.size());
769  for (unsigned int i=0; i<cJ_parameters_outer.size(); i++)
770  {
771  auto cJ_parameters_inner =
772  cJ_parameters_outer[i];
773 
774  for (unsigned int j=0; j<cJ_parameters_inner.size(); j++)
775  {
776  std::string param_name = cJ_parameters_inner[j].getName();
777  Real param_value = cJ_parameters_inner[j].getValue();
778  rb_scm_eval.C_J[i].set_value(param_name, param_value);
779  }
780  }
781  }
782 
783  {
784  auto scm_ub_vectors =
785  rb_scm_evaluation_reader.getScmUbVectors();
786 
787  // The number of UB vectors is the same as the number of C_J values
788  unsigned int n_C_J_values = rb_scm_evaluation_reader.getCJ().size();
789 
790  if (scm_ub_vectors.size() != n_C_J_values*n_A_terms)
791  libmesh_error_msg("Size mismatch in SCB UB vectors");
792 
793  rb_scm_eval.SCM_UB_vectors.resize( n_C_J_values );
794  for (unsigned int i=0; i<n_C_J_values; i++)
795  {
796  rb_scm_eval.SCM_UB_vectors[i].resize(n_A_terms);
797  for (unsigned int j=0; j<n_A_terms; j++)
798  {
799  unsigned int offset = i*n_A_terms + j;
800  rb_scm_eval.SCM_UB_vectors[i][j] = scm_ub_vectors[offset];
801 
802  }
803  }
804  }
805 
806 }
807 #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
808 
809 
810 
811 void load_point(RBData::Point3D::Reader point_reader, Point & point)
812 {
813  point(0) = point_reader.getX();
814 
815  if (LIBMESH_DIM >= 2)
816  point(1) = point_reader.getY();
817 
818  if (LIBMESH_DIM >= 3)
819  point(2) = point_reader.getZ();
820 }
821 
822 
823 void load_elem_into_mesh(RBData::MeshElem::Reader mesh_elem_reader,
824  libMesh::Elem * elem,
826 {
827  auto mesh_elem_point_list = mesh_elem_reader.getPoints();
828  unsigned int n_points = mesh_elem_point_list.size();
829 
830  if (n_points != elem->n_nodes())
831  libmesh_error_msg("Wrong number of nodes for element type");
832 
833  for (unsigned int i=0; i < n_points; ++i)
834  {
835  libMesh::Node * node = new libMesh::Node(mesh_elem_point_list[i].getX(),
836  mesh_elem_point_list[i].getY(),
837  mesh_elem_point_list[i].getZ());
838 
839  mesh.add_node(node);
840 
841  elem->set_node(i) = node;
842  }
843 
844  elem->subdomain_id() = mesh_elem_reader.getSubdomainId();
845 
846  mesh.add_elem(elem);
847 }
848 
849 // ---- Helper functions for adding data to capnp Builders (END) ----
850 
851 } // namespace RBDataSerialization
852 
853 } // namespace libMesh
854 
855 #endif // #if defined(LIBMESH_HAVE_CAPNPROTO)
std::vector< Real > B_min
B_min, B_max define the bounding box.
std::vector< std::vector< std::vector< Number > > > Mq_Mq_representor_innerprods
void load_rb_eim_evaluation_data(RBEIMEvaluation &rb_eim_eval, RBEvaluationReaderNumber &rb_evaluation_reader, RBEIMEvaluationReaderNumber &rb_eim_eval_reader)
Load an EIM RB evaluation from a corresponding reader structure in the buffer.
DenseMatrix< Number > interpolation_matrix
Dense matrix that stores the lower triangular interpolation matrix that can be used.
bool compute_RB_inner_product
Boolean flag to indicate whether we compute the RB_inner_product_matrix.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
RBSCMEvaluation & _rb_scm_eval
The RBSCMEvaluation object we will read into.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:238
virtual void resize_data_structures(const unsigned int Nmax, bool resize_error_bound_data=true)
Resize and clear the data vectors corresponding to the value of Nmax.
virtual Elem * add_elem(Elem *e) libmesh_override
Add elem e to the end of the element array.
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.
void load_rb_evaluation_data(RBEvaluation &rb_evaluation, RBEvaluationReaderNumber &rb_evaluation_reader, bool read_error_bound_data)
Load an RB evaluation from a corresponding reader structure in the buffer.
DenseMatrix< Number > RB_inner_product_matrix
The inner product matrix.
ElemType
Defines an enum for geometric element types.
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
virtual void set_n_basis_functions(unsigned int n_bfs)
Set the number of basis functions.
virtual Node * add_node(Node *n) libmesh_override
Add Node n to the end of the vertex array.
ReplicatedMesh & get_interpolation_points_mesh()
Get a writable reference to the interpolation points mesh.
virtual void clear() libmesh_override
Clear all internal data.
void load_rb_scm_evaluation_data(RBSCMEvaluation &rb_scm_eval, RBData::RBSCMEvaluation::Reader &rb_scm_eval_reader)
Load an SCM RB evaluation from a corresponding reader structure in the buffer.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
This class is part of the rbOOmit framework.
void load_transient_rb_evaluation_data(TransientRBEvaluation &trans_rb_eval, RBEvaluationReaderNumber &rb_evaluation_reader, TransRBEvaluationReaderNumber &trans_rb_eval_reader, bool read_error_bound_data)
Load an RB evaluation from a corresponding reader structure in the buffer.
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
The libMesh namespace provides an interface to certain functionality in the library.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the rb_theta_expansion.
std::vector< Real > B_max
long double max(long double a, double b)
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
std::vector< std::vector< Real > > SCM_UB_vectors
This matrix stores the infimizing vectors y_1( ),...,y_Q_a( ), for each in C_J, which are used in co...
std::vector< Real > C_J_stability_vector
Vector storing the (truth) stability values at the parameters in C_J.
virtual unsigned int n_nodes() const =0
RBSCMEvaluationDeserialization(RBSCMEvaluation &trans_rb_eval)
Initialize a new buffer using the structure from the Cap&#39;n&#39;Proto schema described in rb_data...
std::vector< DenseVector< Number > > RB_Fq_vector
Dense vector for the RHS.
RBEvaluationDeserialization(RBEvaluation &rb_eval)
Initialize a new buffer using the structure from the Cap&#39;n&#39;Proto schema described in rb_data...
std::vector< Real > initial_L2_error_all_N
Vector storing initial L2 error for all 1 <= N <= RB_size.
std::vector< RBParameters > C_J
Vector storing the greedily selected parameters during SCM training.
void load_parameter_ranges(RBParametrized &rb_evaluation, RBData::ParameterRanges::Reader &parameter_ranges, RBData::DiscreteParameterList::Reader &discrete_parameters_list)
Load parameter ranges and discrete parameter values into an RBEvaluation from the corresponding struc...
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...
void read_from_file(const std::string &path, bool read_error_bound_data)
Write the Cap&#39;n&#39;Proto buffer to disk.
RBEIMEvaluation & _rb_eim_eval
The RBEIMEvaluation object we will read into.
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
std::vector< unsigned int > interpolation_points_var
The corresponding list of variables indices at which the interpolation points were identified...
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
void read_from_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
std::vector< Elem * > interpolation_points_elem
The corresponding list of elements at which the interpolation points were identified.
virtual unsigned int get_n_M_terms()
Get Q_m, the number of terms in the affine expansion for the mass operator.
RBEIMEvaluationDeserialization(RBEIMEvaluation &trans_rb_eval)
Initialize a new buffer using the structure from the Cap&#39;n&#39;Proto schema described in rb_data...
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...
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.
This class is part of the rbOOmit framework.
Definition: rb_evaluation.h:50
std::vector< DenseMatrix< Number > > RB_Aq_vector
Dense matrices for the RB computations.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void load_point(RBData::Point3D::Reader point_reader, Point &point)
Helper function that loads point data.
void set_value(const std::string &param_name, Real value)
Set the value of the specified parameter.
Definition: rb_parameters.C:57
This class is part of the rbOOmit framework.
void read_from_file(const std::string &path, bool read_error_bound_data)
Write the Cap&#39;n&#39;Proto buffer to disk.
static const bool value
Definition: xdr_io.C:108
void set_euler_theta(const Real euler_theta_in)
TransientRBEvaluation & _trans_rb_eval
The TransientRBEvaluation object that we will read into.
RBEvaluation & _rb_eval
The RBEvaluation object that we will read into.
void load_elem_into_mesh(RBData::MeshElem::Reader mesh_elem_reader, libMesh::Elem *elem, libMesh::ReplicatedMesh &mesh)
Helper function that loads element data.
std::vector< std::vector< std::vector< std::vector< Number > > > > Aq_Mq_representor_innerprods
void set_n_time_steps(const unsigned int K)
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
TransientRBEvaluationDeserialization(TransientRBEvaluation &trans_rb_eval)
Initialize a new buffer using the structure from the Cap&#39;n&#39;Proto schema described in rb_data...
void initialize_parameters(const RBParameters &mu_min_in, const RBParameters &mu_max_in, const std::map< std::string, std::vector< Real >> &discrete_parameter_values)
Initialize the parameter ranges and set current_parameters.
std::vector< Point > interpolation_points
The list of interpolation points, i.e.
This class is part of the rbOOmit framework.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
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.
void read_from_file(const std::string &path)
Write the Cap&#39;n&#39;Proto buffer to disk.
This class is part of the rbOOmit framework.
std::vector< std::vector< std::vector< Number > > > Aq_Aq_representor_innerprods
uint8_t dof_id_type
Definition: id_types.h:64
std::vector< std::vector< DenseVector< Number > > > RB_output_vectors
The vectors storing the RB output vectors.