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