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 568 : RBEvaluation::RBEvaluation (const Parallel::Communicator & comm_in)
49 : :
50 : ParallelObject(comm_in),
51 536 : evaluate_RB_error_bound(true),
52 536 : compute_RB_inner_product(false),
53 568 : rb_theta_expansion(nullptr)
54 : {
55 568 : }
56 :
57 1608 : 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 284 : void RBEvaluation::set_n_basis_functions(unsigned int n_bfs)
76 : {
77 284 : basis_functions.resize(n_bfs);
78 284 : }
79 :
80 568 : void RBEvaluation::set_rb_theta_expansion(RBThetaExpansion & rb_theta_expansion_in)
81 : {
82 568 : rb_theta_expansion = &rb_theta_expansion_in;
83 568 : }
84 :
85 6156325 : RBThetaExpansion & RBEvaluation::get_rb_theta_expansion()
86 : {
87 6156325 : libmesh_error_msg_if(!is_rb_theta_expansion_initialized(),
88 : "Error: rb_theta_expansion hasn't been initialized yet");
89 :
90 6156325 : return *rb_theta_expansion;
91 : }
92 :
93 2342 : const RBThetaExpansion & RBEvaluation::get_rb_theta_expansion() const
94 : {
95 2342 : libmesh_error_msg_if(!is_rb_theta_expansion_initialized(),
96 : "Error: rb_theta_expansion hasn't been initialized yet");
97 :
98 2342 : return *rb_theta_expansion;
99 : }
100 :
101 6158667 : bool RBEvaluation::is_rb_theta_expansion_initialized() const
102 : {
103 6158667 : if (rb_theta_expansion)
104 : {
105 178454 : return true;
106 : }
107 : else
108 : {
109 0 : return false;
110 : }
111 : }
112 :
113 568 : 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 568 : 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 568 : if (compute_RB_inner_product)
123 136 : RB_inner_product_matrix.resize(Nmax,Nmax);
124 :
125 : // Allocate dense matrices for RB solves
126 568 : RB_Aq_vector.resize(rb_theta_expansion->get_n_A_terms());
127 :
128 2272 : for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
129 : {
130 : // Initialize the memory for the RB matrices
131 1704 : RB_Aq_vector[q].resize(Nmax,Nmax);
132 : }
133 :
134 568 : RB_Fq_vector.resize(rb_theta_expansion->get_n_F_terms());
135 :
136 2556 : for (unsigned int q=0; q<rb_theta_expansion->get_n_F_terms(); q++)
137 : {
138 : // Initialize the memory for the RB vectors
139 1988 : RB_Fq_vector[q].resize(Nmax);
140 : }
141 :
142 :
143 : // Initialize the RB output vectors
144 568 : RB_output_vectors.resize(rb_theta_expansion->get_n_outputs());
145 1704 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
146 : {
147 1168 : RB_output_vectors[n].resize(rb_theta_expansion->get_n_output_terms(n));
148 2272 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
149 : {
150 1168 : RB_output_vectors[n][q_l].resize(Nmax);
151 : }
152 : }
153 :
154 : // Initialize vectors storing output data
155 568 : RB_outputs.resize(rb_theta_expansion->get_n_outputs(), 0.);
156 :
157 :
158 568 : if (resize_error_bound_data)
159 : {
160 : // Initialize vectors for the norms of the Fq representors
161 568 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
162 568 : Fq_representor_innerprods.resize(Q_f_hat);
163 :
164 : // Initialize vectors for the norms of the representors
165 568 : Fq_Aq_representor_innerprods.resize(rb_theta_expansion->get_n_F_terms());
166 2556 : for (unsigned int i=0; i<rb_theta_expansion->get_n_F_terms(); i++)
167 : {
168 2044 : Fq_Aq_representor_innerprods[i].resize(rb_theta_expansion->get_n_A_terms());
169 7952 : for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
170 : {
171 6300 : Fq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
172 : }
173 : }
174 :
175 568 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
176 568 : Aq_Aq_representor_innerprods.resize(Q_a_hat);
177 3976 : for (unsigned int i=0; i<Q_a_hat; i++)
178 : {
179 3504 : Aq_Aq_representor_innerprods[i].resize(Nmax);
180 62904 : for (unsigned int j=0; j<Nmax; j++)
181 : {
182 62856 : Aq_Aq_representor_innerprods[i][j].resize(Nmax, 0.);
183 : }
184 : }
185 :
186 568 : RB_output_error_bounds.resize(rb_theta_expansion->get_n_outputs(), 0.);
187 :
188 : // Resize the output dual norm vectors
189 568 : output_dual_innerprods.resize(rb_theta_expansion->get_n_outputs());
190 1704 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
191 : {
192 1136 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
193 1168 : output_dual_innerprods[n].resize(Q_l_hat);
194 : }
195 :
196 : // Clear and resize the vector of Aq_representors
197 568 : clear_riesz_representors();
198 :
199 568 : Aq_representor.resize(rb_theta_expansion->get_n_A_terms());
200 2272 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
201 : {
202 1752 : Aq_representor[q_a].resize(Nmax);
203 : }
204 : }
205 568 : }
206 :
207 3426803 : NumericVector<Number> & RBEvaluation::get_basis_function(unsigned int i)
208 : {
209 97718 : libmesh_assert_less (i, basis_functions.size());
210 :
211 3524521 : return *(basis_functions[i]);
212 : }
213 :
214 56800 : const NumericVector<Number> & RBEvaluation::get_basis_function(unsigned int i) const
215 : {
216 1600 : libmesh_assert_less (i, basis_functions.size());
217 :
218 58400 : return *(basis_functions[i]);
219 : }
220 :
221 384914 : Real RBEvaluation::rb_solve(unsigned int N)
222 : {
223 384914 : return rb_solve(N, nullptr);
224 : }
225 :
226 384914 : Real RBEvaluation::rb_solve(unsigned int N,
227 : const std::vector<Number> * evaluated_thetas)
228 : {
229 64012 : LOG_SCOPE("rb_solve()", "RBEvaluation");
230 :
231 384914 : 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 384914 : this->check_evaluated_thetas_size(evaluated_thetas);
238 :
239 384914 : const RBParameters & mu = get_parameters();
240 :
241 : // Resize (and clear) the solution vector
242 352908 : RB_solution.resize(N);
243 :
244 : // Assemble the RB system
245 448926 : DenseMatrix<Number> RB_system_matrix(N,N);
246 32006 : RB_system_matrix.zero();
247 :
248 448926 : DenseMatrix<Number> RB_Aq_a;
249 1539656 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
250 : {
251 1250760 : RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);
252 :
253 1154742 : if (evaluated_thetas)
254 0 : RB_system_matrix.add((*evaluated_thetas)[q_a], RB_Aq_a);
255 : else
256 1154742 : RB_system_matrix.add(rb_theta_expansion->eval_A_theta(q_a, mu), RB_Aq_a);
257 : }
258 :
259 : // Assemble the RB rhs
260 384914 : DenseVector<Number> RB_rhs(N);
261 32006 : RB_rhs.zero();
262 :
263 384914 : DenseVector<Number> RB_Fq_f;
264 2570538 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
265 : {
266 2367650 : RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);
267 :
268 2185624 : 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 2185624 : RB_rhs.add(rb_theta_expansion->eval_F_theta(q_f, mu), RB_Fq_f);
272 : }
273 :
274 : // Solve the linear system
275 384914 : if (N > 0)
276 : {
277 359614 : RB_system_matrix.lu_solve(RB_rhs, RB_solution);
278 : }
279 :
280 : // Place to store the output of get_principal_subvector() calls
281 384914 : DenseVector<Number> RB_output_vector_N;
282 :
283 : // Evaluate RB outputs
284 32006 : unsigned int output_counter = 0;
285 484002 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
286 : {
287 107096 : RB_outputs[n] = 0.;
288 198176 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
289 : {
290 115104 : 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 99088 : 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 99088 : 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 99088 : rb_theta_expansion->eval_output_theta(n, q_l, mu);
305 :
306 : // Finally, accumulate the result in RB_outputs[n]
307 91084 : RB_outputs[n] += coeff * dot_prod;
308 :
309 : // Go to next output
310 99088 : output_counter++;
311 : }
312 : }
313 :
314 384914 : if (evaluate_RB_error_bound) // Calculate the error bounds
315 : {
316 : // Evaluate the dual norm of the residual for RB_solution_vector
317 384914 : Real epsilon_N = compute_residual_dual_norm(N, evaluated_thetas);
318 :
319 : // Get lower bound for coercivity constant
320 384914 : 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 384914 : Real abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB);
326 :
327 : // Now compute the output error bounds
328 484002 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
329 99088 : 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 320902 : }
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 384914 : 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 384914 : this->check_evaluated_thetas_size(evaluated_thetas);
363 :
364 : // If evaluated_thetas is provided, then mu is not actually used for anything
365 384914 : const RBParameters & mu = get_parameters();
366 :
367 384914 : const unsigned int n_A_terms = rb_theta_expansion->get_n_A_terms();
368 384914 : 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 11959002 : auto eval_F = [&](unsigned int index)
378 : {
379 11959002 : return (evaluated_thetas) ? (*evaluated_thetas)[index + n_A_terms] : rb_theta_expansion->eval_F_theta(index, mu);
380 384914 : };
381 10021098 : auto eval_A = [&](unsigned int index)
382 : {
383 10021098 : return (evaluated_thetas) ? (*evaluated_thetas)[index] : rb_theta_expansion->eval_A_theta(index, mu);
384 384914 : };
385 :
386 32006 : unsigned int q=0;
387 2570538 : for (unsigned int q_f1=0; q_f1<n_F_terms; q_f1++)
388 : {
389 2185624 : const Number val_q_f1 = eval_F(q_f1);
390 :
391 9773378 : for (unsigned int q_f2=q_f1; q_f2<n_F_terms; q_f2++)
392 : {
393 7587754 : const Number val_q_f2 = eval_F(q_f2);
394 :
395 7587754 : Real delta = (q_f1==q_f2) ? 1. : 2.;
396 7587754 : residual_norm_sq += delta * libmesh_real(val_q_f1 * libmesh_conj(val_q_f2) * Fq_representor_innerprods[q] );
397 :
398 7587754 : q++;
399 : }
400 : }
401 :
402 2570538 : for (unsigned int q_f=0; q_f<n_F_terms; q_f++)
403 : {
404 2185624 : const Number val_q_f = eval_F(q_f);
405 :
406 8742496 : for (unsigned int q_a=0; q_a<n_A_terms; q_a++)
407 : {
408 6556872 : const Number val_q_a = eval_A(q_a);
409 :
410 52649790 : for (unsigned int i=0; i<N; i++)
411 : {
412 3838200 : Real delta = 2.;
413 42255318 : residual_norm_sq +=
414 49931718 : delta * libmesh_real( val_q_f * libmesh_conj(val_q_a) *
415 53769318 : libmesh_conj(RB_solution(i)) * Fq_Aq_representor_innerprods[q_f][q_a][i] );
416 : }
417 : }
418 : }
419 :
420 32006 : q=0;
421 1539656 : for (unsigned int q_a1=0; q_a1<n_A_terms; q_a1++)
422 : {
423 1154742 : const Number val_q_a1 = eval_A(q_a1);
424 :
425 3464226 : for (unsigned int q_a2=q_a1; q_a2<n_A_terms; q_a2++)
426 : {
427 2309484 : const Number val_q_a2 = eval_A(q_a2);
428 :
429 2309484 : Real delta = (q_a1==q_a2) ? 1. : 2.;
430 :
431 18831420 : for (unsigned int i=0; i<N; i++)
432 : {
433 180882852 : for (unsigned int j=0; j<N; j++)
434 : {
435 150693816 : residual_norm_sq +=
436 178038216 : delta * libmesh_real( libmesh_conj(val_q_a1) * val_q_a2 *
437 205377516 : libmesh_conj(RB_solution(i)) * RB_solution(j) * Aq_Aq_representor_innerprods[q][i][j] );
438 : }
439 : }
440 :
441 2309484 : q++;
442 : }
443 : }
444 :
445 384914 : 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 190 : residual_norm_sq = std::abs(residual_norm_sq);
454 : }
455 :
456 416920 : 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 2591914 : Real RBEvaluation::residual_scaling_denom(Real alpha_LB)
468 : {
469 : // Here we implement the residual scaling for a coercive
470 : // problem.
471 2591914 : return alpha_LB;
472 : }
473 :
474 9015368 : 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 9015368 : const RBParameters & mu = this->get_parameters();
482 :
483 : // Index into output_dual_innerprods
484 816816 : unsigned int q=0;
485 18030736 : for (unsigned int q_l1=0; q_l1<rb_theta_expansion->get_n_output_terms(n); q_l1++)
486 : {
487 18030736 : for (unsigned int q_l2=q_l1; q_l2<rb_theta_expansion->get_n_output_terms(n); q_l2++)
488 : {
489 9015368 : Real delta = (q_l1==q_l2) ? 1. : 2.;
490 :
491 : Number val_l1 =
492 9015368 : 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 9015368 : 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 9015368 : rb_theta_expansion->eval_output_theta(n, q_l2, mu);
500 :
501 9007364 : output_bound_sq += delta * libmesh_real(
502 9832184 : libmesh_conj(val_l1) * val_l2 * output_dual_innerprods[n][q]);
503 :
504 9015368 : q++;
505 : }
506 : }
507 :
508 9015368 : return libmesh_real(std::sqrt( output_bound_sq ));
509 : }
510 :
511 568 : void RBEvaluation::clear_riesz_representors()
512 : {
513 552 : Aq_representor.clear();
514 568 : }
515 :
516 284 : 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 284 : unsigned int n_bfs = get_n_basis_functions();
523 :
524 : // The writing mode: ENCODE for binary, WRITE for ASCII
525 284 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
526 :
527 : // The suffix to use for all the files that are written out
528 292 : const std::string suffix = write_binary_data ? ".xdr" : ".dat";
529 :
530 292 : if (this->processor_id() == 0)
531 : {
532 :
533 : // Make a directory to store all the data files
534 48 : 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 56 : std::ostringstream file_name;
543 : {
544 44 : file_name << directory_name << "/n_bfs" << suffix;
545 96 : Xdr n_bfs_out(file_name.str(), mode);
546 8 : n_bfs_out << n_bfs;
547 48 : n_bfs_out.close();
548 40 : }
549 :
550 : // Write out the parameter ranges
551 92 : file_name.str("");
552 44 : 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 92 : file_name.str("");
557 44 : file_name << directory_name << "/discrete_parameter_values" << suffix;
558 8 : std::string discrete_param_file_name = file_name.str();
559 :
560 48 : 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 92 : file_name.str("");
566 44 : file_name << directory_name << "/Fq_innerprods" << suffix;
567 56 : Xdr RB_Fq_innerprods_out(file_name.str(), mode);
568 48 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
569 576 : for (unsigned int i=0; i<Q_f_hat; i++)
570 : {
571 572 : RB_Fq_innerprods_out << Fq_representor_innerprods[i];
572 : }
573 48 : RB_Fq_innerprods_out.close();
574 :
575 : // Write out output data
576 144 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
577 : {
578 184 : file_name.str("");
579 96 : 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 88 : file_name << "_dual_innerprods" << suffix;
587 112 : Xdr output_dual_innerprods_out(file_name.str(), mode);
588 :
589 96 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
590 192 : for (unsigned int q=0; q<Q_l_hat; q++)
591 : {
592 112 : output_dual_innerprods_out << output_dual_innerprods[n][q];
593 : }
594 96 : output_dual_innerprods_out.close();
595 80 : }
596 :
597 :
598 : // Write out output data to multiple files
599 144 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
600 : {
601 192 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
602 : {
603 184 : file_name.str("");
604 96 : 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 96 : 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 120 : Xdr output_n_out(file_name.str(), mode);
618 :
619 1960 : for (unsigned int j=0; j<n_bfs; j++)
620 : {
621 640 : output_n_out << RB_output_vectors[n][q_l](j);
622 : }
623 96 : output_n_out.close();
624 80 : }
625 : }
626 :
627 48 : if (compute_RB_inner_product)
628 : {
629 : // Next write out the inner product matrix
630 21 : file_name.str("");
631 10 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
632 14 : Xdr RB_inner_product_matrix_out(file_name.str(), mode);
633 231 : for (unsigned int i=0; i<n_bfs; i++)
634 : {
635 4620 : 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 11 : RB_inner_product_matrix_out.close();
641 9 : }
642 :
643 : // Next write out the Fq vectors
644 216 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
645 : {
646 322 : file_name.str("");
647 168 : 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 210 : Xdr RB_Fq_f_out(file_name.str(), mode);
655 :
656 2794 : for (unsigned int i=0; i<n_bfs; i++)
657 : {
658 660 : RB_Fq_f_out << RB_Fq_vector[q_f](i);
659 : }
660 168 : RB_Fq_f_out.close();
661 140 : }
662 :
663 : // Next write out the Aq matrices
664 192 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
665 : {
666 276 : file_name.str("");
667 144 : 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 180 : Xdr RB_Aq_a_out(file_name.str(), mode);
675 :
676 2622 : for (unsigned int i=0; i<n_bfs; i++)
677 : {
678 46386 : 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 144 : RB_Aq_a_out.close();
684 120 : }
685 :
686 : // Next write out Fq_Aq representor norm data
687 92 : file_name.str("");
688 44 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
689 60 : Xdr RB_Fq_Aq_innerprods_out(file_name.str(), mode);
690 :
691 216 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
692 : {
693 672 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
694 : {
695 8382 : for (unsigned int i=0; i<n_bfs; i++)
696 : {
697 9858 : RB_Fq_Aq_innerprods_out << Fq_Aq_representor_innerprods[q_f][q_a][i];
698 : }
699 : }
700 : }
701 48 : RB_Fq_Aq_innerprods_out.close();
702 :
703 : // Next write out Aq_Aq representor norm data
704 92 : file_name.str("");
705 44 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
706 56 : Xdr RB_Aq_Aq_innerprods_out(file_name.str(), mode);
707 :
708 48 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
709 336 : for (unsigned int i=0; i<Q_a_hat; i++)
710 : {
711 5244 : for (unsigned int j=0; j<n_bfs; j++)
712 : {
713 92772 : for (unsigned int l=0; l<n_bfs; l++)
714 : {
715 110316 : RB_Aq_Aq_innerprods_out << Aq_Aq_representor_innerprods[i][j][l];
716 : }
717 : }
718 : }
719 48 : RB_Aq_Aq_innerprods_out.close();
720 :
721 : // Also, write out the greedily selected parameters
722 : {
723 92 : file_name.str("");
724 44 : file_name << directory_name << "/greedy_params" << suffix;
725 100 : Xdr greedy_params_out(file_name.str(), mode);
726 :
727 922 : for (const auto & param : greedy_param_list)
728 4549 : for (const auto & pr : param)
729 7350 : 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 3983 : libmesh_error_msg_if(value_vector.size() != 1,
734 : "Error: multi-value RB parameters are not yet supported here.");
735 3675 : Real param_value = value_vector[0];
736 616 : greedy_params_out << param_value;
737 : }
738 48 : greedy_params_out.close();
739 40 : }
740 :
741 80 : }
742 284 : }
743 :
744 284 : 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 284 : XdrMODE mode = read_binary_data ? DECODE : READ;
752 :
753 : // The suffix to use for all the files that are written out
754 292 : const std::string suffix = read_binary_data ? ".xdr" : ".dat";
755 :
756 : // The string stream we'll use to make the file names
757 300 : 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 276 : file_name << directory_name << "/n_bfs" << suffix;
763 560 : assert_file_exists(file_name.str());
764 :
765 568 : Xdr n_bfs_in(file_name.str(), mode);
766 16 : n_bfs_in >> n_bfs;
767 284 : n_bfs_in.close();
768 268 : }
769 284 : resize_data_structures(n_bfs, read_error_bound_data);
770 :
771 : // Read in the parameter ranges
772 560 : file_name.str("");
773 276 : 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 560 : file_name.str("");
778 276 : file_name << directory_name << "/discrete_parameter_values" << suffix;
779 16 : std::string discrete_param_file_name = file_name.str();
780 284 : 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 852 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
786 : {
787 1136 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
788 : {
789 1120 : file_name.str("");
790 568 : 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 568 : 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 1120 : assert_file_exists(file_name.str());
804 :
805 616 : Xdr output_n_in(file_name.str(), mode);
806 :
807 11872 : for (unsigned int j=0; j<n_bfs; j++)
808 : {
809 80 : Number value;
810 640 : output_n_in >> value;
811 11944 : RB_output_vectors[n][q_l](j) = value;
812 : }
813 568 : output_n_in.close();
814 536 : }
815 : }
816 :
817 284 : if (compute_RB_inner_product)
818 : {
819 : // Next read in the inner product matrix
820 138 : file_name.str("");
821 68 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
822 138 : assert_file_exists(file_name.str());
823 :
824 76 : Xdr RB_inner_product_matrix_in(file_name.str(), mode);
825 :
826 1470 : for (unsigned int i=0; i<n_bfs; i++)
827 : {
828 29400 : for (unsigned int j=0; j<n_bfs; j++)
829 : {
830 0 : Number value;
831 1600 : RB_inner_product_matrix_in >> value;
832 28800 : RB_inner_product_matrix(i,j) = value;
833 : }
834 : }
835 70 : RB_inner_product_matrix_in.close();
836 66 : }
837 :
838 : // Next read in the Fq vectors
839 1278 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
840 : {
841 1960 : file_name.str("");
842 994 : 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 1960 : assert_file_exists(file_name.str());
850 :
851 1078 : Xdr RB_Fq_f_in(file_name.str(), mode);
852 :
853 16600 : for (unsigned int i=0; i<n_bfs; i++)
854 : {
855 200 : Number value;
856 880 : RB_Fq_f_in >> value;
857 16046 : RB_Fq_vector[q_f](i) = value;
858 : }
859 994 : RB_Fq_f_in.close();
860 938 : }
861 :
862 : // Next read in the Aq matrices
863 1136 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
864 : {
865 1680 : file_name.str("");
866 852 : 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 1680 : assert_file_exists(file_name.str());
874 :
875 924 : Xdr RB_Aq_a_in(file_name.str(), mode);
876 :
877 15720 : for (unsigned int i=0; i<n_bfs; i++)
878 : {
879 280026 : for (unsigned int j=0; j<n_bfs; j++)
880 : {
881 2550 : Number value;
882 15000 : RB_Aq_a_in >> value;
883 272658 : RB_Aq_vector[q_a](i,j) = value;
884 : }
885 : }
886 852 : RB_Aq_a_in.close();
887 804 : }
888 :
889 :
890 284 : if (read_error_bound_data)
891 : {
892 : // Next read in Fq representor norm data
893 560 : file_name.str("");
894 276 : file_name << directory_name << "/Fq_innerprods" << suffix;
895 560 : assert_file_exists(file_name.str());
896 :
897 300 : Xdr RB_Fq_innerprods_in(file_name.str(), mode);
898 :
899 284 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
900 3408 : for (unsigned int i=0; i<Q_f_hat; i++)
901 : {
902 3212 : RB_Fq_innerprods_in >> Fq_representor_innerprods[i];
903 : }
904 284 : RB_Fq_innerprods_in.close();
905 :
906 : // Read in output data
907 852 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
908 : {
909 1120 : file_name.str("");
910 568 : 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 552 : file_name << "_dual_innerprods" << suffix;
917 1120 : assert_file_exists(file_name.str());
918 :
919 600 : Xdr output_dual_innerprods_in(file_name.str(), mode);
920 :
921 568 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
922 1136 : for (unsigned int q=0; q<Q_l_hat; q++)
923 : {
924 600 : output_dual_innerprods_in >> output_dual_innerprods[n][q];
925 : }
926 568 : output_dual_innerprods_in.close();
927 536 : }
928 :
929 :
930 : // Next read in Fq_Aq representor norm data
931 560 : file_name.str("");
932 276 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
933 560 : assert_file_exists(file_name.str());
934 :
935 308 : Xdr RB_Fq_Aq_innerprods_in(file_name.str(), mode);
936 :
937 1278 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
938 : {
939 3976 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
940 : {
941 49800 : for (unsigned int i=0; i<n_bfs; i++)
942 : {
943 50778 : RB_Fq_Aq_innerprods_in >> Fq_Aq_representor_innerprods[q_f][q_a][i];
944 : }
945 : }
946 : }
947 284 : RB_Fq_Aq_innerprods_in.close();
948 :
949 : // Next read in Aq_Aq representor norm data
950 560 : file_name.str("");
951 276 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
952 560 : assert_file_exists(file_name.str());
953 :
954 300 : Xdr RB_Aq_Aq_innerprods_in(file_name.str(), mode);
955 :
956 284 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
957 1988 : for (unsigned int i=0; i<Q_a_hat; i++)
958 : {
959 31440 : for (unsigned int j=0; j<n_bfs; j++)
960 : {
961 560052 : for (unsigned int l=0; l<n_bfs; l++)
962 : {
963 575316 : RB_Aq_Aq_innerprods_in >> Aq_Aq_representor_innerprods[i][j][l];
964 : }
965 : }
966 : }
967 284 : RB_Aq_Aq_innerprods_in.close();
968 268 : }
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 276 : basis_functions.clear();
974 284 : set_n_basis_functions(n_bfs);
975 552 : }
976 :
977 5032 : void RBEvaluation::assert_file_exists(const std::string & file_name)
978 : {
979 5032 : libmesh_error_msg_if(!std::ifstream(file_name.c_str()), "File missing: " << file_name);
980 5032 : }
981 :
982 284 : 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 5388 : for (std::size_t i=0; i<basis_functions.size(); i++)
990 : {
991 4956 : basis_functions_ptrs.push_back(basis_functions[i].get());
992 : }
993 :
994 560 : write_out_vectors(sys,
995 : basis_functions_ptrs,
996 : directory_name,
997 : "bf",
998 : write_binary_basis_functions);
999 284 : }
1000 :
1001 284 : 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 284 : if (sys.comm().rank() == 0)
1010 : {
1011 : // Make a directory to store all the data files
1012 48 : Utility::mkdir(directory_name.c_str());
1013 : }
1014 :
1015 : // Make sure processors are synced up before we begin
1016 284 : sys.comm().barrier();
1017 :
1018 300 : std::ostringstream file_name;
1019 292 : const std::string basis_function_suffix = (write_binary_vectors ? ".xdr" : ".dat");
1020 :
1021 544 : file_name << directory_name << "/" << data_name << "_data" << basis_function_suffix;
1022 276 : Xdr bf_data(file_name.str(),
1023 584 : write_binary_vectors ? ENCODE : WRITE);
1024 :
1025 : // Following EquationSystems::write(), we should only write the header information
1026 : // if we're proc 0
1027 284 : if (sys.comm().rank() == 0)
1028 : {
1029 96 : std::string version("libMesh-" + libMesh::get_io_compatibility_version());
1030 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1031 4 : version += " with infinite elements";
1032 : #endif
1033 48 : bf_data.data(version ,"# File Format Identifier");
1034 :
1035 48 : 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 284 : 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 5240 : for (const auto & vec : vectors)
1047 4956 : bf_out.push_back(vec);
1048 :
1049 : // for (auto & val : vectors)
1050 : // bf_out.push_back(val);
1051 284 : 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 284 : sys.get_mesh().fix_broken_node_and_element_numbering();
1063 552 : }
1064 :
1065 284 : 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 300 : read_in_vectors(sys,
1072 284 : basis_functions,
1073 : directory_name,
1074 : "bf",
1075 : read_binary_basis_functions);
1076 284 : }
1077 :
1078 284 : 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 284 : vectors_vec.push_back(&vectors);
1086 :
1087 24 : std::vector<std::string> directory_name_vec;
1088 284 : directory_name_vec.push_back(directory_name);
1089 :
1090 16 : std::vector<std::string> data_name_vec;
1091 284 : data_name_vec.push_back(data_name);
1092 :
1093 292 : read_in_vectors_from_multiple_files(sys,
1094 : vectors_vec,
1095 : directory_name_vec,
1096 : data_name_vec,
1097 : read_binary_vectors);
1098 552 : }
1099 :
1100 284 : 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 284 : if (n_files == 0)
1113 0 : return;
1114 :
1115 : // Make sure processors are synced up before we begin
1116 284 : sys.comm().barrier();
1117 :
1118 300 : std::ostringstream file_name;
1119 292 : 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 284 : MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());
1125 :
1126 568 : for (std::size_t data_index=0; data_index<n_directories; data_index++)
1127 : {
1128 292 : std::vector<std::unique_ptr<NumericVector<Number>>> & vectors = *multiple_vectors[data_index];
1129 :
1130 : // Allocate storage for each vector
1131 5240 : 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 4956 : libmesh_error_msg_if(vec, "Non-nullptr vector passed to read_in_vectors_from_multiple_files");
1136 :
1137 9772 : vec = NumericVector<Number>::build(sys.comm());
1138 :
1139 5096 : vec->init (sys.n_dofs(),
1140 : sys.n_local_dofs(),
1141 : false,
1142 280 : PARALLEL);
1143 : }
1144 :
1145 552 : file_name.str("");
1146 16 : file_name << multiple_directory_names[data_index]
1147 16 : << "/" << multiple_data_names[data_index]
1148 544 : << "_data" << basis_function_suffix;
1149 :
1150 : // On processor zero check to be sure the file exists
1151 284 : if (sys.comm().rank() == 0)
1152 : {
1153 : struct stat stat_info;
1154 48 : int stat_result = stat(file_name.str().c_str(), &stat_info);
1155 :
1156 48 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
1157 : }
1158 :
1159 284 : assert_file_exists(file_name.str());
1160 284 : Xdr vector_data(file_name.str(),
1161 584 : 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 284 : vector_data.data(version);
1167 :
1168 292 : const std::string libMesh_label = "libMesh-";
1169 8 : std::string::size_type lm_pos = version.find(libMesh_label);
1170 284 : libmesh_error_msg_if(lm_pos == std::string::npos, "version info missing in Xdr header");
1171 :
1172 300 : std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
1173 284 : int ver_major = 0, ver_minor = 0, ver_patch = 0;
1174 : char dot;
1175 284 : iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
1176 284 : 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 284 : sys.read_header(vector_data, version, /*read_header=*/false, /*read_additional_data=*/false);
1182 268 : }
1183 :
1184 : // Comply with the System::read_serialized_vectors() interface which uses dumb pointers.
1185 16 : std::vector<NumericVector<Number> *> vec_in;
1186 5240 : for (auto & vec : vectors)
1187 4956 : vec_in.push_back(vec.get());
1188 :
1189 8 : sys.read_serialized_vectors (vector_data, vec_in);
1190 268 : }
1191 :
1192 : // Undo the temporary renumbering
1193 284 : sys.get_mesh().fix_broken_node_and_element_numbering();
1194 268 : }
1195 :
1196 769828 : void RBEvaluation::check_evaluated_thetas_size(const std::vector<Number> * evaluated_thetas) const
1197 : {
1198 769828 : 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 769828 : }
1215 :
1216 : } // namespace libMesh
|
