Line data Source code
1 : // rbOOmit: An implementation of the Certified Reduced Basis method.
2 : // Copyright (C) 2009, 2010 David J. Knezevic
3 :
4 : // This file is part of rbOOmit.
5 :
6 : // rbOOmit is free software; you can redistribute it and/or
7 : // modify it under the terms of the GNU Lesser General Public
8 : // License as published by the Free Software Foundation; either
9 : // version 2.1 of the License, or (at your option) any later version.
10 :
11 : // rbOOmit is distributed in the hope that it will be useful,
12 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 : // Lesser General Public License for more details.
15 :
16 : // You should have received a copy of the GNU Lesser General Public
17 : // License along with this library; if not, write to the Free Software
18 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 :
20 : // rbOOmit includes
21 : #include "libmesh/rb_evaluation.h"
22 : #include "libmesh/rb_theta_expansion.h"
23 :
24 : // libMesh includes
25 : #include "libmesh/libmesh_common.h"
26 : #include "libmesh/libmesh_version.h"
27 : #include "libmesh/system.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/libmesh_logging.h"
30 : #include "libmesh/xdr_cxx.h"
31 : #include "libmesh/mesh_tools.h"
32 : #include "libmesh/utility.h"
33 :
34 : // TIMPI includes
35 : #include "timpi/communicator.h"
36 :
37 : // C/C++ includes
38 : #include <sys/types.h>
39 : #include <sys/stat.h>
40 : #include <errno.h>
41 : #include <fstream>
42 : #include <sstream>
43 : #include <iomanip>
44 :
45 : namespace libMesh
46 : {
47 :
48 54 : RBEvaluation::RBEvaluation (const Parallel::Communicator & comm_in)
49 : :
50 : ParallelObject(comm_in),
51 22 : evaluate_RB_error_bound(true),
52 22 : compute_RB_inner_product(false),
53 54 : rb_theta_expansion(nullptr)
54 : {
55 54 : }
56 :
57 66 : RBEvaluation::~RBEvaluation() = default;
58 :
59 0 : void RBEvaluation::clear()
60 : {
61 0 : LOG_SCOPE("clear()", "RBEvaluation");
62 :
63 : // Clear the basis functions
64 0 : basis_functions.clear();
65 0 : set_n_basis_functions(0);
66 :
67 0 : clear_riesz_representors();
68 :
69 : // Clear the Greedy param list
70 0 : for (auto & plist : greedy_param_list)
71 0 : plist.clear();
72 0 : greedy_param_list.clear();
73 0 : }
74 :
75 27 : void RBEvaluation::set_n_basis_functions(unsigned int n_bfs)
76 : {
77 27 : basis_functions.resize(n_bfs);
78 27 : }
79 :
80 54 : void RBEvaluation::set_rb_theta_expansion(RBThetaExpansion & rb_theta_expansion_in)
81 : {
82 54 : rb_theta_expansion = &rb_theta_expansion_in;
83 54 : }
84 :
85 544316 : RBThetaExpansion & RBEvaluation::get_rb_theta_expansion()
86 : {
87 544316 : libmesh_error_msg_if(!is_rb_theta_expansion_initialized(),
88 : "Error: rb_theta_expansion hasn't been initialized yet");
89 :
90 544316 : return *rb_theta_expansion;
91 : }
92 :
93 218 : const RBThetaExpansion & RBEvaluation::get_rb_theta_expansion() const
94 : {
95 218 : libmesh_error_msg_if(!is_rb_theta_expansion_initialized(),
96 : "Error: rb_theta_expansion hasn't been initialized yet");
97 :
98 218 : return *rb_theta_expansion;
99 : }
100 :
101 544534 : bool RBEvaluation::is_rb_theta_expansion_initialized() const
102 : {
103 544534 : if (rb_theta_expansion)
104 : {
105 178454 : return true;
106 : }
107 : else
108 : {
109 0 : return false;
110 : }
111 : }
112 :
113 54 : void RBEvaluation::resize_data_structures(const unsigned int Nmax,
114 : bool resize_error_bound_data)
115 : {
116 32 : LOG_SCOPE("resize_data_structures()", "RBEvaluation");
117 :
118 54 : libmesh_error_msg_if(Nmax < this->get_n_basis_functions(),
119 : "Error: Cannot set Nmax to be less than the current number of basis functions.");
120 :
121 : // Resize/clear inner product matrix
122 54 : if (compute_RB_inner_product)
123 8 : RB_inner_product_matrix.resize(Nmax,Nmax);
124 :
125 : // Allocate dense matrices for RB solves
126 54 : RB_Aq_vector.resize(rb_theta_expansion->get_n_A_terms());
127 :
128 216 : for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
129 : {
130 : // Initialize the memory for the RB matrices
131 162 : RB_Aq_vector[q].resize(Nmax,Nmax);
132 : }
133 :
134 54 : RB_Fq_vector.resize(rb_theta_expansion->get_n_F_terms());
135 :
136 248 : for (unsigned int q=0; q<rb_theta_expansion->get_n_F_terms(); q++)
137 : {
138 : // Initialize the memory for the RB vectors
139 194 : RB_Fq_vector[q].resize(Nmax);
140 : }
141 :
142 :
143 : // Initialize the RB output vectors
144 54 : RB_output_vectors.resize(rb_theta_expansion->get_n_outputs());
145 158 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
146 : {
147 136 : RB_output_vectors[n].resize(rb_theta_expansion->get_n_output_terms(n));
148 208 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
149 : {
150 136 : RB_output_vectors[n][q_l].resize(Nmax);
151 : }
152 : }
153 :
154 : // Initialize vectors storing output data
155 54 : RB_outputs.resize(rb_theta_expansion->get_n_outputs(), 0.);
156 :
157 :
158 54 : if (resize_error_bound_data)
159 : {
160 : // Initialize vectors for the norms of the Fq representors
161 54 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
162 54 : Fq_representor_innerprods.resize(Q_f_hat);
163 :
164 : // Initialize vectors for the norms of the representors
165 54 : Fq_Aq_representor_innerprods.resize(rb_theta_expansion->get_n_F_terms());
166 248 : for (unsigned int i=0; i<rb_theta_expansion->get_n_F_terms(); i++)
167 : {
168 250 : Fq_Aq_representor_innerprods[i].resize(rb_theta_expansion->get_n_A_terms());
169 776 : for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
170 : {
171 918 : Fq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
172 : }
173 : }
174 :
175 54 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
176 54 : Aq_Aq_representor_innerprods.resize(Q_a_hat);
177 378 : for (unsigned int i=0; i<Q_a_hat; i++)
178 : {
179 420 : Aq_Aq_representor_innerprods[i].resize(Nmax);
180 5964 : for (unsigned int j=0; j<Nmax; j++)
181 : {
182 9000 : Aq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
183 : }
184 : }
185 :
186 54 : RB_output_error_bounds.resize(rb_theta_expansion->get_n_outputs(), 0.);
187 :
188 : // Resize the output dual norm vectors
189 54 : output_dual_innerprods.resize(rb_theta_expansion->get_n_outputs());
190 158 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
191 : {
192 104 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
193 136 : output_dual_innerprods[n].resize(Q_l_hat);
194 : }
195 :
196 : // Clear and resize the vector of Aq_representors
197 54 : clear_riesz_representors();
198 :
199 54 : Aq_representor.resize(rb_theta_expansion->get_n_A_terms());
200 216 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
201 : {
202 210 : Aq_representor[q_a].resize(Nmax);
203 : }
204 : }
205 54 : }
206 :
207 299526 : NumericVector<Number> & RBEvaluation::get_basis_function(unsigned int i)
208 : {
209 97718 : libmesh_assert_less (i, basis_functions.size());
210 :
211 397244 : return *(basis_functions[i]);
212 : }
213 :
214 5600 : const NumericVector<Number> & RBEvaluation::get_basis_function(unsigned int i) const
215 : {
216 1600 : libmesh_assert_less (i, basis_functions.size());
217 :
218 7200 : return *(basis_functions[i]);
219 : }
220 :
221 128021 : Real RBEvaluation::rb_solve(unsigned int N)
222 : {
223 128021 : return rb_solve(N, nullptr);
224 : }
225 :
226 128021 : Real RBEvaluation::rb_solve(unsigned int N,
227 : const std::vector<Number> * evaluated_thetas)
228 : {
229 64012 : LOG_SCOPE("rb_solve()", "RBEvaluation");
230 :
231 128021 : libmesh_error_msg_if(N > get_n_basis_functions(),
232 : "ERROR: N cannot be larger than the number of basis functions in rb_solve");
233 :
234 : // In case the theta functions have been pre-evaluated, first check the size for consistency. The
235 : // size of the input "evaluated_thetas" vector must match the sum of the "A", "F", and "output" terms
236 : // in the expansion.
237 128021 : this->check_evaluated_thetas_size(evaluated_thetas);
238 :
239 128021 : const RBParameters & mu = get_parameters();
240 :
241 : // Resize (and clear) the solution vector
242 96015 : RB_solution.resize(N);
243 :
244 : // Assemble the RB system
245 192033 : DenseMatrix<Number> RB_system_matrix(N,N);
246 32006 : RB_system_matrix.zero();
247 :
248 192033 : DenseMatrix<Number> RB_Aq_a;
249 512084 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
250 : {
251 480081 : RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);
252 :
253 384063 : if (evaluated_thetas)
254 0 : RB_system_matrix.add((*evaluated_thetas)[q_a], RB_Aq_a);
255 : else
256 384063 : RB_system_matrix.add(rb_theta_expansion->eval_A_theta(q_a, mu), RB_Aq_a);
257 : }
258 :
259 : // Assemble the RB rhs
260 128021 : DenseVector<Number> RB_rhs(N);
261 32006 : RB_rhs.zero();
262 :
263 128021 : DenseVector<Number> RB_Fq_f;
264 856112 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
265 : {
266 910117 : RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);
267 :
268 728091 : if (evaluated_thetas)
269 0 : RB_rhs.add((*evaluated_thetas)[q_f+rb_theta_expansion->get_n_A_terms()], RB_Fq_f);
270 : else
271 728091 : RB_rhs.add(rb_theta_expansion->eval_F_theta(q_f, mu), RB_Fq_f);
272 : }
273 :
274 : // Solve the linear system
275 128021 : if (N > 0)
276 : {
277 119621 : RB_system_matrix.lu_solve(RB_rhs, RB_solution);
278 : }
279 :
280 : // Place to store the output of get_principal_subvector() calls
281 128021 : DenseVector<Number> RB_output_vector_N;
282 :
283 : // Evaluate RB outputs
284 32006 : unsigned int output_counter = 0;
285 160049 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
286 : {
287 40036 : RB_outputs[n] = 0.;
288 64056 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
289 : {
290 48044 : RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
291 :
292 : // Compute dot product with current output vector and RB_solution
293 32028 : auto dot_prod = RB_output_vector_N.dot(RB_solution);
294 :
295 : // Determine the coefficient depending on whether or not
296 : // pre-evaluated thetas were provided. Note that if
297 : // pre-evaluated thetas were provided, they must come after
298 : // the "A" and "F" thetas and be in "row-major" order. In
299 : // other words, there is a possibly "ragged" 2D array of
300 : // output thetas ordered by output index "n" and term index
301 : // "q_l" which is accessed in row-major order.
302 32028 : auto coeff = evaluated_thetas ?
303 0 : (*evaluated_thetas)[output_counter + rb_theta_expansion->get_n_A_terms() + rb_theta_expansion->get_n_F_terms()] :
304 32028 : rb_theta_expansion->eval_output_theta(n, q_l, mu);
305 :
306 : // Finally, accumulate the result in RB_outputs[n]
307 24024 : RB_outputs[n] += coeff * dot_prod;
308 :
309 : // Go to next output
310 32028 : output_counter++;
311 : }
312 : }
313 :
314 128021 : if (evaluate_RB_error_bound) // Calculate the error bounds
315 : {
316 : // Evaluate the dual norm of the residual for RB_solution_vector
317 128021 : Real epsilon_N = compute_residual_dual_norm(N, evaluated_thetas);
318 :
319 : // Get lower bound for coercivity constant
320 128021 : const Real alpha_LB = get_stability_lower_bound();
321 : // alpha_LB needs to be positive to get a valid error bound
322 32006 : libmesh_assert_greater ( alpha_LB, 0. );
323 :
324 : // Evaluate the (absolute) error bound
325 128021 : Real abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB);
326 :
327 : // Now compute the output error bounds
328 160049 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
329 32028 : RB_output_error_bounds[n] = abs_error_bound * this->eval_output_dual_norm(n, evaluated_thetas);
330 :
331 32006 : return abs_error_bound;
332 : }
333 : else // Don't calculate the error bounds
334 : {
335 : // Just return -1. if we did not compute the error bound
336 0 : return -1.;
337 : }
338 64009 : }
339 :
340 0 : Real RBEvaluation::get_error_bound_normalization()
341 : {
342 : // Normalize the error based on the error bound in the
343 : // case of an empty reduced basis. The error bound is based
344 : // on the residual F - AU, so with an empty basis this gives
345 : // a value based on the norm of F at the current parameters.
346 :
347 0 : Real normalization = rb_solve(0);
348 0 : return normalization;
349 : }
350 :
351 0 : Real RBEvaluation::compute_residual_dual_norm(const unsigned int N)
352 : {
353 0 : return compute_residual_dual_norm(N, nullptr);
354 : }
355 :
356 128021 : Real RBEvaluation::compute_residual_dual_norm(const unsigned int N,
357 : const std::vector<Number> * evaluated_thetas)
358 : {
359 64012 : LOG_SCOPE("compute_residual_dual_norm()", "RBEvaluation");
360 :
361 : // In case the theta functions have been pre-evaluated, first check the size for consistency
362 128021 : this->check_evaluated_thetas_size(evaluated_thetas);
363 :
364 : // If evaluated_thetas is provided, then mu is not actually used for anything
365 128021 : const RBParameters & mu = get_parameters();
366 :
367 128021 : const unsigned int n_A_terms = rb_theta_expansion->get_n_A_terms();
368 128021 : const unsigned int n_F_terms = rb_theta_expansion->get_n_F_terms();
369 :
370 : // Use the stored representor inner product values
371 : // to evaluate the residual norm
372 64009 : Number residual_norm_sq = 0.;
373 :
374 : // Lambdas to help with evaluating F_theta and A_theta functions
375 : // using either the pre-evaluated thetas (if provided) or by calling
376 : // eval_{F,A}_theta()
377 3984483 : auto eval_F = [&](unsigned int index)
378 : {
379 3984483 : return (evaluated_thetas) ? (*evaluated_thetas)[index + n_A_terms] : rb_theta_expansion->eval_F_theta(index, mu);
380 128021 : };
381 3336462 : auto eval_A = [&](unsigned int index)
382 : {
383 3336462 : return (evaluated_thetas) ? (*evaluated_thetas)[index] : rb_theta_expansion->eval_A_theta(index, mu);
384 128021 : };
385 :
386 32006 : unsigned int q=0;
387 856112 : for (unsigned int q_f1=0; q_f1<n_F_terms; q_f1++)
388 : {
389 728091 : const Number val_q_f1 = eval_F(q_f1);
390 :
391 3256392 : for (unsigned int q_f2=q_f1; q_f2<n_F_terms; q_f2++)
392 : {
393 2528301 : const Number val_q_f2 = eval_F(q_f2);
394 :
395 2528301 : Real delta = (q_f1==q_f2) ? 1. : 2.;
396 2528301 : residual_norm_sq += delta * libmesh_real(val_q_f1 * libmesh_conj(val_q_f2) * Fq_representor_innerprods[q] );
397 :
398 2528301 : q++;
399 : }
400 : }
401 :
402 856112 : for (unsigned int q_f=0; q_f<n_F_terms; q_f++)
403 : {
404 728091 : const Number val_q_f = eval_F(q_f);
405 :
406 2912364 : for (unsigned int q_a=0; q_a<n_A_terms; q_a++)
407 : {
408 2184273 : const Number val_q_a = eval_A(q_a);
409 :
410 17536473 : for (unsigned int i=0; i<N; i++)
411 : {
412 3838200 : Real delta = 2.;
413 11514600 : residual_norm_sq +=
414 19191000 : delta * libmesh_real( val_q_f * libmesh_conj(val_q_a) *
415 23028600 : libmesh_conj(RB_solution(i)) * Fq_Aq_representor_innerprods[q_f][q_a][i] );
416 : }
417 : }
418 : }
419 :
420 32006 : q=0;
421 512084 : for (unsigned int q_a1=0; q_a1<n_A_terms; q_a1++)
422 : {
423 384063 : const Number val_q_a1 = eval_A(q_a1);
424 :
425 1152189 : for (unsigned int q_a2=q_a1; q_a2<n_A_terms; q_a2++)
426 : {
427 768126 : const Number val_q_a2 = eval_A(q_a2);
428 :
429 768126 : Real delta = (q_a1==q_a2) ? 1. : 2.;
430 :
431 6266226 : for (unsigned int i=0; i<N; i++)
432 : {
433 60181800 : for (unsigned int j=0; j<N; j++)
434 : {
435 41016600 : residual_norm_sq +=
436 68361000 : delta * libmesh_real( libmesh_conj(val_q_a1) * val_q_a2 *
437 95700300 : libmesh_conj(RB_solution(i)) * RB_solution(j) * Aq_Aq_representor_innerprods[q][i][j] );
438 : }
439 : }
440 :
441 768126 : q++;
442 : }
443 : }
444 :
445 128021 : if (libmesh_real(residual_norm_sq) < 0.)
446 : {
447 : // libMesh::out << "Warning: Square of residual norm is negative "
448 : // << "in RBSystem::compute_residual_dual_norm()" << std::endl;
449 :
450 : // Sometimes this is negative due to rounding error,
451 : // but when this occurs the error is on the order of 1.e-10,
452 : // so shouldn't affect error bound much...
453 177 : residual_norm_sq = std::abs(residual_norm_sq);
454 : }
455 :
456 160027 : return std::sqrt( libmesh_real(residual_norm_sq) );
457 : }
458 :
459 0 : Real RBEvaluation::get_stability_lower_bound()
460 : {
461 : // Return a default value of 1, this function should
462 : // be overloaded to specify a problem-dependent stability
463 : // factor lower bound
464 0 : return 1.;
465 : }
466 :
467 728621 : Real RBEvaluation::residual_scaling_denom(Real alpha_LB)
468 : {
469 : // Here we implement the residual scaling for a coercive
470 : // problem.
471 728621 : return alpha_LB;
472 : }
473 :
474 2458452 : Real RBEvaluation::eval_output_dual_norm(unsigned int n,
475 : const std::vector<Number> * evaluated_thetas)
476 : {
477 : // Return value
478 824820 : Number output_bound_sq = 0.;
479 :
480 : // mu is only used if evaluated_thetas == nullptr
481 2458452 : const RBParameters & mu = this->get_parameters();
482 :
483 : // Index into output_dual_innerprods
484 816816 : unsigned int q=0;
485 4916904 : for (unsigned int q_l1=0; q_l1<rb_theta_expansion->get_n_output_terms(n); q_l1++)
486 : {
487 4916904 : for (unsigned int q_l2=q_l1; q_l2<rb_theta_expansion->get_n_output_terms(n); q_l2++)
488 : {
489 2458452 : Real delta = (q_l1==q_l2) ? 1. : 2.;
490 :
491 : Number val_l1 =
492 2458452 : evaluated_thetas ?
493 0 : (*evaluated_thetas)[rb_theta_expansion->output_index_1D(n, q_l1) + rb_theta_expansion->get_n_A_terms() + rb_theta_expansion->get_n_F_terms()] :
494 2458452 : rb_theta_expansion->eval_output_theta(n, q_l1, mu);
495 :
496 : Number val_l2 =
497 1633632 : evaluated_thetas ?
498 0 : (*evaluated_thetas)[rb_theta_expansion->output_index_1D(n, q_l2) + rb_theta_expansion->get_n_A_terms() + rb_theta_expansion->get_n_F_terms()] :
499 2458452 : rb_theta_expansion->eval_output_theta(n, q_l2, mu);
500 :
501 2450448 : output_bound_sq += delta * libmesh_real(
502 3275268 : libmesh_conj(val_l1) * val_l2 * output_dual_innerprods[n][q]);
503 :
504 2458452 : q++;
505 : }
506 : }
507 :
508 2458452 : return libmesh_real(std::sqrt( output_bound_sq ));
509 : }
510 :
511 54 : void RBEvaluation::clear_riesz_representors()
512 : {
513 38 : Aq_representor.clear();
514 54 : }
515 :
516 27 : void RBEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
517 : const bool write_binary_data)
518 : {
519 16 : LOG_SCOPE("legacy_write_offline_data_to_files()", "RBEvaluation");
520 :
521 : // Get the number of basis functions
522 27 : unsigned int n_bfs = get_n_basis_functions();
523 :
524 : // The writing mode: ENCODE for binary, WRITE for ASCII
525 27 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
526 :
527 : // The suffix to use for all the files that are written out
528 35 : const std::string suffix = write_binary_data ? ".xdr" : ".dat";
529 :
530 35 : if (this->processor_id() == 0)
531 : {
532 :
533 : // Make a directory to store all the data files
534 15 : Utility::mkdir(directory_name.c_str());
535 : // if (mkdir(directory_name.c_str(), 0777) == -1)
536 : // {
537 : // libMesh::out << "In RBEvaluation::write_offline_data_to_files, directory "
538 : // << directory_name << " already exists, overwriting contents." << std::endl;
539 : // }
540 :
541 : // First, write out how many basis functions we have generated
542 23 : std::ostringstream file_name;
543 : {
544 11 : file_name << directory_name << "/n_bfs" << suffix;
545 30 : Xdr n_bfs_out(file_name.str(), mode);
546 8 : n_bfs_out << n_bfs;
547 15 : n_bfs_out.close();
548 7 : }
549 :
550 : // Write out the parameter ranges
551 26 : file_name.str("");
552 11 : file_name << directory_name << "/parameter_ranges" << suffix;
553 8 : std::string continuous_param_file_name = file_name.str();
554 :
555 : // Write out the discrete parameter values
556 26 : file_name.str("");
557 11 : file_name << directory_name << "/discrete_parameter_values" << suffix;
558 8 : std::string discrete_param_file_name = file_name.str();
559 :
560 15 : write_parameter_data_to_files(continuous_param_file_name,
561 : discrete_param_file_name,
562 : write_binary_data);
563 :
564 : // Write out Fq representor norm data
565 26 : file_name.str("");
566 11 : file_name << directory_name << "/Fq_innerprods" << suffix;
567 23 : Xdr RB_Fq_innerprods_out(file_name.str(), mode);
568 15 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
569 190 : for (unsigned int i=0; i<Q_f_hat; i++)
570 : {
571 219 : RB_Fq_innerprods_out << Fq_representor_innerprods[i];
572 : }
573 15 : RB_Fq_innerprods_out.close();
574 :
575 : // Write out output data
576 43 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
577 : {
578 48 : file_name.str("");
579 28 : file_name << directory_name << "/output_";
580 : file_name << std::setw(3)
581 : << std::setprecision(0)
582 8 : << std::setfill('0')
583 8 : << std::right
584 8 : << n;
585 :
586 20 : file_name << "_dual_innerprods" << suffix;
587 44 : Xdr output_dual_innerprods_out(file_name.str(), mode);
588 :
589 28 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
590 56 : for (unsigned int q=0; q<Q_l_hat; q++)
591 : {
592 44 : output_dual_innerprods_out << output_dual_innerprods[n][q];
593 : }
594 28 : output_dual_innerprods_out.close();
595 12 : }
596 :
597 :
598 : // Write out output data to multiple files
599 43 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
600 : {
601 56 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
602 : {
603 48 : file_name.str("");
604 28 : file_name << directory_name << "/output_";
605 : file_name << std::setw(3)
606 : << std::setprecision(0)
607 8 : << std::setfill('0')
608 8 : << std::right
609 8 : << n;
610 28 : file_name << "_";
611 : file_name << std::setw(3)
612 : << std::setprecision(0)
613 8 : << std::setfill('0')
614 8 : << std::right
615 8 : << q_l;
616 8 : file_name << suffix;
617 52 : Xdr output_n_out(file_name.str(), mode);
618 :
619 588 : for (unsigned int j=0; j<n_bfs; j++)
620 : {
621 640 : output_n_out << RB_output_vectors[n][q_l](j);
622 : }
623 28 : output_n_out.close();
624 12 : }
625 : }
626 :
627 15 : if (compute_RB_inner_product)
628 : {
629 : // Next write out the inner product matrix
630 5 : file_name.str("");
631 2 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
632 6 : Xdr RB_inner_product_matrix_out(file_name.str(), mode);
633 63 : for (unsigned int i=0; i<n_bfs; i++)
634 : {
635 1260 : for (unsigned int j=0; j<n_bfs; j++)
636 : {
637 800 : RB_inner_product_matrix_out << RB_inner_product_matrix(i,j);
638 : }
639 : }
640 3 : RB_inner_product_matrix_out.close();
641 1 : }
642 :
643 : // Next write out the Fq vectors
644 70 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
645 : {
646 96 : file_name.str("");
647 55 : file_name << directory_name << "/RB_F_";
648 : file_name << std::setw(3)
649 : << std::setprecision(0)
650 14 : << std::setfill('0')
651 14 : << std::right
652 14 : << q_f;
653 14 : file_name << suffix;
654 97 : Xdr RB_Fq_f_out(file_name.str(), mode);
655 :
656 915 : for (unsigned int i=0; i<n_bfs; i++)
657 : {
658 660 : RB_Fq_f_out << RB_Fq_vector[q_f](i);
659 : }
660 55 : RB_Fq_f_out.close();
661 27 : }
662 :
663 : // Next write out the Aq matrices
664 60 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
665 : {
666 78 : file_name.str("");
667 45 : file_name << directory_name << "/RB_A_";
668 : file_name << std::setw(3)
669 : << std::setprecision(0)
670 12 : << std::setfill('0')
671 12 : << std::right
672 12 : << q_a;
673 12 : file_name << suffix;
674 81 : Xdr RB_Aq_a_out(file_name.str(), mode);
675 :
676 825 : for (unsigned int i=0; i<n_bfs; i++)
677 : {
678 14580 : for (unsigned int j=0; j<n_bfs; j++)
679 : {
680 11250 : RB_Aq_a_out << RB_Aq_vector[q_a](i,j);
681 : }
682 : }
683 45 : RB_Aq_a_out.close();
684 21 : }
685 :
686 : // Next write out Fq_Aq representor norm data
687 26 : file_name.str("");
688 11 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
689 27 : Xdr RB_Fq_Aq_innerprods_out(file_name.str(), mode);
690 :
691 70 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
692 : {
693 220 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
694 : {
695 2745 : for (unsigned int i=0; i<n_bfs; i++)
696 : {
697 4560 : RB_Fq_Aq_innerprods_out << Fq_Aq_representor_innerprods[q_f][q_a][i];
698 : }
699 : }
700 : }
701 15 : RB_Fq_Aq_innerprods_out.close();
702 :
703 : // Next write out Aq_Aq representor norm data
704 26 : file_name.str("");
705 11 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
706 23 : Xdr RB_Aq_Aq_innerprods_out(file_name.str(), mode);
707 :
708 15 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
709 105 : for (unsigned int i=0; i<Q_a_hat; i++)
710 : {
711 1650 : for (unsigned int j=0; j<n_bfs; j++)
712 : {
713 29160 : for (unsigned int l=0; l<n_bfs; l++)
714 : {
715 50100 : RB_Aq_Aq_innerprods_out << Aq_Aq_representor_innerprods[i][j][l];
716 : }
717 : }
718 : }
719 15 : RB_Aq_Aq_innerprods_out.close();
720 :
721 : // Also, write out the greedily selected parameters
722 : {
723 26 : file_name.str("");
724 11 : file_name << directory_name << "/greedy_params" << suffix;
725 34 : Xdr greedy_params_out(file_name.str(), mode);
726 :
727 290 : for (const auto & param : greedy_param_list)
728 1465 : for (const auto & pr : param)
729 2380 : for (const auto & value_vector : pr.second)
730 : {
731 : // Need to make a copy of the value so that it's not const
732 : // Xdr is not templated on const's
733 1498 : libmesh_error_msg_if(value_vector.size() != 1,
734 : "Error: multi-value RB parameters are not yet supported here.");
735 1190 : Real param_value = value_vector[0];
736 616 : greedy_params_out << param_value;
737 : }
738 15 : greedy_params_out.close();
739 7 : }
740 :
741 14 : }
742 27 : }
743 :
744 27 : void RBEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
745 : bool read_error_bound_data,
746 : const bool read_binary_data)
747 : {
748 16 : LOG_SCOPE("legacy_read_offline_data_from_files()", "RBEvaluation");
749 :
750 : // The reading mode: DECODE for binary, READ for ASCII
751 27 : XdrMODE mode = read_binary_data ? DECODE : READ;
752 :
753 : // The suffix to use for all the files that are written out
754 35 : const std::string suffix = read_binary_data ? ".xdr" : ".dat";
755 :
756 : // The string stream we'll use to make the file names
757 43 : std::ostringstream file_name;
758 :
759 : // First, find out how many basis functions we had when Greedy terminated
760 : unsigned int n_bfs;
761 : {
762 19 : file_name << directory_name << "/n_bfs" << suffix;
763 46 : assert_file_exists(file_name.str());
764 :
765 54 : Xdr n_bfs_in(file_name.str(), mode);
766 16 : n_bfs_in >> n_bfs;
767 27 : n_bfs_in.close();
768 11 : }
769 27 : resize_data_structures(n_bfs, read_error_bound_data);
770 :
771 : // Read in the parameter ranges
772 46 : file_name.str("");
773 19 : file_name << directory_name << "/parameter_ranges" << suffix;
774 16 : std::string continuous_param_file_name = file_name.str();
775 :
776 : // Read in the discrete parameter values
777 46 : file_name.str("");
778 19 : file_name << directory_name << "/discrete_parameter_values" << suffix;
779 16 : std::string discrete_param_file_name = file_name.str();
780 27 : read_parameter_data_from_files(continuous_param_file_name,
781 : discrete_param_file_name,
782 : read_binary_data);
783 :
784 : // Read in output data in multiple files
785 79 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
786 : {
787 104 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
788 : {
789 88 : file_name.str("");
790 52 : file_name << directory_name << "/output_";
791 : file_name << std::setw(3)
792 : << std::setprecision(0)
793 16 : << std::setfill('0')
794 16 : << std::right
795 16 : << n;
796 52 : file_name << "_";
797 : file_name << std::setw(3)
798 : << std::setprecision(0)
799 16 : << std::setfill('0')
800 16 : << std::right
801 16 : << q_l;
802 16 : file_name << suffix;
803 88 : assert_file_exists(file_name.str());
804 :
805 100 : Xdr output_n_in(file_name.str(), mode);
806 :
807 1092 : for (unsigned int j=0; j<n_bfs; j++)
808 : {
809 80 : Number value;
810 640 : output_n_in >> value;
811 1680 : RB_output_vectors[n][q_l](j) = value;
812 : }
813 52 : output_n_in.close();
814 20 : }
815 : }
816 :
817 27 : if (compute_RB_inner_product)
818 : {
819 : // Next read in the inner product matrix
820 10 : file_name.str("");
821 4 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
822 10 : assert_file_exists(file_name.str());
823 :
824 12 : Xdr RB_inner_product_matrix_in(file_name.str(), mode);
825 :
826 126 : for (unsigned int i=0; i<n_bfs; i++)
827 : {
828 2520 : for (unsigned int j=0; j<n_bfs; j++)
829 : {
830 0 : Number value;
831 1600 : RB_inner_product_matrix_in >> value;
832 3200 : RB_inner_product_matrix(i,j) = value;
833 : }
834 : }
835 6 : RB_inner_product_matrix_in.close();
836 2 : }
837 :
838 : // Next read in the Fq vectors
839 124 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
840 : {
841 166 : file_name.str("");
842 97 : file_name << directory_name << "/RB_F_";
843 : file_name << std::setw(3)
844 : << std::setprecision(0)
845 28 : << std::setfill('0')
846 28 : << std::right
847 28 : << q_f;
848 28 : file_name << suffix;
849 166 : assert_file_exists(file_name.str());
850 :
851 181 : Xdr RB_Fq_f_in(file_name.str(), mode);
852 :
853 1617 : for (unsigned int i=0; i<n_bfs; i++)
854 : {
855 200 : Number value;
856 880 : RB_Fq_f_in >> value;
857 1960 : RB_Fq_vector[q_f](i) = value;
858 : }
859 97 : RB_Fq_f_in.close();
860 41 : }
861 :
862 : // Next read in the Aq matrices
863 108 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
864 : {
865 138 : file_name.str("");
866 81 : file_name << directory_name << "/RB_A_";
867 : file_name << std::setw(3)
868 : << std::setprecision(0)
869 24 : << std::setfill('0')
870 24 : << std::right
871 24 : << q_a;
872 24 : file_name << suffix;
873 138 : assert_file_exists(file_name.str());
874 :
875 153 : Xdr RB_Aq_a_in(file_name.str(), mode);
876 :
877 1491 : for (unsigned int i=0; i<n_bfs; i++)
878 : {
879 26460 : for (unsigned int j=0; j<n_bfs; j++)
880 : {
881 2550 : Number value;
882 15000 : RB_Aq_a_in >> value;
883 32550 : RB_Aq_vector[q_a](i,j) = value;
884 : }
885 : }
886 81 : RB_Aq_a_in.close();
887 33 : }
888 :
889 :
890 27 : if (read_error_bound_data)
891 : {
892 : // Next read in Fq representor norm data
893 46 : file_name.str("");
894 19 : file_name << directory_name << "/Fq_innerprods" << suffix;
895 46 : assert_file_exists(file_name.str());
896 :
897 43 : Xdr RB_Fq_innerprods_in(file_name.str(), mode);
898 :
899 27 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
900 334 : for (unsigned int i=0; i<Q_f_hat; i++)
901 : {
902 395 : RB_Fq_innerprods_in >> Fq_representor_innerprods[i];
903 : }
904 27 : RB_Fq_innerprods_in.close();
905 :
906 : // Read in output data
907 79 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
908 : {
909 88 : file_name.str("");
910 52 : file_name << directory_name << "/output_";
911 : file_name << std::setw(3)
912 : << std::setprecision(0)
913 16 : << std::setfill('0')
914 16 : << std::right
915 16 : << n;
916 36 : file_name << "_dual_innerprods" << suffix;
917 88 : assert_file_exists(file_name.str());
918 :
919 84 : Xdr output_dual_innerprods_in(file_name.str(), mode);
920 :
921 52 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
922 104 : for (unsigned int q=0; q<Q_l_hat; q++)
923 : {
924 84 : output_dual_innerprods_in >> output_dual_innerprods[n][q];
925 : }
926 52 : output_dual_innerprods_in.close();
927 20 : }
928 :
929 :
930 : // Next read in Fq_Aq representor norm data
931 46 : file_name.str("");
932 19 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
933 46 : assert_file_exists(file_name.str());
934 :
935 51 : Xdr RB_Fq_Aq_innerprods_in(file_name.str(), mode);
936 :
937 124 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
938 : {
939 388 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
940 : {
941 4851 : for (unsigned int i=0; i<n_bfs; i++)
942 : {
943 8520 : RB_Fq_Aq_innerprods_in >> Fq_Aq_representor_innerprods[q_f][q_a][i];
944 : }
945 : }
946 : }
947 27 : RB_Fq_Aq_innerprods_in.close();
948 :
949 : // Next read in Aq_Aq representor norm data
950 46 : file_name.str("");
951 19 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
952 46 : assert_file_exists(file_name.str());
953 :
954 43 : Xdr RB_Aq_Aq_innerprods_in(file_name.str(), mode);
955 :
956 27 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
957 189 : for (unsigned int i=0; i<Q_a_hat; i++)
958 : {
959 2982 : for (unsigned int j=0; j<n_bfs; j++)
960 : {
961 52920 : for (unsigned int l=0; l<n_bfs; l++)
962 : {
963 95100 : RB_Aq_Aq_innerprods_in >> Aq_Aq_representor_innerprods[i][j][l];
964 : }
965 : }
966 : }
967 27 : RB_Aq_Aq_innerprods_in.close();
968 11 : }
969 :
970 : // Resize basis_functions even if we don't read them in so that
971 : // get_n_bfs() returns the correct value. Initialize the pointers
972 : // to nullptr.
973 19 : basis_functions.clear();
974 27 : set_n_basis_functions(n_bfs);
975 38 : }
976 :
977 471 : void RBEvaluation::assert_file_exists(const std::string & file_name)
978 : {
979 471 : libmesh_error_msg_if(!std::ifstream(file_name.c_str()), "File missing: " << file_name);
980 471 : }
981 :
982 27 : void RBEvaluation::write_out_basis_functions(System & sys,
983 : const std::string & directory_name,
984 : const bool write_binary_basis_functions)
985 : {
986 16 : LOG_SCOPE("write_out_basis_functions()", "RBEvaluation");
987 :
988 8 : std::vector<NumericVector<Number>*> basis_functions_ptrs;
989 645 : for (std::size_t i=0; i<basis_functions.size(); i++)
990 : {
991 470 : basis_functions_ptrs.push_back(basis_functions[i].get());
992 : }
993 :
994 46 : write_out_vectors(sys,
995 : basis_functions_ptrs,
996 : directory_name,
997 : "bf",
998 : write_binary_basis_functions);
999 27 : }
1000 :
1001 27 : void RBEvaluation::write_out_vectors(System & sys,
1002 : std::vector<NumericVector<Number>*> & vectors,
1003 : const std::string & directory_name,
1004 : const std::string & data_name,
1005 : const bool write_binary_vectors)
1006 : {
1007 16 : LOG_SCOPE("write_out_vectors()", "RBEvaluation");
1008 :
1009 27 : if (sys.comm().rank() == 0)
1010 : {
1011 : // Make a directory to store all the data files
1012 15 : Utility::mkdir(directory_name.c_str());
1013 : }
1014 :
1015 : // Make sure processors are synced up before we begin
1016 27 : sys.comm().barrier();
1017 :
1018 43 : std::ostringstream file_name;
1019 35 : const std::string basis_function_suffix = (write_binary_vectors ? ".xdr" : ".dat");
1020 :
1021 30 : file_name << directory_name << "/" << data_name << "_data" << basis_function_suffix;
1022 19 : Xdr bf_data(file_name.str(),
1023 70 : write_binary_vectors ? ENCODE : WRITE);
1024 :
1025 : // Following EquationSystems::write(), we should only write the header information
1026 : // if we're proc 0
1027 27 : if (sys.comm().rank() == 0)
1028 : {
1029 30 : std::string version("libMesh-" + libMesh::get_io_compatibility_version());
1030 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1031 4 : version += " with infinite elements";
1032 : #endif
1033 15 : bf_data.data(version ,"# File Format Identifier");
1034 :
1035 15 : sys.write_header(bf_data, /*(unused arg)*/ version, /*write_additional_data=*/false);
1036 : }
1037 :
1038 : // Following EquationSystemsIO::write, we use a temporary numbering (node major)
1039 : // before writing out the data
1040 27 : MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());
1041 :
1042 : // Write all vectors at once.
1043 : {
1044 : // Note the API wants pointers to constant vectors, hence this...
1045 16 : std::vector<const NumericVector<Number> *> bf_out;
1046 497 : for (const auto & vec : vectors)
1047 470 : bf_out.push_back(vec);
1048 :
1049 : // for (auto & val : vectors)
1050 : // bf_out.push_back(val);
1051 27 : sys.write_serialized_vectors (bf_data, bf_out);
1052 : }
1053 :
1054 :
1055 : // set the current version
1056 8 : bf_data.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
1057 : LIBMESH_MINOR_VERSION,
1058 : LIBMESH_MICRO_VERSION));
1059 :
1060 :
1061 : // Undo the temporary renumbering
1062 27 : sys.get_mesh().fix_broken_node_and_element_numbering();
1063 38 : }
1064 :
1065 27 : void RBEvaluation::read_in_basis_functions(System & sys,
1066 : const std::string & directory_name,
1067 : const bool read_binary_basis_functions)
1068 : {
1069 8 : LOG_SCOPE("read_in_basis_functions()", "RBEvaluation");
1070 :
1071 43 : read_in_vectors(sys,
1072 27 : basis_functions,
1073 : directory_name,
1074 : "bf",
1075 : read_binary_basis_functions);
1076 27 : }
1077 :
1078 27 : void RBEvaluation::read_in_vectors(System & sys,
1079 : std::vector<std::unique_ptr<NumericVector<Number>>> & vectors,
1080 : const std::string & directory_name,
1081 : const std::string & data_name,
1082 : const bool read_binary_vectors)
1083 : {
1084 16 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>> *> vectors_vec;
1085 27 : vectors_vec.push_back(&vectors);
1086 :
1087 24 : std::vector<std::string> directory_name_vec;
1088 27 : directory_name_vec.push_back(directory_name);
1089 :
1090 16 : std::vector<std::string> data_name_vec;
1091 27 : data_name_vec.push_back(data_name);
1092 :
1093 35 : read_in_vectors_from_multiple_files(sys,
1094 : vectors_vec,
1095 : directory_name_vec,
1096 : data_name_vec,
1097 : read_binary_vectors);
1098 38 : }
1099 :
1100 27 : void RBEvaluation::read_in_vectors_from_multiple_files(System & sys,
1101 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>> *> multiple_vectors,
1102 : const std::vector<std::string> & multiple_directory_names,
1103 : const std::vector<std::string> & multiple_data_names,
1104 : const bool read_binary_vectors)
1105 : {
1106 8 : LOG_SCOPE("read_in_vectors_from_multiple_files()", "RBEvaluation");
1107 :
1108 16 : std::size_t n_files = multiple_vectors.size();
1109 16 : std::size_t n_directories = multiple_directory_names.size();
1110 8 : libmesh_assert((n_files == n_directories) && (n_files == multiple_data_names.size()));
1111 :
1112 27 : if (n_files == 0)
1113 0 : return;
1114 :
1115 : // Make sure processors are synced up before we begin
1116 27 : sys.comm().barrier();
1117 :
1118 43 : std::ostringstream file_name;
1119 35 : const std::string basis_function_suffix = (read_binary_vectors ? ".xdr" : ".dat");
1120 :
1121 : // Following EquationSystemsIO::read, we use a temporary numbering (node major)
1122 : // before writing out the data. For the sake of efficiency, we do this once for
1123 : // all the vectors that we read in.
1124 27 : MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());
1125 :
1126 54 : for (std::size_t data_index=0; data_index<n_directories; data_index++)
1127 : {
1128 35 : std::vector<std::unique_ptr<NumericVector<Number>>> & vectors = *multiple_vectors[data_index];
1129 :
1130 : // Allocate storage for each vector
1131 497 : for (auto & vec : vectors)
1132 : {
1133 : // vectors should all be nullptr, otherwise we get a memory leak when
1134 : // we create the new vectors in RBEvaluation::read_in_vectors.
1135 470 : libmesh_error_msg_if(vec, "Non-nullptr vector passed to read_in_vectors_from_multiple_files");
1136 :
1137 800 : vec = NumericVector<Number>::build(sys.comm());
1138 :
1139 610 : vec->init (sys.n_dofs(),
1140 : sys.n_local_dofs(),
1141 : false,
1142 280 : PARALLEL);
1143 : }
1144 :
1145 38 : file_name.str("");
1146 16 : file_name << multiple_directory_names[data_index]
1147 16 : << "/" << multiple_data_names[data_index]
1148 30 : << "_data" << basis_function_suffix;
1149 :
1150 : // On processor zero check to be sure the file exists
1151 27 : if (sys.comm().rank() == 0)
1152 : {
1153 : struct stat stat_info;
1154 15 : int stat_result = stat(file_name.str().c_str(), &stat_info);
1155 :
1156 15 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
1157 : }
1158 :
1159 27 : assert_file_exists(file_name.str());
1160 27 : Xdr vector_data(file_name.str(),
1161 70 : read_binary_vectors ? DECODE : READ);
1162 :
1163 : // Read the header data. This block of code is based on EquationSystems::_read_impl.
1164 : {
1165 16 : std::string version;
1166 27 : vector_data.data(version);
1167 :
1168 35 : const std::string libMesh_label = "libMesh-";
1169 8 : std::string::size_type lm_pos = version.find(libMesh_label);
1170 27 : libmesh_error_msg_if(lm_pos == std::string::npos, "version info missing in Xdr header");
1171 :
1172 43 : std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
1173 27 : int ver_major = 0, ver_minor = 0, ver_patch = 0;
1174 : char dot;
1175 27 : iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
1176 27 : vector_data.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));
1177 :
1178 : // Actually read the header data. When we do this, set read_header=false
1179 : // so that we do not reinit sys, since we assume that it has already been
1180 : // set up properly (e.g. the appropriate variables have already been added).
1181 27 : sys.read_header(vector_data, version, /*read_header=*/false, /*read_additional_data=*/false);
1182 11 : }
1183 :
1184 : // Comply with the System::read_serialized_vectors() interface which uses dumb pointers.
1185 16 : std::vector<NumericVector<Number> *> vec_in;
1186 497 : for (auto & vec : vectors)
1187 470 : vec_in.push_back(vec.get());
1188 :
1189 8 : sys.read_serialized_vectors (vector_data, vec_in);
1190 11 : }
1191 :
1192 : // Undo the temporary renumbering
1193 27 : sys.get_mesh().fix_broken_node_and_element_numbering();
1194 11 : }
1195 :
1196 256042 : void RBEvaluation::check_evaluated_thetas_size(const std::vector<Number> * evaluated_thetas) const
1197 : {
1198 256042 : if (rb_theta_expansion && evaluated_thetas)
1199 : {
1200 0 : auto actual_size = evaluated_thetas->size();
1201 : auto expected_size =
1202 0 : rb_theta_expansion->get_n_A_terms() +
1203 0 : rb_theta_expansion->get_n_F_terms() +
1204 0 : rb_theta_expansion->get_total_n_output_terms();
1205 :
1206 0 : libmesh_error_msg_if(actual_size != expected_size,
1207 : "ERROR: Evaluated thetas vector has size " <<
1208 : actual_size << ", but we expected " <<
1209 : rb_theta_expansion->get_n_A_terms() << " A term(s), " <<
1210 : rb_theta_expansion->get_n_F_terms() << " F term(s), and " <<
1211 : rb_theta_expansion->get_total_n_output_terms() << " output term(s), " <<
1212 : "for a total of " << expected_size << " terms.");
1213 : }
1214 256042 : }
1215 :
1216 : } // namespace libMesh
|
