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_eim_evaluation.h"
22 : #include "libmesh/rb_eim_theta.h"
23 : #include "libmesh/rb_parametrized_function.h"
24 : #include "libmesh/rb_evaluation.h"
25 : #include "libmesh/utility.h" // Utility::mkdir
26 :
27 : // libMesh includes
28 : #include "libmesh/xdr_cxx.h"
29 : #include "libmesh/libmesh_logging.h"
30 : #include "libmesh/replicated_mesh.h"
31 : #include "libmesh/elem.h"
32 : #include "libmesh/system.h"
33 : #include "libmesh/equation_systems.h"
34 : #include "libmesh/numeric_vector.h"
35 : #include "libmesh/quadrature.h"
36 : #include "libmesh/boundary_info.h"
37 : #include "timpi/parallel_implementation.h"
38 :
39 : // C++ includes
40 : #include <sstream>
41 : #include <fstream>
42 : #include <numeric> // std::accumulate
43 : #include <iterator> // std::advance
44 :
45 : namespace libMesh
46 : {
47 :
48 0 : RBEIMEvaluation::RBEIMEvaluation(const Parallel::Communicator & comm)
49 : :
50 : ParallelObject(comm),
51 0 : eim_error_indicator_normalization(RBEIMEvaluation::MAX_RHS),
52 0 : limit_eim_error_indicator_to_one(true),
53 0 : _rb_eim_solves_N(0),
54 0 : _preserve_rb_eim_solutions(false),
55 0 : _is_eim_error_indicator_active(false)
56 : {
57 0 : }
58 :
59 0 : RBEIMEvaluation::~RBEIMEvaluation() = default;
60 :
61 0 : void RBEIMEvaluation::clear()
62 : {
63 0 : _vec_eval_input.clear();
64 :
65 0 : _interpolation_points_comp.clear();
66 0 : _interpolation_points_spatial_indices.clear();
67 :
68 :
69 0 : _interpolation_matrix.resize(0,0);
70 :
71 : // Delete any RBTheta objects that were created
72 0 : _rb_eim_theta_objects.clear();
73 0 : }
74 :
75 0 : void RBEIMEvaluation::resize_data_structures(const unsigned int Nmax)
76 : {
77 : // Resize the data structures relevant to the EIM system
78 0 : _vec_eval_input.clear();
79 0 : _vec_eval_input.all_xyz.clear();
80 :
81 0 : _interpolation_points_comp.clear();
82 0 : _interpolation_points_spatial_indices.clear();
83 :
84 0 : _interpolation_matrix.resize(Nmax,Nmax);
85 0 : }
86 :
87 0 : void RBEIMEvaluation::set_parametrized_function(std::unique_ptr<RBParametrizedFunction> pf)
88 : {
89 0 : _parametrized_function = std::move(pf);
90 0 : }
91 :
92 0 : RBParametrizedFunction & RBEIMEvaluation::get_parametrized_function()
93 : {
94 0 : libmesh_error_msg_if(!_parametrized_function, "Parametrized function not initialized yet");
95 :
96 0 : return *_parametrized_function;
97 : }
98 :
99 0 : const RBParametrizedFunction & RBEIMEvaluation::get_parametrized_function() const
100 : {
101 0 : libmesh_error_msg_if(!_parametrized_function, "Parametrized function not initialized yet");
102 :
103 0 : return *_parametrized_function;
104 : }
105 :
106 0 : DenseVector<Number> RBEIMEvaluation::rb_eim_solve(DenseVector<Number> & EIM_rhs)
107 : {
108 0 : LOG_SCOPE("rb_eim_solve()", "RBEIMEvaluation");
109 :
110 0 : libmesh_error_msg_if(EIM_rhs.size() > get_n_basis_functions(),
111 : "Error: N cannot be larger than the number of basis functions in rb_solve");
112 :
113 0 : libmesh_error_msg_if(EIM_rhs.size()==0, "Error: N must be greater than 0 in rb_solve");
114 :
115 0 : const unsigned int N = EIM_rhs.size();
116 0 : DenseVector<Number> rb_eim_solution(N);
117 0 : DenseMatrix<Number> interpolation_matrix_N;
118 0 : _interpolation_matrix.get_principal_submatrix(N, interpolation_matrix_N);
119 :
120 0 : interpolation_matrix_N.lu_solve(EIM_rhs, rb_eim_solution);
121 :
122 0 : return rb_eim_solution;
123 0 : }
124 :
125 0 : void RBEIMEvaluation::rb_eim_solves(const std::vector<RBParameters> & mus,
126 : unsigned int N)
127 : {
128 0 : if (_preserve_rb_eim_solutions)
129 : {
130 : // In this case we preserve _rb_eim_solutions and hence we
131 : // just return immediately so that we skip updating
132 : // _rb_eim_solutions below. This is relevant in cases where
133 : // we set up _rb_eim_solutions elsewhere and we don't want
134 : // to override it.
135 0 : return;
136 : }
137 :
138 0 : libmesh_error_msg_if(N > get_n_basis_functions(),
139 : "Error: N cannot be larger than the number of basis functions in rb_eim_solves");
140 0 : libmesh_error_msg_if(N==0, "Error: N must be greater than 0 in rb_eim_solves");
141 :
142 : // If mus and N are the same as before, then we return early
143 0 : if ((_rb_eim_solves_mus == mus) && (_rb_eim_solves_N == N))
144 0 : return;
145 :
146 0 : LOG_SCOPE("rb_eim_solves()", "RBEIMEvaluation");
147 :
148 0 : _rb_eim_solves_mus = mus;
149 0 : _rb_eim_solves_N = N;
150 :
151 0 : if (get_parametrized_function().is_lookup_table)
152 : {
153 0 : _rb_eim_solutions.resize(mus.size());
154 0 : for (auto mu_index : index_range(mus))
155 : {
156 : Real lookup_table_param =
157 0 : mus[mu_index].get_value(get_parametrized_function().lookup_table_param_name);
158 :
159 : // Cast lookup_table_param to an unsigned integer so that we can use
160 : // it as an index into the EIM rhs values obtained from the lookup table.
161 : unsigned int lookup_table_index =
162 0 : cast_int<unsigned int>(std::round(lookup_table_param));
163 :
164 0 : DenseVector<Number> values;
165 0 : _eim_solutions_for_training_set[lookup_table_index].get_principal_subvector(N, values);
166 0 : _rb_eim_solutions[mu_index] = values;
167 : }
168 :
169 0 : return;
170 : }
171 :
172 : // output all comps indexing is as follows:
173 : // mu index --> interpolation point index --> component index --> value.
174 0 : std::vector<std::vector<std::vector<Number>>> output_all_comps;
175 0 : if (get_parametrized_function().on_mesh_sides())
176 0 : get_parametrized_function().side_vectorized_evaluate(mus, _vec_eval_input, output_all_comps);
177 0 : else if (get_parametrized_function().on_mesh_nodes())
178 0 : get_parametrized_function().node_vectorized_evaluate(mus, _vec_eval_input, output_all_comps);
179 : else
180 0 : get_parametrized_function().vectorized_evaluate(mus, _vec_eval_input, output_all_comps);
181 :
182 : // Previously we did one RB-EIM solve per input mu, but now we do
183 : // one RB-EIM solve per input mu, per sample. In order for this to
184 : // work, we require that all the input mu objects have the same
185 : // number of samples.
186 0 : auto n_samples_0 = mus[0].n_samples();
187 0 : for (const auto & mu : mus)
188 0 : libmesh_error_msg_if(mu.n_samples() != n_samples_0, "All RBParameters objects must have same n_samples()");
189 :
190 : // After we verified that all mus have the same number of samples,
191 : // the total number of RB-EIM solves is simply the number of mus
192 : // times the number of samples.
193 0 : unsigned int num_rb_eim_solves = mus.size() * n_samples_0;
194 :
195 : // A special case is when we are passed a single RBParameters object
196 : // with no parameters stored on it. In this case, we effectively
197 : // have Theta(mu) == const, and therefore we still need to do at
198 : // least one RB-EIM solve. In this case, we require that there is
199 : // only one entry in "mus" for simplicity.
200 0 : if (num_rb_eim_solves == 0 && mus[0].n_parameters() == 0)
201 : {
202 0 : libmesh_error_msg_if(mus.size() != 1, "Must pass in only a single RBParameters object when solving with no parameters.");
203 0 : num_rb_eim_solves = 1;
204 : }
205 :
206 0 : std::vector<std::vector<Number>> evaluated_values_at_interp_points(num_rb_eim_solves);
207 :
208 0 : std::vector<Number> evaluated_values_at_err_indicator_point;
209 0 : if (_is_eim_error_indicator_active)
210 0 : evaluated_values_at_err_indicator_point.resize(num_rb_eim_solves);
211 :
212 : // In this loop, counter goes from 0 to num_rb_eim_solves. The
213 : // purpose of this loop is to strip out the "columns" of the
214 : // output_all_comps array into rows.
215 : {
216 0 : unsigned int counter = 0;
217 0 : for (auto mu_index : index_range(mus))
218 0 : for (auto sample_index : make_range(mus[mu_index].n_samples()))
219 : {
220 : // Ignore compiler warnings about unused loop index
221 0 : libmesh_ignore(sample_index);
222 :
223 0 : evaluated_values_at_interp_points[counter].resize(N);
224 :
225 0 : for (unsigned int interp_pt_index=0; interp_pt_index<N; interp_pt_index++)
226 : {
227 0 : unsigned int comp = _interpolation_points_comp[interp_pt_index];
228 :
229 : // This line of code previously used "mu_index", now we use
230 : // "counter" handle the multi-sample RBParameters case.
231 0 : evaluated_values_at_interp_points[counter][interp_pt_index] =
232 0 : output_all_comps[counter][interp_pt_index][comp];
233 : }
234 :
235 0 : if (_is_eim_error_indicator_active)
236 : {
237 0 : unsigned int comp = _interpolation_points_comp[N];
238 :
239 0 : evaluated_values_at_err_indicator_point[counter] =
240 0 : output_all_comps[counter][N][comp];
241 : }
242 :
243 0 : counter++;
244 : }
245 :
246 : // Throw an error if we didn't do the required number of solves for
247 : // some reason
248 0 : libmesh_error_msg_if(counter != num_rb_eim_solves,
249 : "We should have done " << num_rb_eim_solves <<
250 : " solves, instead we did " << counter);
251 : }
252 :
253 0 : DenseMatrix<Number> interpolation_matrix_N;
254 0 : _interpolation_matrix.get_principal_submatrix(N, interpolation_matrix_N);
255 :
256 : // The number of RB EIM solutions is equal to the size of the
257 : // "evaluated_values_at_interp_points" vector which we determined
258 : // earlier.
259 0 : _rb_eim_solutions.resize(num_rb_eim_solves);
260 0 : if (_is_eim_error_indicator_active)
261 0 : _rb_eim_error_indicators.resize(num_rb_eim_solves);
262 :
263 : {
264 0 : unsigned int counter = 0;
265 0 : for (auto mu_index : index_range(mus))
266 0 : for (auto sample_index : make_range(mus[mu_index].n_samples()))
267 : {
268 : // Ignore compiler warnings about unused loop index
269 0 : libmesh_ignore(sample_index);
270 :
271 0 : DenseVector<Number> EIM_rhs = evaluated_values_at_interp_points[counter];
272 0 : interpolation_matrix_N.lu_solve(EIM_rhs, _rb_eim_solutions[counter]);
273 :
274 : // If we're using the EIM error indicator, then we compute it via the approach
275 : // proposed in Proposition 3.3 of "An empirical interpolation method: application
276 : // to efficient reduced-basis discretization of partial differential equations",
277 : // Barrault et al.
278 0 : if (_is_eim_error_indicator_active)
279 : {
280 0 : Number error_indicator_rhs = evaluated_values_at_err_indicator_point[counter];
281 0 : _rb_eim_error_indicators[counter] =
282 0 : get_eim_error_indicator(
283 0 : error_indicator_rhs, _rb_eim_solutions[counter], EIM_rhs);
284 : }
285 :
286 0 : counter++;
287 : }
288 : }
289 0 : }
290 :
291 0 : void RBEIMEvaluation::initialize_interpolation_points_spatial_indices()
292 : {
293 0 : _interpolation_points_spatial_indices.clear();
294 :
295 0 : get_parametrized_function().get_spatial_indices(_interpolation_points_spatial_indices,
296 0 : _vec_eval_input);
297 0 : }
298 :
299 0 : void RBEIMEvaluation::initialize_param_fn_spatial_indices()
300 : {
301 0 : get_parametrized_function().initialize_spatial_indices(_interpolation_points_spatial_indices,
302 0 : _vec_eval_input);
303 0 : }
304 :
305 0 : unsigned int RBEIMEvaluation::get_n_basis_functions() const
306 : {
307 0 : if (get_parametrized_function().on_mesh_sides())
308 0 : return _local_side_eim_basis_functions.size();
309 0 : else if (get_parametrized_function().on_mesh_nodes())
310 0 : return _local_node_eim_basis_functions.size();
311 : else
312 0 : return _local_eim_basis_functions.size();
313 : }
314 :
315 0 : unsigned int RBEIMEvaluation::get_n_interpolation_points() const
316 : {
317 0 : return _vec_eval_input.all_xyz.size();
318 : }
319 :
320 0 : unsigned int RBEIMEvaluation::get_n_elems() const
321 : {
322 0 : return _vec_eval_input.elem_id_to_local_index.size();
323 : }
324 :
325 0 : unsigned int RBEIMEvaluation::get_n_properties() const
326 : {
327 0 : return _vec_eval_input.rb_property_map.size();
328 : }
329 :
330 0 : void RBEIMEvaluation::set_n_basis_functions(unsigned int n_bfs)
331 : {
332 0 : if (get_parametrized_function().on_mesh_sides())
333 0 : _local_side_eim_basis_functions.resize(n_bfs);
334 0 : else if (get_parametrized_function().on_mesh_nodes())
335 0 : _local_node_eim_basis_functions.resize(n_bfs);
336 : else
337 0 : _local_eim_basis_functions.resize(n_bfs);
338 0 : }
339 :
340 0 : void RBEIMEvaluation::decrement_vector(QpDataMap & v,
341 : const DenseVector<Number> & coeffs)
342 : {
343 0 : LOG_SCOPE("decrement_vector()", "RBEIMEvaluation");
344 :
345 0 : libmesh_error_msg_if(get_n_basis_functions() != coeffs.size(),
346 : "Error: Number of coefficients should match number of basis functions");
347 :
348 0 : for (auto & [elem_id, v_comp_and_qp] : v)
349 : {
350 0 : for (const auto & comp : index_range(v_comp_and_qp))
351 0 : for (unsigned int qp : index_range(v_comp_and_qp[comp]))
352 0 : for (unsigned int i : index_range(_local_eim_basis_functions))
353 : {
354 : // Check that entry (elem_id,comp,qp) exists in _local_eim_basis_functions so that
355 : // we get a clear error message if there is any missing data
356 0 : const auto & basis_comp_and_qp = libmesh_map_find(_local_eim_basis_functions[i], elem_id);
357 :
358 0 : libmesh_error_msg_if(comp >= basis_comp_and_qp.size(), "Error: Invalid comp");
359 0 : libmesh_error_msg_if(qp >= basis_comp_and_qp[comp].size(), "Error: Invalid qp");
360 :
361 0 : v_comp_and_qp[comp][qp] -= coeffs(i) * basis_comp_and_qp[comp][qp];
362 : }
363 : }
364 0 : }
365 :
366 0 : void RBEIMEvaluation::side_decrement_vector(SideQpDataMap & v,
367 : const DenseVector<Number> & coeffs)
368 : {
369 0 : LOG_SCOPE("side_decrement_vector()", "RBEIMEvaluation");
370 :
371 0 : libmesh_error_msg_if(get_n_basis_functions() != coeffs.size(),
372 : "Error: Number of coefficients should match number of basis functions");
373 :
374 0 : for (auto & [elem_and_side, v_comp_and_qp] : v)
375 : {
376 0 : for (const auto & comp : index_range(v_comp_and_qp))
377 0 : for (unsigned int qp : index_range(v_comp_and_qp[comp]))
378 0 : for (unsigned int i : index_range(_local_side_eim_basis_functions))
379 : {
380 : // Check that entry (elem_and_side,comp,qp) exists in _local_side_eim_basis_functions so that
381 : // we get a clear error message if there is any missing data
382 0 : const auto & basis_comp_and_qp = libmesh_map_find(_local_side_eim_basis_functions[i], elem_and_side);
383 :
384 0 : libmesh_error_msg_if(comp >= basis_comp_and_qp.size(), "Error: Invalid comp");
385 0 : libmesh_error_msg_if(qp >= basis_comp_and_qp[comp].size(), "Error: Invalid qp");
386 :
387 0 : v_comp_and_qp[comp][qp] -= coeffs(i) * basis_comp_and_qp[comp][qp];
388 : }
389 : }
390 0 : }
391 :
392 0 : void RBEIMEvaluation::node_decrement_vector(NodeDataMap & v,
393 : const DenseVector<Number> & coeffs)
394 : {
395 0 : LOG_SCOPE("node_decrement_vector()", "RBEIMEvaluation");
396 :
397 0 : libmesh_error_msg_if(get_n_basis_functions() != coeffs.size(),
398 : "Error: Number of coefficients should match number of basis functions");
399 :
400 0 : for (auto & [node_id, v_comps] : v)
401 : {
402 0 : for (const auto & comp : index_range(v_comps))
403 0 : for (unsigned int i : index_range(_local_node_eim_basis_functions))
404 : {
405 : // Check that entry (node_id,comp) exists in _local_node_eim_basis_functions so that
406 : // we get a clear error message if there is any missing data
407 0 : const auto & basis_comp = libmesh_map_find(_local_node_eim_basis_functions[i], node_id);
408 :
409 0 : libmesh_error_msg_if(comp >= basis_comp.size(), "Error: Invalid comp");
410 :
411 0 : v_comps[comp] -= coeffs(i) * basis_comp[comp];
412 : }
413 : }
414 0 : }
415 :
416 0 : void RBEIMEvaluation::initialize_eim_theta_objects()
417 : {
418 : // Initialize the rb_theta objects that access the solution from this rb_eim_evaluation
419 0 : _rb_eim_theta_objects.clear();
420 0 : for (auto i : make_range(get_n_basis_functions()))
421 0 : _rb_eim_theta_objects.emplace_back(build_eim_theta(i));
422 0 : }
423 :
424 0 : std::vector<std::unique_ptr<RBTheta>> & RBEIMEvaluation::get_eim_theta_objects()
425 : {
426 0 : return _rb_eim_theta_objects;
427 : }
428 :
429 0 : std::unique_ptr<RBTheta> RBEIMEvaluation::build_eim_theta(unsigned int index)
430 : {
431 0 : return std::make_unique<RBEIMTheta>(*this, index);
432 : }
433 :
434 0 : void RBEIMEvaluation::get_parametrized_function_values_at_qps(
435 : const QpDataMap & pf,
436 : dof_id_type elem_id,
437 : unsigned int comp,
438 : std::vector<Number> & values)
439 : {
440 0 : LOG_SCOPE("get_parametrized_function_values_at_qps()", "RBEIMConstruction");
441 :
442 0 : values.clear();
443 :
444 0 : if (const auto it = pf.find(elem_id);
445 0 : it != pf.end())
446 : {
447 0 : const auto & comps_and_qps_on_elem = it->second;
448 0 : libmesh_error_msg_if(comp >= comps_and_qps_on_elem.size(),
449 : "Invalid comp index: " << comp);
450 :
451 0 : values = comps_and_qps_on_elem[comp];
452 : }
453 0 : }
454 :
455 0 : void RBEIMEvaluation::get_parametrized_function_side_values_at_qps(
456 : const SideQpDataMap & pf,
457 : dof_id_type elem_id,
458 : unsigned int side_index,
459 : unsigned int comp,
460 : std::vector<Number> & values)
461 : {
462 0 : LOG_SCOPE("get_parametrized_function_side_values_at_qps()", "RBEIMConstruction");
463 :
464 0 : values.clear();
465 :
466 0 : if (const auto it = pf.find(std::make_pair(elem_id, side_index));
467 0 : it != pf.end())
468 : {
469 0 : const auto & comps_and_qps_on_elem = it->second;
470 0 : libmesh_error_msg_if(comp >= comps_and_qps_on_elem.size(),
471 : "Invalid comp index: " << comp);
472 :
473 0 : values = comps_and_qps_on_elem[comp];
474 : }
475 0 : }
476 :
477 0 : Number RBEIMEvaluation::get_parametrized_function_node_local_value(
478 : const NodeDataMap & pf,
479 : dof_id_type node_id,
480 : unsigned int comp)
481 : {
482 0 : LOG_SCOPE("get_parametrized_function_node_local_value()", "RBEIMConstruction");
483 :
484 0 : if (const auto it = pf.find(node_id);
485 0 : it != pf.end())
486 : {
487 0 : const std::vector<Number> & vec = it->second;
488 0 : libmesh_error_msg_if (comp >= vec.size(), "Error: Invalid comp index");
489 0 : return vec[comp];
490 : }
491 : else
492 0 : return 0.;
493 : }
494 :
495 0 : Number RBEIMEvaluation::get_parametrized_function_value(
496 : const Parallel::Communicator & comm,
497 : const QpDataMap & pf,
498 : dof_id_type elem_id,
499 : unsigned int comp,
500 : unsigned int qp)
501 : {
502 0 : std::vector<Number> values;
503 0 : get_parametrized_function_values_at_qps(pf, elem_id, comp, values);
504 :
505 : // In parallel, values should only be non-empty on one processor
506 0 : Number value = 0.;
507 0 : if (!values.empty())
508 : {
509 0 : libmesh_error_msg_if(qp >= values.size(), "Error: Invalid qp index");
510 :
511 0 : value = values[qp];
512 : }
513 0 : comm.sum(value);
514 :
515 0 : return value;
516 : }
517 :
518 0 : Number RBEIMEvaluation::get_parametrized_function_side_value(
519 : const Parallel::Communicator & comm,
520 : const SideQpDataMap & pf,
521 : dof_id_type elem_id,
522 : unsigned int side_index,
523 : unsigned int comp,
524 : unsigned int qp)
525 : {
526 0 : std::vector<Number> values;
527 0 : get_parametrized_function_side_values_at_qps(pf, elem_id, side_index, comp, values);
528 :
529 : // In parallel, values should only be non-empty on one processor
530 0 : Number value = 0.;
531 0 : if (!values.empty())
532 : {
533 0 : libmesh_error_msg_if(qp >= values.size(), "Error: Invalid qp index");
534 :
535 0 : value = values[qp];
536 : }
537 0 : comm.sum(value);
538 :
539 0 : return value;
540 : }
541 :
542 0 : Number RBEIMEvaluation::get_parametrized_function_node_value(
543 : const Parallel::Communicator & comm,
544 : const NodeDataMap & pf,
545 : dof_id_type node_id,
546 : unsigned int comp)
547 : {
548 0 : LOG_SCOPE("get_parametrized_function_node_value()", "RBEIMConstruction");
549 :
550 0 : Number value = get_parametrized_function_node_local_value(pf, node_id, comp);
551 0 : comm.sum(value);
552 :
553 0 : return value;
554 : }
555 :
556 0 : void RBEIMEvaluation::get_eim_basis_function_values_at_qps(unsigned int basis_function_index,
557 : dof_id_type elem_id,
558 : unsigned int comp,
559 : std::vector<Number> & values) const
560 : {
561 0 : libmesh_error_msg_if(basis_function_index >= _local_eim_basis_functions.size(),
562 : "Invalid basis function index: " << basis_function_index);
563 :
564 0 : get_parametrized_function_values_at_qps(
565 0 : _local_eim_basis_functions[basis_function_index],
566 : elem_id,
567 : comp,
568 : values);
569 0 : }
570 :
571 0 : void RBEIMEvaluation::get_eim_basis_function_side_values_at_qps(unsigned int basis_function_index,
572 : dof_id_type elem_id,
573 : unsigned int side_index,
574 : unsigned int comp,
575 : std::vector<Number> & values) const
576 : {
577 0 : libmesh_error_msg_if(basis_function_index >= _local_side_eim_basis_functions.size(),
578 : "Invalid basis function index: " << basis_function_index);
579 :
580 0 : get_parametrized_function_side_values_at_qps(
581 0 : _local_side_eim_basis_functions[basis_function_index],
582 : elem_id,
583 : side_index,
584 : comp,
585 : values);
586 0 : }
587 :
588 0 : Number RBEIMEvaluation::get_eim_basis_function_node_local_value(unsigned int basis_function_index,
589 : dof_id_type node_id,
590 : unsigned int comp) const
591 : {
592 0 : libmesh_error_msg_if(basis_function_index >= _local_node_eim_basis_functions.size(),
593 : "Invalid basis function index: " << basis_function_index);
594 :
595 0 : return get_parametrized_function_node_local_value(
596 0 : _local_node_eim_basis_functions[basis_function_index],
597 : node_id,
598 0 : comp);
599 : }
600 :
601 0 : Number RBEIMEvaluation::get_eim_basis_function_node_value(unsigned int basis_function_index,
602 : dof_id_type node_id,
603 : unsigned int comp) const
604 : {
605 0 : libmesh_error_msg_if(basis_function_index >= _local_node_eim_basis_functions.size(),
606 : "Invalid basis function index: " << basis_function_index);
607 :
608 0 : return get_parametrized_function_node_value(
609 : comm(),
610 0 : _local_node_eim_basis_functions[basis_function_index],
611 : node_id,
612 0 : comp);
613 : }
614 :
615 0 : Number RBEIMEvaluation::get_eim_basis_function_value(unsigned int basis_function_index,
616 : dof_id_type elem_id,
617 : unsigned int comp,
618 : unsigned int qp) const
619 : {
620 0 : libmesh_error_msg_if(basis_function_index >= _local_eim_basis_functions.size(),
621 : "Invalid basis function index: " << basis_function_index);
622 :
623 0 : return get_parametrized_function_value(
624 : comm(),
625 0 : _local_eim_basis_functions[basis_function_index],
626 : elem_id,
627 : comp,
628 0 : qp);
629 : }
630 :
631 0 : Number RBEIMEvaluation::get_eim_basis_function_side_value(unsigned int basis_function_index,
632 : dof_id_type elem_id,
633 : unsigned int side_index,
634 : unsigned int comp,
635 : unsigned int qp) const
636 : {
637 0 : libmesh_error_msg_if(basis_function_index >= _local_side_eim_basis_functions.size(),
638 : "Invalid side basis function index: " << basis_function_index);
639 :
640 0 : return get_parametrized_function_side_value(
641 : comm(),
642 0 : _local_side_eim_basis_functions[basis_function_index],
643 : elem_id,
644 : side_index,
645 : comp,
646 0 : qp);
647 : }
648 :
649 : const RBEIMEvaluation::QpDataMap &
650 0 : RBEIMEvaluation::get_basis_function(unsigned int i) const
651 : {
652 0 : return _local_eim_basis_functions[i];
653 : }
654 :
655 : const RBEIMEvaluation::SideQpDataMap &
656 0 : RBEIMEvaluation::get_side_basis_function(unsigned int i) const
657 : {
658 0 : return _local_side_eim_basis_functions[i];
659 : }
660 :
661 : const RBEIMEvaluation::NodeDataMap &
662 0 : RBEIMEvaluation::get_node_basis_function(unsigned int i) const
663 : {
664 0 : return _local_node_eim_basis_functions[i];
665 : }
666 :
667 0 : void RBEIMEvaluation::set_rb_eim_solutions(const std::vector<DenseVector<Number>> & rb_eim_solutions)
668 : {
669 0 : _rb_eim_solutions = rb_eim_solutions;
670 0 : }
671 :
672 0 : const std::vector<DenseVector<Number>> & RBEIMEvaluation::get_rb_eim_solutions() const
673 : {
674 0 : return _rb_eim_solutions;
675 : }
676 :
677 0 : std::vector<Number> RBEIMEvaluation::get_rb_eim_solutions_entries(unsigned int index) const
678 : {
679 0 : LOG_SCOPE("get_rb_eim_solutions_entries()", "RBEIMEvaluation");
680 :
681 0 : std::vector<Number> rb_eim_solutions_entries(_rb_eim_solutions.size());
682 0 : for (unsigned int mu_index : index_range(_rb_eim_solutions))
683 : {
684 0 : libmesh_error_msg_if(index >= _rb_eim_solutions[mu_index].size(),
685 : "Error: Requested solution index " << index <<
686 : ", but only have " << _rb_eim_solutions[mu_index].size() << " entries.");
687 0 : rb_eim_solutions_entries[mu_index] = _rb_eim_solutions[mu_index](index);
688 : }
689 :
690 0 : return rb_eim_solutions_entries;
691 : }
692 :
693 0 : const std::vector<DenseVector<Number>> & RBEIMEvaluation::get_eim_solutions_for_training_set() const
694 : {
695 0 : return _eim_solutions_for_training_set;
696 : }
697 :
698 0 : std::vector<DenseVector<Number>> & RBEIMEvaluation::get_eim_solutions_for_training_set()
699 : {
700 0 : return _eim_solutions_for_training_set;
701 : }
702 :
703 0 : const std::vector<std::pair<Real,Real>> & RBEIMEvaluation::get_rb_eim_error_indicators() const
704 : {
705 0 : return _rb_eim_error_indicators;
706 : }
707 :
708 0 : void RBEIMEvaluation::add_interpolation_points_xyz(Point p)
709 : {
710 0 : _vec_eval_input.all_xyz.emplace_back(p);
711 0 : }
712 :
713 0 : void RBEIMEvaluation::add_interpolation_points_comp(unsigned int comp)
714 : {
715 0 : _interpolation_points_comp.emplace_back(comp);
716 0 : }
717 :
718 0 : void RBEIMEvaluation::add_interpolation_points_subdomain_id(subdomain_id_type sbd_id)
719 : {
720 0 : _vec_eval_input.sbd_ids.emplace_back(sbd_id);
721 0 : }
722 :
723 0 : void RBEIMEvaluation::add_interpolation_points_boundary_id(boundary_id_type b_id)
724 : {
725 0 : _vec_eval_input.boundary_ids.emplace_back(b_id);
726 0 : }
727 :
728 0 : void RBEIMEvaluation::add_interpolation_points_xyz_perturbations(const std::vector<Point> & perturbs)
729 : {
730 0 : _vec_eval_input.all_xyz_perturb.emplace_back(perturbs);
731 0 : }
732 :
733 0 : void RBEIMEvaluation::add_interpolation_points_elem_id(dof_id_type elem_id)
734 : {
735 0 : _vec_eval_input.elem_ids.emplace_back(elem_id);
736 0 : }
737 :
738 0 : void RBEIMEvaluation::add_interpolation_points_side_index(unsigned int side_index)
739 : {
740 0 : _vec_eval_input.side_indices.emplace_back(side_index);
741 0 : }
742 :
743 0 : void RBEIMEvaluation::add_interpolation_points_node_id(dof_id_type node_id)
744 : {
745 0 : _vec_eval_input.node_ids.emplace_back(node_id);
746 0 : }
747 :
748 0 : void RBEIMEvaluation::add_interpolation_points_qp(unsigned int qp)
749 : {
750 0 : _vec_eval_input.qps.emplace_back(qp);
751 0 : }
752 :
753 0 : void RBEIMEvaluation::add_interpolation_points_elem_type(ElemType elem_type)
754 : {
755 0 : _vec_eval_input.elem_types.emplace_back(elem_type);
756 0 : }
757 :
758 0 : void RBEIMEvaluation::add_interpolation_points_JxW_all_qp(const std::vector<Real> & JxW_all_qp)
759 : {
760 0 : _vec_eval_input.JxW_all_qp.emplace_back(JxW_all_qp);
761 0 : }
762 :
763 0 : void RBEIMEvaluation::add_interpolation_points_phi_i_all_qp(const std::vector<std::vector<Real>> & phi_i_all_qp)
764 : {
765 0 : _vec_eval_input.phi_i_all_qp.emplace_back(phi_i_all_qp);
766 0 : }
767 :
768 0 : void RBEIMEvaluation::add_interpolation_points_phi_i_qp(const std::vector<Real> & phi_i_qp)
769 : {
770 0 : _vec_eval_input.phi_i_qp.emplace_back(phi_i_qp);
771 0 : }
772 :
773 0 : void RBEIMEvaluation::add_elem_center_dxyzdxi(const Point & dxyzdxi)
774 : {
775 0 : _vec_eval_input.dxyzdxi_elem_center.emplace_back(dxyzdxi);
776 0 : }
777 :
778 0 : void RBEIMEvaluation::add_elem_center_dxyzdeta(const Point & dxyzdeta)
779 : {
780 0 : _vec_eval_input.dxyzdeta_elem_center.emplace_back(dxyzdeta);
781 0 : }
782 :
783 0 : void RBEIMEvaluation::add_interpolation_points_qrule_order(Order qrule_order)
784 : {
785 0 : _vec_eval_input.qrule_orders.emplace_back(qrule_order);
786 0 : }
787 :
788 0 : void RBEIMEvaluation::add_interpolation_points_spatial_indices(const std::vector<unsigned int> & spatial_indices)
789 : {
790 0 : _interpolation_points_spatial_indices.emplace_back(spatial_indices);
791 0 : }
792 :
793 0 : void RBEIMEvaluation::add_elem_id_local_index_map_entry(dof_id_type elem_id, unsigned int local_index)
794 : {
795 : // Try to insert object and return an error if object already inserted as duplicated should not happen.
796 0 : bool insert_succeed = _vec_eval_input.elem_id_to_local_index.insert({elem_id, local_index}).second;
797 0 : libmesh_error_msg_if(!insert_succeed, "Entry already added, duplicate detected.");
798 0 : }
799 :
800 0 : void RBEIMEvaluation::add_rb_property_map_entry(std::string & property_name, std::set<dof_id_type> & entity_ids)
801 : {
802 0 : bool insert_succeed = _vec_eval_input.rb_property_map.insert({property_name, entity_ids}).second;
803 0 : libmesh_error_msg_if(!insert_succeed, "Entry already added, duplicate detected.");
804 0 : }
805 :
806 0 : Point RBEIMEvaluation::get_interpolation_points_xyz(unsigned int index) const
807 : {
808 0 : libmesh_error_msg_if(index >= _vec_eval_input.all_xyz.size(), "Error: Invalid index");
809 :
810 0 : return _vec_eval_input.all_xyz[index];
811 : }
812 :
813 0 : unsigned int RBEIMEvaluation::get_interpolation_points_comp(unsigned int index) const
814 : {
815 0 : libmesh_error_msg_if(index >= _interpolation_points_comp.size(), "Error: Invalid index");
816 :
817 0 : return _interpolation_points_comp[index];
818 : }
819 :
820 0 : subdomain_id_type RBEIMEvaluation::get_interpolation_points_subdomain_id(unsigned int index) const
821 : {
822 0 : libmesh_error_msg_if(index >= _vec_eval_input.sbd_ids.size(), "Error: Invalid index");
823 :
824 0 : return _vec_eval_input.sbd_ids[index];
825 : }
826 :
827 0 : boundary_id_type RBEIMEvaluation::get_interpolation_points_boundary_id(unsigned int index) const
828 : {
829 0 : libmesh_error_msg_if(index >= _vec_eval_input.boundary_ids.size(), "Error: Invalid index");
830 :
831 0 : return _vec_eval_input.boundary_ids[index];
832 : }
833 :
834 0 : const std::vector<Point> & RBEIMEvaluation::get_interpolation_points_xyz_perturbations(unsigned int index) const
835 : {
836 0 : libmesh_error_msg_if(index >= _vec_eval_input.all_xyz_perturb.size(), "Error: Invalid index");
837 :
838 0 : return _vec_eval_input.all_xyz_perturb[index];
839 : }
840 :
841 0 : dof_id_type RBEIMEvaluation::get_interpolation_points_elem_id(unsigned int index) const
842 : {
843 0 : libmesh_error_msg_if(index >= _vec_eval_input.elem_ids.size(), "Error: Invalid index");
844 :
845 0 : return _vec_eval_input.elem_ids[index];
846 : }
847 :
848 0 : unsigned int RBEIMEvaluation::get_interpolation_points_side_index(unsigned int index) const
849 : {
850 0 : libmesh_error_msg_if(index >= _vec_eval_input.side_indices.size(), "Error: Invalid index");
851 :
852 0 : return _vec_eval_input.side_indices[index];
853 : }
854 :
855 0 : dof_id_type RBEIMEvaluation::get_interpolation_points_node_id(unsigned int index) const
856 : {
857 0 : libmesh_error_msg_if(index >= _vec_eval_input.node_ids.size(), "Error: Invalid index");
858 :
859 0 : return _vec_eval_input.node_ids[index];
860 : }
861 :
862 0 : unsigned int RBEIMEvaluation::get_interpolation_points_qp(unsigned int index) const
863 : {
864 0 : libmesh_error_msg_if(index >= _vec_eval_input.qps.size(), "Error: Invalid index");
865 :
866 0 : return _vec_eval_input.qps[index];
867 : }
868 :
869 0 : ElemType RBEIMEvaluation::get_interpolation_points_elem_type(unsigned int index) const
870 : {
871 0 : libmesh_error_msg_if(index >= _vec_eval_input.elem_types.size(), "Error: Invalid index");
872 :
873 0 : return _vec_eval_input.elem_types[index];
874 : }
875 :
876 0 : const std::vector<Real> & RBEIMEvaluation::get_interpolation_points_JxW_all_qp(unsigned int index) const
877 : {
878 0 : libmesh_error_msg_if(index >= _vec_eval_input.JxW_all_qp.size(), "Error: Invalid index");
879 :
880 0 : return _vec_eval_input.JxW_all_qp[index];
881 : }
882 :
883 0 : const std::map<dof_id_type, unsigned int> & RBEIMEvaluation::get_elem_id_to_local_index_map() const
884 : {
885 0 : return _vec_eval_input.elem_id_to_local_index;
886 : }
887 :
888 0 : const std::unordered_map<std::string, std::set<dof_id_type>> & RBEIMEvaluation::get_rb_property_map() const
889 : {
890 0 : return _vec_eval_input.rb_property_map;
891 : }
892 :
893 0 : const std::vector<std::vector<Real>> & RBEIMEvaluation::get_interpolation_points_phi_i_all_qp(unsigned int index) const
894 : {
895 0 : libmesh_error_msg_if(index >= _vec_eval_input.phi_i_all_qp.size(), "Error: Invalid index");
896 :
897 0 : return _vec_eval_input.phi_i_all_qp[index];
898 : }
899 :
900 0 : const std::vector<Real> & RBEIMEvaluation::get_interpolation_points_phi_i_qp(unsigned int index) const
901 : {
902 0 : libmesh_error_msg_if(index >= _vec_eval_input.phi_i_qp.size(), "Error: Invalid index");
903 :
904 0 : return _vec_eval_input.phi_i_qp[index];
905 : }
906 :
907 0 : const Point & RBEIMEvaluation::get_elem_center_dxyzdxi(unsigned int index) const
908 : {
909 0 : libmesh_error_msg_if(index >= _vec_eval_input.dxyzdxi_elem_center.size(), "Error: Invalid index");
910 :
911 0 : return _vec_eval_input.dxyzdxi_elem_center[index];
912 : }
913 :
914 0 : const Point & RBEIMEvaluation::get_elem_center_dxyzdeta(unsigned int index) const
915 : {
916 0 : libmesh_error_msg_if(index >= _vec_eval_input.dxyzdeta_elem_center.size(), "Error: Invalid index");
917 :
918 0 : return _vec_eval_input.dxyzdeta_elem_center[index];
919 : }
920 :
921 0 : Order RBEIMEvaluation::get_interpolation_points_qrule_order(unsigned int index) const
922 : {
923 0 : libmesh_error_msg_if(index >= _vec_eval_input.qrule_orders.size(), "Error: Invalid index");
924 :
925 0 : return _vec_eval_input.qrule_orders[index];
926 : }
927 :
928 0 : const std::vector<unsigned int> & RBEIMEvaluation::get_interpolation_points_spatial_indices(unsigned int index) const
929 : {
930 0 : libmesh_error_msg_if(index >= _interpolation_points_spatial_indices.size(), "Error: Invalid index");
931 :
932 0 : return _interpolation_points_spatial_indices[index];
933 : }
934 :
935 0 : unsigned int RBEIMEvaluation::get_n_interpolation_points_spatial_indices() const
936 : {
937 0 : return _interpolation_points_spatial_indices.size();
938 : }
939 :
940 0 : void RBEIMEvaluation::set_interpolation_matrix_entry(unsigned int i, unsigned int j, Number value)
941 : {
942 0 : libmesh_error_msg_if((i >= _interpolation_matrix.m()) || (j >= _interpolation_matrix.n()),
943 : "Error: Invalid matrix indices");
944 :
945 0 : _interpolation_matrix(i,j) = value;
946 0 : }
947 :
948 0 : const DenseMatrix<Number> & RBEIMEvaluation::get_interpolation_matrix() const
949 : {
950 0 : return _interpolation_matrix;
951 : }
952 :
953 0 : void RBEIMEvaluation::set_preserve_rb_eim_solutions(bool preserve_rb_eim_solutions)
954 : {
955 0 : _preserve_rb_eim_solutions = preserve_rb_eim_solutions;
956 0 : }
957 :
958 0 : bool RBEIMEvaluation::get_preserve_rb_eim_solutions() const
959 : {
960 0 : return _preserve_rb_eim_solutions;
961 : }
962 :
963 0 : void RBEIMEvaluation::add_basis_function(
964 : const QpDataMap & bf)
965 : {
966 0 : _local_eim_basis_functions.emplace_back(bf);
967 0 : }
968 :
969 0 : void RBEIMEvaluation::add_interpolation_data(
970 : Point p,
971 : unsigned int comp,
972 : dof_id_type elem_id,
973 : subdomain_id_type subdomain_id,
974 : unsigned int qp,
975 : const std::vector<Point> & perturbs,
976 : const std::vector<Real> & phi_i_qp,
977 : ElemType elem_type,
978 : const std::vector<Real> & JxW_all_qp,
979 : const std::vector<std::vector<Real>> & phi_i_all_qp,
980 : Order qrule_order,
981 : const Point & dxyzdxi_elem_center,
982 : const Point & dxyzdeta_elem_center)
983 : {
984 0 : _vec_eval_input.all_xyz.emplace_back(p);
985 0 : _interpolation_points_comp.emplace_back(comp);
986 0 : _vec_eval_input.elem_ids.emplace_back(elem_id);
987 0 : _vec_eval_input.sbd_ids.emplace_back(subdomain_id);
988 0 : _vec_eval_input.qps.emplace_back(qp);
989 0 : _vec_eval_input.all_xyz_perturb.emplace_back(perturbs);
990 0 : _vec_eval_input.phi_i_qp.emplace_back(phi_i_qp);
991 0 : _vec_eval_input.elem_types.emplace_back(elem_type);
992 :
993 : // The following quantities are indexed by elem id to prevent duplicated data
994 : // If an entry is already present then we should not need to add that data again as it is element
995 : // based so it should be identical between 2 points if they refer to the same element id.
996 0 : if (_vec_eval_input.elem_id_to_local_index.count(elem_id) == 0)
997 : {
998 0 : unsigned int local_index = _vec_eval_input.JxW_all_qp.size();
999 : // Maybe check that getting local index from a a different source is still ok.
1000 0 : _vec_eval_input.elem_id_to_local_index[elem_id] = local_index;
1001 0 : _vec_eval_input.JxW_all_qp.emplace_back(JxW_all_qp);
1002 0 : _vec_eval_input.phi_i_all_qp.emplace_back(phi_i_all_qp);
1003 0 : _vec_eval_input.qrule_orders.emplace_back(qrule_order);
1004 0 : _vec_eval_input.dxyzdxi_elem_center.emplace_back(dxyzdxi_elem_center);
1005 0 : _vec_eval_input.dxyzdeta_elem_center.emplace_back(dxyzdeta_elem_center);
1006 : }
1007 :
1008 : // add_interpolation_data_to_property_map is a virtual function that aims to add only the property data
1009 : // required by the corresponding rb_parametrized_function.
1010 : // Right now, only elem_ids are supported but if subdomain_ids need to be stored for a property
1011 : // it can be added to the list of arguments and used in the overridden subclass function, the map type
1012 : // might also required an update to support various id types.
1013 0 : get_parametrized_function().add_interpolation_data_to_rb_property_map(this->comm(),
1014 0 : _vec_eval_input.rb_property_map,
1015 0 : elem_id);
1016 0 : }
1017 :
1018 0 : void RBEIMEvaluation::add_side_basis_function(
1019 : const SideQpDataMap & side_bf)
1020 : {
1021 0 : _local_side_eim_basis_functions.emplace_back(side_bf);
1022 0 : }
1023 :
1024 0 : void RBEIMEvaluation::add_side_interpolation_data(
1025 : Point p,
1026 : unsigned int comp,
1027 : dof_id_type elem_id,
1028 : unsigned int side_index,
1029 : subdomain_id_type subdomain_id,
1030 : boundary_id_type boundary_id,
1031 : unsigned int qp,
1032 : const std::vector<Point> & perturbs,
1033 : const std::vector<Real> & phi_i_qp)
1034 : {
1035 0 : _vec_eval_input.all_xyz.emplace_back(p);
1036 0 : _interpolation_points_comp.emplace_back(comp);
1037 0 : _vec_eval_input.elem_ids.emplace_back(elem_id);
1038 0 : _vec_eval_input.side_indices.emplace_back(side_index);
1039 0 : _vec_eval_input.sbd_ids.emplace_back(subdomain_id);
1040 0 : _vec_eval_input.boundary_ids.emplace_back(boundary_id);
1041 0 : _vec_eval_input.qps.emplace_back(qp);
1042 0 : _vec_eval_input.all_xyz_perturb.emplace_back(perturbs);
1043 0 : _vec_eval_input.phi_i_qp.emplace_back(phi_i_qp);
1044 :
1045 : // Add dummy values for the other properties, which are unused in the
1046 : // node case.
1047 0 : _vec_eval_input.elem_types.emplace_back(INVALID_ELEM);
1048 0 : _vec_eval_input.JxW_all_qp.emplace_back();
1049 0 : _vec_eval_input.phi_i_all_qp.emplace_back();
1050 0 : _vec_eval_input.qrule_orders.emplace_back(INVALID_ORDER);
1051 0 : _vec_eval_input.dxyzdxi_elem_center.emplace_back();
1052 0 : _vec_eval_input.dxyzdeta_elem_center.emplace_back();
1053 0 : }
1054 :
1055 0 : void RBEIMEvaluation::add_node_basis_function(
1056 : const NodeDataMap & node_bf)
1057 : {
1058 0 : _local_node_eim_basis_functions.emplace_back(node_bf);
1059 0 : }
1060 :
1061 0 : void RBEIMEvaluation::add_node_interpolation_data(
1062 : Point p,
1063 : unsigned int comp,
1064 : dof_id_type node_id,
1065 : boundary_id_type boundary_id)
1066 : {
1067 0 : _vec_eval_input.all_xyz.emplace_back(p);
1068 0 : _interpolation_points_comp.emplace_back(comp);
1069 0 : _vec_eval_input.node_ids.emplace_back(node_id);
1070 0 : _vec_eval_input.boundary_ids.emplace_back(boundary_id);
1071 :
1072 : // Add dummy values for the other properties, which are unused in the
1073 : // node case.
1074 0 : _vec_eval_input.elem_ids.emplace_back(0);
1075 0 : _vec_eval_input.side_indices.emplace_back(0);
1076 0 : _vec_eval_input.sbd_ids.emplace_back(0);
1077 0 : _vec_eval_input.qps.emplace_back(0);
1078 0 : _vec_eval_input.all_xyz_perturb.emplace_back();
1079 0 : _vec_eval_input.phi_i_qp.emplace_back();
1080 0 : _vec_eval_input.elem_types.emplace_back(INVALID_ELEM);
1081 0 : _vec_eval_input.JxW_all_qp.emplace_back();
1082 0 : _vec_eval_input.phi_i_all_qp.emplace_back();
1083 0 : _vec_eval_input.qrule_orders.emplace_back(INVALID_ORDER);
1084 0 : _vec_eval_input.dxyzdxi_elem_center.emplace_back();
1085 0 : _vec_eval_input.dxyzdeta_elem_center.emplace_back();
1086 0 : }
1087 :
1088 0 : void RBEIMEvaluation::
1089 : write_out_basis_functions(const std::string & directory_name,
1090 : bool write_binary_basis_functions)
1091 : {
1092 0 : LOG_SCOPE("write_out_basis_functions()", "RBEIMEvaluation");
1093 :
1094 0 : if (get_parametrized_function().on_mesh_sides())
1095 0 : write_out_side_basis_functions(directory_name, write_binary_basis_functions);
1096 0 : else if (get_parametrized_function().on_mesh_nodes())
1097 0 : write_out_node_basis_functions(directory_name, write_binary_basis_functions);
1098 : else
1099 0 : write_out_interior_basis_functions(directory_name, write_binary_basis_functions);
1100 0 : }
1101 :
1102 0 : void RBEIMEvaluation::
1103 : write_out_interior_basis_functions(const std::string & directory_name,
1104 : bool write_binary_basis_functions)
1105 : {
1106 0 : LOG_SCOPE("write_out_interior_basis_functions()", "RBEIMEvaluation");
1107 :
1108 : // Quick return if there is no work to do. Note: make sure all procs
1109 : // agree there is no work to do.
1110 0 : bool is_empty = _local_eim_basis_functions.empty();
1111 0 : libmesh_assert(this->comm().verify(is_empty));
1112 :
1113 0 : if (is_empty)
1114 0 : return;
1115 :
1116 : // Gather basis function data from other procs, storing it in
1117 : // _local_eim_basis_functions, so that we can then print everything
1118 : // from processor 0.
1119 0 : this->gather_bfs();
1120 :
1121 : // Write values from processor 0 only.
1122 0 : if (this->processor_id() == 0)
1123 : {
1124 0 : std::vector<unsigned int> n_qp_per_elem;
1125 : auto interior_basis_function_sizes =
1126 0 : get_interior_basis_function_sizes(n_qp_per_elem);
1127 :
1128 : // Make a directory to store all the data files
1129 0 : Utility::mkdir(directory_name.c_str());
1130 :
1131 : // Create filename
1132 0 : std::ostringstream file_name;
1133 0 : const std::string basis_function_suffix = (write_binary_basis_functions ? ".xdr" : ".dat");
1134 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1135 :
1136 : // Create XDR writer object
1137 0 : Xdr xdr(file_name.str(), write_binary_basis_functions ? ENCODE : WRITE);
1138 :
1139 : // Write number of basis functions to file. Note: the
1140 : // Xdr::data() function takes non-const references, so you can't
1141 : // pass e.g. vec.size() to that interface.
1142 0 : auto n_bf = libmesh_map_find(interior_basis_function_sizes,"n_bf");
1143 0 : xdr.data(n_bf, "# Number of basis functions");
1144 :
1145 : // We assume that each basis function has data for the same
1146 : // number of elements as basis function 0, which is equal to the
1147 : // size of the map.
1148 0 : auto n_elem = libmesh_map_find(interior_basis_function_sizes,"n_elem");
1149 0 : xdr.data(n_elem, "# Number of elements");
1150 :
1151 : // We assume that each element has the same number of variables,
1152 : // and we get the number of vars from the first element of the
1153 : // first basis function.
1154 0 : auto n_vars = libmesh_map_find(interior_basis_function_sizes,"n_vars");
1155 0 : xdr.data(n_vars, "# Number of variables");
1156 :
1157 : // We assume that the list of elements for each basis function
1158 : // is the same as basis function 0. We also assume that all vars
1159 : // have the same number of qps.
1160 0 : xdr.data(n_qp_per_elem, "# Number of QPs per Elem");
1161 :
1162 : // The total amount of qp data for each var is the sum of the
1163 : // entries in the "n_qp_per_elem" array.
1164 0 : auto n_qp_data = libmesh_map_find(interior_basis_function_sizes,"n_qp_data");
1165 :
1166 : // Now we construct a vector for each basis function, for each
1167 : // variable which is ordered according to:
1168 : // [ [qp vals for Elem 0], [qp vals for Elem 1], ... [qp vals for Elem N] ]
1169 : // and write it to file.
1170 0 : for (auto bf_index : index_range(_local_eim_basis_functions))
1171 : {
1172 0 : auto qp_data = get_interior_basis_function_as_vec_helper(n_vars, n_qp_data, bf_index);
1173 :
1174 : // Write all the var values for this bf
1175 0 : for (auto var : index_range(qp_data))
1176 0 : xdr.data_stream(qp_data[var].data(), qp_data[var].size(), /*line_break=*/qp_data[var].size());
1177 0 : }
1178 0 : }
1179 : }
1180 :
1181 0 : std::map<std::string,std::size_t> RBEIMEvaluation::
1182 : get_interior_basis_function_sizes(std::vector<unsigned int> & n_qp_per_elem)
1183 : {
1184 0 : std::map<std::string,std::size_t> interior_basis_function_sizes;
1185 :
1186 : // Write number of basis functions to file. Note: the
1187 : // Xdr::data() function takes non-const references, so you can't
1188 : // pass e.g. vec.size() to that interface.
1189 0 : auto n_bf = _local_eim_basis_functions.size();
1190 0 : interior_basis_function_sizes["n_bf"] = n_bf;
1191 :
1192 : // We assume that each basis function has data for the same
1193 : // number of elements as basis function 0, which is equal to the
1194 : // size of the map.
1195 0 : auto n_elem = _local_eim_basis_functions[0].size();
1196 0 : interior_basis_function_sizes["n_elem"] = n_elem;
1197 :
1198 : // We assume that each element has the same number of variables,
1199 : // and we get the number of vars from the first element of the
1200 : // first basis function.
1201 0 : auto n_vars = _local_eim_basis_functions[0].begin()->second.size();
1202 0 : interior_basis_function_sizes["n_vars"] = n_vars;
1203 :
1204 : // We assume that the list of elements for each basis function
1205 : // is the same as basis function 0. We also assume that all vars
1206 : // have the same number of qps.
1207 0 : n_qp_per_elem.clear();
1208 0 : n_qp_per_elem.reserve(n_elem);
1209 0 : dof_id_type expected_elem_id = 0;
1210 0 : for (const auto & [actual_elem_id, array] : _local_eim_basis_functions[0])
1211 : {
1212 : // Note: Currently we require that the Elems are numbered
1213 : // contiguously from [0..n_elem). This allows us to avoid
1214 : // writing the Elem ids to the Xdr file, but if we need to
1215 : // generalize this assumption later, we can.
1216 0 : libmesh_error_msg_if(actual_elem_id != expected_elem_id++,
1217 : "RBEIMEvaluation currently assumes a contiguous Elem numbering starting from 0.");
1218 :
1219 : // array[n_vars][n_qp] per Elem. We get the number of QPs
1220 : // for variable 0, assuming they are all the same.
1221 0 : n_qp_per_elem.push_back(array[0].size());
1222 : }
1223 :
1224 : // The total amount of qp data for each var is the sum of the
1225 : // entries in the "n_qp_per_elem" array.
1226 : auto n_qp_data =
1227 0 : std::accumulate(n_qp_per_elem.begin(),
1228 : n_qp_per_elem.end(),
1229 : 0u);
1230 0 : interior_basis_function_sizes["n_qp_data"] = n_qp_data;
1231 :
1232 0 : return interior_basis_function_sizes;
1233 : }
1234 :
1235 0 : std::vector<std::vector<Number>> RBEIMEvaluation::
1236 : get_interior_basis_function_as_vec_helper(
1237 : unsigned int n_vars,
1238 : unsigned int n_qp_data,
1239 : unsigned int bf_index)
1240 : {
1241 0 : LOG_SCOPE("get_interior_basis_function_as_vec_helper()", "RBEIMEvaluation");
1242 :
1243 0 : std::vector<std::vector<Number>> qp_data(n_vars);
1244 :
1245 : // Reserve enough capacity in qp_data in order to do the insertions below
1246 : // without further memory allocation.
1247 0 : for (auto var : index_range(qp_data))
1248 0 : qp_data[var].reserve(n_qp_data);
1249 :
1250 : // Now we construct a vector for each basis function, for each
1251 : // variable which is ordered according to:
1252 : // [ [qp vals for Elem 0], [qp vals for Elem 1], ... [qp vals for Elem N] ]
1253 : // and write it to file.
1254 :
1255 0 : libmesh_error_msg_if(bf_index >= _local_eim_basis_functions.size(), "bf_index not valid");
1256 0 : for (const auto & pr : _local_eim_basis_functions[bf_index])
1257 : {
1258 : // array[n_vars][n_qp] per Elem
1259 0 : const auto & array = pr.second;
1260 0 : for (auto var : index_range(array))
1261 : {
1262 : // Insert all qp values for this var
1263 0 : qp_data[var].insert(/*insert at*/qp_data[var].end(),
1264 0 : /*data start*/array[var].begin(),
1265 0 : /*data end*/array[var].end());
1266 : }
1267 : }
1268 :
1269 0 : return qp_data;
1270 0 : }
1271 :
1272 0 : std::vector<std::vector<std::vector<Number>>> RBEIMEvaluation::
1273 : get_interior_basis_functions_as_vecs()
1274 : {
1275 0 : LOG_SCOPE("get_interior_basis_function_as_vec()", "RBEIMEvaluation");
1276 :
1277 0 : std::vector<std::vector<std::vector<Number>>> interior_basis_functions;
1278 :
1279 : // Quick return if there is no work to do. Note: make sure all procs
1280 : // agree there is no work to do.
1281 0 : bool is_empty = _local_eim_basis_functions.empty();
1282 0 : libmesh_assert(this->comm().verify(is_empty));
1283 :
1284 0 : if (is_empty)
1285 0 : return interior_basis_functions;
1286 :
1287 : // Gather basis function data from other procs, storing it in
1288 : // _local_eim_basis_functions, so that we can then print everything
1289 : // from processor 0.
1290 0 : this->gather_bfs();
1291 :
1292 0 : if (this->processor_id() == 0)
1293 : {
1294 0 : std::vector<unsigned int> n_qp_per_elem;
1295 : std::map<std::string,std::size_t> interior_basis_function_sizes =
1296 0 : get_interior_basis_function_sizes(n_qp_per_elem);
1297 :
1298 0 : for (auto bf_index : index_range(_local_eim_basis_functions))
1299 0 : interior_basis_functions.emplace_back(
1300 0 : get_interior_basis_function_as_vec_helper(
1301 0 : libmesh_map_find(interior_basis_function_sizes,"n_vars"),
1302 0 : libmesh_map_find(interior_basis_function_sizes,"n_qp_data"),
1303 0 : bf_index));
1304 : }
1305 :
1306 0 : return interior_basis_functions;
1307 0 : }
1308 :
1309 0 : void RBEIMEvaluation::
1310 : write_out_side_basis_functions(const std::string & directory_name,
1311 : bool write_binary_basis_functions)
1312 : {
1313 0 : LOG_SCOPE("write_out_side_basis_functions()", "RBEIMEvaluation");
1314 :
1315 : // Quick return if there is no work to do. Note: make sure all procs
1316 : // agree there is no work to do.
1317 0 : bool is_empty = _local_side_eim_basis_functions.empty();
1318 0 : libmesh_assert(this->comm().verify(is_empty));
1319 :
1320 0 : if (is_empty)
1321 0 : return;
1322 :
1323 : // Gather basis function data from other procs, storing it in
1324 : // _local_side_eim_basis_functions, so that we can then print everything
1325 : // from processor 0.
1326 0 : this->side_gather_bfs();
1327 :
1328 : // Write values from processor 0 only.
1329 0 : if (this->processor_id() == 0)
1330 : {
1331 : // Make a directory to store all the data files
1332 0 : Utility::mkdir(directory_name.c_str());
1333 :
1334 : // Create filename
1335 0 : std::ostringstream file_name;
1336 0 : const std::string basis_function_suffix = (write_binary_basis_functions ? ".xdr" : ".dat");
1337 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1338 :
1339 : // Create XDR writer object
1340 0 : Xdr xdr(file_name.str(), write_binary_basis_functions ? ENCODE : WRITE);
1341 :
1342 : // Write number of basis functions to file. Note: the
1343 : // Xdr::data() function takes non-const references, so you can't
1344 : // pass e.g. vec.size() to that interface.
1345 0 : auto n_bf = _local_side_eim_basis_functions.size();
1346 0 : xdr.data(n_bf, "# Number of basis functions");
1347 :
1348 : // We assume that each basis function has data for the same
1349 : // number of (elem,side) pairs as basis function 0, which is equal to the
1350 : // size of the map.
1351 0 : auto n_elem = _local_side_eim_basis_functions[0].size();
1352 0 : xdr.data(n_elem, "# Number of (elem,side) pairs");
1353 :
1354 : // We assume that each element has the same number of variables,
1355 : // and we get the number of vars from the first element of the
1356 : // first basis function.
1357 0 : auto n_vars = _local_side_eim_basis_functions[0].begin()->second.size();
1358 0 : xdr.data(n_vars, "# Number of variables");
1359 :
1360 : // We write out the following arrays:
1361 : // - element IDs
1362 : // - side indices
1363 : // - n_qp_per_elem_side
1364 0 : std::vector<unsigned int> n_qp_per_elem_side;
1365 0 : std::vector<unsigned int> elem_ids;
1366 0 : std::vector<unsigned int> side_indices;
1367 0 : elem_ids.reserve(n_elem);
1368 0 : side_indices.reserve(n_elem);
1369 0 : n_qp_per_elem_side.reserve(n_elem);
1370 0 : for (const auto & [elem_side_pair, array] : _local_side_eim_basis_functions[0])
1371 : {
1372 0 : elem_ids.push_back(elem_side_pair.first);
1373 0 : side_indices.push_back(elem_side_pair.second);
1374 :
1375 : // array[n_vars][n_qp] per Elem. We get the number of QPs
1376 : // for variable 0, assuming they are all the same.
1377 0 : n_qp_per_elem_side.push_back(array[0].size());
1378 : }
1379 0 : xdr.data(elem_ids, "# Elem IDs");
1380 0 : xdr.data(side_indices, "# Side indices");
1381 0 : xdr.data(n_qp_per_elem_side, "# Number of QPs per Elem");
1382 :
1383 : // The total amount of qp data for each var is the sum of the
1384 : // entries in the "n_qp_per_elem" array.
1385 : auto n_qp_data =
1386 0 : std::accumulate(n_qp_per_elem_side.begin(),
1387 : n_qp_per_elem_side.end(),
1388 : 0u);
1389 :
1390 : // Reserve space to store contiguous vectors of qp data for each var
1391 0 : std::vector<std::vector<Number>> qp_data(n_vars);
1392 0 : for (auto var : index_range(qp_data))
1393 0 : qp_data[var].reserve(n_qp_data);
1394 :
1395 : // Now we construct a vector for each basis function, for each
1396 : // variable which is ordered according to:
1397 : // [ [qp vals for Elem 0], [qp vals for Elem 1], ... [qp vals for Elem N] ]
1398 : // and write it to file.
1399 0 : for (auto bf : index_range(_local_side_eim_basis_functions))
1400 : {
1401 : // Clear any data from previous bf
1402 0 : for (auto var : index_range(qp_data))
1403 0 : qp_data[var].clear();
1404 :
1405 0 : for (const auto & pr : _local_side_eim_basis_functions[bf])
1406 : {
1407 : // array[n_vars][n_qp] per Elem
1408 0 : const auto & array = pr.second;
1409 0 : for (auto var : index_range(array))
1410 : {
1411 : // Insert all qp values for this var
1412 0 : qp_data[var].insert(/*insert at*/qp_data[var].end(),
1413 0 : /*data start*/array[var].begin(),
1414 0 : /*data end*/array[var].end());
1415 : }
1416 : }
1417 :
1418 : // Write all the var values for this bf
1419 0 : for (auto var : index_range(qp_data))
1420 0 : xdr.data_stream(qp_data[var].data(), qp_data[var].size(), /*line_break=*/qp_data[var].size());
1421 : }
1422 0 : }
1423 : }
1424 :
1425 0 : void RBEIMEvaluation::
1426 : write_out_node_basis_functions(const std::string & directory_name,
1427 : bool write_binary_basis_functions)
1428 : {
1429 0 : LOG_SCOPE("write_out_node_basis_functions()", "RBEIMEvaluation");
1430 :
1431 : // Quick return if there is no work to do. Note: make sure all procs
1432 : // agree there is no work to do.
1433 0 : bool is_empty = _local_node_eim_basis_functions.empty();
1434 0 : libmesh_assert(this->comm().verify(is_empty));
1435 :
1436 0 : if (is_empty)
1437 0 : return;
1438 :
1439 : // Gather basis function data from other procs, storing it in
1440 : // _local_node_eim_basis_functions, so that we can then print everything
1441 : // from processor 0.
1442 0 : this->node_gather_bfs();
1443 :
1444 : // Write values from processor 0 only.
1445 0 : if (this->processor_id() == 0)
1446 : {
1447 : // Make a directory to store all the data files
1448 0 : Utility::mkdir(directory_name.c_str());
1449 :
1450 : // Create filename
1451 0 : std::ostringstream file_name;
1452 0 : const std::string basis_function_suffix = (write_binary_basis_functions ? ".xdr" : ".dat");
1453 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1454 :
1455 : // Create XDR writer object
1456 0 : Xdr xdr(file_name.str(), write_binary_basis_functions ? ENCODE : WRITE);
1457 :
1458 : // Write number of basis functions to file. Note: the
1459 : // Xdr::data() function takes non-const references, so you can't
1460 : // pass e.g. vec.size() to that interface.
1461 0 : auto n_bf = _local_node_eim_basis_functions.size();
1462 0 : xdr.data(n_bf, "# Number of basis functions");
1463 :
1464 : // We assume that each basis function has data for the same
1465 : // number of elements as basis function 0, which is equal to the
1466 : // size of the map.
1467 0 : auto n_node = _local_node_eim_basis_functions[0].size();
1468 0 : xdr.data(n_node, "# Number of nodes");
1469 :
1470 : // We assume that each element has the same number of variables,
1471 : // and we get the number of vars from the first element of the
1472 : // first basis function.
1473 0 : auto n_vars = _local_node_eim_basis_functions[0].begin()->second.size();
1474 0 : xdr.data(n_vars, "# Number of variables");
1475 :
1476 : // We write out the following arrays:
1477 : // - node IDs
1478 0 : std::vector<unsigned int> node_ids;
1479 0 : node_ids.reserve(n_node);
1480 0 : for (const auto & pr : _local_node_eim_basis_functions[0])
1481 : {
1482 0 : node_ids.push_back(pr.first);
1483 : }
1484 0 : xdr.data(node_ids, "# Node IDs");
1485 :
1486 : // Now we construct a vector for each basis function, for each
1487 : // variable which is ordered according to:
1488 : // [ [val for Node 0], [val for Node 1], ... [val for Node N] ]
1489 : // and write it to file.
1490 :
1491 0 : std::vector<std::vector<Number>> var_data(n_vars);
1492 0 : for (unsigned int var=0; var<n_vars; var++)
1493 0 : var_data[var].resize(n_node);
1494 :
1495 0 : for (auto bf : index_range(_local_node_eim_basis_functions))
1496 : {
1497 0 : unsigned int node_counter = 0;
1498 0 : for (const auto & pr : _local_node_eim_basis_functions[bf])
1499 : {
1500 : // array[n_vars] per Node
1501 0 : const auto & array = pr.second;
1502 0 : for (auto var : index_range(array))
1503 : {
1504 : // Based on the error check above, we know that node_id is numbered
1505 : // contiguously from [0..nodes], so we can use it as the vector
1506 : // index here.
1507 0 : var_data[var][node_counter] = array[var];
1508 : }
1509 :
1510 0 : node_counter++;
1511 : }
1512 :
1513 : // Write all the var values for this bf
1514 0 : for (auto var : index_range(var_data))
1515 0 : xdr.data_stream(var_data[var].data(), var_data[var].size(), /*line_break=*/var_data[var].size());
1516 : }
1517 0 : }
1518 : }
1519 :
1520 0 : void RBEIMEvaluation::
1521 : read_in_basis_functions(const System & sys,
1522 : const std::string & directory_name,
1523 : bool read_binary_basis_functions)
1524 : {
1525 0 : LOG_SCOPE("read_in_basis_functions()", "RBEIMEvaluation");
1526 :
1527 : // Return early without reading in anything if there are no basis functions
1528 0 : if (get_n_basis_functions() == 0)
1529 0 : return;
1530 :
1531 0 : if (get_parametrized_function().on_mesh_sides())
1532 0 : read_in_side_basis_functions(sys, directory_name, read_binary_basis_functions);
1533 0 : else if (get_parametrized_function().on_mesh_nodes())
1534 0 : read_in_node_basis_functions(sys, directory_name, read_binary_basis_functions);
1535 : else
1536 0 : read_in_interior_basis_functions(sys, directory_name, read_binary_basis_functions);
1537 : }
1538 :
1539 0 : void RBEIMEvaluation::
1540 : read_in_interior_basis_functions(const System & sys,
1541 : const std::string & directory_name,
1542 : bool read_binary_basis_functions)
1543 : {
1544 0 : LOG_SCOPE("read_in_interior_basis_functions()", "RBEIMEvaluation");
1545 :
1546 : // Read values on processor 0 only.
1547 0 : if (sys.comm().rank() == 0)
1548 : {
1549 : // Create filename
1550 0 : std::ostringstream file_name;
1551 0 : const std::string basis_function_suffix = (read_binary_basis_functions ? ".xdr" : ".dat");
1552 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1553 :
1554 : // Create XDR reader object
1555 0 : Xdr xdr(file_name.str(), read_binary_basis_functions ? DECODE : READ);
1556 :
1557 : // Read in the number of basis functions. The comment parameter
1558 : // is ignored when reading.
1559 : std::size_t n_bf;
1560 0 : xdr.data(n_bf);
1561 :
1562 : // Read in the number of elements
1563 : std::size_t n_elem;
1564 0 : xdr.data(n_elem);
1565 :
1566 : // Read in the number of variables.
1567 : std::size_t n_vars;
1568 0 : xdr.data(n_vars);
1569 :
1570 : // Read in vector containing the number of QPs per elem. We can
1571 : // create this vector with the required size or let it be read
1572 : // from the file and sized for us.
1573 0 : std::vector<unsigned int> n_qp_per_elem(n_elem);
1574 0 : xdr.data(n_qp_per_elem);
1575 :
1576 : // The total amount of qp data for each var is the sum of the
1577 : // entries in the "n_qp_per_elem" array.
1578 : auto n_qp_data =
1579 0 : std::accumulate(n_qp_per_elem.begin(),
1580 : n_qp_per_elem.end(),
1581 : 0u);
1582 :
1583 : // Allocate space to store all required basis functions,
1584 : // clearing any data that may have been there previously.
1585 : //
1586 : // TODO: Do we need to also write out/read in Elem ids?
1587 : // Or can we assume they will always be contiguously
1588 : // numbered (at least on proc 0)?
1589 0 : _local_eim_basis_functions.clear();
1590 0 : _local_eim_basis_functions.resize(n_bf);
1591 0 : for (auto i : index_range(_local_eim_basis_functions))
1592 0 : for (std::size_t elem_id=0; elem_id<n_elem; ++elem_id)
1593 : {
1594 0 : auto & array = _local_eim_basis_functions[i][elem_id];
1595 0 : array.resize(n_vars);
1596 : }
1597 :
1598 : // Allocate temporary storage for one var's worth of qp data.
1599 0 : std::vector<Number> qp_data;
1600 :
1601 : // Read in data for each basis function
1602 0 : for (auto i : index_range(_local_eim_basis_functions))
1603 : {
1604 : // Reference to the data map for the current basis function.
1605 0 : auto & bf_map = _local_eim_basis_functions[i];
1606 :
1607 0 : for (std::size_t var=0; var<n_vars; ++var)
1608 : {
1609 0 : qp_data.clear();
1610 0 : qp_data.resize(n_qp_data);
1611 :
1612 : // Read data using data_stream() since that is
1613 : // (currently) how we write it out. The "line_break"
1614 : // parameter of data_stream() is ignored while reading.
1615 0 : xdr.data_stream(qp_data.data(), qp_data.size());
1616 :
1617 : // Iterate over the qp_data vector, filling in the
1618 : // "small" vectors for each Elem.
1619 0 : auto cursor = qp_data.begin();
1620 0 : for (std::size_t elem_id=0; elem_id<n_elem; ++elem_id)
1621 : {
1622 : // Get reference to the [n_vars][n_qp] array for
1623 : // this Elem. We assign() into the vector of
1624 : // quadrature point values, which allocates space if
1625 : // it doesn't already exist.
1626 0 : auto & array = bf_map[elem_id];
1627 0 : array[var].assign(cursor, cursor + n_qp_per_elem[elem_id]);
1628 0 : std::advance(cursor, n_qp_per_elem[elem_id]);
1629 : }
1630 : } // end for (var)
1631 : } // end for (i)
1632 0 : } // end if processor 0
1633 :
1634 : // Distribute the basis function information to the processors that require it
1635 0 : this->distribute_bfs(sys);
1636 0 : }
1637 :
1638 0 : void RBEIMEvaluation::
1639 : read_in_side_basis_functions(const System & sys,
1640 : const std::string & directory_name,
1641 : bool read_binary_basis_functions)
1642 : {
1643 0 : LOG_SCOPE("read_in_basis_functions()", "RBEIMEvaluation");
1644 :
1645 : // Read values on processor 0 only.
1646 0 : if (sys.comm().rank() == 0)
1647 : {
1648 : // Create filename
1649 0 : std::ostringstream file_name;
1650 0 : const std::string basis_function_suffix = (read_binary_basis_functions ? ".xdr" : ".dat");
1651 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1652 :
1653 : // Create XDR reader object
1654 0 : Xdr xdr(file_name.str(), read_binary_basis_functions ? DECODE : READ);
1655 :
1656 : // Read in the number of basis functions. The comment parameter
1657 : // is ignored when reading.
1658 : std::size_t n_bf;
1659 0 : xdr.data(n_bf);
1660 :
1661 : // Read in the number of elements
1662 : std::size_t n_elem_side;
1663 0 : xdr.data(n_elem_side);
1664 :
1665 : // Read in the number of variables.
1666 : std::size_t n_vars;
1667 0 : xdr.data(n_vars);
1668 :
1669 0 : std::vector<unsigned int> elem_ids(n_elem_side);
1670 0 : xdr.data(elem_ids);
1671 0 : std::vector<unsigned int> side_indices(n_elem_side);
1672 0 : xdr.data(side_indices);
1673 :
1674 : // Read in vector containing the number of QPs per elem. We can
1675 : // create this vector with the required size or let it be read
1676 : // from the file and sized for us.
1677 0 : std::vector<unsigned int> n_qp_per_elem_side(n_elem_side);
1678 0 : xdr.data(n_qp_per_elem_side);
1679 :
1680 : // The total amount of qp data for each var is the sum of the
1681 : // entries in the "n_qp_per_elem" array.
1682 : auto n_qp_data =
1683 0 : std::accumulate(n_qp_per_elem_side.begin(),
1684 : n_qp_per_elem_side.end(),
1685 : 0u);
1686 :
1687 : // Allocate space to store all required basis functions,
1688 : // clearing any data that may have been there previously.
1689 0 : _local_side_eim_basis_functions.clear();
1690 0 : _local_side_eim_basis_functions.resize(n_bf);
1691 0 : for (auto i : index_range(_local_side_eim_basis_functions))
1692 0 : for (std::size_t elem_side_idx=0; elem_side_idx<n_elem_side; ++elem_side_idx)
1693 : {
1694 0 : unsigned int elem_id = elem_ids[elem_side_idx];
1695 0 : unsigned int side_index = side_indices[elem_side_idx];
1696 0 : auto elem_side_pair = std::make_pair(elem_id, side_index);
1697 :
1698 0 : auto & array = _local_side_eim_basis_functions[i][elem_side_pair];
1699 0 : array.resize(n_vars);
1700 : }
1701 :
1702 : // Allocate temporary storage for one var's worth of qp data.
1703 0 : std::vector<Number> qp_data;
1704 :
1705 : // Read in data for each basis function
1706 0 : for (auto i : index_range(_local_side_eim_basis_functions))
1707 : {
1708 : // Reference to the data map for the current basis function.
1709 0 : auto & bf_map = _local_side_eim_basis_functions[i];
1710 :
1711 0 : for (std::size_t var=0; var<n_vars; ++var)
1712 : {
1713 0 : qp_data.clear();
1714 0 : qp_data.resize(n_qp_data);
1715 :
1716 : // Read data using data_stream() since that is
1717 : // (currently) how we write it out. The "line_break"
1718 : // parameter of data_stream() is ignored while reading.
1719 0 : xdr.data_stream(qp_data.data(), qp_data.size());
1720 :
1721 : // Iterate over the qp_data vector, filling in the
1722 : // "small" vectors for each Elem.
1723 0 : auto cursor = qp_data.begin();
1724 0 : for (std::size_t elem_side_idx=0; elem_side_idx<n_elem_side; ++elem_side_idx)
1725 : {
1726 0 : unsigned int elem_id = elem_ids[elem_side_idx];
1727 0 : unsigned int side_index = side_indices[elem_side_idx];
1728 0 : auto elem_side_pair = std::make_pair(elem_id, side_index);
1729 :
1730 : // Get reference to the [n_vars][n_qp] array for
1731 : // this Elem. We assign() into the vector of
1732 : // quadrature point values, which allocates space if
1733 : // it doesn't already exist.
1734 0 : auto & array = bf_map[elem_side_pair];
1735 0 : array[var].assign(cursor, cursor + n_qp_per_elem_side[elem_side_idx]);
1736 0 : std::advance(cursor, n_qp_per_elem_side[elem_side_idx]);
1737 : }
1738 : } // end for (var)
1739 : } // end for (i)
1740 0 : } // end if processor 0
1741 :
1742 : // Distribute the basis function information to the processors that require it
1743 0 : this->side_distribute_bfs(sys);
1744 0 : }
1745 :
1746 0 : void RBEIMEvaluation::
1747 : read_in_node_basis_functions(const System & sys,
1748 : const std::string & directory_name,
1749 : bool read_binary_basis_functions)
1750 : {
1751 0 : LOG_SCOPE("read_in_node_basis_functions()", "RBEIMEvaluation");
1752 :
1753 : // Read values on processor 0 only.
1754 0 : if (sys.comm().rank() == 0)
1755 : {
1756 : // Create filename
1757 0 : std::ostringstream file_name;
1758 0 : const std::string basis_function_suffix = (read_binary_basis_functions ? ".xdr" : ".dat");
1759 0 : file_name << directory_name << "/" << "bf_data" << basis_function_suffix;
1760 :
1761 : // Create XDR reader object
1762 0 : Xdr xdr(file_name.str(), read_binary_basis_functions ? DECODE : READ);
1763 :
1764 : // Read in the number of basis functions. The comment parameter
1765 : // is ignored when reading.
1766 : std::size_t n_bf;
1767 0 : xdr.data(n_bf);
1768 :
1769 : // Read in the number of nodes
1770 : std::size_t n_node;
1771 0 : xdr.data(n_node);
1772 :
1773 : // Read in the number of variables.
1774 : std::size_t n_vars;
1775 0 : xdr.data(n_vars);
1776 :
1777 0 : std::vector<unsigned int> node_ids(n_node);
1778 0 : xdr.data(node_ids);
1779 :
1780 : // Allocate space to store all required basis functions,
1781 : // clearing any data that may have been there previously.
1782 : //
1783 : // TODO: Do we need to also write out/read in Node ids?
1784 : // Or can we assume they will always be contiguously
1785 : // numbered (at least on proc 0)?
1786 0 : _local_node_eim_basis_functions.clear();
1787 0 : _local_node_eim_basis_functions.resize(n_bf);
1788 0 : for (auto i : index_range(_local_node_eim_basis_functions))
1789 0 : for (auto node_id : node_ids)
1790 : {
1791 0 : auto & array = _local_node_eim_basis_functions[i][node_id];
1792 0 : array.resize(n_vars);
1793 : }
1794 :
1795 : // Read data into node_value from xdr
1796 0 : std::vector<Number> node_value;
1797 :
1798 : // Read in data for each basis function
1799 0 : for (auto i : index_range(_local_node_eim_basis_functions))
1800 : {
1801 : // Reference to the data map for the current basis function.
1802 0 : auto & bf_map = _local_node_eim_basis_functions[i];
1803 :
1804 0 : for (std::size_t var=0; var<n_vars; ++var)
1805 : {
1806 0 : node_value.clear();
1807 0 : node_value.resize(n_node);
1808 :
1809 : // Read data using data_stream() since that is
1810 : // (currently) how we write it out. The "line_break"
1811 : // parameter of data_stream() is ignored while reading.
1812 0 : xdr.data_stream(node_value.data(), node_value.size());
1813 :
1814 0 : for (unsigned int node_counter=0; node_counter<n_node; node_counter++)
1815 : {
1816 0 : auto & array = bf_map[node_ids[node_counter]];
1817 0 : array[var] = node_value[node_counter];
1818 : }
1819 : } // end for (var)
1820 : } // end for (i)
1821 0 : } // end if processor 0
1822 :
1823 : // Distribute the basis function information to the processors that require it
1824 0 : this->node_distribute_bfs(sys);
1825 0 : }
1826 :
1827 0 : void RBEIMEvaluation::print_local_eim_basis_functions() const
1828 : {
1829 0 : for (auto bf : index_range(_local_eim_basis_functions))
1830 : {
1831 0 : libMesh::out << "Interior basis function " << bf << std::endl;
1832 0 : for (const auto & [elem_id, array] : _local_eim_basis_functions[bf])
1833 : {
1834 0 : libMesh::out << "Elem " << elem_id << std::endl;
1835 0 : for (auto var : index_range(array))
1836 : {
1837 0 : libMesh::out << "Variable " << var << std::endl;
1838 0 : for (auto qp : index_range(array[var]))
1839 0 : libMesh::out << array[var][qp] << " ";
1840 0 : libMesh::out << std::endl;
1841 : }
1842 : }
1843 : }
1844 :
1845 0 : for (auto bf : index_range(_local_side_eim_basis_functions))
1846 : {
1847 0 : libMesh::out << "Side basis function " << bf << std::endl;
1848 0 : for (const auto & [pr, array] : _local_side_eim_basis_functions[bf])
1849 : {
1850 0 : const auto & elem_id = pr.first;
1851 0 : const auto & side_index = pr.second;
1852 0 : libMesh::out << "Elem " << elem_id << ", Side " << side_index << std::endl;
1853 0 : for (auto var : index_range(array))
1854 : {
1855 0 : libMesh::out << "Variable " << var << std::endl;
1856 0 : for (auto qp : index_range(array[var]))
1857 0 : libMesh::out << array[var][qp] << " ";
1858 0 : libMesh::out << std::endl;
1859 : }
1860 : }
1861 : }
1862 :
1863 0 : for (auto bf : index_range(_local_node_eim_basis_functions))
1864 : {
1865 0 : libMesh::out << "Node basis function " << bf << std::endl;
1866 0 : for (const auto & [node_id, array] : _local_node_eim_basis_functions[bf])
1867 : {
1868 0 : libMesh::out << "Node " << node_id << std::endl;
1869 0 : for (auto var : index_range(array))
1870 : {
1871 0 : libMesh::out << "Variable " << var << ": " << array[var] << std::endl;
1872 : }
1873 0 : libMesh::out << std::endl;
1874 : }
1875 : }
1876 0 : }
1877 :
1878 0 : void RBEIMEvaluation::gather_bfs()
1879 : {
1880 : // We need to gather _local_eim_basis_functions data from other
1881 : // procs for printing.
1882 : //
1883 : // Ideally, this could be accomplished by simply calling:
1884 : // this->comm().gather(/*root_id=*/0, _local_eim_basis_functions);
1885 : //
1886 : // but the data structure seems to be too complicated for this to
1887 : // work automatically. (I get some error about the function called
1888 : // being "private within this context".) Therefore, we have to
1889 : // gather the information manually.
1890 :
1891 : // So we can avoid calling this many times below
1892 0 : auto n_procs = this->n_processors();
1893 :
1894 : // In serial there's nothing to gather
1895 0 : if (n_procs == 1)
1896 0 : return;
1897 :
1898 : // Current assumption is that the number of basis functions stored on
1899 : // each processor is the same, the only thing that differs is the number
1900 : // of elements, so make sure that is the case now.
1901 0 : auto n_bf = _local_eim_basis_functions.size();
1902 0 : libmesh_assert(this->comm().verify(n_bf));
1903 :
1904 : // This function should never be called if there are no basis
1905 : // functions, so if it was, something went wrong.
1906 0 : libmesh_error_msg_if(!n_bf, "RBEIMEvaluation::gather_bfs() should not be called with 0 basis functions.");
1907 :
1908 : // The number of variables should be the same on all processors
1909 : // and we can get this from _local_eim_basis_functions. However,
1910 : // it may be that some processors have no local elements, so on
1911 : // those processors we cannot look up the size from
1912 : // _local_eim_basis_functions. As a result we use comm().max(n_vars)
1913 : // to make sure all processors agree on the final value.
1914 : std::size_t n_vars =
1915 0 : _local_eim_basis_functions[0].empty() ? 0 : _local_eim_basis_functions[0].begin()->second.size();
1916 0 : this->comm().max(n_vars);
1917 :
1918 : // Gather list of Elem ids stored on each processor to proc 0. We
1919 : // use basis function 0 as an example and assume all the basis
1920 : // functions are distributed similarly.
1921 0 : std::vector<dof_id_type> elem_ids;
1922 0 : elem_ids.reserve(_local_eim_basis_functions[0].size());
1923 0 : for (const auto & pr : _local_eim_basis_functions[0])
1924 0 : elem_ids.push_back(pr.first);
1925 0 : this->comm().gather(/*root_id=*/0, elem_ids);
1926 :
1927 : // Store the number of qps per Elem on this processor. Again, use
1928 : // basis function 0 (and variable 0) to get this information, then
1929 : // apply it to all basis functions.
1930 0 : std::vector<unsigned int> n_qp_per_elem;
1931 0 : n_qp_per_elem.reserve(_local_eim_basis_functions[0].size());
1932 0 : for (const auto & pr : _local_eim_basis_functions[0])
1933 : {
1934 : // array[n_vars][n_qp] per Elem. We get the number of QPs
1935 : // for variable 0, assuming they are all the same.
1936 0 : const auto & array = pr.second;
1937 0 : n_qp_per_elem.push_back(array[0].size());
1938 : }
1939 :
1940 : // Before gathering, compute the total amount of local qp data for
1941 : // each var, which is the sum of the entries in the "n_qp_per_elem" array.
1942 : // This will be used to reserve space in a vector below.
1943 : auto n_local_qp_data =
1944 0 : std::accumulate(n_qp_per_elem.begin(),
1945 : n_qp_per_elem.end(),
1946 : 0u);
1947 :
1948 : // Gather the number of qps per Elem for each processor onto processor 0.
1949 0 : this->comm().gather(/*root_id=*/0, n_qp_per_elem);
1950 :
1951 : // Sanity check: On processor 0, this checks that we have gathered the same number
1952 : // of elem ids and qp counts.
1953 0 : libmesh_error_msg_if(elem_ids.size() != n_qp_per_elem.size(),
1954 : "Must gather same number of Elem ids as qps per Elem.");
1955 :
1956 : // Reserve space to store contiguous vectors of qp data for each var
1957 0 : std::vector<std::vector<Number>> gathered_qp_data(n_vars);
1958 0 : for (auto var : index_range(gathered_qp_data))
1959 0 : gathered_qp_data[var].reserve(n_local_qp_data);
1960 :
1961 : // Now we construct a vector for each basis function, for each
1962 : // variable, which is ordered according to:
1963 : // [ [qp vals for Elem 0], [qp vals for Elem 1], ... [qp vals for Elem N] ]
1964 : // and gather it to processor 0.
1965 0 : for (auto bf : index_range(_local_eim_basis_functions))
1966 : {
1967 : // Clear any data from previous bf
1968 0 : for (auto var : index_range(gathered_qp_data))
1969 0 : gathered_qp_data[var].clear();
1970 :
1971 0 : for (const auto & pr : _local_eim_basis_functions[bf])
1972 : {
1973 : // array[n_vars][n_qp] per Elem
1974 0 : const auto & array = pr.second;
1975 0 : for (auto var : index_range(array))
1976 : {
1977 : // Insert all qp values for this var
1978 0 : gathered_qp_data[var].insert(/*insert at*/gathered_qp_data[var].end(),
1979 0 : /*data start*/array[var].begin(),
1980 0 : /*data end*/array[var].end());
1981 : }
1982 : }
1983 :
1984 : // Reference to the data map for the current basis function.
1985 0 : auto & bf_map = _local_eim_basis_functions[bf];
1986 :
1987 0 : for (auto var : index_range(gathered_qp_data))
1988 : {
1989 : // For each var, gather gathered_qp_data[var] onto processor
1990 : // 0. There apparently is not a gather overload for
1991 : // vector-of-vectors...
1992 0 : this->comm().gather(/*root_id=*/0, gathered_qp_data[var]);
1993 :
1994 : // On processor 0, iterate over the gathered_qp_data[var]
1995 : // vector we just gathered, filling in the "small" vectors
1996 : // for each Elem. Note: here we ignore the fact that we
1997 : // already have the data on processor 0 and just overwrite
1998 : // it, this makes the indexing logic a bit simpler.
1999 0 : if (this->processor_id() == 0)
2000 : {
2001 0 : auto cursor = gathered_qp_data[var].begin();
2002 0 : for (auto i : index_range(elem_ids))
2003 : {
2004 0 : auto elem_id = elem_ids[i];
2005 0 : auto n_qp_this_elem = n_qp_per_elem[i];
2006 :
2007 : // Get reference to the [n_vars][n_qp] array for
2008 : // this Elem. We assign() into the vector of
2009 : // quadrature point values, which allocates space if
2010 : // it doesn't already exist.
2011 0 : auto & array = bf_map[elem_id];
2012 :
2013 : // Possibly allocate space if this is data for a new
2014 : // element we haven't seen before.
2015 0 : if (array.empty())
2016 0 : array.resize(n_vars);
2017 :
2018 0 : array[var].assign(cursor, cursor + n_qp_this_elem);
2019 0 : std::advance(cursor, n_qp_this_elem);
2020 : }
2021 : }
2022 : }
2023 : } // end loop over basis functions
2024 0 : }
2025 :
2026 0 : void RBEIMEvaluation::side_gather_bfs()
2027 : {
2028 : // We need to gather _local_side_eim_basis_functions data from other
2029 : // procs for printing.
2030 : //
2031 : // Ideally, this could be accomplished by simply calling:
2032 : // this->comm().gather(/*root_id=*/0, _local_side_eim_basis_functions);
2033 : //
2034 : // but the data structure seems to be too complicated for this to
2035 : // work automatically. (I get some error about the function called
2036 : // being "private within this context".) Therefore, we have to
2037 : // gather the information manually.
2038 :
2039 : // So we can avoid calling this many times below
2040 0 : auto n_procs = this->n_processors();
2041 :
2042 : // In serial there's nothing to gather
2043 0 : if (n_procs == 1)
2044 0 : return;
2045 :
2046 : // Current assumption is that the number of basis functions stored on
2047 : // each processor is the same, the only thing that differs is the number
2048 : // of elements, so make sure that is the case now.
2049 0 : auto n_bf = _local_side_eim_basis_functions.size();
2050 0 : libmesh_assert(this->comm().verify(n_bf));
2051 :
2052 : // This function should never be called if there are no basis
2053 : // functions, so if it was, something went wrong.
2054 0 : libmesh_error_msg_if(!n_bf, "SideRBEIMEvaluation::gather_bfs() should not be called with 0 basis functions.");
2055 :
2056 : // The number of variables should be the same on all processors
2057 : // and we can get this from _local_side_eim_basis_functions. However,
2058 : // it may be that some processors have no local elements, so on
2059 : // those processors we cannot look up the size from
2060 : // _local_side_eim_basis_functions. As a result we use comm().max(n_vars)
2061 : // to make sure all processors agree on the final value.
2062 : std::size_t n_vars =
2063 0 : _local_side_eim_basis_functions[0].empty() ? 0 : _local_side_eim_basis_functions[0].begin()->second.size();
2064 0 : this->comm().max(n_vars);
2065 :
2066 : // Gather list of (elem,side) pairs stored on each processor to proc 0. We
2067 : // use basis function 0 as an example and assume all the basis
2068 : // functions are distributed similarly.
2069 0 : std::vector<std::pair<dof_id_type,unsigned int>> elem_side_pairs;
2070 0 : elem_side_pairs.reserve(_local_side_eim_basis_functions[0].size());
2071 0 : for (const auto & pr : _local_side_eim_basis_functions[0])
2072 0 : elem_side_pairs.push_back(pr.first);
2073 0 : this->comm().gather(/*root_id=*/0, elem_side_pairs);
2074 :
2075 : // Store the number of qps per Elem on this processor. Again, use
2076 : // basis function 0 (and variable 0) to get this information, then
2077 : // apply it to all basis functions.
2078 0 : std::vector<unsigned int> n_qp_per_elem_side;
2079 0 : n_qp_per_elem_side.reserve(_local_side_eim_basis_functions[0].size());
2080 0 : for (const auto & pr : _local_side_eim_basis_functions[0])
2081 : {
2082 : // array[n_vars][n_qp] per (elem,side). We get the number of QPs
2083 : // for variable 0, assuming they are all the same.
2084 0 : const auto & array = pr.second;
2085 0 : n_qp_per_elem_side.push_back(array[0].size());
2086 : }
2087 :
2088 : // Before gathering, compute the total amount of local qp data for
2089 : // each var, which is the sum of the entries in the "n_qp_per_elem_side" array.
2090 : // This will be used to reserve space in a vector below.
2091 : auto n_local_qp_data =
2092 0 : std::accumulate(n_qp_per_elem_side.begin(),
2093 : n_qp_per_elem_side.end(),
2094 : 0u);
2095 :
2096 : // Gather the number of qps per Elem for each processor onto processor 0.
2097 0 : this->comm().gather(/*root_id=*/0, n_qp_per_elem_side);
2098 :
2099 : // Sanity check: On processor 0, this checks that we have gathered the same number
2100 : // of (elem,side) pairs and qp counts.
2101 0 : libmesh_error_msg_if(elem_side_pairs.size() != n_qp_per_elem_side.size(),
2102 : "Must gather same number of Elem ids as qps per Elem.");
2103 :
2104 : // Reserve space to store contiguous vectors of qp data for each var
2105 0 : std::vector<std::vector<Number>> gathered_qp_data(n_vars);
2106 0 : for (auto var : index_range(gathered_qp_data))
2107 0 : gathered_qp_data[var].reserve(n_local_qp_data);
2108 :
2109 : // Now we construct a vector for each basis function, for each
2110 : // variable, which is ordered according to:
2111 : // [ [qp vals for Elem 0], [qp vals for Elem 1], ... [qp vals for Elem N] ]
2112 : // and gather it to processor 0.
2113 0 : for (auto bf : index_range(_local_side_eim_basis_functions))
2114 : {
2115 : // Clear any data from previous bf
2116 0 : for (auto var : index_range(gathered_qp_data))
2117 0 : gathered_qp_data[var].clear();
2118 :
2119 0 : for (const auto & pr : _local_side_eim_basis_functions[bf])
2120 : {
2121 : // array[n_vars][n_qp] per (elem,side) pair
2122 0 : const auto & array = pr.second;
2123 0 : for (auto var : index_range(array))
2124 : {
2125 : // Insert all qp values for this var
2126 0 : gathered_qp_data[var].insert(/*insert at*/gathered_qp_data[var].end(),
2127 0 : /*data start*/array[var].begin(),
2128 0 : /*data end*/array[var].end());
2129 : }
2130 : }
2131 :
2132 : // Reference to the data map for the current basis function.
2133 0 : auto & bf_map = _local_side_eim_basis_functions[bf];
2134 :
2135 0 : for (auto var : index_range(gathered_qp_data))
2136 : {
2137 : // For each var, gather gathered_qp_data[var] onto processor
2138 : // 0. There apparently is not a gather overload for
2139 : // vector-of-vectors...
2140 0 : this->comm().gather(/*root_id=*/0, gathered_qp_data[var]);
2141 :
2142 : // On processor 0, iterate over the gathered_qp_data[var]
2143 : // vector we just gathered, filling in the "small" vectors
2144 : // for each Elem. Note: here we ignore the fact that we
2145 : // already have the data on processor 0 and just overwrite
2146 : // it, this makes the indexing logic a bit simpler.
2147 0 : if (this->processor_id() == 0)
2148 : {
2149 0 : auto cursor = gathered_qp_data[var].begin();
2150 0 : for (auto i : index_range(elem_side_pairs))
2151 : {
2152 0 : auto elem_side_pair = elem_side_pairs[i];
2153 0 : auto n_qp_this_elem_side = n_qp_per_elem_side[i];
2154 :
2155 : // Get reference to the [n_vars][n_qp] array for
2156 : // this Elem. We assign() into the vector of
2157 : // quadrature point values, which allocates space if
2158 : // it doesn't already exist.
2159 0 : auto & array = bf_map[elem_side_pair];
2160 :
2161 : // Possibly allocate space if this is data for a new
2162 : // element we haven't seen before.
2163 0 : if (array.empty())
2164 0 : array.resize(n_vars);
2165 :
2166 0 : array[var].assign(cursor, cursor + n_qp_this_elem_side);
2167 0 : std::advance(cursor, n_qp_this_elem_side);
2168 : }
2169 : }
2170 : }
2171 : } // end loop over basis functions
2172 0 : }
2173 :
2174 0 : void RBEIMEvaluation::node_gather_bfs()
2175 : {
2176 : // We need to gather _local_node_eim_basis_functions data from other
2177 : // procs for printing.
2178 : //
2179 : // Ideally, this could be accomplished by simply calling:
2180 : // this->comm().gather(/*root_id=*/0, _local_node_eim_basis_functions);
2181 : //
2182 : // but the data structure seems to be too complicated for this to
2183 : // work automatically. (I get some error about the function called
2184 : // being "private within this context".) Therefore, we have to
2185 : // gather the information manually.
2186 :
2187 : // So we can avoid calling this many times below
2188 0 : auto n_procs = this->n_processors();
2189 :
2190 : // In serial there's nothing to gather
2191 0 : if (n_procs == 1)
2192 0 : return;
2193 :
2194 : // Current assumption is that the number of basis functions stored on
2195 : // each processor is the same, the only thing that differs is the number
2196 : // of elements, so make sure that is the case now.
2197 0 : auto n_bf = _local_node_eim_basis_functions.size();
2198 0 : libmesh_assert(this->comm().verify(n_bf));
2199 :
2200 : // This function should never be called if there are no basis
2201 : // functions, so if it was, something went wrong.
2202 0 : libmesh_error_msg_if(!n_bf, "RBEIMEvaluation::gather_bfs() should not be called with 0 basis functions.");
2203 :
2204 : // The number of variables should be the same on all processors
2205 : // and we can get this from _local_eim_basis_functions. However,
2206 : // it may be that some processors have no local elements, so on
2207 : // those processors we cannot look up the size from
2208 : // _local_eim_basis_functions. As a result we use comm().max(n_vars)
2209 : // to make sure all processors agree on the final value.
2210 : std::size_t n_vars =
2211 0 : _local_node_eim_basis_functions[0].empty() ? 0 : _local_node_eim_basis_functions[0].begin()->second.size();
2212 0 : this->comm().max(n_vars);
2213 :
2214 : // Gather list of Node ids stored on each processor to proc 0. We
2215 : // use basis function 0 as an example and assume all the basis
2216 : // functions are distributed similarly.
2217 0 : unsigned int n_local_nodes = _local_node_eim_basis_functions[0].size();
2218 0 : std::vector<dof_id_type> node_ids;
2219 0 : node_ids.reserve(n_local_nodes);
2220 0 : for (const auto & pr : _local_node_eim_basis_functions[0])
2221 0 : node_ids.push_back(pr.first);
2222 0 : this->comm().gather(/*root_id=*/0, node_ids);
2223 :
2224 : // Now we construct a vector for each basis function, for each
2225 : // variable, which is ordered according to:
2226 : // [ [val for Node 0], [val for Node 1], ... [val for Node N] ]
2227 : // and gather it to processor 0.
2228 0 : std::vector<std::vector<Number>> gathered_node_data(n_vars);
2229 0 : for (auto bf : index_range(_local_node_eim_basis_functions))
2230 : {
2231 : // Clear any data from previous bf
2232 0 : for (auto var : index_range(gathered_node_data))
2233 : {
2234 0 : gathered_node_data[var].clear();
2235 0 : gathered_node_data[var].resize(n_local_nodes);
2236 : }
2237 :
2238 0 : unsigned int local_node_idx = 0;
2239 0 : for (const auto & pr : _local_node_eim_basis_functions[bf])
2240 : {
2241 : // array[n_vars] per Node
2242 0 : const auto & array = pr.second;
2243 0 : for (auto var : index_range(array))
2244 : {
2245 0 : gathered_node_data[var][local_node_idx] = array[var];
2246 : }
2247 :
2248 0 : local_node_idx++;
2249 : }
2250 :
2251 : // Reference to the data map for the current basis function.
2252 0 : auto & bf_map = _local_node_eim_basis_functions[bf];
2253 :
2254 0 : for (auto var : index_range(gathered_node_data))
2255 : {
2256 : // For each var, gather gathered_qp_data[var] onto processor
2257 : // 0. There apparently is not a gather overload for
2258 : // vector-of-vectors...
2259 0 : this->comm().gather(/*root_id=*/0, gathered_node_data[var]);
2260 :
2261 : // On processor 0, iterate over the gathered_qp_data[var]
2262 : // vector we just gathered, filling in the "small" vectors
2263 : // for each Elem. Note: here we ignore the fact that we
2264 : // already have the data on processor 0 and just overwrite
2265 : // it, this makes the indexing logic a bit simpler.
2266 0 : if (this->processor_id() == 0)
2267 : {
2268 0 : auto cursor = gathered_node_data[var].begin();
2269 0 : for (auto i : index_range(node_ids))
2270 : {
2271 0 : auto node_id = node_ids[i];
2272 :
2273 : // Get reference to the [n_vars] array for
2274 : // this Node. We assign() into the vector of
2275 : // node values, which allocates space if
2276 : // it doesn't already exist.
2277 0 : auto & array = bf_map[node_id];
2278 :
2279 : // Possibly allocate space if this is data for a new
2280 : // node we haven't seen before.
2281 0 : if (array.empty())
2282 0 : array.resize(n_vars);
2283 :
2284 : // There is only one value per variable per node, so
2285 : // we set the value by de-referencing cursor, and
2286 : // then advance the cursor by 1.
2287 0 : array[var] = *cursor;
2288 0 : std::advance(cursor, 1);
2289 : }
2290 : }
2291 : }
2292 : } // end loop over basis functions
2293 0 : }
2294 :
2295 :
2296 :
2297 0 : void RBEIMEvaluation::distribute_bfs(const System & sys)
2298 : {
2299 : // So we can avoid calling these many times below
2300 0 : auto n_procs = sys.comm().size();
2301 0 : auto rank = sys.comm().rank();
2302 :
2303 : // In serial there's nothing to distribute
2304 0 : if (n_procs == 1)
2305 0 : return;
2306 :
2307 : // Broadcast the number of basis functions from proc 0. After
2308 : // distributing, all procs should have the same number of basis
2309 : // functions.
2310 0 : auto n_bf = _local_eim_basis_functions.size();
2311 0 : sys.comm().broadcast(n_bf);
2312 :
2313 : // Allocate enough space to store n_bf basis functions on non-zero ranks
2314 0 : if (rank != 0)
2315 0 : _local_eim_basis_functions.resize(n_bf);
2316 :
2317 : // Broadcast the number of variables from proc 0. After
2318 : // distributing, all procs should have the same number of variables.
2319 0 : auto n_vars = _local_eim_basis_functions[0].begin()->second.size();
2320 0 : sys.comm().broadcast(n_vars);
2321 :
2322 : // Construct lists of elem ids owned by different processors
2323 0 : const MeshBase & mesh = sys.get_mesh();
2324 :
2325 0 : std::vector<dof_id_type> gathered_local_elem_ids;
2326 0 : gathered_local_elem_ids.reserve(mesh.n_elem());
2327 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
2328 0 : gathered_local_elem_ids.push_back(elem->id());
2329 :
2330 : // I _think_ the local elem ids are likely to already be sorted in
2331 : // ascending order, since that is how they are stored on the Mesh,
2332 : // but we can always just guarantee this to be on the safe side as
2333 : // well.
2334 0 : std::sort(gathered_local_elem_ids.begin(), gathered_local_elem_ids.end());
2335 :
2336 : // Gather the number of local elems from all procs to proc 0
2337 0 : auto n_local_elems = gathered_local_elem_ids.size();
2338 0 : std::vector<std::size_t> gathered_n_local_elems = {n_local_elems};
2339 0 : sys.comm().gather(/*root_id=*/0, gathered_n_local_elems);
2340 :
2341 : // Gather the elem ids owned by each processor onto processor 0.
2342 0 : sys.comm().gather(/*root_id=*/0, gathered_local_elem_ids);
2343 :
2344 : // Construct vectors of "start" and "one-past-the-end" indices into
2345 : // the gathered_local_elem_ids vector for each proc. Only valid on
2346 : // processor 0.
2347 0 : std::vector<std::size_t> start_elem_ids_index, end_elem_ids_index;
2348 :
2349 0 : if (rank == 0)
2350 : {
2351 0 : start_elem_ids_index.resize(n_procs);
2352 0 : start_elem_ids_index[0] = 0;
2353 0 : for (processor_id_type p=1; p<n_procs; ++p)
2354 0 : start_elem_ids_index[p] = start_elem_ids_index[p-1] + gathered_n_local_elems[p-1];
2355 :
2356 0 : end_elem_ids_index.resize(n_procs);
2357 0 : end_elem_ids_index[n_procs - 1] = gathered_local_elem_ids.size();
2358 0 : for (processor_id_type p=0; p<n_procs - 1; ++p)
2359 0 : end_elem_ids_index[p] = start_elem_ids_index[p+1];
2360 : }
2361 :
2362 : // On processor 0, using basis function 0 and variable 0, prepare a
2363 : // vector with the number of qps per Elem. Then scatter this vector
2364 : // out to the processors that require it. The order of this vector
2365 : // matches the gathered_local_elem_ids ordering. The counts will be
2366 : // gathered_n_local_elems, since there will be one qp count per Elem.
2367 0 : std::vector<unsigned int> n_qp_per_elem_data;
2368 :
2369 : // On rank 0, the "counts" vector holds the number of floating point values that
2370 : // are to be scattered to each proc. It is only required on proc 0.
2371 0 : std::vector<int> counts;
2372 :
2373 0 : if (rank == 0)
2374 : {
2375 0 : n_qp_per_elem_data.reserve(gathered_local_elem_ids.size());
2376 0 : counts.resize(n_procs);
2377 :
2378 0 : auto & bf_map = _local_eim_basis_functions[0];
2379 :
2380 0 : for (processor_id_type p=0; p<n_procs; ++p)
2381 : {
2382 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2383 : {
2384 0 : auto elem_id = gathered_local_elem_ids[e];
2385 :
2386 : // Get reference to array[n_vars][n_qp] for current Elem.
2387 : // Throws an error if the required elem_id is not found.
2388 0 : const auto & array = libmesh_map_find(bf_map, elem_id);
2389 :
2390 0 : auto n_qps = array[0].size();
2391 :
2392 : // We use var==0 to set the number of qps for all vars
2393 0 : n_qp_per_elem_data.push_back(n_qps);
2394 :
2395 : // Accumulate the count for this proc
2396 0 : counts[p] += n_qps;
2397 : } // end for (e)
2398 : } // end for proc_id
2399 : } // if (rank == 0)
2400 :
2401 : // Now scatter the n_qp_per_elem_data to all procs (must call the
2402 : // scatter on all procs, it is a collective).
2403 : {
2404 0 : std::vector<unsigned int> recv;
2405 0 : std::vector<int> tmp(gathered_n_local_elems.begin(), gathered_n_local_elems.end());
2406 0 : sys.comm().scatter(n_qp_per_elem_data, tmp, recv, /*root_id=*/0);
2407 :
2408 : // Now swap n_qp_per_elem_data and recv. All processors now have a
2409 : // vector of length n_local_elems containing the number of
2410 : // quadarature points per Elem.
2411 0 : n_qp_per_elem_data.swap(recv);
2412 : }
2413 :
2414 : // For each basis function and each variable, build a vector
2415 : // of qp data in the Elem ordering given by the
2416 : // gathered_local_elem_ids, then call
2417 : //
2418 : // sys.comm().scatter(data, counts, recv, /*root_id=*/0);
2419 0 : std::vector<std::vector<Number>> qp_data(n_vars);
2420 0 : if (rank == 0)
2421 : {
2422 : // The total amount of qp data is given by summing the entries
2423 : // of the "counts" vector.
2424 : auto n_qp_data =
2425 0 : std::accumulate(counts.begin(), counts.end(), 0u);
2426 :
2427 : // On processor 0, reserve enough space to hold all the qp
2428 : // data for a single basis function for each var.
2429 0 : for (auto var : index_range(qp_data))
2430 0 : qp_data[var].reserve(n_qp_data);
2431 : }
2432 :
2433 : // The recv_qp_data vector will be used on the receiving end of all
2434 : // the scatters below.
2435 0 : std::vector<Number> recv_qp_data;
2436 :
2437 : // Loop from 0..n_bf on _all_ procs, since the scatters inside this
2438 : // loop are collective.
2439 0 : for (auto bf : make_range(n_bf))
2440 : {
2441 : // Prepare data for scattering (only on proc 0)
2442 0 : if (rank == 0)
2443 : {
2444 : // Reference to the data map for the current basis function.
2445 0 : auto & bf_map = _local_eim_basis_functions[bf];
2446 :
2447 : // Clear any data from previous bf
2448 0 : for (auto var : index_range(qp_data))
2449 0 : qp_data[var].clear();
2450 :
2451 0 : for (processor_id_type p=0; p<n_procs; ++p)
2452 : {
2453 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2454 : {
2455 0 : auto elem_id = gathered_local_elem_ids[e];
2456 :
2457 : // Get reference to array[n_vars][n_qp] for current Elem.
2458 : // Throws an error if the required elem_id is not found.
2459 0 : const auto & array = libmesh_map_find(bf_map, elem_id);
2460 :
2461 0 : for (auto var : index_range(array))
2462 : {
2463 : // Insert all qp values for this var
2464 0 : qp_data[var].insert(/*insert at*/qp_data[var].end(),
2465 0 : /*data start*/array[var].begin(),
2466 0 : /*data end*/array[var].end());
2467 : } // end for (var)
2468 : } // end for (e)
2469 : } // end for proc_id
2470 : } // end if rank==0
2471 :
2472 : // Perform the scatters (all procs)
2473 0 : for (auto var : make_range(n_vars))
2474 : {
2475 : // Do the scatter for the current var
2476 0 : sys.comm().scatter(qp_data[var], counts, recv_qp_data, /*root_id=*/0);
2477 :
2478 0 : if (rank != 0)
2479 : {
2480 : // Store the scattered data we received in _local_eim_basis_functions[bf]
2481 0 : auto & bf_map = _local_eim_basis_functions[bf];
2482 0 : auto cursor = recv_qp_data.begin();
2483 :
2484 0 : for (auto i : index_range(gathered_local_elem_ids))
2485 : {
2486 0 : auto elem_id = gathered_local_elem_ids[i];
2487 0 : auto n_qp_this_elem = n_qp_per_elem_data[i];
2488 0 : auto & array = bf_map[elem_id];
2489 :
2490 : // Create space to store the data if it doesn't already exist.
2491 0 : if (array.empty())
2492 0 : array.resize(n_vars);
2493 :
2494 0 : array[var].assign(cursor, cursor + n_qp_this_elem);
2495 0 : std::advance(cursor, n_qp_this_elem);
2496 : }
2497 : } // if (rank != 0)
2498 : } // end for (var)
2499 : } // end for (bf)
2500 :
2501 : // Now that the scattering is done, delete non-local Elem
2502 : // information from processor 0's _local_eim_basis_functions data
2503 : // structure.
2504 0 : if (rank == 0)
2505 : {
2506 0 : for (processor_id_type p=1; p<n_procs; ++p)
2507 : {
2508 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2509 : {
2510 0 : auto elem_id = gathered_local_elem_ids[e];
2511 :
2512 : // Delete this Elem's information from every basis function.
2513 0 : for (auto & bf_map : _local_eim_basis_functions)
2514 0 : bf_map.erase(elem_id);
2515 : } // end for (e)
2516 : } // end for proc_id
2517 : } // if (rank == 0)
2518 0 : }
2519 :
2520 0 : void RBEIMEvaluation::side_distribute_bfs(const System & sys)
2521 : {
2522 : // So we can avoid calling these many times below
2523 0 : auto n_procs = sys.comm().size();
2524 0 : auto rank = sys.comm().rank();
2525 :
2526 : // In serial there's nothing to distribute
2527 0 : if (n_procs == 1)
2528 0 : return;
2529 :
2530 : // Broadcast the number of basis functions from proc 0. After
2531 : // distributing, all procs should have the same number of basis
2532 : // functions.
2533 0 : auto n_bf = _local_side_eim_basis_functions.size();
2534 0 : sys.comm().broadcast(n_bf);
2535 :
2536 : // Allocate enough space to store n_bf basis functions on non-zero ranks
2537 0 : if (rank != 0)
2538 0 : _local_side_eim_basis_functions.resize(n_bf);
2539 :
2540 : // Broadcast the number of variables from proc 0. After
2541 : // distributing, all procs should have the same number of variables.
2542 0 : auto n_vars = _local_side_eim_basis_functions[0].begin()->second.size();
2543 0 : sys.comm().broadcast(n_vars);
2544 :
2545 : const std::set<boundary_id_type> & parametrized_function_boundary_ids =
2546 0 : get_parametrized_function().get_parametrized_function_boundary_ids();
2547 :
2548 : // Construct lists of elem ids owned by different processors
2549 0 : const MeshBase & mesh = sys.get_mesh();
2550 :
2551 : // BoundaryInfo and related data structures
2552 0 : const auto & binfo = mesh.get_boundary_info();
2553 0 : std::vector<boundary_id_type> side_boundary_ids;
2554 :
2555 0 : std::vector<dof_id_type> gathered_local_elem_ids;
2556 0 : std::vector<dof_id_type> gathered_local_side_indices;
2557 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
2558 : {
2559 0 : for (unsigned int side = 0; side != elem->n_sides(); ++side)
2560 : {
2561 : // skip non-boundary elements
2562 0 : if (!elem->neighbor_ptr(side))
2563 : {
2564 0 : binfo.boundary_ids(elem, side, side_boundary_ids);
2565 :
2566 0 : bool has_side_boundary_id = false;
2567 0 : for (boundary_id_type side_boundary_id : side_boundary_ids)
2568 0 : if (parametrized_function_boundary_ids.count(side_boundary_id))
2569 : {
2570 0 : has_side_boundary_id = true;
2571 0 : break;
2572 : }
2573 :
2574 0 : if (has_side_boundary_id)
2575 : {
2576 0 : gathered_local_elem_ids.push_back(elem->id());
2577 0 : gathered_local_side_indices.push_back(side);
2578 : }
2579 : }
2580 : }
2581 :
2582 : // In the case of 2D elements, we also check the shellfaces
2583 0 : if (elem->dim() == 2)
2584 0 : for (unsigned int shellface_index=0; shellface_index<2; shellface_index++)
2585 : {
2586 0 : binfo.shellface_boundary_ids(elem, shellface_index, side_boundary_ids);
2587 :
2588 0 : bool has_side_boundary_id = false;
2589 0 : for (boundary_id_type side_boundary_id : side_boundary_ids)
2590 0 : if (parametrized_function_boundary_ids.count(side_boundary_id))
2591 : {
2592 0 : has_side_boundary_id = true;
2593 0 : break;
2594 : }
2595 :
2596 0 : if (has_side_boundary_id)
2597 : {
2598 : // We use shellface_index as the side_index since shellface boundary conditions
2599 : // are stored separately from side boundary conditions in BoundaryInfo.
2600 0 : gathered_local_elem_ids.push_back(elem->id());
2601 0 : gathered_local_side_indices.push_back(shellface_index);
2602 : }
2603 : }
2604 0 : }
2605 :
2606 : // Gather the number of local elems from all procs to proc 0
2607 0 : auto n_local_elems = gathered_local_elem_ids.size();
2608 0 : std::vector<std::size_t> gathered_n_local_elems = {n_local_elems};
2609 0 : sys.comm().gather(/*root_id=*/0, gathered_n_local_elems);
2610 :
2611 : // Gather the (elem,side) owned by each processor onto processor 0.
2612 0 : sys.comm().gather(/*root_id=*/0, gathered_local_elem_ids);
2613 0 : sys.comm().gather(/*root_id=*/0, gathered_local_side_indices);
2614 :
2615 : // Construct vectors of "start" and "one-past-the-end" indices into
2616 : // the gathered_local_elem_ids vector for each proc. Only valid on
2617 : // processor 0.
2618 0 : std::vector<std::size_t> start_elem_ids_index, end_elem_ids_index;
2619 :
2620 0 : if (rank == 0)
2621 : {
2622 0 : start_elem_ids_index.resize(n_procs);
2623 0 : start_elem_ids_index[0] = 0;
2624 0 : for (processor_id_type p=1; p<n_procs; ++p)
2625 0 : start_elem_ids_index[p] = start_elem_ids_index[p-1] + gathered_n_local_elems[p-1];
2626 :
2627 0 : end_elem_ids_index.resize(n_procs);
2628 0 : end_elem_ids_index[n_procs - 1] = gathered_local_elem_ids.size();
2629 0 : for (processor_id_type p=0; p<n_procs - 1; ++p)
2630 0 : end_elem_ids_index[p] = start_elem_ids_index[p+1];
2631 : }
2632 :
2633 : // On processor 0, using basis function 0 and variable 0, prepare a
2634 : // vector with the number of qps per Elem. Then scatter this vector
2635 : // out to the processors that require it. The order of this vector
2636 : // matches the gathered_local_elem_ids ordering. The counts will be
2637 : // gathered_n_local_elems, since there will be one qp count per Elem.
2638 0 : std::vector<unsigned int> n_qp_per_elem_data;
2639 :
2640 : // On rank 0, the "counts" vector holds the number of floating point values that
2641 : // are to be scattered to each proc. It is only required on proc 0.
2642 0 : std::vector<int> counts;
2643 :
2644 0 : if (rank == 0)
2645 : {
2646 0 : n_qp_per_elem_data.reserve(gathered_local_elem_ids.size());
2647 0 : counts.resize(n_procs);
2648 :
2649 0 : auto & bf_map = _local_side_eim_basis_functions[0];
2650 :
2651 0 : for (processor_id_type p=0; p<n_procs; ++p)
2652 : {
2653 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2654 : {
2655 0 : auto elem_id = gathered_local_elem_ids[e];
2656 0 : auto side_index = gathered_local_side_indices[e];
2657 :
2658 : // Get reference to array[n_vars][n_qp] for current Elem.
2659 : // Throws an error if the required elem_id is not found.
2660 0 : const auto & array = libmesh_map_find(bf_map, std::make_pair(elem_id, side_index));
2661 :
2662 0 : auto n_qps = array[0].size();
2663 :
2664 : // We use var==0 to set the number of qps for all vars
2665 0 : n_qp_per_elem_data.push_back(n_qps);
2666 :
2667 : // Accumulate the count for this proc
2668 0 : counts[p] += n_qps;
2669 : } // end for (e)
2670 : } // end for proc_id
2671 : } // if (rank == 0)
2672 :
2673 : // Now scatter the n_qp_per_elem_data to all procs (must call the
2674 : // scatter on all procs, it is a collective).
2675 : {
2676 0 : std::vector<unsigned int> recv;
2677 0 : std::vector<int> tmp(gathered_n_local_elems.begin(), gathered_n_local_elems.end());
2678 0 : sys.comm().scatter(n_qp_per_elem_data, tmp, recv, /*root_id=*/0);
2679 :
2680 : // Now swap n_qp_per_elem_data and recv. All processors now have a
2681 : // vector of length n_local_elems containing the number of
2682 : // quadarature points per Elem.
2683 0 : n_qp_per_elem_data.swap(recv);
2684 : }
2685 :
2686 : // For each basis function and each variable, build a vector
2687 : // of qp data in the Elem ordering given by the
2688 : // gathered_local_elem_ids, then call
2689 : //
2690 : // sys.comm().scatter(data, counts, recv, /*root_id=*/0);
2691 0 : std::vector<std::vector<Number>> qp_data(n_vars);
2692 0 : if (rank == 0)
2693 : {
2694 : // The total amount of qp data is given by summing the entries
2695 : // of the "counts" vector.
2696 : auto n_qp_data =
2697 0 : std::accumulate(counts.begin(), counts.end(), 0u);
2698 :
2699 : // On processor 0, reserve enough space to hold all the qp
2700 : // data for a single basis function for each var.
2701 0 : for (auto var : index_range(qp_data))
2702 0 : qp_data[var].reserve(n_qp_data);
2703 : }
2704 :
2705 : // The recv_qp_data vector will be used on the receiving end of all
2706 : // the scatters below.
2707 0 : std::vector<Number> recv_qp_data;
2708 :
2709 : // Loop from 0..n_bf on _all_ procs, since the scatters inside this
2710 : // loop are collective.
2711 0 : for (auto bf : make_range(n_bf))
2712 : {
2713 : // Prepare data for scattering (only on proc 0)
2714 0 : if (rank == 0)
2715 : {
2716 : // Reference to the data map for the current basis function.
2717 0 : auto & bf_map = _local_side_eim_basis_functions[bf];
2718 :
2719 : // Clear any data from previous bf
2720 0 : for (auto var : index_range(qp_data))
2721 0 : qp_data[var].clear();
2722 :
2723 0 : for (processor_id_type p=0; p<n_procs; ++p)
2724 : {
2725 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2726 : {
2727 0 : auto elem_id = gathered_local_elem_ids[e];
2728 0 : auto side_index = gathered_local_side_indices[e];
2729 :
2730 : // Get reference to array[n_vars][n_qp] for current Elem.
2731 : // Throws an error if the required (elem,side) is not found.
2732 0 : const auto & array = libmesh_map_find(bf_map, std::make_pair(elem_id, side_index));
2733 :
2734 0 : for (auto var : index_range(array))
2735 : {
2736 : // Insert all qp values for this var
2737 0 : qp_data[var].insert(/*insert at*/qp_data[var].end(),
2738 0 : /*data start*/array[var].begin(),
2739 0 : /*data end*/array[var].end());
2740 : } // end for (var)
2741 : } // end for (e)
2742 : } // end for proc_id
2743 : } // end if rank==0
2744 :
2745 : // Perform the scatters (all procs)
2746 0 : for (auto var : make_range(n_vars))
2747 : {
2748 : // Do the scatter for the current var
2749 0 : sys.comm().scatter(qp_data[var], counts, recv_qp_data, /*root_id=*/0);
2750 :
2751 0 : if (rank != 0)
2752 : {
2753 : // Store the scattered data we received in _local_side_eim_basis_functions[bf]
2754 0 : auto & bf_map = _local_side_eim_basis_functions[bf];
2755 0 : auto cursor = recv_qp_data.begin();
2756 :
2757 0 : for (auto i : index_range(gathered_local_elem_ids))
2758 : {
2759 0 : auto elem_id = gathered_local_elem_ids[i];
2760 0 : auto side_index = gathered_local_side_indices[i];
2761 0 : auto n_qp_this_elem = n_qp_per_elem_data[i];
2762 0 : auto & array = bf_map[std::make_pair(elem_id, side_index)];
2763 :
2764 : // Create space to store the data if it doesn't already exist.
2765 0 : if (array.empty())
2766 0 : array.resize(n_vars);
2767 :
2768 0 : array[var].assign(cursor, cursor + n_qp_this_elem);
2769 0 : std::advance(cursor, n_qp_this_elem);
2770 : }
2771 : } // if (rank != 0)
2772 : } // end for (var)
2773 : } // end for (bf)
2774 :
2775 : // Now that the scattering is done, delete non-local Elem
2776 : // information from processor 0's _local_side_eim_basis_functions data
2777 : // structure.
2778 0 : if (rank == 0)
2779 : {
2780 0 : for (processor_id_type p=1; p<n_procs; ++p)
2781 : {
2782 0 : for (auto e : make_range(start_elem_ids_index[p], end_elem_ids_index[p]))
2783 : {
2784 0 : auto elem_id = gathered_local_elem_ids[e];
2785 0 : auto side_index = gathered_local_side_indices[e];
2786 :
2787 : // Delete this Elem's information from every basis function.
2788 0 : for (auto & bf_map : _local_side_eim_basis_functions)
2789 0 : bf_map.erase(std::make_pair(elem_id, side_index));
2790 : } // end for (e)
2791 : } // end for proc_id
2792 : } // if (rank == 0)
2793 0 : }
2794 :
2795 0 : void RBEIMEvaluation::node_distribute_bfs(const System & sys)
2796 : {
2797 : // So we can avoid calling these many times below
2798 0 : auto n_procs = sys.comm().size();
2799 0 : auto rank = sys.comm().rank();
2800 :
2801 : // In serial there's nothing to distribute
2802 0 : if (n_procs == 1)
2803 0 : return;
2804 :
2805 : // Broadcast the number of basis functions from proc 0. After
2806 : // distributing, all procs should have the same number of basis
2807 : // functions.
2808 0 : auto n_bf = _local_node_eim_basis_functions.size();
2809 0 : sys.comm().broadcast(n_bf);
2810 :
2811 : // Allocate enough space to store n_bf basis functions on non-zero ranks
2812 0 : if (rank != 0)
2813 0 : _local_node_eim_basis_functions.resize(n_bf);
2814 :
2815 : // Broadcast the number of variables from proc 0. After
2816 : // distributing, all procs should have the same number of variables.
2817 0 : auto n_vars = _local_node_eim_basis_functions[0].begin()->second.size();
2818 0 : sys.comm().broadcast(n_vars);
2819 :
2820 : // Construct lists of elem ids owned by different processors
2821 0 : const MeshBase & mesh = sys.get_mesh();
2822 :
2823 0 : std::vector<dof_id_type> gathered_local_node_ids;
2824 : {
2825 : const std::set<boundary_id_type> & parametrized_function_boundary_ids =
2826 0 : get_parametrized_function().get_parametrized_function_boundary_ids();
2827 :
2828 0 : const auto & binfo = mesh.get_boundary_info();
2829 :
2830 : // Make a set with all the nodes that have nodesets. Use
2831 : // a set so that we don't have any duplicate entries. We
2832 : // deal with duplicate entries below by getting all boundary
2833 : // IDs on each node.
2834 0 : std::set<dof_id_type> nodes_with_nodesets;
2835 0 : for (const auto & t : binfo.build_node_list())
2836 0 : nodes_with_nodesets.insert(std::get<0>(t));
2837 :
2838 : // To be filled in by BoundaryInfo calls in loop below
2839 0 : std::vector<boundary_id_type> node_boundary_ids;
2840 :
2841 0 : for(dof_id_type node_id : nodes_with_nodesets)
2842 : {
2843 0 : const Node * node = mesh.node_ptr(node_id);
2844 :
2845 0 : if (node->processor_id() != mesh.comm().rank())
2846 0 : continue;
2847 :
2848 0 : binfo.boundary_ids(node, node_boundary_ids);
2849 :
2850 0 : bool has_node_boundary_id = false;
2851 0 : for(boundary_id_type node_boundary_id : node_boundary_ids)
2852 0 : if(parametrized_function_boundary_ids.count(node_boundary_id))
2853 : {
2854 0 : has_node_boundary_id = true;
2855 0 : break;
2856 : }
2857 :
2858 0 : if(has_node_boundary_id)
2859 : {
2860 0 : gathered_local_node_ids.push_back(node_id);
2861 : }
2862 : }
2863 : }
2864 :
2865 : // I _think_ the local node ids are likely to already be sorted in
2866 : // ascending order, since that is how they are stored on the Mesh,
2867 : // but we can always just guarantee this to be on the safe side as
2868 : // well.
2869 0 : std::sort(gathered_local_node_ids.begin(), gathered_local_node_ids.end());
2870 :
2871 : // Gather the number of local nodes from all procs to proc 0
2872 0 : auto n_local_nodes = gathered_local_node_ids.size();
2873 0 : std::vector<std::size_t> gathered_n_local_nodes = {n_local_nodes};
2874 0 : sys.comm().gather(/*root_id=*/0, gathered_n_local_nodes);
2875 :
2876 : // Gather the node ids owned by each processor onto processor 0.
2877 0 : sys.comm().gather(/*root_id=*/0, gathered_local_node_ids);
2878 :
2879 : // Construct vectors of "start" and "one-past-the-end" indices into
2880 : // the gathered_local_node_ids vector for each proc. Only valid on
2881 : // processor 0.
2882 0 : std::vector<std::size_t> start_node_ids_index, end_node_ids_index;
2883 :
2884 0 : if (rank == 0)
2885 : {
2886 0 : start_node_ids_index.resize(n_procs);
2887 0 : start_node_ids_index[0] = 0;
2888 0 : for (processor_id_type p=1; p<n_procs; ++p)
2889 0 : start_node_ids_index[p] = start_node_ids_index[p-1] + gathered_n_local_nodes[p-1];
2890 :
2891 0 : end_node_ids_index.resize(n_procs);
2892 0 : end_node_ids_index[n_procs - 1] = gathered_local_node_ids.size();
2893 0 : for (processor_id_type p=0; p<n_procs - 1; ++p)
2894 0 : end_node_ids_index[p] = start_node_ids_index[p+1];
2895 : }
2896 :
2897 : // On processor 0, using basis function 0 and variable 0, prepare a
2898 : // vector with the nodes. Then scatter this vector
2899 : // out to the processors that require it. The order of this vector
2900 : // matches the gathered_local_node_ids ordering. The counts will be
2901 : // gathered_n_local_nodes.
2902 :
2903 : // On rank 0, the "counts" vector holds the number of floating point values that
2904 : // are to be scattered to each proc. It is only required on proc 0.
2905 0 : std::vector<int> counts;
2906 :
2907 0 : if (rank == 0)
2908 : {
2909 0 : counts.resize(n_procs);
2910 :
2911 0 : for (processor_id_type p=0; p<n_procs; ++p)
2912 : {
2913 0 : auto node_ids_range = (end_node_ids_index[p] - start_node_ids_index[p]);
2914 :
2915 : // Accumulate the count for this proc
2916 0 : counts[p] += node_ids_range;
2917 : } // end for proc_id
2918 : } // if (rank == 0)
2919 :
2920 : // The recv_node_data vector will be used on the receiving end of all
2921 : // the scatters below.
2922 0 : std::vector<Number> recv_node_data;
2923 :
2924 : // For each basis function and each variable, build a vector
2925 : // data in the Node ordering given by the
2926 : // gathered_local_node_ids, then call
2927 : //
2928 : // sys.comm().scatter(data, counts, recv, /*root_id=*/0);
2929 0 : std::vector<std::vector<Number>> node_data(n_vars);
2930 :
2931 : // We also reserve space in node_data, since we will push_back into it below.
2932 0 : int count_sum = std::accumulate(counts.begin(), counts.end(), 0);
2933 0 : for (auto var : index_range(node_data))
2934 0 : node_data[var].reserve(count_sum);
2935 :
2936 : // Loop from 0..n_bf on _all_ procs, since the scatters inside this
2937 : // loop are collective.
2938 0 : for (auto bf : make_range(n_bf))
2939 : {
2940 : // Prepare data for scattering (only on proc 0)
2941 0 : if (rank == 0)
2942 : {
2943 : // Reference to the data map for the current basis function.
2944 0 : auto & bf_map = _local_node_eim_basis_functions[bf];
2945 :
2946 : // Clear any data from previous bf (this does not change the capacity
2947 : // that was reserved above).
2948 0 : for (auto var : index_range(node_data))
2949 0 : node_data[var].clear();
2950 :
2951 0 : for (processor_id_type p=0; p<n_procs; ++p)
2952 : {
2953 0 : for (auto n : make_range(start_node_ids_index[p], end_node_ids_index[p]))
2954 : {
2955 0 : auto node_id = gathered_local_node_ids[n];
2956 :
2957 : // Get reference to array[n_vars] for current Node.
2958 : // Throws an error if the required node_id is not found.
2959 0 : const auto & array = libmesh_map_find(bf_map, node_id);
2960 :
2961 0 : for (auto var : index_range(array))
2962 0 : node_data[var].push_back(array[var]);
2963 : } // end for (n)
2964 : } // end for proc_id
2965 : } // end if rank==0
2966 :
2967 : // Perform the scatters (all procs)
2968 0 : for (auto var : make_range(n_vars))
2969 : {
2970 : // Do the scatter for the current var
2971 0 : sys.comm().scatter(node_data[var], counts, recv_node_data, /*root_id=*/0);
2972 :
2973 0 : if (rank != 0)
2974 : {
2975 : // Store the scattered data we received in _local_eim_basis_functions[bf]
2976 0 : auto & bf_map = _local_node_eim_basis_functions[bf];
2977 0 : auto cursor = recv_node_data.begin();
2978 :
2979 0 : for (auto i : index_range(gathered_local_node_ids))
2980 : {
2981 0 : auto node_id = gathered_local_node_ids[i];
2982 0 : auto & array = bf_map[node_id];
2983 :
2984 : // Create space to store the data if it doesn't already exist.
2985 0 : if (array.empty())
2986 0 : array.resize(n_vars);
2987 :
2988 : // There is only one value per variable per node, so
2989 : // we set the value by de-referencing cursor, and
2990 : // then advance the cursor by 1.
2991 0 : array[var] = *cursor;
2992 0 : std::advance(cursor, 1);
2993 : }
2994 : } // if (rank != 0)
2995 : } // end for (var)
2996 : } // end for (bf)
2997 :
2998 : // Now that the scattering is done, delete non-local Elem
2999 : // information from processor 0's _local_eim_basis_functions data
3000 : // structure.
3001 0 : if (rank == 0)
3002 : {
3003 0 : for (processor_id_type p=1; p<n_procs; ++p)
3004 : {
3005 0 : for (auto n : make_range(start_node_ids_index[p], end_node_ids_index[p]))
3006 : {
3007 0 : auto node_id = gathered_local_node_ids[n];
3008 :
3009 : // Delete this Node's information from every basis function.
3010 0 : for (auto & bf_map : _local_node_eim_basis_functions)
3011 0 : bf_map.erase(node_id);
3012 : } // end for (n)
3013 : } // end for proc_id
3014 : } // if (rank == 0)
3015 0 : }
3016 :
3017 0 : void RBEIMEvaluation::project_qp_data_vector_onto_system(System & /*sys*/,
3018 : const std::vector<Number> & /*bf_data*/,
3019 : const EIMVarGroupPlottingInfo & /*eim_vargroup*/,
3020 : const std::map<std::string,std::string> & /*extra_options*/)
3021 : {
3022 : // No-op by default, implement in subclasses if needed
3023 0 : }
3024 :
3025 0 : const std::vector<EIMVarGroupPlottingInfo> & RBEIMEvaluation::get_eim_vars_to_project_and_write() const
3026 : {
3027 0 : return _eim_vars_to_project_and_write;
3028 : }
3029 :
3030 0 : const std::set<unsigned int> & RBEIMEvaluation::scale_components_in_enrichment() const
3031 : {
3032 0 : return _scale_components_in_enrichment;
3033 : }
3034 :
3035 0 : bool RBEIMEvaluation::use_eim_error_indicator() const
3036 : {
3037 : // Return false by default, but we override this in subclasses
3038 : // for cases where we want to use the error indicator.
3039 0 : return false;
3040 : }
3041 :
3042 0 : void RBEIMEvaluation::set_eim_error_indicator_active(bool is_active)
3043 : {
3044 : // We skip setting _is_eim_error_indicator_active in the case that
3045 : // we have no parameters, since we do not use the EIM error indicator
3046 : // in that case. We also check if the number of interpolation points
3047 : // is larger than the number of EIM basis functions, since that is
3048 : // also always the case when the error indicator is active.
3049 0 : if ((get_n_params() > 0) && (get_n_interpolation_points() > get_n_basis_functions()))
3050 0 : _is_eim_error_indicator_active = (is_active && use_eim_error_indicator());
3051 0 : }
3052 :
3053 0 : std::pair<Real,Real> RBEIMEvaluation::get_eim_error_indicator(
3054 : Number error_indicator_rhs,
3055 : const DenseVector<Number> & eim_solution,
3056 : const DenseVector<Number> & eim_rhs)
3057 : {
3058 0 : DenseVector<Number> coeffs;
3059 0 : _error_indicator_interpolation_row.get_principal_subvector(eim_solution.size(), coeffs);
3060 :
3061 0 : Number EIM_val_at_error_indicator_pt = coeffs.dot(eim_solution);
3062 : Real error_indicator_val =
3063 0 : std::real(error_indicator_rhs - EIM_val_at_error_indicator_pt);
3064 :
3065 0 : Real normalization = 0.;
3066 0 : if (eim_error_indicator_normalization == RESIDUAL_SUM)
3067 : {
3068 : // This normalization is based on the sum of terms from the "EIM residual" calculation
3069 : // used in the calculation of error_indicator_val. This ensures that the error indicator
3070 : // will always be less than or equal to one, which is a useful property for an error
3071 : // indicator.
3072 0 : normalization =
3073 0 : std::abs(error_indicator_rhs) + std::abs(EIM_val_at_error_indicator_pt);
3074 : }
3075 0 : else if (eim_error_indicator_normalization == RESIDUAL_RHS)
3076 : {
3077 : // Normalize with respect to the right-hand side from the "EIM residual" calculation.
3078 0 : normalization = std::abs(error_indicator_rhs);
3079 : }
3080 0 : else if (eim_error_indicator_normalization == MAX_RHS)
3081 : {
3082 : // Normalize the error indicator based on the max-norm of the EIM RHS vector.
3083 : // This approach handles the case where different EIM variables have different
3084 : // magnitudes well, i.e. if error_indicator_val is based on a
3085 : // "small magnitude" variable, then by normalizing based on the entire
3086 : // RHS vector (which will typically include values from multiple different
3087 : // EIM variables) we will effectively scale down the error indicator
3088 : // corresponding to small variables, which is typically what we want.
3089 0 : normalization = std::max(eim_rhs.linfty_norm(), std::abs(error_indicator_rhs));
3090 : }
3091 : else
3092 : {
3093 0 : libmesh_error_msg("unsupported eim_error_indicator_normalization");
3094 : }
3095 :
3096 : // We avoid NaNs by setting normalization to 1 in the case that it is exactly 0.
3097 : // But we return the "original normalization" as well (as opposed to the modified
3098 : // normalization) since that can be useful information since it can indicate that
3099 : // the EIM approximation was identically zero, for example.
3100 0 : Real orig_normalization = normalization;
3101 0 : if (normalization == 0.)
3102 0 : normalization = 1.;
3103 :
3104 : // Return the relative error indicator, and the normalization that we used. By returning
3105 : // the normalization, we can subsequently recover the absolute error indicator if
3106 : // desired.
3107 : //
3108 : // We also optionally clamp the relative error indicator to 1.0 since this typically
3109 : // indicators 100% error and hence it may be the maximum value that we want to see.
3110 0 : Real rel_error_indicator = std::abs(error_indicator_val) / normalization;
3111 0 : if (limit_eim_error_indicator_to_one && (rel_error_indicator > 1.0))
3112 0 : rel_error_indicator = 1.0;
3113 :
3114 0 : return std::make_pair(rel_error_indicator, orig_normalization);
3115 : }
3116 :
3117 0 : const VectorizedEvalInput & RBEIMEvaluation::get_vec_eval_input() const
3118 : {
3119 0 : return _vec_eval_input;
3120 : }
3121 :
3122 0 : void RBEIMEvaluation::initialize_rb_property_map()
3123 : {
3124 0 : const auto & rb_property_map = get_parametrized_function().get_rb_property_map();
3125 : // Initialize rb_eim_eval VectorizedEvaluateInput from the one in rb_parametrized_function with
3126 : // empty sets as it will be filled by subclasses virtual functions.
3127 0 : for (const auto & [key, val] : rb_property_map)
3128 0 : _vec_eval_input.rb_property_map[key] = {};
3129 0 : }
3130 :
3131 0 : const DenseVector<Number> & RBEIMEvaluation::get_error_indicator_interpolation_row() const
3132 : {
3133 0 : return _error_indicator_interpolation_row;
3134 : }
3135 :
3136 0 : void RBEIMEvaluation::set_error_indicator_interpolation_row(const DenseVector<Number> & extra_point_row)
3137 : {
3138 0 : _error_indicator_interpolation_row = extra_point_row;
3139 0 : }
3140 :
3141 : } // namespace libMesh
|
