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_construction_base.h"
22 : #include "libmesh/rb_construction.h"
23 : #include "libmesh/rb_assembly_expansion.h"
24 : #include "libmesh/rb_evaluation.h"
25 : #include "libmesh/elem_assembly.h"
26 :
27 : // LibMesh includes
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/sparse_matrix.h"
30 : #include "libmesh/dof_map.h"
31 : #include "libmesh/libmesh_logging.h"
32 : #include "libmesh/equation_systems.h"
33 : #include "libmesh/exodusII_io.h"
34 : #include "libmesh/gmv_io.h"
35 : #include "libmesh/linear_solver.h"
36 : #include "libmesh/getpot.h"
37 : #include "libmesh/int_range.h"
38 : #include "libmesh/mesh_base.h"
39 : #include "libmesh/parallel.h"
40 : #include "libmesh/xdr_cxx.h"
41 : #include "libmesh/timestamp.h"
42 : #include "libmesh/petsc_linear_solver.h"
43 : #include "libmesh/dg_fem_context.h"
44 : #include "libmesh/dirichlet_boundaries.h"
45 : #include "libmesh/zero_function.h"
46 : #include "libmesh/coupling_matrix.h"
47 : #include "libmesh/face_tri3_subdivision.h"
48 : #include "libmesh/quadrature.h"
49 : #include "libmesh/utility.h"
50 :
51 : // C++ includes
52 : #include <sys/types.h>
53 : #include <sys/stat.h>
54 : #include <errno.h>
55 : #include <fstream>
56 : #include <sstream>
57 : #include <limits>
58 : #include <stdlib.h> // mkstemps on Linux
59 : #ifdef LIBMESH_HAVE_UNISTD_H
60 : #include <unistd.h> // mkstemps on MacOS
61 : #endif
62 :
63 : namespace libMesh
64 : {
65 :
66 568 : RBConstruction::RBConstruction (EquationSystems & es,
67 : const std::string & name_in,
68 568 : const unsigned int number_in)
69 : : Parent(es, name_in, number_in),
70 536 : inner_product_solver(LinearSolver<Number>::build(es.comm())),
71 536 : extra_linear_solver(nullptr),
72 536 : inner_product_matrix(SparseMatrix<Number>::build(es.comm())),
73 536 : skip_residual_in_train_reduced_basis(false),
74 536 : exit_on_repeated_greedy_parameters(true),
75 536 : impose_internal_fluxes(false),
76 536 : skip_degenerate_sides(true),
77 536 : compute_RB_inner_product(false),
78 536 : store_non_dirichlet_operators(false),
79 536 : store_untransformed_basis(false),
80 536 : use_empty_rb_solve_in_greedy(true),
81 536 : Fq_representor_innerprods_computed(false),
82 536 : Nmax(0),
83 536 : delta_N(1),
84 536 : output_dual_innerprods_computed(false),
85 536 : assert_convergence(true),
86 536 : rb_eval(nullptr),
87 536 : inner_product_assembly(nullptr),
88 536 : use_energy_inner_product(false),
89 536 : rel_training_tolerance(1.e-4),
90 536 : abs_training_tolerance(1.e-12),
91 536 : normalize_rb_bound_in_greedy(false),
92 536 : RB_training_type("Greedy"),
93 536 : _preevaluate_thetas_flag(false),
94 1104 : _preevaluate_thetas_completed(false)
95 : {
96 : // set assemble_before_solve flag to false
97 : // so that we control matrix assembly.
98 568 : assemble_before_solve = false;
99 568 : }
100 :
101 2144 : RBConstruction::~RBConstruction () = default;
102 :
103 0 : void RBConstruction::clear()
104 : {
105 0 : LOG_SCOPE("clear()", "RBConstruction");
106 :
107 0 : Parent::clear();
108 :
109 0 : Aq_vector.clear();
110 0 : Fq_vector.clear();
111 0 : outputs_vector.clear();
112 :
113 0 : if (store_non_dirichlet_operators)
114 : {
115 0 : non_dirichlet_Aq_vector.clear();
116 0 : non_dirichlet_Fq_vector.clear();
117 0 : non_dirichlet_outputs_vector.clear();
118 : }
119 :
120 : // Also delete the Fq representors
121 0 : Fq_representor.clear();
122 :
123 : // Set Fq_representor_innerprods_computed flag to false now
124 : // that we've cleared the Fq representors
125 0 : Fq_representor_innerprods_computed = false;
126 0 : }
127 :
128 568 : std::string RBConstruction::system_type () const
129 : {
130 568 : return "RBConstruction";
131 : }
132 :
133 160467 : void RBConstruction::solve_for_matrix_and_rhs(LinearSolver<Number> & input_solver,
134 : SparseMatrix<Number> & input_matrix,
135 : NumericVector<Number> & input_rhs)
136 : {
137 : // This is similar to LinearImplicitSysmte::solve()
138 :
139 : // Get a reference to the EquationSystems
140 : const EquationSystems & es =
141 9156 : this->get_equation_systems();
142 :
143 : // If the linear solver hasn't been initialized, we do so here.
144 160467 : input_solver.init();
145 :
146 : // Get the user-specifiied linear solver tolerance
147 : const double tol =
148 160467 : double(es.parameters.get<Real>("linear solver tolerance"));
149 :
150 : // Get the user-specified maximum # of linear solver iterations
151 : const unsigned int maxits =
152 160467 : es.parameters.get<unsigned int>("linear solver maximum iterations");
153 :
154 : // It's good practice to clear the solution vector first since it can
155 : // affect convergence of iterative solvers
156 160467 : solution->zero();
157 :
158 : // Solve the linear system.
159 : // Store the number of linear iterations required to
160 : // solve and the final residual.
161 160467 : std::tie(_n_linear_iterations, _final_linear_residual) =
162 169623 : input_solver.solve (input_matrix, *solution, input_rhs, tol, maxits);
163 :
164 160467 : get_dof_map().enforce_constraints_exactly(*this);
165 :
166 : // Update the system after the solve
167 160467 : this->update();
168 160467 : }
169 :
170 568 : void RBConstruction::set_rb_evaluation(RBEvaluation & rb_eval_in)
171 : {
172 568 : rb_eval = &rb_eval_in;
173 568 : }
174 :
175 30451994 : RBEvaluation & RBConstruction::get_rb_evaluation()
176 : {
177 30451994 : libmesh_error_msg_if(!rb_eval, "Error: RBEvaluation object hasn't been initialized yet");
178 :
179 30451994 : return *rb_eval;
180 : }
181 :
182 90453 : const RBEvaluation & RBConstruction::get_rb_evaluation() const
183 : {
184 90453 : libmesh_error_msg_if(!rb_eval, "Error: RBEvaluation object hasn't been initialized yet");
185 :
186 90453 : return *rb_eval;
187 : }
188 :
189 354 : bool RBConstruction::is_rb_eval_initialized() const
190 : {
191 354 : return (rb_eval != nullptr);
192 : }
193 :
194 6111903 : RBThetaExpansion & RBConstruction::get_rb_theta_expansion()
195 : {
196 6111903 : return get_rb_evaluation().get_rb_theta_expansion();
197 : }
198 :
199 2342 : const RBThetaExpansion & RBConstruction::get_rb_theta_expansion() const
200 : {
201 2342 : return get_rb_evaluation().get_rb_theta_expansion();
202 : }
203 :
204 284 : void RBConstruction::process_parameters_file (const std::string & parameters_filename)
205 : {
206 : // First read in data from input_filename
207 584 : GetPot infile(parameters_filename);
208 :
209 284 : const unsigned int n_training_samples = infile("n_training_samples",0);
210 284 : const bool deterministic_training = infile("deterministic_training",false);
211 284 : unsigned int training_parameters_random_seed_in =
212 : static_cast<unsigned int>(-1);
213 284 : training_parameters_random_seed_in = infile("training_parameters_random_seed",
214 : training_parameters_random_seed_in);
215 284 : const bool quiet_mode_in = infile("quiet_mode", quiet_mode);
216 284 : const unsigned int Nmax_in = infile("Nmax", Nmax);
217 292 : const Real rel_training_tolerance_in = infile("rel_training_tolerance",
218 284 : rel_training_tolerance);
219 292 : const Real abs_training_tolerance_in = infile("abs_training_tolerance",
220 284 : abs_training_tolerance);
221 :
222 : // Initialize value to false, let the input file value override.
223 292 : const bool normalize_rb_bound_in_greedy_in = infile("normalize_rb_bound_in_greedy",
224 284 : false);
225 :
226 292 : const std::string RB_training_type_in = infile("RB_training_type", "Greedy");
227 :
228 : // Read in the parameters from the input file too
229 276 : unsigned int n_continuous_parameters = infile.vector_variable_size("parameter_names");
230 300 : RBParameters mu_min_in;
231 300 : RBParameters mu_max_in;
232 1563 : for (unsigned int i=0; i<n_continuous_parameters; i++)
233 : {
234 : // Read in the parameter names
235 1315 : std::string param_name = infile("parameter_names", "NONE", i);
236 :
237 : {
238 1279 : Real min_val = infile(param_name, 0., 0);
239 1279 : mu_min_in.set_value(param_name, min_val);
240 : }
241 :
242 : {
243 1279 : Real max_val = infile(param_name, 0., 1);
244 1279 : mu_max_in.set_value(param_name, max_val);
245 : }
246 : }
247 :
248 16 : std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
249 :
250 276 : unsigned int n_discrete_parameters = infile.vector_variable_size("discrete_parameter_names");
251 284 : for (unsigned int i=0; i<n_discrete_parameters; i++)
252 : {
253 0 : std::string param_name = infile("discrete_parameter_names", "NONE", i);
254 :
255 0 : unsigned int n_vals_for_param = infile.vector_variable_size(param_name);
256 0 : std::vector<Real> vals_for_param(n_vals_for_param);
257 0 : for (auto j : make_range(vals_for_param.size()))
258 0 : vals_for_param[j] = infile(param_name, 0., j);
259 :
260 0 : discrete_parameter_values_in[param_name] = vals_for_param;
261 : }
262 :
263 16 : std::map<std::string,bool> log_scaling_in;
264 : // For now, just set all entries to false.
265 : // TODO: Implement a decent way to specify log-scaling true/false
266 : // in the input text file
267 1563 : for (const auto & pr : mu_min_in)
268 1279 : log_scaling_in[pr.first] = false;
269 :
270 : // Set the parameters that have been read in
271 284 : set_rb_construction_parameters(n_training_samples,
272 : deterministic_training,
273 : static_cast<int>(training_parameters_random_seed_in),
274 : quiet_mode_in,
275 : Nmax_in,
276 : rel_training_tolerance_in,
277 : abs_training_tolerance_in,
278 : normalize_rb_bound_in_greedy_in,
279 : RB_training_type_in,
280 : mu_min_in,
281 : mu_max_in,
282 : discrete_parameter_values_in,
283 : log_scaling_in);
284 552 : }
285 :
286 284 : void RBConstruction::set_rb_construction_parameters(
287 : unsigned int n_training_samples_in,
288 : bool deterministic_training_in,
289 : int training_parameters_random_seed_in,
290 : bool quiet_mode_in,
291 : unsigned int Nmax_in,
292 : Real rel_training_tolerance_in,
293 : Real abs_training_tolerance_in,
294 : bool normalize_rb_bound_in_greedy_in,
295 : const std::string & RB_training_type_in,
296 : const RBParameters & mu_min_in,
297 : const RBParameters & mu_max_in,
298 : const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
299 : const std::map<std::string,bool> & log_scaling_in,
300 : std::map<std::string, std::vector<RBParameter>> * training_sample_list)
301 : {
302 : // Read in training_parameters_random_seed value. This is used to
303 : // seed the RNG when picking the training parameters. By default the
304 : // value is -1, which means use std::time to seed the RNG.
305 284 : set_training_random_seed(training_parameters_random_seed_in);
306 :
307 : // Set quiet mode
308 8 : set_quiet_mode(quiet_mode_in);
309 :
310 : // Initialize RB parameters
311 284 : set_Nmax(Nmax_in);
312 :
313 8 : set_rel_training_tolerance(rel_training_tolerance_in);
314 8 : set_abs_training_tolerance(abs_training_tolerance_in);
315 :
316 8 : set_normalize_rb_bound_in_greedy(normalize_rb_bound_in_greedy_in);
317 :
318 284 : set_RB_training_type(RB_training_type_in);
319 :
320 : // Initialize the parameter ranges and the parameters themselves
321 284 : initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_in);
322 :
323 8 : bool updated_deterministic_training = deterministic_training_in;
324 284 : if (training_sample_list && (this->get_parameters_min().n_parameters() > 3))
325 : {
326 : // In this case we force deterministic_training to be false because
327 : // a) deterministic training samples are not currrently supported with
328 : // more than 3 parameters, and
329 : // b) we will overwrite the training samples anyway in the call to
330 : // load_training_set() below, so we do not want to generate an
331 : // error due to deterministic training sample generation when
332 : // the samples will be overwritten anyway.
333 0 : updated_deterministic_training = false;
334 : }
335 :
336 284 : initialize_training_parameters(this->get_parameters_min(),
337 : this->get_parameters_max(),
338 : n_training_samples_in,
339 : log_scaling_in,
340 16 : updated_deterministic_training); // use deterministic parameters
341 :
342 284 : if (training_sample_list)
343 : {
344 : // Note that we must call initialize_training_parameters() before
345 : // load_training_set() in order to initialize the parameter vectors.
346 0 : load_training_set(*training_sample_list);
347 : }
348 284 : }
349 :
350 284 : void RBConstruction::print_info() const
351 : {
352 : // Print out info that describes the current setup
353 8 : libMesh::out << std::endl << "RBConstruction parameters:" << std::endl;
354 8 : libMesh::out << "system name: " << this->name() << std::endl;
355 16 : libMesh::out << "Nmax: " << Nmax << std::endl;
356 16 : libMesh::out << "Greedy relative error tolerance: " << get_rel_training_tolerance() << std::endl;
357 16 : libMesh::out << "Greedy absolute error tolerance: " << get_abs_training_tolerance() << std::endl;
358 16 : libMesh::out << "Do we normalize RB error bound in greedy? " << get_normalize_rb_bound_in_greedy() << std::endl;
359 284 : libMesh::out << "RB training type: " << get_RB_training_type() << std::endl;
360 284 : if (is_rb_eval_initialized())
361 : {
362 284 : libMesh::out << "Aq operators attached: " << get_rb_theta_expansion().get_n_A_terms() << std::endl;
363 284 : libMesh::out << "Fq functions attached: " << get_rb_theta_expansion().get_n_F_terms() << std::endl;
364 284 : libMesh::out << "n_outputs: " << get_rb_theta_expansion().get_n_outputs() << std::endl;
365 852 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
366 568 : libMesh::out << "output " << n << ", Q_l = " << get_rb_theta_expansion().get_n_output_terms(n) << std::endl;
367 : }
368 : else
369 : {
370 0 : libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
371 : }
372 284 : libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
373 1563 : for (const auto & pr : get_parameters())
374 1279 : if (!is_discrete_parameter(pr.first))
375 : {
376 36 : libMesh::out << "Parameter " << pr.first
377 1279 : << ": Min = " << get_parameter_min(pr.first)
378 1279 : << ", Max = " << get_parameter_max(pr.first) << std::endl;
379 : }
380 :
381 284 : print_discrete_parameter_values();
382 284 : libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
383 16 : libMesh::out << "quiet mode? " << is_quiet() << std::endl;
384 8 : libMesh::out << std::endl;
385 284 : }
386 :
387 71 : void RBConstruction::print_basis_function_orthogonality() const
388 : {
389 73 : std::unique_ptr<NumericVector<Number>> temp = solution->clone();
390 :
391 1491 : for (unsigned int i=0; i<get_rb_evaluation().get_n_basis_functions(); i++)
392 : {
393 29820 : for (unsigned int j=0; j<get_rb_evaluation().get_n_basis_functions(); j++)
394 : {
395 28400 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
396 28400 : Number value = temp->dot( get_rb_evaluation().get_basis_function(i) );
397 :
398 800 : libMesh::out << value << " ";
399 : }
400 40 : libMesh::out << std::endl;
401 : }
402 2 : libMesh::out << std::endl;
403 71 : }
404 :
405 568 : void RBConstruction::set_rb_assembly_expansion(RBAssemblyExpansion & rb_assembly_expansion_in)
406 : {
407 568 : rb_assembly_expansion = &rb_assembly_expansion_in;
408 568 : }
409 :
410 112 : RBAssemblyExpansion & RBConstruction::get_rb_assembly_expansion()
411 : {
412 112 : libmesh_error_msg_if(!rb_assembly_expansion, "Error: RBAssemblyExpansion object hasn't been initialized yet");
413 :
414 112 : return *rb_assembly_expansion;
415 : }
416 :
417 568 : void RBConstruction::set_inner_product_assembly(ElemAssembly & inner_product_assembly_in)
418 : {
419 568 : use_energy_inner_product = false;
420 568 : inner_product_assembly = &inner_product_assembly_in;
421 568 : }
422 :
423 0 : ElemAssembly & RBConstruction::get_inner_product_assembly()
424 : {
425 0 : libmesh_error_msg_if(use_energy_inner_product,
426 : "Error: inner_product_assembly not available since we're using energy inner-product");
427 :
428 0 : libmesh_error_msg_if(!inner_product_assembly,
429 : "Error: inner_product_assembly hasn't been initialized yet");
430 :
431 0 : return *inner_product_assembly;
432 : }
433 :
434 0 : void RBConstruction::set_energy_inner_product(const std::vector<Number> & energy_inner_product_coeffs_in)
435 : {
436 0 : use_energy_inner_product = true;
437 0 : energy_inner_product_coeffs = energy_inner_product_coeffs_in;
438 0 : }
439 :
440 0 : void RBConstruction::zero_constrained_dofs_on_vector(NumericVector<Number> & vector) const
441 : {
442 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
443 0 : const DofMap & dof_map = get_dof_map();
444 :
445 0 : for (dof_id_type i=dof_map.first_dof(); i<dof_map.end_dof(); i++)
446 : {
447 0 : if (get_dof_map().is_constrained_dof(i))
448 : {
449 0 : vector.set(i, 0.);
450 : }
451 : }
452 : #endif
453 :
454 0 : vector.close();
455 0 : }
456 :
457 4956 : bool RBConstruction::check_if_zero_truth_solve() const
458 : {
459 4956 : return (solution->l2_norm() == 0.);
460 : }
461 :
462 284 : void RBConstruction::initialize_rb_construction(bool skip_matrix_assembly,
463 : bool skip_vector_assembly)
464 : {
465 284 : if (!skip_matrix_assembly && !skip_vector_assembly)
466 : {
467 : // Check that the theta and assembly objects are consistently sized
468 8 : libmesh_assert_equal_to (get_rb_theta_expansion().get_n_A_terms(), get_rb_assembly_expansion().get_n_A_terms());
469 8 : libmesh_assert_equal_to (get_rb_theta_expansion().get_n_F_terms(), get_rb_assembly_expansion().get_n_F_terms());
470 8 : libmesh_assert_equal_to (get_rb_theta_expansion().get_n_outputs(), get_rb_assembly_expansion().get_n_outputs());
471 852 : for (unsigned int i=0; i<get_rb_theta_expansion().get_n_outputs(); i++)
472 : {
473 16 : libmesh_assert_equal_to (get_rb_theta_expansion().get_n_output_terms(i),
474 : get_rb_assembly_expansion().get_n_output_terms(i));
475 : }
476 : }
477 :
478 : // Perform the initialization
479 284 : allocate_data_structures();
480 284 : assemble_affine_expansion(skip_matrix_assembly, skip_vector_assembly);
481 :
482 : // inner_product_solver performs solves with the same matrix every time
483 : // hence we can set reuse_preconditioner(true).
484 284 : inner_product_solver->reuse_preconditioner(true);
485 :
486 : // The primary solver is used for truth solves and other solves that
487 : // require different matrices, so set reuse_preconditioner(false).
488 284 : get_linear_solver()->reuse_preconditioner(false);
489 :
490 284 : }
491 :
492 284 : void RBConstruction::assemble_affine_expansion(bool skip_matrix_assembly,
493 : bool skip_vector_assembly)
494 : {
495 284 : if (!skip_matrix_assembly)
496 : {
497 : // Assemble and store all of the matrices
498 284 : this->assemble_misc_matrices();
499 284 : this->assemble_all_affine_operators();
500 : }
501 :
502 284 : if (!skip_vector_assembly)
503 : {
504 : // Assemble and store all of the vectors
505 284 : this->assemble_all_affine_vectors();
506 284 : this->assemble_all_output_vectors();
507 : }
508 284 : }
509 :
510 284 : void RBConstruction::allocate_data_structures()
511 : {
512 : // Resize vectors for storing mesh-dependent data
513 284 : Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
514 284 : Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
515 :
516 : // Resize the Fq_representors and initialize each to nullptr.
517 : // These are basis independent and hence stored here, whereas
518 : // the Aq_representors are stored in RBEvaluation
519 284 : Fq_representor.resize(get_rb_theta_expansion().get_n_F_terms());
520 :
521 : // Initialize vectors for the inner products of the Fq representors
522 : // These are basis independent and therefore stored here.
523 284 : unsigned int Q_f_hat = get_rb_theta_expansion().get_n_F_terms()*(get_rb_theta_expansion().get_n_F_terms()+1)/2;
524 284 : Fq_representor_innerprods.resize(Q_f_hat);
525 :
526 : // Resize the output vectors
527 284 : outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
528 852 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
529 584 : outputs_vector[n].resize( get_rb_theta_expansion().get_n_output_terms(n) );
530 :
531 : // Resize the output dual norm vectors
532 284 : output_dual_innerprods.resize(get_rb_theta_expansion().get_n_outputs());
533 852 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
534 : {
535 568 : unsigned int Q_l_hat = get_rb_theta_expansion().get_n_output_terms(n)*(get_rb_theta_expansion().get_n_output_terms(n)+1)/2;
536 584 : output_dual_innerprods[n].resize(Q_l_hat);
537 : }
538 :
539 : {
540 8 : DofMap & dof_map = this->get_dof_map();
541 :
542 284 : dof_map.attach_matrix(*inner_product_matrix);
543 284 : inner_product_matrix->init();
544 284 : inner_product_matrix->zero();
545 :
546 284 : if(store_non_dirichlet_operators)
547 : {
548 : // We also need a non-Dirichlet inner-product matrix
549 0 : non_dirichlet_inner_product_matrix = SparseMatrix<Number>::build(this->comm());
550 0 : dof_map.attach_matrix(*non_dirichlet_inner_product_matrix);
551 0 : non_dirichlet_inner_product_matrix->init();
552 0 : non_dirichlet_inner_product_matrix->zero();
553 : }
554 :
555 1136 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
556 : {
557 : // Initialize the memory for the matrices
558 852 : Aq_vector[q] = SparseMatrix<Number>::build(this->comm());
559 876 : dof_map.attach_matrix(*Aq_vector[q]);
560 876 : Aq_vector[q]->init();
561 876 : Aq_vector[q]->zero();
562 : }
563 :
564 : // We also need to initialize a second set of non-Dirichlet operators
565 284 : if (store_non_dirichlet_operators)
566 : {
567 0 : non_dirichlet_Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
568 0 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
569 : {
570 : // Initialize the memory for the matrices
571 0 : non_dirichlet_Aq_vector[q] = SparseMatrix<Number>::build(this->comm());
572 0 : dof_map.attach_matrix(*non_dirichlet_Aq_vector[q]);
573 0 : non_dirichlet_Aq_vector[q]->init();
574 0 : non_dirichlet_Aq_vector[q]->zero();
575 : }
576 : }
577 : }
578 :
579 : // Initialize the vectors
580 1278 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
581 : {
582 : // Initialize the memory for the vectors
583 994 : Fq_vector[q] = NumericVector<Number>::build(this->comm());
584 1022 : Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
585 : }
586 :
587 : // We also need to initialize a second set of non-Dirichlet operators
588 284 : if (store_non_dirichlet_operators)
589 : {
590 0 : non_dirichlet_Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
591 0 : for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
592 : {
593 : // Initialize the memory for the vectors
594 0 : non_dirichlet_Fq_vector[q] = NumericVector<Number>::build(this->comm());
595 0 : non_dirichlet_Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
596 : }
597 : }
598 :
599 852 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
600 1136 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
601 : {
602 : // Initialize the memory for the truth output vectors
603 568 : outputs_vector[n][q_l] = NumericVector<Number>::build(this->comm());
604 600 : outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
605 : }
606 :
607 284 : if (store_non_dirichlet_operators)
608 : {
609 0 : non_dirichlet_outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
610 0 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
611 : {
612 0 : non_dirichlet_outputs_vector[n].resize( get_rb_theta_expansion().get_n_output_terms(n) );
613 0 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
614 : {
615 : // Initialize the memory for the truth output vectors
616 0 : non_dirichlet_outputs_vector[n][q_l] = NumericVector<Number>::build(this->comm());
617 0 : non_dirichlet_outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
618 : }
619 : }
620 : }
621 :
622 : // Resize truth_outputs vector
623 284 : truth_outputs.resize(this->get_rb_theta_expansion().get_n_outputs());
624 284 : }
625 :
626 2883 : std::unique_ptr<DGFEMContext> RBConstruction::build_context ()
627 : {
628 2883 : return std::make_unique<DGFEMContext>(*this);
629 : }
630 :
631 2883 : void RBConstruction::add_scaled_matrix_and_vector(Number scalar,
632 : ElemAssembly * elem_assembly,
633 : SparseMatrix<Number> * input_matrix,
634 : NumericVector<Number> * input_vector,
635 : bool symmetrize,
636 : bool apply_dof_constraints)
637 : {
638 80 : LOG_SCOPE("add_scaled_matrix_and_vector()", "RBConstruction");
639 :
640 2883 : bool assemble_matrix = (input_matrix != nullptr);
641 80 : bool assemble_vector = (input_vector != nullptr);
642 :
643 2883 : if (!assemble_matrix && !assemble_vector)
644 0 : return;
645 :
646 160 : const MeshBase & mesh = this->get_mesh();
647 :
648 : // First add any node-based terms (e.g. point loads)
649 :
650 : // Make a std::set of all the nodes that are in 1 or more
651 : // nodesets. We only want to call get_nodal_values() once per Node
652 : // per ElemAssembly object, regardless of how many nodesets it
653 : // appears in.
654 160 : std::set<dof_id_type> nodes_with_nodesets;
655 779730 : for (const auto & t : mesh.get_boundary_info().build_node_list())
656 715147 : nodes_with_nodesets.insert(std::get<0>(t));
657 :
658 : // It's possible for the node assembly loop below to throw an
659 : // exception on one (or some subset) of the processors. In that
660 : // case, we stop assembling on the processor(s) that threw and let
661 : // the other processors finish assembly. Then we synchronize to
662 : // check whether an exception was thrown on _any_ processor, and if
663 : // so, re-throw it on _all_ processors.
664 2883 : int nodal_assembly_threw = 0;
665 :
666 : libmesh_try
667 : {
668 :
669 702238 : for (const auto & id : nodes_with_nodesets)
670 : {
671 699355 : const Node & node = mesh.node_ref(id);
672 :
673 : // If node is on this processor, then all dofs on node are too
674 : // so we can do the add below safely
675 699355 : if (node.processor_id() == this->comm().rank())
676 : {
677 : // Get the values to add to the rhs vector
678 54800 : std::vector<dof_id_type> nodal_dof_indices;
679 387900 : DenseMatrix<Number> nodal_matrix;
680 333100 : DenseVector<Number> nodal_rhs;
681 333100 : elem_assembly->get_nodal_values(nodal_dof_indices,
682 : nodal_matrix,
683 : nodal_rhs,
684 : *this,
685 54800 : node);
686 :
687 : // Perform any required user-defined postprocessing on
688 : // the matrix and rhs.
689 : //
690 : // TODO: We need to postprocess node matrices and vectors
691 : // in some cases (e.g. when rotations are applied to
692 : // nodes), but since we don't have a FEMContext at this
693 : // point we would need to have a different interface
694 : // taking the DenseMatrix, DenseVector, and probably the
695 : // current node that we are on...
696 : // this->post_process_elem_matrix_and_vector(nodal_matrix, nodal_rhs);
697 :
698 333100 : if (!nodal_dof_indices.empty())
699 : {
700 0 : if (apply_dof_constraints)
701 : {
702 : // Apply constraints, e.g. Dirichlet and periodic constraints
703 0 : this->get_dof_map().constrain_element_matrix_and_vector(
704 : nodal_matrix,
705 : nodal_rhs,
706 : nodal_dof_indices,
707 : /*asymmetric_constraint_rows*/ false);
708 : }
709 :
710 : // Scale and add to global matrix and/or vector
711 0 : nodal_matrix *= scalar;
712 0 : nodal_rhs *= scalar;
713 :
714 0 : if (assemble_vector)
715 0 : input_vector->add_vector(nodal_rhs, nodal_dof_indices);
716 :
717 0 : if (assemble_matrix)
718 0 : input_matrix->add_matrix(nodal_matrix, nodal_dof_indices);
719 : }
720 278300 : }
721 : }
722 : } // libmesh_try
723 0 : libmesh_catch(...)
724 : {
725 0 : nodal_assembly_threw = 1;
726 0 : }
727 :
728 : // Check for exceptions on any procs and if there is one, re-throw
729 : // it on all procs.
730 2883 : this->comm().max(nodal_assembly_threw);
731 :
732 2883 : if (nodal_assembly_threw)
733 0 : libmesh_error_msg("Error during assembly in RBConstruction::add_scaled_matrix_and_vector()");
734 :
735 2963 : std::unique_ptr<DGFEMContext> c = this->build_context();
736 80 : DGFEMContext & context = cast_ref<DGFEMContext &>(*c);
737 :
738 2883 : this->init_context(context);
739 :
740 : // It's possible for the assembly loop below to throw an exception
741 : // on one (or some subset) of the processors. This can happen when
742 : // e.g. the mesh contains one or more elements with negative
743 : // Jacobian. In that case, we stop assembling on the processor(s)
744 : // that threw and let the other processors finish assembly. Then we
745 : // synchronize to check whether an exception was thrown on _any_
746 : // processor, and if so, re-throw it on _all_ processors. This way,
747 : // we make it easier for callers to handle exceptions in parallel
748 : // during assembly: they can simply assume that the code either
749 : // throws on all procs or on no procs.
750 2883 : int assembly_threw = 0;
751 :
752 : libmesh_try
753 : {
754 :
755 837686 : for (const auto & elem : mesh.active_local_element_ptr_range())
756 : {
757 451375 : const ElemType elemtype = elem->type();
758 :
759 451375 : if(elemtype == NODEELEM)
760 : {
761 : // We assume that we do not perform any assembly directly on
762 : // NodeElems, so we skip the assembly calls.
763 :
764 : // However, in a spline basis with Dirichlet constraints on
765 : // spline nodes, a constrained matrix has to take those
766 : // nodes into account.
767 0 : if (!apply_dof_constraints)
768 0 : continue;
769 : }
770 :
771 : // Subdivision elements need special care:
772 : // - skip ghost elements
773 : // - init special quadrature rule
774 380625 : std::unique_ptr<QBase> qrule;
775 451375 : if (elemtype == TRI3SUBDIVISION)
776 : {
777 0 : const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
778 0 : if (gh_elem->is_ghost())
779 0 : continue ;
780 : // A Gauss quadrature rule for numerical integration.
781 : // For subdivision shell elements, a single Gauss point per
782 : // element is sufficient, hence we use extraorder = 0.
783 0 : const int extraorder = 0;
784 0 : FEBase * elem_fe = nullptr;
785 0 : context.get_element_fe( 0, elem_fe );
786 :
787 0 : qrule = elem_fe->get_fe_type().default_quadrature_rule (2, extraorder);
788 :
789 : // Tell the finite element object to use our quadrature rule.
790 0 : elem_fe->attach_quadrature_rule (qrule.get());
791 : }
792 :
793 451375 : context.pre_fe_reinit(*this, elem);
794 :
795 : // Do nothing in case there are no dof_indices on the current element
796 451375 : if ( context.get_dof_indices().empty() )
797 0 : continue;
798 :
799 451375 : context.elem_fe_reinit();
800 :
801 451375 : if (elemtype != NODEELEM)
802 : {
803 451375 : elem_assembly->interior_assembly(context);
804 :
805 451375 : const unsigned char n_sides = context.get_elem().n_sides();
806 2400875 : for (context.side = 0; context.side != n_sides; ++context.side)
807 : {
808 : // May not need to apply fluxes on non-boundary elements
809 2103000 : if ((context.get_elem().neighbor_ptr(context.get_side()) != nullptr) && !impose_internal_fluxes)
810 1684120 : continue;
811 :
812 : // skip degenerate sides with zero area
813 131440 : if( (context.get_elem().side_ptr(context.get_side())->volume() <= 0.) && skip_degenerate_sides)
814 0 : continue;
815 :
816 121660 : context.side_fe_reinit();
817 121660 : elem_assembly->boundary_assembly(context);
818 :
819 121660 : if (context.dg_terms_are_active())
820 : {
821 0 : input_matrix->add_matrix (context.get_elem_elem_jacobian(),
822 0 : context.get_dof_indices(),
823 0 : context.get_dof_indices());
824 :
825 0 : input_matrix->add_matrix (context.get_elem_neighbor_jacobian(),
826 0 : context.get_dof_indices(),
827 0 : context.get_neighbor_dof_indices());
828 :
829 0 : input_matrix->add_matrix (context.get_neighbor_elem_jacobian(),
830 : context.get_neighbor_dof_indices(),
831 0 : context.get_dof_indices());
832 :
833 0 : input_matrix->add_matrix (context.get_neighbor_neighbor_jacobian(),
834 : context.get_neighbor_dof_indices(),
835 0 : context.get_neighbor_dof_indices());
836 : }
837 : }
838 : }
839 :
840 : // Do any required user post-processing before symmetrizing and/or applying
841 : // constraints.
842 : //
843 : // We only do this if apply_dof_constraints is true because we want to be
844 : // able to set apply_dof_constraints=false in order to obtain a matrix
845 : // A with no dof constraints or dof transformations, as opposed to C^T A C,
846 : // which includes constraints and/or dof transformations. Here C refers to
847 : // the matrix that imposes dof constraints and transformations on the
848 : // solution u.
849 : //
850 : // Matrices such as A are what we store in our "non_dirichlet" operators, and
851 : // they are useful for computing terms such as (C u_i)^T A (C u_j) (e.g. see
852 : // update_RB_system_matrices()), where C u is the result of a "truth_solve",
853 : // which includes calls to both enforce_constraints_exactly() and
854 : // post_process_truth_solution(). If we use C^T A C to compute these terms then
855 : // we would "double apply" the matrix C, which can give incorrect results.
856 451375 : if (apply_dof_constraints)
857 451375 : this->post_process_elem_matrix_and_vector(context);
858 :
859 : // Need to symmetrize before imposing
860 : // periodic constraints
861 451375 : if (assemble_matrix && symmetrize)
862 : {
863 28125 : DenseMatrix<Number> Ke_transpose;
864 28125 : context.get_elem_jacobian().get_transpose(Ke_transpose);
865 28125 : context.get_elem_jacobian() += Ke_transpose;
866 0 : context.get_elem_jacobian() *= 0.5;
867 28125 : }
868 :
869 : // As discussed above, we can set apply_dof_constraints=false to
870 : // get A instead of C^T A C
871 451375 : if (apply_dof_constraints)
872 : {
873 : // Apply constraints, e.g. Dirichlet and periodic constraints
874 35375 : this->get_dof_map().constrain_element_matrix_and_vector
875 451375 : (context.get_elem_jacobian(),
876 : context.get_elem_residual(),
877 : context.get_dof_indices(),
878 : /*asymmetric_constraint_rows*/ false );
879 : }
880 :
881 : // Scale and add to global matrix and/or vector
882 63875 : context.get_elem_jacobian() *= scalar;
883 63875 : context.get_elem_residual() *= scalar;
884 :
885 451375 : if (assemble_matrix)
886 : {
887 :
888 220675 : CouplingMatrix * coupling_matrix = get_dof_map()._dof_coupling;
889 220675 : if (!coupling_matrix)
890 : {
891 : // If we haven't defined a _dof_coupling matrix then just add
892 : // the whole matrix
893 220675 : input_matrix->add_matrix (context.get_elem_jacobian(),
894 32300 : context.get_dof_indices() );
895 : }
896 : else
897 : {
898 : // Otherwise we should only add the relevant submatrices
899 0 : for (unsigned int var1=0; var1<n_vars(); var1++)
900 : {
901 0 : ConstCouplingRow ccr(var1, *coupling_matrix);
902 0 : for (const auto & var2 : ccr)
903 : {
904 0 : unsigned int sub_m = context.get_elem_jacobian( var1, var2 ).m();
905 0 : unsigned int sub_n = context.get_elem_jacobian( var1, var2 ).n();
906 0 : DenseMatrix<Number> sub_jac(sub_m, sub_n);
907 0 : for (unsigned int row=0; row<sub_m; row++)
908 0 : for (unsigned int col=0; col<sub_n; col++)
909 : {
910 0 : sub_jac(row,col) = context.get_elem_jacobian( var1, var2 ).el(row,col);
911 : }
912 0 : input_matrix->add_matrix (sub_jac,
913 0 : context.get_dof_indices(var1),
914 0 : context.get_dof_indices(var2) );
915 0 : }
916 : }
917 : }
918 :
919 : }
920 :
921 451375 : if (assemble_vector)
922 211475 : input_vector->add_vector (context.get_elem_residual(),
923 19225 : context.get_dof_indices() );
924 383348 : } // end for (elem)
925 :
926 : } // libmesh_try
927 0 : libmesh_catch(...)
928 : {
929 0 : assembly_threw = 1;
930 0 : }
931 :
932 : // Note: regardless of whether any procs threw during assembly (and
933 : // thus didn't finish assembling), we should not leave the matrix
934 : // and vector in an inconsistent state, since it may be possible to
935 : // recover from the exception. Therefore, we close them now. The
936 : // assumption here is that the nature of the exception does not
937 : // prevent the matrix and vector from still being assembled (albeit
938 : // with incomplete data).
939 2883 : if (assemble_matrix)
940 1321 : input_matrix->close();
941 2883 : if (assemble_vector)
942 1562 : input_vector->close();
943 :
944 : // Check for exceptions on any procs and if there is one, re-throw
945 : // it on all procs.
946 2883 : this->comm().max(assembly_threw);
947 :
948 2883 : if (assembly_threw)
949 0 : libmesh_error_msg("Error during assembly in RBConstruction::add_scaled_matrix_and_vector()");
950 2723 : }
951 :
952 0 : void RBConstruction::set_context_solution_vec(NumericVector<Number> & vec)
953 : {
954 : // Set current_local_solution = vec so that we can access
955 : // vec from DGFEMContext during assembly
956 : vec.localize
957 0 : (*current_local_solution, this->get_dof_map().get_send_list());
958 0 : }
959 :
960 3556 : void RBConstruction::truth_assembly()
961 : {
962 200 : LOG_SCOPE("truth_assembly()", "RBConstruction");
963 :
964 3556 : const RBParameters & mu = get_parameters();
965 :
966 3556 : this->matrix->zero();
967 3556 : this->rhs->zero();
968 :
969 3556 : this->matrix->close();
970 3556 : this->rhs->close();
971 :
972 : {
973 : // We should have already assembled the matrices
974 : // and vectors in the affine expansion, so
975 : // just use them
976 :
977 14224 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
978 : {
979 10668 : matrix->add(get_rb_theta_expansion().eval_A_theta(q_a, mu), *get_Aq(q_a));
980 : }
981 :
982 3656 : std::unique_ptr<NumericVector<Number>> temp_vec = NumericVector<Number>::build(this->comm());
983 3556 : temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
984 17762 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
985 : {
986 14206 : *temp_vec = *get_Fq(q_f);
987 14206 : temp_vec->scale( get_rb_theta_expansion().eval_F_theta(q_f, mu) );
988 14606 : rhs->add(*temp_vec);
989 : }
990 3356 : }
991 :
992 3556 : this->matrix->close();
993 3556 : this->rhs->close();
994 3556 : }
995 :
996 284 : void RBConstruction::assemble_inner_product_matrix(SparseMatrix<Number> * input_matrix,
997 : bool apply_dof_constraints)
998 : {
999 284 : input_matrix->zero();
1000 :
1001 284 : if(!use_energy_inner_product)
1002 : {
1003 284 : add_scaled_matrix_and_vector(1.,
1004 : inner_product_assembly,
1005 : input_matrix,
1006 : nullptr,
1007 : false, /* symmetrize */
1008 : apply_dof_constraints);
1009 : }
1010 : else
1011 : {
1012 0 : libmesh_error_msg_if(energy_inner_product_coeffs.size() != get_rb_theta_expansion().get_n_A_terms(),
1013 : "Error: invalid number of entries in energy_inner_product_coeffs.");
1014 :
1015 : // We symmetrize below so that we may use the energy inner-product even in cases
1016 : // where the A_q are not symmetric.
1017 0 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1018 : {
1019 0 : add_scaled_matrix_and_vector(energy_inner_product_coeffs[q_a],
1020 0 : &rb_assembly_expansion->get_A_assembly(q_a),
1021 : input_matrix,
1022 : nullptr,
1023 : true, /* symmetrize */
1024 : apply_dof_constraints);
1025 : }
1026 : }
1027 284 : }
1028 :
1029 852 : void RBConstruction::assemble_Aq_matrix(unsigned int q,
1030 : SparseMatrix<Number> * input_matrix,
1031 : bool apply_dof_constraints)
1032 : {
1033 852 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
1034 : "Error: We must have q < Q_a in assemble_Aq_matrix.");
1035 :
1036 852 : input_matrix->zero();
1037 :
1038 876 : add_scaled_matrix_and_vector(1.,
1039 852 : &rb_assembly_expansion->get_A_assembly(q),
1040 : input_matrix,
1041 : nullptr,
1042 : false, /* symmetrize */
1043 : apply_dof_constraints);
1044 852 : }
1045 :
1046 45 : void RBConstruction::add_scaled_Aq(Number scalar,
1047 : unsigned int q_a,
1048 : SparseMatrix<Number> * input_matrix,
1049 : bool symmetrize)
1050 : {
1051 0 : LOG_SCOPE("add_scaled_Aq()", "RBConstruction");
1052 :
1053 45 : libmesh_error_msg_if(q_a >= get_rb_theta_expansion().get_n_A_terms(),
1054 : "Error: We must have q < Q_a in add_scaled_Aq.");
1055 :
1056 45 : if (!symmetrize)
1057 : {
1058 0 : input_matrix->add(scalar, *get_Aq(q_a));
1059 0 : input_matrix->close();
1060 : }
1061 : else
1062 : {
1063 45 : add_scaled_matrix_and_vector(scalar,
1064 45 : &rb_assembly_expansion->get_A_assembly(q_a),
1065 : input_matrix,
1066 : nullptr,
1067 : symmetrize);
1068 : }
1069 45 : }
1070 :
1071 284 : void RBConstruction::assemble_misc_matrices()
1072 : {
1073 8 : libMesh::out << "Assembling inner product matrix" << std::endl;
1074 284 : assemble_inner_product_matrix(inner_product_matrix.get());
1075 :
1076 284 : if (store_non_dirichlet_operators)
1077 : {
1078 0 : libMesh::out << "Assembling non-Dirichlet inner product matrix" << std::endl;
1079 0 : assemble_inner_product_matrix(non_dirichlet_inner_product_matrix.get(), false);
1080 : }
1081 284 : }
1082 :
1083 284 : void RBConstruction::assemble_all_affine_operators()
1084 : {
1085 1136 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1086 : {
1087 852 : libMesh::out << "Assembling affine operator " << (q_a+1) << " of "
1088 852 : << get_rb_theta_expansion().get_n_A_terms() << std::endl;
1089 852 : assemble_Aq_matrix(q_a, get_Aq(q_a));
1090 : }
1091 :
1092 284 : if (store_non_dirichlet_operators)
1093 : {
1094 0 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1095 : {
1096 0 : libMesh::out << "Assembling non-Dirichlet affine operator " << (q_a+1) << " of "
1097 0 : << get_rb_theta_expansion().get_n_A_terms() << std::endl;
1098 0 : assemble_Aq_matrix(q_a, get_non_dirichlet_Aq(q_a), false);
1099 : }
1100 : }
1101 284 : }
1102 :
1103 284 : void RBConstruction::assemble_all_affine_vectors()
1104 : {
1105 1278 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1106 : {
1107 994 : libMesh::out << "Assembling affine vector " << (q_f+1) << " of "
1108 994 : << get_rb_theta_expansion().get_n_F_terms() << std::endl;
1109 994 : assemble_Fq_vector(q_f, get_Fq(q_f));
1110 : }
1111 :
1112 284 : if (store_non_dirichlet_operators)
1113 : {
1114 0 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1115 : {
1116 0 : libMesh::out << "Assembling non-Dirichlet affine vector " << (q_f+1) << " of "
1117 0 : << get_rb_theta_expansion().get_n_F_terms() << std::endl;
1118 0 : assemble_Fq_vector(q_f, get_non_dirichlet_Fq(q_f), false);
1119 : }
1120 : }
1121 :
1122 284 : }
1123 :
1124 994 : void RBConstruction::assemble_Fq_vector(unsigned int q,
1125 : NumericVector<Number> * input_vector,
1126 : bool apply_dof_constraints)
1127 : {
1128 994 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_F_terms(),
1129 : "Error: We must have q < Q_f in assemble_Fq_vector.");
1130 :
1131 994 : input_vector->zero();
1132 :
1133 1022 : add_scaled_matrix_and_vector(1.,
1134 994 : &rb_assembly_expansion->get_F_assembly(q),
1135 : nullptr,
1136 : input_vector,
1137 : false, /* symmetrize */
1138 : apply_dof_constraints /* apply_dof_constraints */);
1139 994 : }
1140 :
1141 284 : void RBConstruction::assemble_all_output_vectors()
1142 : {
1143 852 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1144 1136 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1145 : {
1146 1120 : libMesh::out << "Assembling output vector, (" << (n+1) << "," << (q_l+1)
1147 568 : << ") of (" << get_rb_theta_expansion().get_n_outputs()
1148 568 : << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1149 16 : << std::endl;
1150 568 : get_output_vector(n, q_l)->zero();
1151 568 : add_scaled_matrix_and_vector(1., &rb_assembly_expansion->get_output_assembly(n,q_l),
1152 : nullptr,
1153 : get_output_vector(n,q_l),
1154 : false, /* symmetrize */
1155 : true /* apply_dof_constraints */);
1156 : }
1157 :
1158 284 : if (store_non_dirichlet_operators)
1159 : {
1160 0 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1161 0 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1162 : {
1163 0 : libMesh::out << "Assembling non-Dirichlet output vector, (" << (n+1) << "," << (q_l+1)
1164 0 : << ") of (" << get_rb_theta_expansion().get_n_outputs()
1165 0 : << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1166 0 : << std::endl;
1167 0 : get_non_dirichlet_output_vector(n, q_l)->zero();
1168 0 : add_scaled_matrix_and_vector(1., &rb_assembly_expansion->get_output_assembly(n,q_l),
1169 : nullptr,
1170 : get_non_dirichlet_output_vector(n,q_l),
1171 : false, /* symmetrize */
1172 : false /* apply_dof_constraints */);
1173 : }
1174 : }
1175 284 : }
1176 :
1177 284 : Real RBConstruction::train_reduced_basis(const bool resize_rb_eval_data)
1178 : {
1179 560 : if(get_RB_training_type() == "Greedy")
1180 : {
1181 284 : return train_reduced_basis_with_greedy(resize_rb_eval_data);
1182 : }
1183 0 : else if (get_RB_training_type() == "POD")
1184 : {
1185 0 : train_reduced_basis_with_POD();
1186 0 : return 0.;
1187 : }
1188 : else
1189 : {
1190 0 : libmesh_error_msg("RB training type not recognized: " + get_RB_training_type());
1191 : }
1192 :
1193 : return 0.;
1194 : }
1195 :
1196 284 : Real RBConstruction::train_reduced_basis_with_greedy(const bool resize_rb_eval_data)
1197 : {
1198 16 : LOG_SCOPE("train_reduced_basis_with_greedy()", "RBConstruction");
1199 :
1200 8 : int count = 0;
1201 :
1202 284 : RBEvaluation & rbe = get_rb_evaluation();
1203 :
1204 : // initialize rbe's parameters
1205 284 : rbe.initialize_parameters(*this);
1206 :
1207 : // possibly resize data structures according to Nmax
1208 284 : if (resize_rb_eval_data)
1209 284 : rbe.resize_data_structures(get_Nmax());
1210 :
1211 : // Clear the Greedy param list
1212 284 : for (auto & plist : rbe.greedy_param_list)
1213 0 : plist.clear();
1214 :
1215 276 : rbe.greedy_param_list.clear();
1216 :
1217 8 : Real training_greedy_error = 0.;
1218 :
1219 :
1220 : // If we are continuing from a previous training run,
1221 : // we might already be at the max number of basis functions.
1222 : // If so, we can just return.
1223 284 : if (rbe.get_n_basis_functions() >= get_Nmax())
1224 : {
1225 0 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1226 0 : << get_Nmax() << std::endl;
1227 0 : return 0.;
1228 : }
1229 :
1230 : // Optionally pre-evaluate the theta functions on the entire (local) training parameter set.
1231 284 : if (get_preevaluate_thetas_flag())
1232 0 : preevaluate_thetas();
1233 :
1234 284 : if(!skip_residual_in_train_reduced_basis)
1235 : {
1236 : // Compute the dual norms of the outputs if we haven't already done so.
1237 284 : compute_output_dual_innerprods();
1238 :
1239 : // Compute the Fq Riesz representor dual norms if we haven't already done so.
1240 284 : compute_Fq_representor_innerprods();
1241 : }
1242 :
1243 8 : libMesh::out << std::endl << "---- Performing Greedy basis enrichment ----" << std::endl;
1244 8 : Real initial_greedy_error = 0.;
1245 8 : bool initial_greedy_error_initialized = false;
1246 : while (true)
1247 : {
1248 140 : libMesh::out << std::endl << "---- Basis dimension: "
1249 4957 : << rbe.get_n_basis_functions() << " ----" << std::endl;
1250 :
1251 4957 : if (count > 0 || (count==0 && use_empty_rb_solve_in_greedy))
1252 : {
1253 140 : libMesh::out << "Performing RB solves on training set" << std::endl;
1254 4957 : training_greedy_error = compute_max_error_bound();
1255 :
1256 140 : libMesh::out << "Maximum error bound is " << training_greedy_error << std::endl << std::endl;
1257 :
1258 : // record the initial error
1259 4957 : if (!initial_greedy_error_initialized)
1260 : {
1261 8 : initial_greedy_error = training_greedy_error;
1262 8 : initial_greedy_error_initialized = true;
1263 : }
1264 :
1265 : // Break out of training phase if we have reached Nmax
1266 : // or if the training_tolerance is satisfied.
1267 4957 : if (greedy_termination_test(training_greedy_error, initial_greedy_error, count))
1268 0 : break;
1269 : }
1270 :
1271 140 : libMesh::out << "Performing truth solve at parameter:" << std::endl;
1272 4956 : print_parameters();
1273 :
1274 : // Update the list of Greedily selected parameters
1275 4956 : this->update_greedy_param_list();
1276 :
1277 : // Perform an Offline truth solve for the current parameter
1278 4956 : truth_solve(-1);
1279 :
1280 4956 : if (check_if_zero_truth_solve())
1281 : {
1282 0 : libMesh::out << "Zero basis function encountered hence ending basis enrichment" << std::endl;
1283 0 : break;
1284 : }
1285 :
1286 : // Add orthogonal part of the snapshot to the RB space
1287 140 : libMesh::out << "Enriching the RB space" << std::endl;
1288 4956 : enrich_RB_space();
1289 :
1290 4956 : update_system();
1291 :
1292 : // Check if we've reached Nmax now. We do this before calling
1293 : // update_residual_terms() since we can skip that step if we've
1294 : // already reached Nmax.
1295 4956 : if (rbe.get_n_basis_functions() >= this->get_Nmax())
1296 : {
1297 8 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1298 16 : << get_Nmax() << std::endl;
1299 8 : break;
1300 : }
1301 :
1302 4673 : if(!skip_residual_in_train_reduced_basis)
1303 : {
1304 4673 : update_residual_terms();
1305 : }
1306 :
1307 : // Increment counter
1308 4673 : count++;
1309 : }
1310 284 : this->update_greedy_param_list();
1311 :
1312 8 : return training_greedy_error;
1313 : }
1314 :
1315 0 : void RBConstruction::enrich_basis_from_rhs_terms(const bool resize_rb_eval_data)
1316 : {
1317 0 : LOG_SCOPE("enrich_basis_from_rhs_terms()", "RBConstruction");
1318 :
1319 : // initialize rb_eval's parameters
1320 0 : get_rb_evaluation().initialize_parameters(*this);
1321 :
1322 : // possibly resize data structures according to Nmax
1323 0 : if (resize_rb_eval_data)
1324 : {
1325 0 : get_rb_evaluation().resize_data_structures(get_Nmax());
1326 : }
1327 :
1328 0 : libMesh::out << std::endl << "---- Enriching basis from rhs terms ----" << std::endl;
1329 :
1330 0 : truth_assembly();
1331 :
1332 0 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1333 : {
1334 0 : libMesh::out << std::endl << "Performing truth solve with rhs from rhs term " << q_f << std::endl;
1335 :
1336 0 : *rhs = *get_Fq(q_f);
1337 :
1338 0 : if (rhs->l2_norm() == 0)
1339 : {
1340 : // Skip enrichment if the rhs is zero
1341 0 : continue;
1342 : }
1343 :
1344 : // truth_assembly assembles into matrix and rhs, so use those for the solve
1345 0 : if (extra_linear_solver)
1346 : {
1347 : // If extra_linear_solver has been initialized, then we use it for the
1348 : // truth solves.
1349 0 : solve_for_matrix_and_rhs(*extra_linear_solver, *matrix, *rhs);
1350 :
1351 0 : if (assert_convergence)
1352 0 : check_convergence(*extra_linear_solver);
1353 : }
1354 : else
1355 : {
1356 0 : solve_for_matrix_and_rhs(*get_linear_solver(), *matrix, *rhs);
1357 :
1358 0 : if (assert_convergence)
1359 0 : check_convergence(*get_linear_solver());
1360 : }
1361 :
1362 : // Debugging: enable this code to print the rhs that was used in
1363 : // the most recent truth solve to a uniquely-named file.
1364 : #if 0
1365 : {
1366 : char temp_file[] = "truth_rhs_XXXXXX.dat";
1367 : int fd = mkstemps(temp_file, 4);
1368 : if (fd != -1)
1369 : {
1370 : libMesh::out << "Writing truth system rhs to file: " << temp_file << std::endl;
1371 : rhs->print_matlab(std::string(temp_file));
1372 : }
1373 : }
1374 : #endif // 0
1375 :
1376 : // Debugging: enable this code to print the most recent truth
1377 : // solution to a uniquely-named file.
1378 : #ifdef LIBMESH_HAVE_EXODUS_API
1379 : #if 0
1380 : {
1381 : // Note: mkstemps creates a file and returns an open file descriptor to it.
1382 : // The filename is created from a template which must have 6 'X' characters followed
1383 : // by a suffix having the specified length (in this case 4, for ".exo").
1384 : char temp_file[] = "truth_XXXXXX.exo";
1385 : int fd = mkstemps(temp_file, 4);
1386 : if (fd != -1)
1387 : {
1388 : libMesh::out << "Writing truth solution to file: " << temp_file << std::endl;
1389 : ExodusII_IO exo_io(this->get_mesh());
1390 : std::set<std::string> system_names = {this->name()};
1391 : exo_io.write_equation_systems(std::string(temp_file),
1392 : this->get_equation_systems(), &system_names);
1393 : }
1394 : }
1395 : #endif // 0
1396 : #endif // LIBMESH_HAVE_EXODUS_API
1397 :
1398 : // Call user-defined post-processing routines on the truth solution.
1399 0 : post_process_truth_solution();
1400 :
1401 : // Add orthogonal part of the snapshot to the RB space
1402 0 : libMesh::out << "Enriching the RB space" << std::endl;
1403 0 : enrich_RB_space();
1404 :
1405 0 : update_system();
1406 : }
1407 0 : }
1408 :
1409 0 : void RBConstruction::train_reduced_basis_with_POD()
1410 : {
1411 : // We need to use the same training set on all processes so that
1412 : // the truth solves below work correctly in parallel.
1413 0 : libmesh_error_msg_if(!serial_training_set, "We must use a serial training set with POD");
1414 0 : libmesh_error_msg_if(get_rb_evaluation().get_n_basis_functions() > 0, "Basis should not already be initialized");
1415 :
1416 0 : get_rb_evaluation().initialize_parameters(*this);
1417 0 : get_rb_evaluation().resize_data_structures(get_Nmax());
1418 :
1419 : // Storage for the POD snapshots
1420 0 : unsigned int n_snapshots = get_n_training_samples();
1421 :
1422 0 : if (get_n_params() == 0)
1423 : {
1424 : // In this case we should have generated an empty training set
1425 : // so assert this
1426 0 : libmesh_assert(n_snapshots == 0);
1427 :
1428 : // If we have no parameters, then we should do exactly one "truth solve"
1429 0 : n_snapshots = 1;
1430 : }
1431 :
1432 0 : std::vector<std::unique_ptr<NumericVector<Number>>> POD_snapshots(n_snapshots);
1433 0 : for (unsigned int i=0; i<n_snapshots; i++)
1434 : {
1435 0 : POD_snapshots[i] = NumericVector<Number>::build(this->comm());
1436 0 : POD_snapshots[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1437 : }
1438 :
1439 : // We use the same training set on all processes
1440 0 : libMesh::out << std::endl;
1441 0 : for (unsigned int i=0; i<n_snapshots; i++)
1442 : {
1443 0 : if (get_n_params() > 0)
1444 : {
1445 0 : set_params_from_training_set(i);
1446 : }
1447 :
1448 0 : libMesh::out << "Truth solve " << (i+1) << " of " << n_snapshots << std::endl;
1449 :
1450 0 : truth_solve(-1);
1451 :
1452 0 : *POD_snapshots[i] = *solution;
1453 : }
1454 0 : libMesh::out << std::endl;
1455 :
1456 0 : if (_normalize_solution_snapshots)
1457 : {
1458 0 : libMesh::out << "Normalizing solution snapshots" << std::endl;
1459 0 : for (unsigned int i=0; i<n_snapshots; i++)
1460 : {
1461 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(
1462 0 : *inner_product_storage_vector, *POD_snapshots[i]);
1463 0 : Real norm = std::sqrt(std::real(POD_snapshots[i]->dot(*inner_product_storage_vector)));
1464 :
1465 0 : if (norm > 0.)
1466 0 : POD_snapshots[i]->scale(1./norm);
1467 : }
1468 : }
1469 :
1470 : // Set up the "correlation matrix"
1471 0 : DenseMatrix<Number> correlation_matrix(n_snapshots,n_snapshots);
1472 0 : for (unsigned int i=0; i<n_snapshots; i++)
1473 : {
1474 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(
1475 0 : *inner_product_storage_vector, *POD_snapshots[i]);
1476 :
1477 0 : for (unsigned int j=0; j<=i; j++)
1478 : {
1479 0 : Number inner_prod = (POD_snapshots[j]->dot(*inner_product_storage_vector));
1480 :
1481 0 : correlation_matrix(i,j) = inner_prod;
1482 0 : if(i != j)
1483 : {
1484 0 : correlation_matrix(j,i) = libmesh_conj(inner_prod);
1485 : }
1486 : }
1487 : }
1488 :
1489 : // compute SVD of correlation matrix
1490 0 : DenseVector<Real> sigma( n_snapshots );
1491 0 : DenseMatrix<Number> U( n_snapshots, n_snapshots );
1492 0 : DenseMatrix<Number> VT( n_snapshots, n_snapshots );
1493 0 : correlation_matrix.svd(sigma, U, VT );
1494 :
1495 0 : if (sigma(0) == 0.)
1496 0 : return;
1497 :
1498 : // Add dominant vectors from the POD as basis functions.
1499 0 : unsigned int j = 0;
1500 : while (true)
1501 : {
1502 0 : if (j >= get_Nmax() || j >= n_snapshots)
1503 : {
1504 0 : libMesh::out << "Maximum number of basis functions (" << j << ") reached." << std::endl;
1505 0 : break;
1506 : }
1507 :
1508 : // The "energy" error in the POD approximation is determined by the first omitted
1509 : // singular value, i.e. sigma(j). We normalize by sigma(0), which gives the total
1510 : // "energy", in order to obtain a relative error.
1511 0 : const Real rel_err = std::sqrt(sigma(j)) / std::sqrt(sigma(0));
1512 :
1513 0 : libMesh::out << "Number of basis functions: " << j
1514 0 : << ", POD error norm: " << rel_err << std::endl;
1515 :
1516 0 : if (rel_err < this->rel_training_tolerance)
1517 : {
1518 0 : libMesh::out << "Training tolerance reached." << std::endl;
1519 0 : break;
1520 : }
1521 :
1522 0 : std::unique_ptr< NumericVector<Number> > v = POD_snapshots[j]->zero_clone();
1523 0 : for ( unsigned int i=0; i<n_snapshots; ++i )
1524 : {
1525 0 : v->add( U.el(i, j), *POD_snapshots[i] );
1526 : }
1527 :
1528 0 : Real norm_v = std::sqrt(sigma(j));
1529 0 : v->scale( 1./norm_v );
1530 :
1531 0 : get_rb_evaluation().basis_functions.emplace_back( std::move(v) );
1532 :
1533 0 : j++;
1534 0 : }
1535 0 : libMesh::out << std::endl;
1536 :
1537 0 : this->delta_N = get_rb_evaluation().get_n_basis_functions();
1538 0 : update_system();
1539 :
1540 : // We now compute all terms required to evaluate the RB error indicator.
1541 : // Unlike in the case of the RB Greedy algorithm, for the POD approach
1542 : // we do not need this data in order to compute the basis. However, we
1543 : // do need this data in order to evaluate error indicator quantities in
1544 : // the Online stage, so we compute it now so that it can be saved in
1545 : // the training data.
1546 0 : recompute_all_residual_terms();
1547 0 : }
1548 :
1549 4957 : bool RBConstruction::greedy_termination_test(Real abs_greedy_error,
1550 : Real initial_error,
1551 : int)
1552 : {
1553 4957 : if (abs_greedy_error < this->abs_training_tolerance)
1554 : {
1555 0 : libMesh::out << "Absolute error tolerance reached." << std::endl;
1556 0 : return true;
1557 : }
1558 :
1559 4957 : Real rel_greedy_error = abs_greedy_error/initial_error;
1560 4957 : if (rel_greedy_error < this->rel_training_tolerance)
1561 : {
1562 0 : libMesh::out << "Relative error tolerance reached." << std::endl;
1563 1 : return true;
1564 : }
1565 :
1566 4956 : RBEvaluation & rbe = get_rb_evaluation();
1567 :
1568 4956 : if (rbe.get_n_basis_functions() >= this->get_Nmax())
1569 : {
1570 0 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1571 0 : << get_Nmax() << std::endl;
1572 0 : return true;
1573 : }
1574 :
1575 4956 : if (exit_on_repeated_greedy_parameters)
1576 31971 : for (auto & plist : rbe.greedy_param_list)
1577 28415 : if (plist == get_parameters())
1578 : {
1579 0 : libMesh::out << "Exiting greedy because the same parameters were selected twice" << std::endl;
1580 0 : return true;
1581 : }
1582 :
1583 140 : return false;
1584 : }
1585 :
1586 5240 : void RBConstruction::update_greedy_param_list()
1587 : {
1588 5240 : get_rb_evaluation().greedy_param_list.push_back( get_parameters() );
1589 5240 : }
1590 :
1591 0 : const RBParameters & RBConstruction::get_greedy_parameter(unsigned int i)
1592 : {
1593 0 : libmesh_error_msg_if(i >= get_rb_evaluation().greedy_param_list.size(),
1594 : "Error: Argument in RBConstruction::get_greedy_parameter is too large.");
1595 :
1596 0 : return get_rb_evaluation().greedy_param_list[i];
1597 : }
1598 :
1599 3556 : Real RBConstruction::truth_solve(int plot_solution)
1600 : {
1601 200 : LOG_SCOPE("truth_solve()", "RBConstruction");
1602 :
1603 3556 : if(store_untransformed_basis && !_untransformed_solution)
1604 : {
1605 0 : _untransformed_solution = solution->zero_clone();
1606 : }
1607 :
1608 3556 : truth_assembly();
1609 :
1610 : // truth_assembly assembles into matrix and rhs, so use those for the solve
1611 3556 : if (extra_linear_solver)
1612 : {
1613 : // If extra_linear_solver has been initialized, then we use it for the
1614 : // truth solves.
1615 0 : solve_for_matrix_and_rhs(*extra_linear_solver, *matrix, *rhs);
1616 :
1617 0 : if (assert_convergence)
1618 0 : check_convergence(*extra_linear_solver);
1619 : }
1620 : else
1621 : {
1622 3556 : solve_for_matrix_and_rhs(*get_linear_solver(), *matrix, *rhs);
1623 :
1624 3556 : if (assert_convergence)
1625 3556 : check_convergence(*get_linear_solver());
1626 : }
1627 :
1628 3556 : if (store_untransformed_basis)
1629 : {
1630 0 : *_untransformed_solution = *solution;
1631 : }
1632 :
1633 : // Call user-defined post-processing routines on the truth solution.
1634 3556 : post_process_truth_solution();
1635 :
1636 3556 : const RBParameters & mu = get_parameters();
1637 :
1638 9260 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1639 : {
1640 5864 : truth_outputs[n] = 0.;
1641 11408 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1642 11488 : truth_outputs[n] += get_rb_theta_expansion().eval_output_theta(n, q_l, mu)*
1643 5784 : get_output_vector(n,q_l)->dot(*solution);
1644 : }
1645 :
1646 : #ifdef LIBMESH_HAVE_EXODUS_API
1647 3556 : if (plot_solution > 0)
1648 : {
1649 0 : ExodusII_IO exo_io(this->get_mesh());
1650 0 : std::set<std::string> system_names = {this->name()};
1651 0 : exo_io.write_equation_systems("truth.exo", this->get_equation_systems(), &system_names);
1652 0 : }
1653 : #else
1654 : libmesh_ignore(plot_solution);
1655 : #endif
1656 :
1657 : // Get the X norm of the truth solution
1658 : // Useful for normalizing our true error data
1659 3556 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *solution);
1660 3656 : Number truth_X_norm = std::sqrt(inner_product_storage_vector->dot(*solution));
1661 :
1662 3656 : return libmesh_real(truth_X_norm);
1663 0 : }
1664 :
1665 284 : bool RBConstruction::is_serial_training_type(const std::string & RB_training_type_in)
1666 : {
1667 284 : return (RB_training_type_in == "POD");
1668 : }
1669 :
1670 284 : void RBConstruction::set_RB_training_type(const std::string & RB_training_type_in)
1671 : {
1672 284 : this->RB_training_type = RB_training_type_in;
1673 :
1674 284 : if(is_serial_training_type(RB_training_type_in))
1675 : {
1676 : // We need to use a serial training set (so that the training
1677 : // set is the same on all processes) if we're using POD
1678 0 : this->serial_training_set = true;
1679 : }
1680 284 : }
1681 :
1682 638 : const std::string & RBConstruction::get_RB_training_type() const
1683 : {
1684 638 : return RB_training_type;
1685 : }
1686 :
1687 284 : void RBConstruction::set_Nmax(unsigned int Nmax_in)
1688 : {
1689 284 : this->Nmax = Nmax_in;
1690 284 : }
1691 :
1692 142 : void RBConstruction::load_basis_function(unsigned int i)
1693 : {
1694 8 : LOG_SCOPE("load_basis_function()", "RBConstruction");
1695 :
1696 4 : libmesh_assert_less (i, get_rb_evaluation().get_n_basis_functions());
1697 :
1698 142 : *solution = get_rb_evaluation().get_basis_function(i);
1699 :
1700 142 : this->update();
1701 142 : }
1702 :
1703 3556 : void RBConstruction::enrich_RB_space()
1704 : {
1705 200 : LOG_SCOPE("enrich_RB_space()", "RBConstruction");
1706 :
1707 3656 : auto new_bf = solution->clone();
1708 :
1709 3556 : std::unique_ptr<NumericVector<Number>> new_untransformed_bf;
1710 3556 : if (store_untransformed_basis)
1711 : {
1712 0 : new_untransformed_bf = _untransformed_solution->clone();
1713 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
1714 : }
1715 :
1716 31971 : for (unsigned int index=0; index<get_rb_evaluation().get_n_basis_functions(); index++)
1717 : {
1718 28415 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *new_bf);
1719 :
1720 : Number scalar =
1721 28415 : inner_product_storage_vector->dot(get_rb_evaluation().get_basis_function(index));
1722 28415 : new_bf->add(-scalar, get_rb_evaluation().get_basis_function(index));
1723 :
1724 28415 : if (store_untransformed_basis)
1725 : {
1726 0 : new_untransformed_bf->add(-scalar, *_untransformed_basis_functions[index]);
1727 : }
1728 : }
1729 :
1730 : // Normalize new_bf
1731 3556 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *new_bf);
1732 3656 : Number new_bf_norm = std::sqrt( inner_product_storage_vector->dot(*new_bf) );
1733 :
1734 3506 : if (new_bf_norm == 0.)
1735 : {
1736 0 : new_bf->zero(); // avoid potential nan's
1737 :
1738 0 : if (store_untransformed_basis)
1739 : {
1740 0 : new_untransformed_bf->zero();
1741 : }
1742 : }
1743 : else
1744 : {
1745 3556 : new_bf->scale(1./new_bf_norm);
1746 :
1747 3556 : if (store_untransformed_basis)
1748 : {
1749 0 : new_untransformed_bf->scale(1./new_bf_norm);
1750 : }
1751 : }
1752 :
1753 : // load the new basis function into the basis_functions vector.
1754 3556 : get_rb_evaluation().basis_functions.emplace_back( std::move(new_bf) );
1755 :
1756 3556 : if (store_untransformed_basis)
1757 : {
1758 0 : _untransformed_basis_functions.emplace_back( std::move(new_untransformed_bf) );
1759 : }
1760 3556 : }
1761 :
1762 4956 : void RBConstruction::update_system()
1763 : {
1764 140 : libMesh::out << "Updating RB matrices" << std::endl;
1765 4956 : update_RB_system_matrices();
1766 4956 : }
1767 :
1768 406700 : Real RBConstruction::get_RB_error_bound()
1769 : {
1770 406700 : get_rb_evaluation().set_parameters( get_parameters() );
1771 :
1772 34000 : Real error_bound = 0.;
1773 406700 : if (get_preevaluate_thetas_flag())
1774 : {
1775 : // Obtain the pre-evaluated theta functions from the current training parameter index
1776 0 : const auto & evaluated_thetas = get_evaluated_thetas(get_current_training_parameter_index());
1777 0 : error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions(),
1778 0 : &evaluated_thetas);
1779 : }
1780 : else
1781 406700 : error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions());
1782 :
1783 :
1784 406700 : if (normalize_rb_bound_in_greedy)
1785 : {
1786 0 : Real error_bound_normalization = get_rb_evaluation().get_error_bound_normalization();
1787 :
1788 0 : if ((error_bound < abs_training_tolerance) ||
1789 0 : (error_bound_normalization < abs_training_tolerance))
1790 : {
1791 : // We don't want to normalize this error bound if the bound or the
1792 : // normalization value are below the absolute tolerance. Hence do nothing
1793 : // in this case.
1794 : }
1795 : else
1796 0 : error_bound /= error_bound_normalization;
1797 : }
1798 :
1799 406700 : return error_bound;
1800 : }
1801 :
1802 0 : void RBConstruction::recompute_all_residual_terms(bool compute_inner_products)
1803 : {
1804 : // Compute the basis independent terms
1805 0 : Fq_representor_innerprods_computed = false;
1806 0 : compute_Fq_representor_innerprods(compute_inner_products);
1807 :
1808 : // and all the basis dependent terms
1809 0 : unsigned int saved_delta_N = delta_N;
1810 0 : delta_N = get_rb_evaluation().get_n_basis_functions();
1811 :
1812 0 : update_residual_terms(compute_inner_products);
1813 :
1814 0 : delta_N = saved_delta_N;
1815 0 : }
1816 :
1817 4957 : Real RBConstruction::compute_max_error_bound()
1818 : {
1819 280 : LOG_SCOPE("compute_max_error_bound()", "RBConstruction");
1820 :
1821 : // Treat the case with no parameters in a special way
1822 4957 : if (get_n_params() == 0)
1823 : {
1824 : Real max_val;
1825 : if (std::numeric_limits<Real>::has_infinity)
1826 : {
1827 0 : max_val = std::numeric_limits<Real>::infinity();
1828 : }
1829 : else
1830 : {
1831 : max_val = std::numeric_limits<Real>::max();
1832 : }
1833 :
1834 : // Make sure we do at least one solve, but otherwise return a zero error bound
1835 : // when we have no parameters
1836 0 : return (get_rb_evaluation().get_n_basis_functions() == 0) ? max_val : 0.;
1837 : }
1838 :
1839 4957 : training_error_bounds.resize(this->get_local_n_training_samples());
1840 :
1841 : // keep track of the maximum error
1842 140 : unsigned int max_err_index = 0;
1843 140 : Real max_err = 0.;
1844 :
1845 4957 : numeric_index_type first_index = get_first_local_training_index();
1846 411657 : for (unsigned int i=0; i<get_local_n_training_samples(); i++)
1847 : {
1848 : // Load training parameter i, this is only loaded
1849 : // locally since the RB solves are local.
1850 406700 : set_params_from_training_set( first_index+i );
1851 :
1852 : // In case we pre-evaluate the theta functions,
1853 : // also keep track of the current training parameter index.
1854 406700 : if (get_preevaluate_thetas_flag())
1855 0 : set_current_training_parameter_index(first_index+i);
1856 :
1857 :
1858 406700 : training_error_bounds[i] = get_RB_error_bound();
1859 :
1860 406700 : if (training_error_bounds[i] > max_err)
1861 : {
1862 970 : max_err_index = i;
1863 970 : max_err = training_error_bounds[i];
1864 : }
1865 : }
1866 :
1867 4957 : std::pair<numeric_index_type, Real> error_pair(first_index+max_err_index, max_err);
1868 4957 : get_global_max_error_pair(this->comm(),error_pair);
1869 :
1870 : // If we have a serial training set (i.e. a training set that is the same on all processors)
1871 : // just set the parameters on all processors
1872 4957 : if (serial_training_set)
1873 : {
1874 0 : set_params_from_training_set( error_pair.first );
1875 : }
1876 : // otherwise, broadcast the parameter that produced the maximum error
1877 : else
1878 : {
1879 4957 : unsigned int root_id=0;
1880 7876 : if ((get_first_local_training_index() <= error_pair.first) &&
1881 2919 : (error_pair.first < get_last_local_training_index()))
1882 : {
1883 827 : set_params_from_training_set( error_pair.first );
1884 897 : root_id = this->processor_id();
1885 : }
1886 :
1887 4957 : this->comm().sum(root_id); // root_id is only non-zero on one processor
1888 4957 : broadcast_parameters(root_id);
1889 : }
1890 :
1891 4957 : return error_pair.second;
1892 : }
1893 :
1894 4956 : void RBConstruction::update_RB_system_matrices()
1895 : {
1896 280 : LOG_SCOPE("update_RB_system_matrices()", "RBConstruction");
1897 :
1898 4956 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1899 :
1900 5096 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
1901 4956 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1902 :
1903 20562 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1904 : {
1905 31212 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1906 : {
1907 15606 : get_rb_evaluation().RB_Fq_vector[q_f](i) = get_non_dirichlet_Fq_if_avail(q_f)->dot(get_rb_evaluation().get_basis_function(i));
1908 : }
1909 : }
1910 :
1911 9912 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1912 : {
1913 16260 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1914 22608 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1915 : {
1916 11304 : get_rb_evaluation().RB_output_vectors[n][q_l](i) =
1917 11304 : get_output_vector(n,q_l)->dot(get_rb_evaluation().get_basis_function(i));
1918 : }
1919 :
1920 51627 : for (unsigned int j=0; j<RB_size; j++)
1921 : {
1922 1320 : Number value = 0.;
1923 :
1924 46671 : if (compute_RB_inner_product)
1925 : {
1926 : // Compute reduced inner_product_matrix
1927 14700 : temp->zero();
1928 14700 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1929 :
1930 14700 : value = temp->dot( get_rb_evaluation().get_basis_function(i) );
1931 14700 : get_rb_evaluation().RB_inner_product_matrix(i,j) = value;
1932 14700 : if (i!=j)
1933 : {
1934 : // The inner product matrix is assumed
1935 : // to be hermitian
1936 13300 : get_rb_evaluation().RB_inner_product_matrix(j,i) = libmesh_conj(value);
1937 : }
1938 : }
1939 :
1940 186684 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1941 : {
1942 : // Compute reduced Aq matrix
1943 140013 : temp->zero();
1944 140013 : get_non_dirichlet_Aq_if_avail(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1945 :
1946 140013 : value = (*temp).dot(get_rb_evaluation().get_basis_function(i));
1947 140013 : get_rb_evaluation().RB_Aq_vector[q_a](i,j) = value;
1948 :
1949 140013 : if (i!=j)
1950 : {
1951 125145 : temp->zero();
1952 125145 : get_non_dirichlet_Aq_if_avail(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(i));
1953 :
1954 125145 : value = (*temp).dot(get_rb_evaluation().get_basis_function(j));
1955 125145 : get_rb_evaluation().RB_Aq_vector[q_a](j,i) = value;
1956 : }
1957 : }
1958 : }
1959 : }
1960 4956 : }
1961 :
1962 :
1963 4673 : void RBConstruction::update_residual_terms(bool compute_inner_products)
1964 : {
1965 264 : LOG_SCOPE("update_residual_terms()", "RBConstruction");
1966 :
1967 132 : libMesh::out << "Updating RB residual terms" << std::endl;
1968 :
1969 4673 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1970 :
1971 132 : if (store_untransformed_basis)
1972 : {
1973 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
1974 : }
1975 :
1976 18692 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1977 : {
1978 28038 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1979 : {
1980 : // Initialize the vector in which we'll store the representor
1981 14019 : if (!get_rb_evaluation().Aq_representor[q_a][i])
1982 : {
1983 27246 : get_rb_evaluation().Aq_representor[q_a][i] = NumericVector<Number>::build(this->comm());
1984 14019 : get_rb_evaluation().Aq_representor[q_a][i]->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1985 : }
1986 :
1987 396 : libmesh_assert(get_rb_evaluation().Aq_representor[q_a][i]->size() == this->n_dofs() &&
1988 : get_rb_evaluation().Aq_representor[q_a][i]->local_size() == this->n_local_dofs() );
1989 :
1990 14019 : rhs->zero();
1991 14019 : if (!store_untransformed_basis)
1992 : {
1993 14019 : get_Aq(q_a)->vector_mult(*rhs, get_rb_evaluation().get_basis_function(i));
1994 : }
1995 : else
1996 : {
1997 0 : get_Aq(q_a)->vector_mult(*rhs, *_untransformed_basis_functions[i]);
1998 : }
1999 14019 : rhs->scale(-1.);
2000 :
2001 14019 : if (!is_quiet())
2002 : {
2003 0 : libMesh::out << "Starting solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
2004 0 : << Utility::get_timestamp() << std::endl;
2005 : }
2006 :
2007 14415 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2008 :
2009 14019 : if (assert_convergence)
2010 14019 : check_convergence(*inner_product_solver);
2011 :
2012 14019 : if (!is_quiet())
2013 : {
2014 0 : libMesh::out << "Finished solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
2015 0 : << Utility::get_timestamp() << std::endl;
2016 0 : libMesh::out << this->n_linear_iterations() << " iterations, final residual "
2017 0 : << this->final_linear_residual() << std::endl;
2018 : }
2019 :
2020 : // Store the representor
2021 14019 : *get_rb_evaluation().Aq_representor[q_a][i] = *solution;
2022 : }
2023 : }
2024 :
2025 : // Now compute and store the inner products (if requested)
2026 4673 : if (compute_inner_products)
2027 : {
2028 :
2029 19286 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2030 : {
2031 14613 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,*Fq_representor[q_f]);
2032 :
2033 58452 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
2034 : {
2035 87678 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
2036 : {
2037 43839 : get_rb_evaluation().Fq_Aq_representor_innerprods[q_f][q_a][i] =
2038 43839 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a][i]);
2039 : }
2040 : }
2041 : }
2042 :
2043 132 : unsigned int q=0;
2044 18692 : for (unsigned int q_a1=0; q_a1<get_rb_theta_expansion().get_n_A_terms(); q_a1++)
2045 : {
2046 42057 : for (unsigned int q_a2=q_a1; q_a2<get_rb_theta_expansion().get_n_A_terms(); q_a2++)
2047 : {
2048 56076 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
2049 : {
2050 278364 : for (unsigned int j=0; j<RB_size; j++)
2051 : {
2052 250326 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][j]);
2053 250326 : get_rb_evaluation().Aq_Aq_representor_innerprods[q][i][j] =
2054 250326 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][i]);
2055 :
2056 250326 : if (i != j)
2057 : {
2058 222288 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][i]);
2059 222288 : get_rb_evaluation().Aq_Aq_representor_innerprods[q][j][i] =
2060 222288 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][j]);
2061 : }
2062 : }
2063 : }
2064 28038 : q++;
2065 : }
2066 : }
2067 : } // end if (compute_inner_products)
2068 4673 : }
2069 :
2070 576 : SparseMatrix<Number> & RBConstruction::get_matrix_for_output_dual_solves()
2071 : {
2072 576 : return *inner_product_matrix;
2073 : }
2074 :
2075 284 : void RBConstruction::compute_output_dual_innerprods()
2076 : {
2077 : // Skip calculations if we've already computed the output dual norms
2078 284 : if (!output_dual_innerprods_computed)
2079 : {
2080 : // Short circuit if we don't have any outputs
2081 284 : if (get_rb_theta_expansion().get_n_outputs() == 0)
2082 : {
2083 142 : output_dual_innerprods_computed = true;
2084 142 : return;
2085 : }
2086 :
2087 : // Only log if we get to here
2088 8 : LOG_SCOPE("compute_output_dual_innerprods()", "RBConstruction");
2089 :
2090 4 : libMesh::out << "Compute output dual inner products" << std::endl;
2091 :
2092 : // Find out the largest value of Q_l
2093 4 : unsigned int max_Q_l = 0;
2094 710 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2095 568 : max_Q_l = (get_rb_theta_expansion().get_n_output_terms(n) > max_Q_l) ? get_rb_theta_expansion().get_n_output_terms(n) : max_Q_l;
2096 :
2097 142 : std::vector<std::unique_ptr<NumericVector<Number>>> L_q_representor(max_Q_l);
2098 284 : for (unsigned int q=0; q<max_Q_l; q++)
2099 : {
2100 142 : L_q_representor[q] = NumericVector<Number>::build(this->comm());
2101 146 : L_q_representor[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2102 : }
2103 :
2104 710 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2105 : {
2106 1136 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
2107 : {
2108 568 : rhs->zero();
2109 568 : rhs->add(1., *get_output_vector(n,q_l));
2110 :
2111 568 : if (!is_quiet())
2112 0 : libMesh::out << "Starting solve n=" << n << ", q_l=" << q_l
2113 0 : << " in RBConstruction::compute_output_dual_innerprods() at "
2114 0 : << Utility::get_timestamp() << std::endl;
2115 :
2116 : // Use the main linear solver here instead of the inner_product solver, since
2117 : // get_matrix_for_output_dual_solves() may not return the inner product matrix.
2118 568 : solve_for_matrix_and_rhs(*get_linear_solver(), get_matrix_for_output_dual_solves(), *rhs);
2119 :
2120 : // We possibly perform multiple solves here with the same matrix, hence
2121 : // set reuse_preconditioner(true) (and set it back to false again below
2122 : // at the end of this function).
2123 568 : linear_solver->reuse_preconditioner(true);
2124 :
2125 568 : if (assert_convergence)
2126 568 : check_convergence(*get_linear_solver());
2127 :
2128 568 : if (!is_quiet())
2129 : {
2130 0 : libMesh::out << "Finished solve n=" << n << ", q_l=" << q_l
2131 0 : << " in RBConstruction::compute_output_dual_innerprods() at "
2132 0 : << Utility::get_timestamp() << std::endl;
2133 :
2134 0 : libMesh::out << this->n_linear_iterations()
2135 0 : << " iterations, final residual "
2136 0 : << this->final_linear_residual() << std::endl;
2137 : }
2138 :
2139 600 : *L_q_representor[q_l] = *solution;
2140 : }
2141 :
2142 16 : unsigned int q=0;
2143 1136 : for (unsigned int q_l1=0; q_l1<get_rb_theta_expansion().get_n_output_terms(n); q_l1++)
2144 : {
2145 568 : get_matrix_for_output_dual_solves().vector_mult(*inner_product_storage_vector, *L_q_representor[q_l1]);
2146 :
2147 1136 : for (unsigned int q_l2=q_l1; q_l2<get_rb_theta_expansion().get_n_output_terms(n); q_l2++)
2148 : {
2149 600 : output_dual_innerprods[n][q] = L_q_representor[q_l2]->dot(*inner_product_storage_vector);
2150 64 : libMesh::out << "output_dual_innerprods[" << n << "][" << q << "] = " << output_dual_innerprods[n][q] << std::endl;
2151 :
2152 568 : q++;
2153 : }
2154 : }
2155 : }
2156 :
2157 : // We may not need to use linear_solver again (e.g. this would happen if we use
2158 : // extra_linear_solver for the truth_solves). As a result, let's clear linear_solver
2159 : // to release any memory it may be taking up. If we do need it again, it will
2160 : // be initialized when necessary.
2161 142 : linear_solver->clear();
2162 142 : linear_solver->reuse_preconditioner(false);
2163 :
2164 142 : output_dual_innerprods_computed = true;
2165 134 : }
2166 :
2167 142 : get_rb_evaluation().output_dual_innerprods = output_dual_innerprods;
2168 : }
2169 :
2170 284 : void RBConstruction::compute_Fq_representor_innerprods(bool compute_inner_products)
2171 : {
2172 :
2173 : // Skip calculations if we've already computed the Fq_representors
2174 284 : if (!Fq_representor_innerprods_computed)
2175 : {
2176 : // Only log if we get to here
2177 8 : LOG_SCOPE("compute_Fq_representor_innerprods()", "RBConstruction");
2178 :
2179 1278 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2180 : {
2181 1022 : if (!Fq_representor[q_f])
2182 : {
2183 1932 : Fq_representor[q_f] = NumericVector<Number>::build(this->comm());
2184 1022 : Fq_representor[q_f]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2185 : }
2186 :
2187 28 : libmesh_assert(Fq_representor[q_f]->size() == this->n_dofs() &&
2188 : Fq_representor[q_f]->local_size() == this->n_local_dofs() );
2189 :
2190 994 : rhs->zero();
2191 994 : rhs->add(1., *get_Fq(q_f));
2192 :
2193 994 : if (!is_quiet())
2194 0 : libMesh::out << "Starting solve q_f=" << q_f
2195 0 : << " in RBConstruction::update_residual_terms() at "
2196 0 : << Utility::get_timestamp() << std::endl;
2197 :
2198 1022 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2199 :
2200 994 : if (assert_convergence)
2201 994 : check_convergence(*inner_product_solver);
2202 :
2203 994 : if (!is_quiet())
2204 : {
2205 0 : libMesh::out << "Finished solve q_f=" << q_f
2206 0 : << " in RBConstruction::update_residual_terms() at "
2207 0 : << Utility::get_timestamp() << std::endl;
2208 :
2209 0 : libMesh::out << this->n_linear_iterations()
2210 0 : << " iterations, final residual "
2211 0 : << this->final_linear_residual() << std::endl;
2212 : }
2213 :
2214 1050 : *Fq_representor[q_f] = *solution;
2215 : }
2216 :
2217 284 : if (compute_inner_products)
2218 : {
2219 8 : unsigned int q=0;
2220 :
2221 1278 : for (unsigned int q_f1=0; q_f1<get_rb_theta_expansion().get_n_F_terms(); q_f1++)
2222 : {
2223 994 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *Fq_representor[q_f1]);
2224 :
2225 4118 : for (unsigned int q_f2=q_f1; q_f2<get_rb_theta_expansion().get_n_F_terms(); q_f2++)
2226 : {
2227 3212 : Fq_representor_innerprods[q] = inner_product_storage_vector->dot(*Fq_representor[q_f2]);
2228 :
2229 3124 : q++;
2230 : }
2231 : }
2232 : } // end if (compute_inner_products)
2233 :
2234 284 : Fq_representor_innerprods_computed = true;
2235 : }
2236 :
2237 284 : get_rb_evaluation().Fq_representor_innerprods = Fq_representor_innerprods;
2238 284 : }
2239 :
2240 214 : void RBConstruction::load_rb_solution()
2241 : {
2242 12 : LOG_SCOPE("load_rb_solution()", "RBConstruction");
2243 :
2244 214 : solution->zero();
2245 :
2246 214 : libmesh_error_msg_if(get_rb_evaluation().RB_solution.size() > get_rb_evaluation().get_n_basis_functions(),
2247 : "ERROR: System contains " << get_rb_evaluation().get_n_basis_functions() << " basis functions."
2248 : << " RB_solution vector contains " << get_rb_evaluation().RB_solution.size() << " entries."
2249 : << " RB_solution in RBConstruction::load_rb_solution is too long!");
2250 :
2251 3776 : for (auto i : make_range(get_rb_evaluation().RB_solution.size()))
2252 3556 : solution->add(get_rb_evaluation().RB_solution(i), get_rb_evaluation().get_basis_function(i));
2253 :
2254 214 : update();
2255 214 : }
2256 :
2257 : // The slow (but simple, non-error prone) way to compute the residual dual norm
2258 : // Useful for error checking
2259 0 : Real RBConstruction::compute_residual_dual_norm_slow(const unsigned int N)
2260 : {
2261 0 : LOG_SCOPE("compute_residual_dual_norm_slow()", "RBConstruction");
2262 :
2263 : // Put the residual in rhs in order to compute the norm of the Riesz representor
2264 : // Note that this only works in serial since otherwise each processor will
2265 : // have a different parameter value during the Greedy training.
2266 :
2267 0 : std::unique_ptr<NumericVector<Number>> RB_sol = NumericVector<Number>::build(comm());
2268 0 : RB_sol->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2269 :
2270 0 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(comm());
2271 0 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2272 :
2273 0 : if (store_untransformed_basis)
2274 : {
2275 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
2276 : }
2277 :
2278 0 : for (unsigned int i=0; i<N; i++)
2279 : {
2280 0 : if (!store_untransformed_basis)
2281 : {
2282 0 : RB_sol->add(get_rb_evaluation().RB_solution(i), get_rb_evaluation().get_basis_function(i));
2283 : }
2284 : else
2285 : {
2286 0 : RB_sol->add(get_rb_evaluation().RB_solution(i), *_untransformed_basis_functions[i]);
2287 : }
2288 : }
2289 :
2290 0 : this->truth_assembly();
2291 0 : matrix->vector_mult(*temp, *RB_sol);
2292 0 : rhs->add(-1., *temp);
2293 :
2294 : // Then solve to get the Reisz representor
2295 0 : matrix->zero();
2296 0 : matrix->add(1., *inner_product_matrix);
2297 :
2298 0 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2299 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *solution);
2300 0 : Number slow_residual_norm_sq = solution->dot(*inner_product_storage_vector);
2301 :
2302 0 : return std::sqrt( libmesh_real(slow_residual_norm_sq) );
2303 0 : }
2304 :
2305 1059389 : SparseMatrix<Number> * RBConstruction::get_inner_product_matrix()
2306 : {
2307 1059389 : return inner_product_matrix.get();
2308 : }
2309 :
2310 28400 : const SparseMatrix<Number> * RBConstruction::get_inner_product_matrix() const
2311 : {
2312 28400 : return inner_product_matrix.get();
2313 : }
2314 :
2315 0 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix()
2316 : {
2317 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2318 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_inner_product_matrix.");
2319 :
2320 0 : return non_dirichlet_inner_product_matrix.get();
2321 : }
2322 :
2323 0 : const SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix() const
2324 : {
2325 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2326 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_inner_product_matrix.");
2327 :
2328 0 : return non_dirichlet_inner_product_matrix.get();
2329 : }
2330 :
2331 1059388 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix_if_avail()
2332 : {
2333 1059388 : if (store_non_dirichlet_operators)
2334 : {
2335 0 : return get_non_dirichlet_inner_product_matrix();
2336 : }
2337 :
2338 1059388 : return get_inner_product_matrix();
2339 : }
2340 :
2341 28400 : const SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix_if_avail() const
2342 : {
2343 28400 : if (store_non_dirichlet_operators)
2344 : {
2345 0 : return get_non_dirichlet_inner_product_matrix();
2346 : }
2347 :
2348 28400 : return get_inner_product_matrix();
2349 : }
2350 :
2351 1130697 : SparseMatrix<Number> * RBConstruction::get_Aq(unsigned int q)
2352 : {
2353 1130697 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
2354 : "Error: We must have q < Q_a in get_Aq.");
2355 :
2356 1162917 : return Aq_vector[q].get();
2357 : }
2358 :
2359 0 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_Aq(unsigned int q)
2360 : {
2361 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2362 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Aq.");
2363 :
2364 0 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
2365 : "Error: We must have q < Q_a in get_Aq.");
2366 :
2367 0 : return non_dirichlet_Aq_vector[q].get();
2368 : }
2369 :
2370 265158 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_Aq_if_avail(unsigned int q)
2371 : {
2372 265158 : if (store_non_dirichlet_operators)
2373 : {
2374 0 : return get_non_dirichlet_Aq(q);
2375 : }
2376 :
2377 265158 : return get_Aq(q);
2378 : }
2379 :
2380 171800 : NumericVector<Number> * RBConstruction::get_Fq(unsigned int q)
2381 : {
2382 171800 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_F_terms(),
2383 : "Error: We must have q < Q_f in get_Fq.");
2384 :
2385 176696 : return Fq_vector[q].get();
2386 : }
2387 :
2388 0 : NumericVector<Number> * RBConstruction::get_non_dirichlet_Fq(unsigned int q)
2389 : {
2390 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2391 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Fq.");
2392 :
2393 0 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_F_terms(),
2394 : "Error: We must have q < Q_f in get_Fq.");
2395 :
2396 0 : return non_dirichlet_Fq_vector[q].get();
2397 : }
2398 :
2399 15606 : NumericVector<Number> * RBConstruction::get_non_dirichlet_Fq_if_avail(unsigned int q)
2400 : {
2401 15606 : if (store_non_dirichlet_operators)
2402 : {
2403 0 : return get_non_dirichlet_Fq(q);
2404 : }
2405 :
2406 15606 : return get_Fq(q);
2407 : }
2408 :
2409 584312 : NumericVector<Number> * RBConstruction::get_output_vector(unsigned int n, unsigned int q_l)
2410 : {
2411 584312 : libmesh_error_msg_if((n >= get_rb_theta_expansion().get_n_outputs()) ||
2412 : (q_l >= get_rb_theta_expansion().get_n_output_terms(n)),
2413 : "Error: We must have n < n_outputs and "
2414 : "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_output_vector.");
2415 :
2416 617688 : return outputs_vector[n][q_l].get();
2417 : }
2418 :
2419 0 : NumericVector<Number> * RBConstruction::get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l)
2420 : {
2421 0 : libmesh_error_msg_if((n >= get_rb_theta_expansion().get_n_outputs()) ||
2422 : (q_l >= get_rb_theta_expansion().get_n_output_terms(n)),
2423 : "Error: We must have n < n_outputs and "
2424 : "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_non_dirichlet_output_vector.");
2425 :
2426 0 : return non_dirichlet_outputs_vector[n][q_l].get();
2427 : }
2428 :
2429 0 : void RBConstruction::get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices)
2430 : {
2431 0 : all_matrices.clear();
2432 :
2433 0 : all_matrices["inner_product"] = get_inner_product_matrix();
2434 :
2435 0 : if (store_non_dirichlet_operators)
2436 : {
2437 0 : all_matrices["non_dirichlet_inner_product"] = get_non_dirichlet_inner_product_matrix();
2438 : }
2439 :
2440 0 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
2441 : {
2442 0 : std::stringstream matrix_name;
2443 0 : matrix_name << "A" << q_a;
2444 0 : all_matrices[matrix_name.str()] = get_Aq(q_a);
2445 :
2446 0 : if (store_non_dirichlet_operators)
2447 : {
2448 0 : matrix_name << "_non_dirichlet";
2449 0 : all_matrices[matrix_name.str()] = get_non_dirichlet_Aq(q_a);
2450 : }
2451 0 : }
2452 0 : }
2453 :
2454 0 : void RBConstruction::get_all_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors)
2455 : {
2456 0 : all_vectors.clear();
2457 :
2458 0 : get_output_vectors(all_vectors);
2459 :
2460 0 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2461 : {
2462 0 : std::stringstream F_vector_name;
2463 0 : F_vector_name << "F" << q_f;
2464 0 : all_vectors[F_vector_name.str()] = get_Fq(q_f);
2465 :
2466 0 : if (store_non_dirichlet_operators)
2467 : {
2468 0 : F_vector_name << "_non_dirichlet";
2469 0 : all_vectors[F_vector_name.str()] = get_non_dirichlet_Fq(q_f);
2470 : }
2471 0 : }
2472 0 : }
2473 :
2474 0 : void RBConstruction::get_output_vectors(std::map<std::string, NumericVector<Number> *> & output_vectors)
2475 : {
2476 0 : output_vectors.clear();
2477 :
2478 0 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2479 0 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
2480 : {
2481 0 : std::stringstream output_name;
2482 0 : output_name << "output_" << n << "_"<< q_l;
2483 0 : output_vectors[output_name.str()] = get_output_vector(n,q_l);
2484 :
2485 0 : if (store_non_dirichlet_operators)
2486 : {
2487 0 : output_name << "_non_dirichlet";
2488 0 : output_vectors[output_name.str()] = get_non_dirichlet_output_vector(n,q_l);
2489 : }
2490 0 : }
2491 0 : }
2492 :
2493 : #ifdef LIBMESH_ENABLE_DIRICHLET
2494 :
2495 568 : std::unique_ptr<DirichletBoundary> RBConstruction::build_zero_dirichlet_boundary_object()
2496 : {
2497 32 : ZeroFunction<> zf;
2498 :
2499 32 : std::set<boundary_id_type> dirichlet_ids;
2500 16 : std::vector<unsigned int> variables;
2501 :
2502 : // The DirichletBoundary constructor clones zf, so it's OK that zf is only in local scope
2503 1136 : return std::make_unique<DirichletBoundary>(dirichlet_ids, variables, &zf);
2504 : }
2505 :
2506 : #endif
2507 :
2508 0 : void RBConstruction::write_riesz_representors_to_files(const std::string & riesz_representors_dir,
2509 : const bool write_binary_residual_representors)
2510 : {
2511 0 : LOG_SCOPE("write_riesz_representors_to_files()", "RBConstruction");
2512 :
2513 : // Write out Riesz representors. These are useful to have when restarting,
2514 : // so you don't have to recompute them all over again.
2515 :
2516 : // First we write out the Fq representors, these are independent of an RBEvaluation object.
2517 0 : libMesh::out << "Writing out the Fq_representors..." << std::endl;
2518 :
2519 0 : std::ostringstream file_name;
2520 0 : const std::string riesz_representor_suffix = (write_binary_residual_representors ? ".xdr" : ".dat");
2521 : struct stat stat_info;
2522 :
2523 : // Residual representors written out to their own separate directory
2524 0 : if ( this->processor_id() == 0)
2525 0 : if ( Utility::mkdir(riesz_representors_dir.c_str()) != 0)
2526 0 : libMesh::out << "Skipping creating residual_representors directory: " << strerror(errno) << std::endl;
2527 :
2528 0 : for (auto i : index_range(Fq_representor))
2529 : {
2530 0 : if (Fq_representor[i])
2531 : {
2532 0 : file_name.str(""); // reset filename
2533 0 : file_name << riesz_representors_dir << "/Fq_representor" << i << riesz_representor_suffix;
2534 :
2535 : // Check to see if file exists, if so, don't overwrite it, we assume it was
2536 : // there from a previous call to this function. Note: if stat returns zero
2537 : // it means it successfully got the file attributes (and therefore the file
2538 : // exists). Because of the following factors:
2539 : // 1.) write_serialized_data takes longer for proc 0 than it does for others,
2540 : // so processors can get out of sync.
2541 : // 2.) The constructor for Xdr opens a (0 length) file on *all* processors,
2542 : // there are typically hundreds of 0-length files created during this loop,
2543 : // and that screws up checking for the existence of files. One way to stay
2544 : // in sync is to but a barrier at each iteration of the loop -- not sure how
2545 : // bad this will affect performance, but it can't be much worse than serialized
2546 : // I/O already is :)
2547 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2548 :
2549 0 : if ( (stat_result != 0) || // file definitely doesn't already exist
2550 0 : (stat_info.st_size == 0)) // file exists, but has zero length (can happen if another proc already opened it!)
2551 : {
2552 : // No need to copy!
2553 : // *solution = *(Fq_representor[i]);
2554 : // std::swap doesn't work on pointers
2555 : //std::swap(solution.get(), Fq_representor[i]);
2556 0 : Fq_representor[i]->swap(*solution);
2557 :
2558 0 : Xdr fqr_data(file_name.str(),
2559 0 : write_binary_residual_representors ? ENCODE : WRITE);
2560 :
2561 0 : write_serialized_data(fqr_data, false);
2562 :
2563 : // Synchronize before moving on
2564 0 : this->comm().barrier();
2565 :
2566 : // Swap back.
2567 0 : Fq_representor[i]->swap(*solution);
2568 :
2569 : // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2570 : // for the system call, be sure to do it only on one processor, etc.
2571 0 : }
2572 : }
2573 : }
2574 :
2575 :
2576 : // Next, write out the Aq representors associated with rb_eval.
2577 0 : libMesh::out << "Writing out the Aq_representors..." << std::endl;
2578 :
2579 0 : const unsigned int jstop = get_rb_evaluation().get_n_basis_functions();
2580 0 : const unsigned int jstart = jstop-get_delta_N();
2581 :
2582 0 : RBEvaluation & rbe = get_rb_evaluation();
2583 0 : for (auto i : index_range(rbe.Aq_representor))
2584 0 : for (unsigned int j=jstart; j<jstop; ++j)
2585 : {
2586 0 : libMesh::out << "Writing out Aq_representor[" << i << "][" << j << "]..." << std::endl;
2587 0 : libmesh_assert(get_rb_evaluation().Aq_representor[i][j]);
2588 :
2589 0 : file_name.str(""); // reset filename
2590 : file_name << riesz_representors_dir
2591 0 : << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2592 :
2593 : {
2594 : // No need to copy! Use swap instead.
2595 : // *solution = *(Aq_representor[i][j]);
2596 0 : rbe.Aq_representor[i][j]->swap(*solution);
2597 :
2598 0 : Xdr aqr_data(file_name.str(),
2599 0 : write_binary_residual_representors ? ENCODE : WRITE);
2600 :
2601 0 : write_serialized_data(aqr_data, false);
2602 :
2603 : // Synchronize before moving on
2604 0 : this->comm().barrier();
2605 :
2606 : // Swap back.
2607 0 : rbe.Aq_representor[i][j]->swap(*solution);
2608 :
2609 : // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2610 : // for the system call, be sure to do it only on one processor, etc.
2611 0 : }
2612 : }
2613 0 : }
2614 :
2615 :
2616 :
2617 0 : void RBConstruction::read_riesz_representors_from_files(const std::string & riesz_representors_dir,
2618 : const bool read_binary_residual_representors)
2619 : {
2620 0 : LOG_SCOPE("read_riesz_representors_from_files()", "RBConstruction");
2621 :
2622 0 : libMesh::out << "Reading in the Fq_representors..." << std::endl;
2623 :
2624 0 : const std::string riesz_representor_suffix = (read_binary_residual_representors ? ".xdr" : ".dat");
2625 0 : std::ostringstream file_name;
2626 : struct stat stat_info;
2627 :
2628 : // Read in the Fq_representors. There should be Q_f of these. FIXME:
2629 : // should we be worried about leaks here?
2630 0 : for (const auto & rep : Fq_representor)
2631 0 : libmesh_error_msg_if(rep, "Error, must delete existing Fq_representor before reading in from file.");
2632 :
2633 0 : for (auto i : index_range(Fq_representor))
2634 : {
2635 0 : file_name.str(""); // reset filename
2636 : file_name << riesz_representors_dir
2637 0 : << "/Fq_representor" << i << riesz_representor_suffix;
2638 :
2639 : // On processor zero check to be sure the file exists
2640 0 : if (this->processor_id() == 0)
2641 : {
2642 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2643 :
2644 0 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
2645 : }
2646 :
2647 0 : Xdr fqr_data(file_name.str(),
2648 0 : read_binary_residual_representors ? DECODE : READ);
2649 :
2650 0 : read_serialized_data(fqr_data, false);
2651 :
2652 0 : Fq_representor[i] = NumericVector<Number>::build(this->comm());
2653 0 : Fq_representor[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2654 :
2655 : // No need to copy, just swap
2656 : // *Fq_representor[i] = *solution;
2657 0 : Fq_representor[i]->swap(*solution);
2658 0 : }
2659 :
2660 : // Alert the update_residual_terms() function that we don't need to recompute
2661 : // the Fq_representors as we have already read them in from file!
2662 0 : Fq_representor_innerprods_computed = true;
2663 :
2664 :
2665 0 : libMesh::out << "Reading in the Aq_representors..." << std::endl;
2666 :
2667 : // Read in the Aq representors. The class makes room for [Q_a][Nmax] of these. We are going to
2668 : // read in [Q_a][get_rb_evaluation().get_n_basis_functions()]. FIXME:
2669 : // should we be worried about leaks in the locations where we're about to fill entries?
2670 0 : RBEvaluation & rbe = get_rb_evaluation();
2671 0 : for (const auto & row : rbe.Aq_representor)
2672 0 : for (const auto & rep : row)
2673 0 : libmesh_error_msg_if(rep, "Error, must delete existing Aq_representor before reading in from file.");
2674 :
2675 : // Now ready to read them in from file!
2676 0 : for (auto i : index_range(rbe.Aq_representor))
2677 0 : for (unsigned int j=0; j<rbe.get_n_basis_functions(); ++j)
2678 : {
2679 0 : file_name.str(""); // reset filename
2680 : file_name << riesz_representors_dir
2681 0 : << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2682 :
2683 : // On processor zero check to be sure the file exists
2684 0 : if (this->processor_id() == 0)
2685 : {
2686 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2687 :
2688 0 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
2689 : }
2690 :
2691 0 : Xdr aqr_data(file_name.str(), read_binary_residual_representors ? DECODE : READ);
2692 :
2693 0 : read_serialized_data(aqr_data, false);
2694 :
2695 0 : rbe.Aq_representor[i][j] = NumericVector<Number>::build(this->comm());
2696 0 : rbe.Aq_representor[i][j]->init (n_dofs(), n_local_dofs(),
2697 0 : false, PARALLEL);
2698 :
2699 : // No need to copy, just swap
2700 0 : rbe.Aq_representor[i][j]->swap(*solution);
2701 0 : }
2702 0 : }
2703 :
2704 160467 : void RBConstruction::check_convergence(LinearSolver<Number> & input_solver)
2705 : {
2706 : libMesh::LinearConvergenceReason conv_flag;
2707 :
2708 160467 : conv_flag = input_solver.get_converged_reason();
2709 :
2710 160467 : libmesh_error_msg_if(conv_flag < 0, "Convergence error. Error id: " << conv_flag);
2711 160467 : }
2712 :
2713 0 : bool RBConstruction::get_convergence_assertion_flag() const
2714 : {
2715 0 : return assert_convergence;
2716 : }
2717 :
2718 0 : void RBConstruction::set_convergence_assertion_flag(bool flag)
2719 : {
2720 0 : assert_convergence = flag;
2721 0 : }
2722 :
2723 813684 : bool RBConstruction::get_preevaluate_thetas_flag() const
2724 : {
2725 813684 : return _preevaluate_thetas_flag;
2726 : }
2727 :
2728 0 : void RBConstruction::set_preevaluate_thetas_flag(bool flag)
2729 : {
2730 0 : _preevaluate_thetas_flag = flag;
2731 0 : }
2732 :
2733 0 : unsigned int RBConstruction::get_current_training_parameter_index() const
2734 : {
2735 0 : return _current_training_parameter_index;
2736 : }
2737 :
2738 0 : void RBConstruction::set_current_training_parameter_index(unsigned int index)
2739 : {
2740 0 : _current_training_parameter_index = index;
2741 0 : }
2742 :
2743 : const std::vector<Number> &
2744 0 : RBConstruction::get_evaluated_thetas(unsigned int training_parameter_index) const
2745 : {
2746 0 : const numeric_index_type first_index = get_first_local_training_index();
2747 0 : libmesh_assert(training_parameter_index >= first_index);
2748 :
2749 0 : const numeric_index_type local_index = training_parameter_index - first_index;
2750 0 : libmesh_assert(local_index < _evaluated_thetas.size());
2751 :
2752 0 : return _evaluated_thetas[local_index];
2753 : }
2754 :
2755 0 : void RBConstruction::preevaluate_thetas()
2756 : {
2757 0 : LOG_SCOPE("preevaluate_thetas()", "RBConstruction");
2758 :
2759 0 : _evaluated_thetas.resize(get_local_n_training_samples());
2760 :
2761 : // Early return if we've already preevaluated thetas.
2762 0 : if (_preevaluate_thetas_completed)
2763 0 : return;
2764 :
2765 0 : if ( get_local_n_training_samples() == 0 )
2766 0 : return;
2767 :
2768 0 : auto & rb_theta_expansion = get_rb_evaluation().get_rb_theta_expansion();
2769 0 : const unsigned int n_A_terms = rb_theta_expansion.get_n_A_terms();
2770 0 : const unsigned int n_F_terms = rb_theta_expansion.get_n_F_terms();
2771 0 : const unsigned int n_outputs = rb_theta_expansion.get_total_n_output_terms();
2772 :
2773 : // Collect all training parameters
2774 : // TODO: Here instead of using a vector of RBParameters objects,
2775 : // we could use a single RBParameters object with multiple samples.
2776 : // This would save memory over the current approach, but that may
2777 : // not be a big deal in practice unless the number of training samples
2778 : // is very large for some reason.
2779 0 : std::vector<RBParameters> mus(get_local_n_training_samples());
2780 0 : const numeric_index_type first_index = get_first_local_training_index();
2781 0 : for (unsigned int i=0; i<get_local_n_training_samples(); i++)
2782 : {
2783 : // Load training parameter i, this is only loaded
2784 : // locally since the RB solves are local.
2785 0 : set_params_from_training_set( first_index+i );
2786 0 : mus[i] = get_parameters();
2787 0 : _evaluated_thetas[i].resize(n_A_terms + n_F_terms + n_outputs);
2788 : }
2789 :
2790 : // Evaluate thetas for all training parameters simultaneously
2791 0 : for (unsigned int q_a=0; q_a<n_A_terms; q_a++)
2792 : {
2793 0 : const auto A_vals = rb_theta_expansion.eval_A_theta(q_a, mus);
2794 0 : for (auto i : make_range(get_local_n_training_samples()))
2795 0 : _evaluated_thetas[i][q_a] = A_vals[i];
2796 : }
2797 :
2798 0 : for (unsigned int q_f=0; q_f<n_F_terms; q_f++)
2799 : {
2800 0 : const auto F_vals = rb_theta_expansion.eval_F_theta(q_f, mus);
2801 0 : for (auto i : make_range(get_local_n_training_samples()))
2802 0 : _evaluated_thetas[i][n_A_terms + q_f] = F_vals[i];
2803 : }
2804 :
2805 : {
2806 0 : unsigned int output_counter = 0;
2807 0 : for (unsigned int n=0; n<rb_theta_expansion.get_n_outputs(); n++)
2808 0 : for (unsigned int q_l=0; q_l<rb_theta_expansion.get_n_output_terms(n); q_l++)
2809 : {
2810 : // Evaluate the current output functional term for all
2811 : // training parameters simultaneously.
2812 0 : const auto output_vals = rb_theta_expansion.eval_output_theta(n, q_l, mus);
2813 :
2814 : // TODO: the size of _evaluated_thetas is currently assumed to be
2815 : // the same as get_local_n_training_samples(), but this won't be
2816 : // the case if we use RBParameters objects that have multiple samples.
2817 : // So just make sure that's the case for now.
2818 0 : libmesh_error_msg_if(output_vals.size() != get_local_n_training_samples(),
2819 : "We currently only support single-sample RBParameters "
2820 : "objects during the training stage.");
2821 :
2822 0 : for (auto i : index_range(output_vals))
2823 0 : _evaluated_thetas[i][n_A_terms + n_F_terms + output_counter] = output_vals[i];
2824 :
2825 : // Go to next output term
2826 0 : output_counter++;
2827 : }
2828 : }
2829 :
2830 0 : _preevaluate_thetas_completed = true;
2831 0 : }
2832 :
2833 0 : void RBConstruction::reset_preevaluate_thetas_completed()
2834 : {
2835 0 : _preevaluate_thetas_completed = false;
2836 0 : }
2837 :
2838 : } // namespace libMesh
|
