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 133 : 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 : #ifdef LIBMESH_ENABLE_DEPRECATED
475 0 : Real RBEvaluation::eval_output_dual_norm(unsigned int n, const RBParameters & /*mu*/)
476 : {
477 : libmesh_deprecated();
478 :
479 : // Call non-deprecated version of this function, ignoring input mu
480 0 : return this->eval_output_dual_norm(n, nullptr);
481 : }
482 : #endif // LIBMESH_ENABLE_DEPRECATED
483 :
484 9015368 : Real RBEvaluation::eval_output_dual_norm(unsigned int n,
485 : const std::vector<Number> * evaluated_thetas)
486 : {
487 : // Return value
488 824820 : Number output_bound_sq = 0.;
489 :
490 : // mu is only used if evaluated_thetas == nullptr
491 9015368 : const RBParameters & mu = this->get_parameters();
492 :
493 : // Index into output_dual_innerprods
494 816816 : unsigned int q=0;
495 18030736 : for (unsigned int q_l1=0; q_l1<rb_theta_expansion->get_n_output_terms(n); q_l1++)
496 : {
497 18030736 : for (unsigned int q_l2=q_l1; q_l2<rb_theta_expansion->get_n_output_terms(n); q_l2++)
498 : {
499 9015368 : Real delta = (q_l1==q_l2) ? 1. : 2.;
500 :
501 : Number val_l1 =
502 9015368 : evaluated_thetas ?
503 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()] :
504 9015368 : rb_theta_expansion->eval_output_theta(n, q_l1, mu);
505 :
506 : Number val_l2 =
507 1633632 : evaluated_thetas ?
508 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()] :
509 9015368 : rb_theta_expansion->eval_output_theta(n, q_l2, mu);
510 :
511 9007364 : output_bound_sq += delta * libmesh_real(
512 9832184 : libmesh_conj(val_l1) * val_l2 * output_dual_innerprods[n][q]);
513 :
514 9015368 : q++;
515 : }
516 : }
517 :
518 9015368 : return libmesh_real(std::sqrt( output_bound_sq ));
519 : }
520 :
521 568 : void RBEvaluation::clear_riesz_representors()
522 : {
523 552 : Aq_representor.clear();
524 568 : }
525 :
526 284 : void RBEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
527 : const bool write_binary_data)
528 : {
529 16 : LOG_SCOPE("legacy_write_offline_data_to_files()", "RBEvaluation");
530 :
531 : // Get the number of basis functions
532 284 : unsigned int n_bfs = get_n_basis_functions();
533 :
534 : // The writing mode: ENCODE for binary, WRITE for ASCII
535 284 : XdrMODE mode = write_binary_data ? ENCODE : WRITE;
536 :
537 : // The suffix to use for all the files that are written out
538 292 : const std::string suffix = write_binary_data ? ".xdr" : ".dat";
539 :
540 292 : if (this->processor_id() == 0)
541 : {
542 :
543 : // Make a directory to store all the data files
544 48 : Utility::mkdir(directory_name.c_str());
545 : // if (mkdir(directory_name.c_str(), 0777) == -1)
546 : // {
547 : // libMesh::out << "In RBEvaluation::write_offline_data_to_files, directory "
548 : // << directory_name << " already exists, overwriting contents." << std::endl;
549 : // }
550 :
551 : // First, write out how many basis functions we have generated
552 56 : std::ostringstream file_name;
553 : {
554 44 : file_name << directory_name << "/n_bfs" << suffix;
555 96 : Xdr n_bfs_out(file_name.str(), mode);
556 8 : n_bfs_out << n_bfs;
557 48 : n_bfs_out.close();
558 40 : }
559 :
560 : // Write out the parameter ranges
561 92 : file_name.str("");
562 44 : file_name << directory_name << "/parameter_ranges" << suffix;
563 8 : std::string continuous_param_file_name = file_name.str();
564 :
565 : // Write out the discrete parameter values
566 92 : file_name.str("");
567 44 : file_name << directory_name << "/discrete_parameter_values" << suffix;
568 8 : std::string discrete_param_file_name = file_name.str();
569 :
570 48 : write_parameter_data_to_files(continuous_param_file_name,
571 : discrete_param_file_name,
572 : write_binary_data);
573 :
574 : // Write out Fq representor norm data
575 92 : file_name.str("");
576 44 : file_name << directory_name << "/Fq_innerprods" << suffix;
577 56 : Xdr RB_Fq_innerprods_out(file_name.str(), mode);
578 48 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
579 576 : for (unsigned int i=0; i<Q_f_hat; i++)
580 : {
581 572 : RB_Fq_innerprods_out << Fq_representor_innerprods[i];
582 : }
583 48 : RB_Fq_innerprods_out.close();
584 :
585 : // Write out output data
586 144 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
587 : {
588 184 : file_name.str("");
589 96 : file_name << directory_name << "/output_";
590 : file_name << std::setw(3)
591 : << std::setprecision(0)
592 8 : << std::setfill('0')
593 8 : << std::right
594 8 : << n;
595 :
596 88 : file_name << "_dual_innerprods" << suffix;
597 112 : Xdr output_dual_innerprods_out(file_name.str(), mode);
598 :
599 96 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
600 192 : for (unsigned int q=0; q<Q_l_hat; q++)
601 : {
602 112 : output_dual_innerprods_out << output_dual_innerprods[n][q];
603 : }
604 96 : output_dual_innerprods_out.close();
605 80 : }
606 :
607 :
608 : // Write out output data to multiple files
609 144 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
610 : {
611 192 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
612 : {
613 184 : file_name.str("");
614 96 : file_name << directory_name << "/output_";
615 : file_name << std::setw(3)
616 : << std::setprecision(0)
617 8 : << std::setfill('0')
618 8 : << std::right
619 8 : << n;
620 96 : file_name << "_";
621 : file_name << std::setw(3)
622 : << std::setprecision(0)
623 8 : << std::setfill('0')
624 8 : << std::right
625 8 : << q_l;
626 8 : file_name << suffix;
627 120 : Xdr output_n_out(file_name.str(), mode);
628 :
629 1960 : for (unsigned int j=0; j<n_bfs; j++)
630 : {
631 640 : output_n_out << RB_output_vectors[n][q_l](j);
632 : }
633 96 : output_n_out.close();
634 80 : }
635 : }
636 :
637 48 : if (compute_RB_inner_product)
638 : {
639 : // Next write out the inner product matrix
640 21 : file_name.str("");
641 10 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
642 14 : Xdr RB_inner_product_matrix_out(file_name.str(), mode);
643 231 : for (unsigned int i=0; i<n_bfs; i++)
644 : {
645 4620 : for (unsigned int j=0; j<n_bfs; j++)
646 : {
647 800 : RB_inner_product_matrix_out << RB_inner_product_matrix(i,j);
648 : }
649 : }
650 11 : RB_inner_product_matrix_out.close();
651 9 : }
652 :
653 : // Next write out the Fq vectors
654 216 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
655 : {
656 322 : file_name.str("");
657 168 : file_name << directory_name << "/RB_F_";
658 : file_name << std::setw(3)
659 : << std::setprecision(0)
660 14 : << std::setfill('0')
661 14 : << std::right
662 14 : << q_f;
663 14 : file_name << suffix;
664 210 : Xdr RB_Fq_f_out(file_name.str(), mode);
665 :
666 2794 : for (unsigned int i=0; i<n_bfs; i++)
667 : {
668 660 : RB_Fq_f_out << RB_Fq_vector[q_f](i);
669 : }
670 168 : RB_Fq_f_out.close();
671 140 : }
672 :
673 : // Next write out the Aq matrices
674 192 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
675 : {
676 276 : file_name.str("");
677 144 : file_name << directory_name << "/RB_A_";
678 : file_name << std::setw(3)
679 : << std::setprecision(0)
680 12 : << std::setfill('0')
681 12 : << std::right
682 12 : << q_a;
683 12 : file_name << suffix;
684 180 : Xdr RB_Aq_a_out(file_name.str(), mode);
685 :
686 2622 : for (unsigned int i=0; i<n_bfs; i++)
687 : {
688 46386 : for (unsigned int j=0; j<n_bfs; j++)
689 : {
690 11250 : RB_Aq_a_out << RB_Aq_vector[q_a](i,j);
691 : }
692 : }
693 144 : RB_Aq_a_out.close();
694 120 : }
695 :
696 : // Next write out Fq_Aq representor norm data
697 92 : file_name.str("");
698 44 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
699 60 : Xdr RB_Fq_Aq_innerprods_out(file_name.str(), mode);
700 :
701 216 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
702 : {
703 672 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
704 : {
705 8382 : for (unsigned int i=0; i<n_bfs; i++)
706 : {
707 9858 : RB_Fq_Aq_innerprods_out << Fq_Aq_representor_innerprods[q_f][q_a][i];
708 : }
709 : }
710 : }
711 48 : RB_Fq_Aq_innerprods_out.close();
712 :
713 : // Next write out Aq_Aq representor norm data
714 92 : file_name.str("");
715 44 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
716 56 : Xdr RB_Aq_Aq_innerprods_out(file_name.str(), mode);
717 :
718 48 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
719 336 : for (unsigned int i=0; i<Q_a_hat; i++)
720 : {
721 5244 : for (unsigned int j=0; j<n_bfs; j++)
722 : {
723 92772 : for (unsigned int l=0; l<n_bfs; l++)
724 : {
725 110316 : RB_Aq_Aq_innerprods_out << Aq_Aq_representor_innerprods[i][j][l];
726 : }
727 : }
728 : }
729 48 : RB_Aq_Aq_innerprods_out.close();
730 :
731 : // Also, write out the greedily selected parameters
732 : {
733 92 : file_name.str("");
734 44 : file_name << directory_name << "/greedy_params" << suffix;
735 100 : Xdr greedy_params_out(file_name.str(), mode);
736 :
737 922 : for (const auto & param : greedy_param_list)
738 4549 : for (const auto & pr : param)
739 7350 : for (const auto & value_vector : pr.second)
740 : {
741 : // Need to make a copy of the value so that it's not const
742 : // Xdr is not templated on const's
743 3983 : libmesh_error_msg_if(value_vector.size() != 1,
744 : "Error: multi-value RB parameters are not yet supported here.");
745 3675 : Real param_value = value_vector[0];
746 616 : greedy_params_out << param_value;
747 : }
748 48 : greedy_params_out.close();
749 40 : }
750 :
751 80 : }
752 284 : }
753 :
754 284 : void RBEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
755 : bool read_error_bound_data,
756 : const bool read_binary_data)
757 : {
758 16 : LOG_SCOPE("legacy_read_offline_data_from_files()", "RBEvaluation");
759 :
760 : // The reading mode: DECODE for binary, READ for ASCII
761 284 : XdrMODE mode = read_binary_data ? DECODE : READ;
762 :
763 : // The suffix to use for all the files that are written out
764 292 : const std::string suffix = read_binary_data ? ".xdr" : ".dat";
765 :
766 : // The string stream we'll use to make the file names
767 300 : std::ostringstream file_name;
768 :
769 : // First, find out how many basis functions we had when Greedy terminated
770 : unsigned int n_bfs;
771 : {
772 276 : file_name << directory_name << "/n_bfs" << suffix;
773 560 : assert_file_exists(file_name.str());
774 :
775 568 : Xdr n_bfs_in(file_name.str(), mode);
776 16 : n_bfs_in >> n_bfs;
777 284 : n_bfs_in.close();
778 268 : }
779 284 : resize_data_structures(n_bfs, read_error_bound_data);
780 :
781 : // Read in the parameter ranges
782 560 : file_name.str("");
783 276 : file_name << directory_name << "/parameter_ranges" << suffix;
784 16 : std::string continuous_param_file_name = file_name.str();
785 :
786 : // Read in the discrete parameter values
787 560 : file_name.str("");
788 276 : file_name << directory_name << "/discrete_parameter_values" << suffix;
789 16 : std::string discrete_param_file_name = file_name.str();
790 284 : read_parameter_data_from_files(continuous_param_file_name,
791 : discrete_param_file_name,
792 : read_binary_data);
793 :
794 : // Read in output data in multiple files
795 852 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
796 : {
797 1136 : for (unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
798 : {
799 1120 : file_name.str("");
800 568 : file_name << directory_name << "/output_";
801 : file_name << std::setw(3)
802 : << std::setprecision(0)
803 16 : << std::setfill('0')
804 16 : << std::right
805 16 : << n;
806 568 : file_name << "_";
807 : file_name << std::setw(3)
808 : << std::setprecision(0)
809 16 : << std::setfill('0')
810 16 : << std::right
811 16 : << q_l;
812 16 : file_name << suffix;
813 1120 : assert_file_exists(file_name.str());
814 :
815 616 : Xdr output_n_in(file_name.str(), mode);
816 :
817 11872 : for (unsigned int j=0; j<n_bfs; j++)
818 : {
819 80 : Number value;
820 640 : output_n_in >> value;
821 11944 : RB_output_vectors[n][q_l](j) = value;
822 : }
823 568 : output_n_in.close();
824 536 : }
825 : }
826 :
827 284 : if (compute_RB_inner_product)
828 : {
829 : // Next read in the inner product matrix
830 138 : file_name.str("");
831 68 : file_name << directory_name << "/RB_inner_product_matrix" << suffix;
832 138 : assert_file_exists(file_name.str());
833 :
834 76 : Xdr RB_inner_product_matrix_in(file_name.str(), mode);
835 :
836 1470 : for (unsigned int i=0; i<n_bfs; i++)
837 : {
838 29400 : for (unsigned int j=0; j<n_bfs; j++)
839 : {
840 0 : Number value;
841 1600 : RB_inner_product_matrix_in >> value;
842 28800 : RB_inner_product_matrix(i,j) = value;
843 : }
844 : }
845 70 : RB_inner_product_matrix_in.close();
846 66 : }
847 :
848 : // Next read in the Fq vectors
849 1278 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
850 : {
851 1960 : file_name.str("");
852 994 : file_name << directory_name << "/RB_F_";
853 : file_name << std::setw(3)
854 : << std::setprecision(0)
855 28 : << std::setfill('0')
856 28 : << std::right
857 28 : << q_f;
858 28 : file_name << suffix;
859 1960 : assert_file_exists(file_name.str());
860 :
861 1078 : Xdr RB_Fq_f_in(file_name.str(), mode);
862 :
863 16600 : for (unsigned int i=0; i<n_bfs; i++)
864 : {
865 200 : Number value;
866 880 : RB_Fq_f_in >> value;
867 16046 : RB_Fq_vector[q_f](i) = value;
868 : }
869 994 : RB_Fq_f_in.close();
870 938 : }
871 :
872 : // Next read in the Aq matrices
873 1136 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
874 : {
875 1680 : file_name.str("");
876 852 : file_name << directory_name << "/RB_A_";
877 : file_name << std::setw(3)
878 : << std::setprecision(0)
879 24 : << std::setfill('0')
880 24 : << std::right
881 24 : << q_a;
882 24 : file_name << suffix;
883 1680 : assert_file_exists(file_name.str());
884 :
885 924 : Xdr RB_Aq_a_in(file_name.str(), mode);
886 :
887 15720 : for (unsigned int i=0; i<n_bfs; i++)
888 : {
889 280026 : for (unsigned int j=0; j<n_bfs; j++)
890 : {
891 2550 : Number value;
892 15000 : RB_Aq_a_in >> value;
893 272658 : RB_Aq_vector[q_a](i,j) = value;
894 : }
895 : }
896 852 : RB_Aq_a_in.close();
897 804 : }
898 :
899 :
900 284 : if (read_error_bound_data)
901 : {
902 : // Next read in Fq representor norm data
903 560 : file_name.str("");
904 276 : file_name << directory_name << "/Fq_innerprods" << suffix;
905 560 : assert_file_exists(file_name.str());
906 :
907 300 : Xdr RB_Fq_innerprods_in(file_name.str(), mode);
908 :
909 284 : unsigned int Q_f_hat = rb_theta_expansion->get_n_F_terms()*(rb_theta_expansion->get_n_F_terms()+1)/2;
910 3408 : for (unsigned int i=0; i<Q_f_hat; i++)
911 : {
912 3212 : RB_Fq_innerprods_in >> Fq_representor_innerprods[i];
913 : }
914 284 : RB_Fq_innerprods_in.close();
915 :
916 : // Read in output data
917 852 : for (unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
918 : {
919 1120 : file_name.str("");
920 568 : file_name << directory_name << "/output_";
921 : file_name << std::setw(3)
922 : << std::setprecision(0)
923 16 : << std::setfill('0')
924 16 : << std::right
925 16 : << n;
926 552 : file_name << "_dual_innerprods" << suffix;
927 1120 : assert_file_exists(file_name.str());
928 :
929 600 : Xdr output_dual_innerprods_in(file_name.str(), mode);
930 :
931 568 : unsigned int Q_l_hat = rb_theta_expansion->get_n_output_terms(n)*(rb_theta_expansion->get_n_output_terms(n)+1)/2;
932 1136 : for (unsigned int q=0; q<Q_l_hat; q++)
933 : {
934 600 : output_dual_innerprods_in >> output_dual_innerprods[n][q];
935 : }
936 568 : output_dual_innerprods_in.close();
937 536 : }
938 :
939 :
940 : // Next read in Fq_Aq representor norm data
941 560 : file_name.str("");
942 276 : file_name << directory_name << "/Fq_Aq_innerprods" << suffix;
943 560 : assert_file_exists(file_name.str());
944 :
945 308 : Xdr RB_Fq_Aq_innerprods_in(file_name.str(), mode);
946 :
947 1278 : for (unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
948 : {
949 3976 : for (unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
950 : {
951 49800 : for (unsigned int i=0; i<n_bfs; i++)
952 : {
953 50778 : RB_Fq_Aq_innerprods_in >> Fq_Aq_representor_innerprods[q_f][q_a][i];
954 : }
955 : }
956 : }
957 284 : RB_Fq_Aq_innerprods_in.close();
958 :
959 : // Next read in Aq_Aq representor norm data
960 560 : file_name.str("");
961 276 : file_name << directory_name << "/Aq_Aq_innerprods" << suffix;
962 560 : assert_file_exists(file_name.str());
963 :
964 300 : Xdr RB_Aq_Aq_innerprods_in(file_name.str(), mode);
965 :
966 284 : unsigned int Q_a_hat = rb_theta_expansion->get_n_A_terms()*(rb_theta_expansion->get_n_A_terms()+1)/2;
967 1988 : for (unsigned int i=0; i<Q_a_hat; i++)
968 : {
969 31440 : for (unsigned int j=0; j<n_bfs; j++)
970 : {
971 560052 : for (unsigned int l=0; l<n_bfs; l++)
972 : {
973 575316 : RB_Aq_Aq_innerprods_in >> Aq_Aq_representor_innerprods[i][j][l];
974 : }
975 : }
976 : }
977 284 : RB_Aq_Aq_innerprods_in.close();
978 268 : }
979 :
980 : // Resize basis_functions even if we don't read them in so that
981 : // get_n_bfs() returns the correct value. Initialize the pointers
982 : // to nullptr.
983 276 : basis_functions.clear();
984 284 : set_n_basis_functions(n_bfs);
985 552 : }
986 :
987 5032 : void RBEvaluation::assert_file_exists(const std::string & file_name)
988 : {
989 5032 : libmesh_error_msg_if(!std::ifstream(file_name.c_str()), "File missing: " << file_name);
990 5032 : }
991 :
992 284 : void RBEvaluation::write_out_basis_functions(System & sys,
993 : const std::string & directory_name,
994 : const bool write_binary_basis_functions)
995 : {
996 16 : LOG_SCOPE("write_out_basis_functions()", "RBEvaluation");
997 :
998 8 : std::vector<NumericVector<Number>*> basis_functions_ptrs;
999 5388 : for (std::size_t i=0; i<basis_functions.size(); i++)
1000 : {
1001 4956 : basis_functions_ptrs.push_back(basis_functions[i].get());
1002 : }
1003 :
1004 560 : write_out_vectors(sys,
1005 : basis_functions_ptrs,
1006 : directory_name,
1007 : "bf",
1008 : write_binary_basis_functions);
1009 284 : }
1010 :
1011 284 : void RBEvaluation::write_out_vectors(System & sys,
1012 : std::vector<NumericVector<Number>*> & vectors,
1013 : const std::string & directory_name,
1014 : const std::string & data_name,
1015 : const bool write_binary_vectors)
1016 : {
1017 16 : LOG_SCOPE("write_out_vectors()", "RBEvaluation");
1018 :
1019 284 : if (sys.comm().rank() == 0)
1020 : {
1021 : // Make a directory to store all the data files
1022 48 : Utility::mkdir(directory_name.c_str());
1023 : }
1024 :
1025 : // Make sure processors are synced up before we begin
1026 284 : sys.comm().barrier();
1027 :
1028 300 : std::ostringstream file_name;
1029 292 : const std::string basis_function_suffix = (write_binary_vectors ? ".xdr" : ".dat");
1030 :
1031 544 : file_name << directory_name << "/" << data_name << "_data" << basis_function_suffix;
1032 276 : Xdr bf_data(file_name.str(),
1033 584 : write_binary_vectors ? ENCODE : WRITE);
1034 :
1035 : // Following EquationSystems::write(), we should only write the header information
1036 : // if we're proc 0
1037 284 : if (sys.comm().rank() == 0)
1038 : {
1039 96 : std::string version("libMesh-" + libMesh::get_io_compatibility_version());
1040 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1041 4 : version += " with infinite elements";
1042 : #endif
1043 48 : bf_data.data(version ,"# File Format Identifier");
1044 :
1045 48 : sys.write_header(bf_data, /*(unused arg)*/ version, /*write_additional_data=*/false);
1046 : }
1047 :
1048 : // Following EquationSystemsIO::write, we use a temporary numbering (node major)
1049 : // before writing out the data
1050 284 : MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());
1051 :
1052 : // Write all vectors at once.
1053 : {
1054 : // Note the API wants pointers to constant vectors, hence this...
1055 16 : std::vector<const NumericVector<Number> *> bf_out;
1056 5240 : for (const auto & vec : vectors)
1057 4956 : bf_out.push_back(vec);
1058 :
1059 : // for (auto & val : vectors)
1060 : // bf_out.push_back(val);
1061 284 : sys.write_serialized_vectors (bf_data, bf_out);
1062 : }
1063 :
1064 :
1065 : // set the current version
1066 8 : bf_data.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
1067 : LIBMESH_MINOR_VERSION,
1068 : LIBMESH_MICRO_VERSION));
1069 :
1070 :
1071 : // Undo the temporary renumbering
1072 284 : sys.get_mesh().fix_broken_node_and_element_numbering();
1073 552 : }
1074 :
1075 284 : void RBEvaluation::read_in_basis_functions(System & sys,
1076 : const std::string & directory_name,
1077 : const bool read_binary_basis_functions)
1078 : {
1079 8 : LOG_SCOPE("read_in_basis_functions()", "RBEvaluation");
1080 :
1081 300 : read_in_vectors(sys,
1082 284 : basis_functions,
1083 : directory_name,
1084 : "bf",
1085 : read_binary_basis_functions);
1086 284 : }
1087 :
1088 284 : void RBEvaluation::read_in_vectors(System & sys,
1089 : std::vector<std::unique_ptr<NumericVector<Number>>> & vectors,
1090 : const std::string & directory_name,
1091 : const std::string & data_name,
1092 : const bool read_binary_vectors)
1093 : {
1094 16 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>> *> vectors_vec;
1095 284 : vectors_vec.push_back(&vectors);
1096 :
1097 24 : std::vector<std::string> directory_name_vec;
1098 284 : directory_name_vec.push_back(directory_name);
1099 :
1100 16 : std::vector<std::string> data_name_vec;
1101 284 : data_name_vec.push_back(data_name);
1102 :
1103 292 : read_in_vectors_from_multiple_files(sys,
1104 : vectors_vec,
1105 : directory_name_vec,
1106 : data_name_vec,
1107 : read_binary_vectors);
1108 552 : }
1109 :
1110 284 : void RBEvaluation::read_in_vectors_from_multiple_files(System & sys,
1111 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>> *> multiple_vectors,
1112 : const std::vector<std::string> & multiple_directory_names,
1113 : const std::vector<std::string> & multiple_data_names,
1114 : const bool read_binary_vectors)
1115 : {
1116 8 : LOG_SCOPE("read_in_vectors_from_multiple_files()", "RBEvaluation");
1117 :
1118 16 : std::size_t n_files = multiple_vectors.size();
1119 16 : std::size_t n_directories = multiple_directory_names.size();
1120 8 : libmesh_assert((n_files == n_directories) && (n_files == multiple_data_names.size()));
1121 :
1122 284 : if (n_files == 0)
1123 0 : return;
1124 :
1125 : // Make sure processors are synced up before we begin
1126 284 : sys.comm().barrier();
1127 :
1128 300 : std::ostringstream file_name;
1129 292 : const std::string basis_function_suffix = (read_binary_vectors ? ".xdr" : ".dat");
1130 :
1131 : // Following EquationSystemsIO::read, we use a temporary numbering (node major)
1132 : // before writing out the data. For the sake of efficiency, we do this once for
1133 : // all the vectors that we read in.
1134 284 : MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());
1135 :
1136 568 : for (std::size_t data_index=0; data_index<n_directories; data_index++)
1137 : {
1138 292 : std::vector<std::unique_ptr<NumericVector<Number>>> & vectors = *multiple_vectors[data_index];
1139 :
1140 : // Allocate storage for each vector
1141 5240 : for (auto & vec : vectors)
1142 : {
1143 : // vectors should all be nullptr, otherwise we get a memory leak when
1144 : // we create the new vectors in RBEvaluation::read_in_vectors.
1145 4956 : libmesh_error_msg_if(vec, "Non-nullptr vector passed to read_in_vectors_from_multiple_files");
1146 :
1147 9772 : vec = NumericVector<Number>::build(sys.comm());
1148 :
1149 5096 : vec->init (sys.n_dofs(),
1150 : sys.n_local_dofs(),
1151 : false,
1152 280 : PARALLEL);
1153 : }
1154 :
1155 552 : file_name.str("");
1156 16 : file_name << multiple_directory_names[data_index]
1157 16 : << "/" << multiple_data_names[data_index]
1158 544 : << "_data" << basis_function_suffix;
1159 :
1160 : // On processor zero check to be sure the file exists
1161 284 : if (sys.comm().rank() == 0)
1162 : {
1163 : struct stat stat_info;
1164 48 : int stat_result = stat(file_name.str().c_str(), &stat_info);
1165 :
1166 48 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
1167 : }
1168 :
1169 284 : assert_file_exists(file_name.str());
1170 284 : Xdr vector_data(file_name.str(),
1171 584 : read_binary_vectors ? DECODE : READ);
1172 :
1173 : // Read the header data. This block of code is based on EquationSystems::_read_impl.
1174 : {
1175 16 : std::string version;
1176 284 : vector_data.data(version);
1177 :
1178 292 : const std::string libMesh_label = "libMesh-";
1179 8 : std::string::size_type lm_pos = version.find(libMesh_label);
1180 284 : libmesh_error_msg_if(lm_pos == std::string::npos, "version info missing in Xdr header");
1181 :
1182 300 : std::istringstream iss(version.substr(lm_pos + libMesh_label.size()));
1183 284 : int ver_major = 0, ver_minor = 0, ver_patch = 0;
1184 : char dot;
1185 284 : iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
1186 284 : vector_data.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));
1187 :
1188 : // Actually read the header data. When we do this, set read_header=false
1189 : // so that we do not reinit sys, since we assume that it has already been
1190 : // set up properly (e.g. the appropriate variables have already been added).
1191 284 : sys.read_header(vector_data, version, /*read_header=*/false, /*read_additional_data=*/false);
1192 268 : }
1193 :
1194 : // Comply with the System::read_serialized_vectors() interface which uses dumb pointers.
1195 16 : std::vector<NumericVector<Number> *> vec_in;
1196 5240 : for (auto & vec : vectors)
1197 4956 : vec_in.push_back(vec.get());
1198 :
1199 8 : sys.read_serialized_vectors (vector_data, vec_in);
1200 268 : }
1201 :
1202 : // Undo the temporary renumbering
1203 284 : sys.get_mesh().fix_broken_node_and_element_numbering();
1204 268 : }
1205 :
1206 769828 : void RBEvaluation::check_evaluated_thetas_size(const std::vector<Number> * evaluated_thetas) const
1207 : {
1208 769828 : if (rb_theta_expansion && evaluated_thetas)
1209 : {
1210 0 : auto actual_size = evaluated_thetas->size();
1211 : auto expected_size =
1212 0 : rb_theta_expansion->get_n_A_terms() +
1213 0 : rb_theta_expansion->get_n_F_terms() +
1214 0 : rb_theta_expansion->get_total_n_output_terms();
1215 :
1216 0 : libmesh_error_msg_if(actual_size != expected_size,
1217 : "ERROR: Evaluated thetas vector has size " <<
1218 : actual_size << ", but we expected " <<
1219 : rb_theta_expansion->get_n_A_terms() << " A term(s), " <<
1220 : rb_theta_expansion->get_n_F_terms() << " F term(s), and " <<
1221 : rb_theta_expansion->get_total_n_output_terms() << " output term(s), " <<
1222 : "for a total of " << expected_size << " terms.");
1223 : }
1224 769828 : }
1225 :
1226 : } // namespace libMesh
|
