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 160 : 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 35375 : context.get_elem_jacobian() *= scalar;
893 35375 : 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 300 : 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 24 : LOG_SCOPE("train_reduced_basis_with_greedy()", "RBConstruction");
1215 :
1216 : // Make sure we clear the "untransformed basis functions" if we're training
1217 : // a new RB basis. We need to do this here since we store the untransformed
1218 : // basis functions in RBConstruction, as opposed to the standard basis functions,
1219 : // which are stored in RBEvaluation.
1220 284 : if (store_untransformed_basis && (get_rb_evaluation().get_n_basis_functions() == 0))
1221 0 : _untransformed_basis_functions.clear();
1222 :
1223 8 : int count = 0;
1224 :
1225 284 : RBEvaluation & rbe = get_rb_evaluation();
1226 :
1227 : // initialize rbe's parameters
1228 284 : rbe.initialize_parameters(*this);
1229 :
1230 : // possibly resize data structures according to Nmax
1231 284 : if (resize_rb_eval_data)
1232 284 : rbe.resize_data_structures(get_Nmax());
1233 :
1234 : // Clear the Greedy param list
1235 284 : for (auto & plist : rbe.greedy_param_list)
1236 0 : plist.clear();
1237 :
1238 276 : rbe.greedy_param_list.clear();
1239 :
1240 8 : Real training_greedy_error = 0.;
1241 :
1242 :
1243 : // If we are continuing from a previous training run,
1244 : // we might already be at the max number of basis functions.
1245 : // If so, we can just return.
1246 284 : if (rbe.get_n_basis_functions() >= get_Nmax())
1247 : {
1248 0 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1249 0 : << get_Nmax() << std::endl;
1250 0 : return 0.;
1251 : }
1252 :
1253 : // Optionally pre-evaluate the theta functions on the entire (local) training parameter set.
1254 284 : if (get_preevaluate_thetas_flag())
1255 0 : preevaluate_thetas();
1256 :
1257 284 : if(!skip_residual_in_train_reduced_basis)
1258 : {
1259 : // Compute the dual norms of the outputs if we haven't already done so.
1260 284 : compute_output_dual_innerprods();
1261 :
1262 : // Compute the Fq Riesz representor dual norms if we haven't already done so.
1263 284 : compute_Fq_representor_innerprods();
1264 : }
1265 :
1266 8 : libMesh::out << std::endl << "---- Performing Greedy basis enrichment ----" << std::endl;
1267 8 : Real initial_greedy_error = 0.;
1268 8 : bool initial_greedy_error_initialized = false;
1269 : while (true)
1270 : {
1271 140 : libMesh::out << std::endl << "---- Basis dimension: "
1272 4957 : << rbe.get_n_basis_functions() << " ----" << std::endl;
1273 :
1274 4957 : if (count > 0 || (count==0 && use_empty_rb_solve_in_greedy))
1275 : {
1276 140 : libMesh::out << "Performing RB solves on training set" << std::endl;
1277 4957 : training_greedy_error = compute_max_error_bound();
1278 :
1279 140 : libMesh::out << "Maximum error bound is " << training_greedy_error << std::endl << std::endl;
1280 :
1281 : // record the initial error
1282 4957 : if (!initial_greedy_error_initialized)
1283 : {
1284 8 : initial_greedy_error = training_greedy_error;
1285 8 : initial_greedy_error_initialized = true;
1286 : }
1287 :
1288 : // Break out of training phase if we have reached Nmax
1289 : // or if the training_tolerance is satisfied.
1290 4957 : if (greedy_termination_test(training_greedy_error, initial_greedy_error, count))
1291 0 : break;
1292 : }
1293 :
1294 140 : libMesh::out << "Performing truth solve at parameter:" << std::endl;
1295 4956 : print_parameters();
1296 :
1297 : // Update the list of Greedily selected parameters
1298 4956 : this->update_greedy_param_list();
1299 :
1300 : // Perform an Offline truth solve for the current parameter
1301 4956 : truth_solve(-1);
1302 :
1303 4956 : if (check_if_zero_truth_solve())
1304 : {
1305 0 : libMesh::out << "Zero basis function encountered hence ending basis enrichment" << std::endl;
1306 0 : break;
1307 : }
1308 :
1309 : // Add orthogonal part of the snapshot to the RB space
1310 140 : libMesh::out << "Enriching the RB space" << std::endl;
1311 4956 : enrich_RB_space();
1312 :
1313 4956 : update_system();
1314 :
1315 : // Check if we've reached Nmax now. We do this before calling
1316 : // update_residual_terms() since we can skip that step if we've
1317 : // already reached Nmax.
1318 4956 : if (rbe.get_n_basis_functions() >= this->get_Nmax())
1319 : {
1320 8 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1321 16 : << get_Nmax() << std::endl;
1322 8 : break;
1323 : }
1324 :
1325 4673 : if(!skip_residual_in_train_reduced_basis)
1326 : {
1327 4673 : update_residual_terms();
1328 : }
1329 :
1330 : // Increment counter
1331 4673 : count++;
1332 : }
1333 284 : this->update_greedy_param_list();
1334 :
1335 8 : return training_greedy_error;
1336 : }
1337 :
1338 0 : void RBConstruction::enrich_basis_from_rhs_terms(const bool resize_rb_eval_data)
1339 : {
1340 0 : LOG_SCOPE("enrich_basis_from_rhs_terms()", "RBConstruction");
1341 :
1342 : // initialize rb_eval's parameters
1343 0 : get_rb_evaluation().initialize_parameters(*this);
1344 :
1345 : // possibly resize data structures according to Nmax
1346 0 : if (resize_rb_eval_data)
1347 : {
1348 0 : get_rb_evaluation().resize_data_structures(get_Nmax());
1349 : }
1350 :
1351 0 : libMesh::out << std::endl << "---- Enriching basis from rhs terms ----" << std::endl;
1352 :
1353 0 : truth_assembly();
1354 :
1355 0 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1356 : {
1357 0 : libMesh::out << std::endl << "Performing truth solve with rhs from rhs term " << q_f << std::endl;
1358 :
1359 0 : *rhs = *get_Fq(q_f);
1360 :
1361 0 : if (rhs->l2_norm() == 0)
1362 : {
1363 : // Skip enrichment if the rhs is zero
1364 0 : continue;
1365 : }
1366 :
1367 : // truth_assembly assembles into matrix and rhs, so use those for the solve
1368 0 : if (extra_linear_solver)
1369 : {
1370 : // If extra_linear_solver has been initialized, then we use it for the
1371 : // truth solves.
1372 0 : solve_for_matrix_and_rhs(*extra_linear_solver, *matrix, *rhs);
1373 :
1374 0 : if (assert_convergence)
1375 0 : check_convergence(*extra_linear_solver);
1376 : }
1377 : else
1378 : {
1379 0 : solve_for_matrix_and_rhs(*get_linear_solver(), *matrix, *rhs);
1380 :
1381 0 : if (assert_convergence)
1382 0 : check_convergence(*get_linear_solver());
1383 : }
1384 :
1385 : // Debugging: enable this code to print the rhs that was used in
1386 : // the most recent truth solve to a uniquely-named file.
1387 : #if 0
1388 : {
1389 : char temp_file[] = "truth_rhs_XXXXXX.dat";
1390 : int fd = mkstemps(temp_file, 4);
1391 : if (fd != -1)
1392 : {
1393 : libMesh::out << "Writing truth system rhs to file: " << temp_file << std::endl;
1394 : rhs->print_matlab(std::string(temp_file));
1395 : }
1396 : }
1397 : #endif // 0
1398 :
1399 : // Debugging: enable this code to print the most recent truth
1400 : // solution to a uniquely-named file.
1401 : #ifdef LIBMESH_HAVE_EXODUS_API
1402 : #if 0
1403 : {
1404 : // Note: mkstemps creates a file and returns an open file descriptor to it.
1405 : // The filename is created from a template which must have 6 'X' characters followed
1406 : // by a suffix having the specified length (in this case 4, for ".exo").
1407 : char temp_file[] = "truth_XXXXXX.exo";
1408 : int fd = mkstemps(temp_file, 4);
1409 : if (fd != -1)
1410 : {
1411 : libMesh::out << "Writing truth solution to file: " << temp_file << std::endl;
1412 : ExodusII_IO exo_io(this->get_mesh());
1413 : std::set<std::string> system_names = {this->name()};
1414 : exo_io.write_equation_systems(std::string(temp_file),
1415 : this->get_equation_systems(), &system_names);
1416 : }
1417 : }
1418 : #endif // 0
1419 : #endif // LIBMESH_HAVE_EXODUS_API
1420 :
1421 : // Call user-defined post-processing routines on the truth solution.
1422 0 : post_process_truth_solution();
1423 :
1424 : // Add orthogonal part of the snapshot to the RB space
1425 0 : libMesh::out << "Enriching the RB space" << std::endl;
1426 0 : enrich_RB_space();
1427 :
1428 0 : update_system();
1429 : }
1430 0 : }
1431 :
1432 0 : void RBConstruction::train_reduced_basis_with_POD()
1433 : {
1434 : // Make sure we clear the "untransformed basis functions" if we're training
1435 : // a new RB basis. We need to do this here since we store the untransformed
1436 : // basis functions in RBConstruction, as opposed to the standard basis functions,
1437 : // which are stored in RBEvaluation.
1438 0 : if (store_untransformed_basis && (get_rb_evaluation().get_n_basis_functions() == 0))
1439 0 : _untransformed_basis_functions.clear();
1440 :
1441 : // We need to use the same training set on all processes so that
1442 : // the truth solves below work correctly in parallel.
1443 0 : libmesh_error_msg_if(!serial_training_set, "We must use a serial training set with POD");
1444 0 : libmesh_error_msg_if(get_rb_evaluation().get_n_basis_functions() > 0, "Basis should not already be initialized");
1445 :
1446 0 : get_rb_evaluation().initialize_parameters(*this);
1447 0 : get_rb_evaluation().resize_data_structures(get_Nmax());
1448 :
1449 : // Storage for the POD snapshots
1450 0 : unsigned int n_snapshots = get_n_training_samples();
1451 :
1452 0 : if (get_n_params() == 0)
1453 : {
1454 : // In this case we should have generated an empty training set
1455 : // so assert this
1456 0 : libmesh_assert(n_snapshots == 0);
1457 :
1458 : // If we have no parameters, then we should do exactly one "truth solve"
1459 0 : n_snapshots = 1;
1460 : }
1461 :
1462 0 : std::vector<std::unique_ptr<NumericVector<Number>>> POD_snapshots(n_snapshots);
1463 0 : for (unsigned int i=0; i<n_snapshots; i++)
1464 : {
1465 0 : POD_snapshots[i] = NumericVector<Number>::build(this->comm());
1466 0 : POD_snapshots[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1467 : }
1468 :
1469 : // We use the same training set on all processes
1470 0 : libMesh::out << std::endl;
1471 0 : for (unsigned int i=0; i<n_snapshots; i++)
1472 : {
1473 0 : if (get_n_params() > 0)
1474 : {
1475 0 : set_params_from_training_set(i);
1476 : }
1477 :
1478 0 : libMesh::out << "Truth solve " << (i+1) << " of " << n_snapshots << std::endl;
1479 :
1480 0 : truth_solve(-1);
1481 :
1482 0 : *POD_snapshots[i] = *solution;
1483 : }
1484 0 : libMesh::out << std::endl;
1485 :
1486 0 : if (_normalize_solution_snapshots)
1487 : {
1488 0 : libMesh::out << "Normalizing solution snapshots" << std::endl;
1489 0 : for (unsigned int i=0; i<n_snapshots; i++)
1490 : {
1491 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(
1492 0 : *inner_product_storage_vector, *POD_snapshots[i]);
1493 0 : Real norm = std::sqrt(std::real(POD_snapshots[i]->dot(*inner_product_storage_vector)));
1494 :
1495 0 : if (norm > 0.)
1496 0 : POD_snapshots[i]->scale(1./norm);
1497 : }
1498 : }
1499 :
1500 : // Set up the "correlation matrix"
1501 0 : DenseMatrix<Number> correlation_matrix(n_snapshots,n_snapshots);
1502 0 : for (unsigned int i=0; i<n_snapshots; i++)
1503 : {
1504 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(
1505 0 : *inner_product_storage_vector, *POD_snapshots[i]);
1506 :
1507 0 : for (unsigned int j=0; j<=i; j++)
1508 : {
1509 0 : Number inner_prod = (POD_snapshots[j]->dot(*inner_product_storage_vector));
1510 :
1511 0 : correlation_matrix(i,j) = inner_prod;
1512 0 : if(i != j)
1513 : {
1514 0 : correlation_matrix(j,i) = libmesh_conj(inner_prod);
1515 : }
1516 : }
1517 : }
1518 :
1519 : // compute SVD of correlation matrix
1520 0 : DenseVector<Real> sigma( n_snapshots );
1521 0 : DenseMatrix<Number> U( n_snapshots, n_snapshots );
1522 0 : DenseMatrix<Number> VT( n_snapshots, n_snapshots );
1523 0 : correlation_matrix.svd(sigma, U, VT );
1524 :
1525 0 : if (sigma(0) == 0.)
1526 0 : return;
1527 :
1528 : // Add dominant vectors from the POD as basis functions.
1529 0 : unsigned int j = 0;
1530 : while (true)
1531 : {
1532 0 : if (j >= get_Nmax() || j >= n_snapshots)
1533 : {
1534 0 : libMesh::out << "Maximum number of basis functions (" << j << ") reached." << std::endl;
1535 0 : break;
1536 : }
1537 :
1538 : // The "energy" error in the POD approximation is determined by the first omitted
1539 : // singular value, i.e. sigma(j). We normalize by sigma(0), which gives the total
1540 : // "energy", in order to obtain a relative error.
1541 0 : const Real rel_err = std::sqrt(sigma(j)) / std::sqrt(sigma(0));
1542 :
1543 0 : libMesh::out << "Number of basis functions: " << j
1544 0 : << ", POD error norm: " << rel_err << std::endl;
1545 :
1546 0 : if (rel_err < this->rel_training_tolerance)
1547 : {
1548 0 : libMesh::out << "Training tolerance reached." << std::endl;
1549 0 : break;
1550 : }
1551 :
1552 0 : std::unique_ptr< NumericVector<Number> > v = POD_snapshots[j]->zero_clone();
1553 0 : for ( unsigned int i=0; i<n_snapshots; ++i )
1554 : {
1555 0 : v->add( U.el(i, j), *POD_snapshots[i] );
1556 : }
1557 :
1558 0 : Real norm_v = std::sqrt(sigma(j));
1559 0 : v->scale( 1./norm_v );
1560 :
1561 0 : get_rb_evaluation().basis_functions.emplace_back( std::move(v) );
1562 :
1563 0 : j++;
1564 0 : }
1565 0 : libMesh::out << std::endl;
1566 :
1567 0 : this->delta_N = get_rb_evaluation().get_n_basis_functions();
1568 0 : update_system();
1569 :
1570 : // We now compute all terms required to evaluate the RB error indicator.
1571 : // Unlike in the case of the RB Greedy algorithm, for the POD approach
1572 : // we do not need this data in order to compute the basis. However, we
1573 : // do need this data in order to evaluate error indicator quantities in
1574 : // the Online stage, so we compute it now so that it can be saved in
1575 : // the training data.
1576 0 : recompute_all_residual_terms();
1577 0 : }
1578 :
1579 4957 : bool RBConstruction::greedy_termination_test(Real abs_greedy_error,
1580 : Real initial_error,
1581 : int)
1582 : {
1583 4957 : if (abs_greedy_error < this->abs_training_tolerance)
1584 : {
1585 0 : libMesh::out << "Absolute error tolerance reached." << std::endl;
1586 0 : return true;
1587 : }
1588 :
1589 4957 : Real rel_greedy_error = abs_greedy_error/initial_error;
1590 4957 : if (rel_greedy_error < this->rel_training_tolerance)
1591 : {
1592 0 : libMesh::out << "Relative error tolerance reached." << std::endl;
1593 1 : return true;
1594 : }
1595 :
1596 4956 : RBEvaluation & rbe = get_rb_evaluation();
1597 :
1598 4956 : if (rbe.get_n_basis_functions() >= this->get_Nmax())
1599 : {
1600 0 : libMesh::out << "Maximum number of basis functions reached: Nmax = "
1601 0 : << get_Nmax() << std::endl;
1602 0 : return true;
1603 : }
1604 :
1605 4956 : if (exit_on_repeated_greedy_parameters)
1606 31971 : for (auto & plist : rbe.greedy_param_list)
1607 28415 : if (plist == get_parameters())
1608 : {
1609 0 : libMesh::out << "Exiting greedy because the same parameters were selected twice" << std::endl;
1610 0 : return true;
1611 : }
1612 :
1613 140 : return false;
1614 : }
1615 :
1616 5240 : void RBConstruction::update_greedy_param_list()
1617 : {
1618 5240 : get_rb_evaluation().greedy_param_list.push_back( get_parameters() );
1619 5240 : }
1620 :
1621 0 : const RBParameters & RBConstruction::get_greedy_parameter(unsigned int i)
1622 : {
1623 0 : libmesh_error_msg_if(i >= get_rb_evaluation().greedy_param_list.size(),
1624 : "Error: Argument in RBConstruction::get_greedy_parameter is too large.");
1625 :
1626 0 : return get_rb_evaluation().greedy_param_list[i];
1627 : }
1628 :
1629 3556 : Real RBConstruction::truth_solve(int plot_solution)
1630 : {
1631 200 : LOG_SCOPE("truth_solve()", "RBConstruction");
1632 :
1633 3556 : if(store_untransformed_basis && !_untransformed_solution)
1634 : {
1635 0 : _untransformed_solution = solution->zero_clone();
1636 : }
1637 :
1638 3556 : truth_assembly();
1639 :
1640 : // truth_assembly assembles into matrix and rhs, so use those for the solve
1641 3556 : if (extra_linear_solver)
1642 : {
1643 : // If extra_linear_solver has been initialized, then we use it for the
1644 : // truth solves.
1645 0 : solve_for_matrix_and_rhs(*extra_linear_solver, *matrix, *rhs);
1646 :
1647 0 : if (assert_convergence)
1648 0 : check_convergence(*extra_linear_solver);
1649 : }
1650 : else
1651 : {
1652 3556 : solve_for_matrix_and_rhs(*get_linear_solver(), *matrix, *rhs);
1653 :
1654 3556 : if (assert_convergence)
1655 3556 : check_convergence(*get_linear_solver());
1656 : }
1657 :
1658 3556 : if (store_untransformed_basis)
1659 : {
1660 0 : *_untransformed_solution = *solution;
1661 : }
1662 :
1663 : // Call user-defined post-processing routines on the truth solution.
1664 3556 : post_process_truth_solution();
1665 :
1666 3556 : const RBParameters & mu = get_parameters();
1667 :
1668 9260 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1669 : {
1670 5864 : truth_outputs[n] = 0.;
1671 11408 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1672 11488 : truth_outputs[n] += get_rb_theta_expansion().eval_output_theta(n, q_l, mu)*
1673 5784 : get_output_vector(n,q_l)->dot(*solution);
1674 : }
1675 :
1676 : #ifdef LIBMESH_HAVE_EXODUS_API
1677 3556 : if (plot_solution > 0)
1678 : {
1679 0 : ExodusII_IO exo_io(this->get_mesh());
1680 0 : std::set<std::string> system_names = {this->name()};
1681 0 : exo_io.write_equation_systems("truth.exo", this->get_equation_systems(), &system_names);
1682 0 : }
1683 : #else
1684 : libmesh_ignore(plot_solution);
1685 : #endif
1686 :
1687 : // Get the X norm of the truth solution
1688 : // Useful for normalizing our true error data
1689 3556 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *solution);
1690 3656 : Number truth_X_norm = std::sqrt(inner_product_storage_vector->dot(*solution));
1691 :
1692 3656 : return libmesh_real(truth_X_norm);
1693 0 : }
1694 :
1695 284 : bool RBConstruction::is_serial_training_type(const std::string & RB_training_type_in)
1696 : {
1697 284 : return (RB_training_type_in == "POD");
1698 : }
1699 :
1700 284 : void RBConstruction::set_RB_training_type(const std::string & RB_training_type_in)
1701 : {
1702 284 : this->RB_training_type = RB_training_type_in;
1703 :
1704 284 : if(is_serial_training_type(RB_training_type_in))
1705 : {
1706 : // We need to use a serial training set (so that the training
1707 : // set is the same on all processes) if we're using POD
1708 0 : this->serial_training_set = true;
1709 : }
1710 284 : }
1711 :
1712 638 : const std::string & RBConstruction::get_RB_training_type() const
1713 : {
1714 638 : return RB_training_type;
1715 : }
1716 :
1717 284 : void RBConstruction::set_Nmax(unsigned int Nmax_in)
1718 : {
1719 284 : this->Nmax = Nmax_in;
1720 284 : }
1721 :
1722 142 : void RBConstruction::load_basis_function(unsigned int i)
1723 : {
1724 12 : LOG_SCOPE("load_basis_function()", "RBConstruction");
1725 :
1726 4 : libmesh_assert_less (i, get_rb_evaluation().get_n_basis_functions());
1727 :
1728 142 : *solution = get_rb_evaluation().get_basis_function(i);
1729 :
1730 142 : this->update();
1731 142 : }
1732 :
1733 3556 : void RBConstruction::enrich_RB_space()
1734 : {
1735 300 : LOG_SCOPE("enrich_RB_space()", "RBConstruction");
1736 :
1737 3656 : auto new_bf = solution->clone();
1738 :
1739 3556 : std::unique_ptr<NumericVector<Number>> new_untransformed_bf;
1740 3556 : if (store_untransformed_basis)
1741 : {
1742 0 : new_untransformed_bf = _untransformed_solution->clone();
1743 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
1744 : }
1745 :
1746 31971 : for (unsigned int index=0; index<get_rb_evaluation().get_n_basis_functions(); index++)
1747 : {
1748 28415 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *new_bf);
1749 :
1750 : Number scalar =
1751 28415 : inner_product_storage_vector->dot(get_rb_evaluation().get_basis_function(index));
1752 28415 : new_bf->add(-scalar, get_rb_evaluation().get_basis_function(index));
1753 :
1754 28415 : if (store_untransformed_basis)
1755 : {
1756 0 : new_untransformed_bf->add(-scalar, *_untransformed_basis_functions[index]);
1757 : }
1758 : }
1759 :
1760 : // Normalize new_bf
1761 3556 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *new_bf);
1762 3656 : Number new_bf_norm = std::sqrt( inner_product_storage_vector->dot(*new_bf) );
1763 :
1764 3506 : if (new_bf_norm == 0.)
1765 : {
1766 0 : new_bf->zero(); // avoid potential nan's
1767 :
1768 0 : if (store_untransformed_basis)
1769 : {
1770 0 : new_untransformed_bf->zero();
1771 : }
1772 : }
1773 : else
1774 : {
1775 3556 : new_bf->scale(1./new_bf_norm);
1776 :
1777 3556 : if (store_untransformed_basis)
1778 : {
1779 0 : new_untransformed_bf->scale(1./new_bf_norm);
1780 : }
1781 : }
1782 :
1783 : // load the new basis function into the basis_functions vector.
1784 3556 : get_rb_evaluation().basis_functions.emplace_back( std::move(new_bf) );
1785 :
1786 3556 : if (store_untransformed_basis)
1787 : {
1788 0 : _untransformed_basis_functions.emplace_back( std::move(new_untransformed_bf) );
1789 : }
1790 3556 : }
1791 :
1792 4956 : void RBConstruction::update_system()
1793 : {
1794 140 : libMesh::out << "Updating RB matrices" << std::endl;
1795 4956 : update_RB_system_matrices();
1796 4956 : }
1797 :
1798 406700 : Real RBConstruction::get_RB_error_bound()
1799 : {
1800 406700 : get_rb_evaluation().set_parameters( get_parameters() );
1801 :
1802 34000 : Real error_bound = 0.;
1803 406700 : if (get_preevaluate_thetas_flag())
1804 : {
1805 : // Obtain the pre-evaluated theta functions from the current training parameter index
1806 0 : const auto & evaluated_thetas = get_evaluated_thetas(get_current_training_parameter_index());
1807 0 : error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions(),
1808 0 : &evaluated_thetas);
1809 : }
1810 : else
1811 406700 : error_bound = get_rb_evaluation().rb_solve(get_rb_evaluation().get_n_basis_functions());
1812 :
1813 :
1814 406700 : if (normalize_rb_bound_in_greedy)
1815 : {
1816 0 : Real error_bound_normalization = get_rb_evaluation().get_error_bound_normalization();
1817 :
1818 0 : if ((error_bound < abs_training_tolerance) ||
1819 0 : (error_bound_normalization < abs_training_tolerance))
1820 : {
1821 : // We don't want to normalize this error bound if the bound or the
1822 : // normalization value are below the absolute tolerance. Hence do nothing
1823 : // in this case.
1824 : }
1825 : else
1826 0 : error_bound /= error_bound_normalization;
1827 : }
1828 :
1829 406700 : return error_bound;
1830 : }
1831 :
1832 0 : void RBConstruction::recompute_all_residual_terms(bool compute_inner_products)
1833 : {
1834 : // Compute the basis independent terms
1835 0 : Fq_representor_innerprods_computed = false;
1836 0 : compute_Fq_representor_innerprods(compute_inner_products);
1837 :
1838 : // and all the basis dependent terms
1839 0 : unsigned int saved_delta_N = delta_N;
1840 0 : delta_N = get_rb_evaluation().get_n_basis_functions();
1841 :
1842 0 : update_residual_terms(compute_inner_products);
1843 :
1844 0 : delta_N = saved_delta_N;
1845 0 : }
1846 :
1847 4957 : Real RBConstruction::compute_max_error_bound()
1848 : {
1849 420 : LOG_SCOPE("compute_max_error_bound()", "RBConstruction");
1850 :
1851 : // Treat the case with no parameters in a special way
1852 4957 : if (get_n_params() == 0)
1853 : {
1854 : Real max_val;
1855 : if (std::numeric_limits<Real>::has_infinity)
1856 : {
1857 0 : max_val = std::numeric_limits<Real>::infinity();
1858 : }
1859 : else
1860 : {
1861 : max_val = std::numeric_limits<Real>::max();
1862 : }
1863 :
1864 : // Make sure we do at least one solve, but otherwise return a zero error bound
1865 : // when we have no parameters
1866 0 : return (get_rb_evaluation().get_n_basis_functions() == 0) ? max_val : 0.;
1867 : }
1868 :
1869 4957 : training_error_bounds.resize(this->get_local_n_training_samples());
1870 :
1871 : // keep track of the maximum error
1872 140 : unsigned int max_err_index = 0;
1873 140 : Real max_err = 0.;
1874 :
1875 4957 : numeric_index_type first_index = get_first_local_training_index();
1876 411657 : for (unsigned int i=0; i<get_local_n_training_samples(); i++)
1877 : {
1878 : // Load training parameter i, this is only loaded
1879 : // locally since the RB solves are local.
1880 406700 : set_params_from_training_set( first_index+i );
1881 :
1882 : // In case we pre-evaluate the theta functions,
1883 : // also keep track of the current training parameter index.
1884 406700 : if (get_preevaluate_thetas_flag())
1885 0 : set_current_training_parameter_index(first_index+i);
1886 :
1887 :
1888 406700 : training_error_bounds[i] = get_RB_error_bound();
1889 :
1890 406700 : if (training_error_bounds[i] > max_err)
1891 : {
1892 980 : max_err_index = i;
1893 980 : max_err = training_error_bounds[i];
1894 : }
1895 : }
1896 :
1897 4957 : std::pair<numeric_index_type, Real> error_pair(first_index+max_err_index, max_err);
1898 4957 : get_global_max_error_pair(this->comm(),error_pair);
1899 :
1900 : // If we have a serial training set (i.e. a training set that is the same on all processors)
1901 : // just set the parameters on all processors
1902 4957 : if (serial_training_set)
1903 : {
1904 0 : set_params_from_training_set( error_pair.first );
1905 : }
1906 : // otherwise, broadcast the parameter that produced the maximum error
1907 : else
1908 : {
1909 4957 : unsigned int root_id=0;
1910 7883 : if ((get_first_local_training_index() <= error_pair.first) &&
1911 2926 : (error_pair.first < get_last_local_training_index()))
1912 : {
1913 827 : set_params_from_training_set( error_pair.first );
1914 897 : root_id = this->processor_id();
1915 : }
1916 :
1917 4957 : this->comm().sum(root_id); // root_id is only non-zero on one processor
1918 4957 : broadcast_parameters(root_id);
1919 : }
1920 :
1921 4957 : return error_pair.second;
1922 : }
1923 :
1924 4956 : void RBConstruction::update_RB_system_matrices()
1925 : {
1926 420 : LOG_SCOPE("update_RB_system_matrices()", "RBConstruction");
1927 :
1928 4956 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
1929 :
1930 5096 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
1931 4956 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1932 :
1933 20562 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1934 : {
1935 31212 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1936 : {
1937 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));
1938 : }
1939 : }
1940 :
1941 9912 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
1942 : {
1943 16260 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1944 22608 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1945 : {
1946 11304 : get_rb_evaluation().RB_output_vectors[n][q_l](i) =
1947 11304 : get_output_vector(n,q_l)->dot(get_rb_evaluation().get_basis_function(i));
1948 : }
1949 :
1950 51627 : for (unsigned int j=0; j<RB_size; j++)
1951 : {
1952 1320 : Number value = 0.;
1953 :
1954 46671 : if (compute_RB_inner_product)
1955 : {
1956 : // Compute reduced inner_product_matrix
1957 14700 : temp->zero();
1958 14700 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1959 :
1960 14700 : value = temp->dot( get_rb_evaluation().get_basis_function(i) );
1961 14700 : get_rb_evaluation().RB_inner_product_matrix(i,j) = value;
1962 14700 : if (i!=j)
1963 : {
1964 : // The inner product matrix is assumed
1965 : // to be hermitian
1966 13300 : get_rb_evaluation().RB_inner_product_matrix(j,i) = libmesh_conj(value);
1967 : }
1968 : }
1969 :
1970 186684 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
1971 : {
1972 : // Compute reduced Aq matrix
1973 140013 : temp->zero();
1974 140013 : get_non_dirichlet_Aq_if_avail(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(j));
1975 :
1976 140013 : value = (*temp).dot(get_rb_evaluation().get_basis_function(i));
1977 140013 : get_rb_evaluation().RB_Aq_vector[q_a](i,j) = value;
1978 :
1979 140013 : if (i!=j)
1980 : {
1981 125145 : temp->zero();
1982 125145 : get_non_dirichlet_Aq_if_avail(q_a)->vector_mult(*temp, get_rb_evaluation().get_basis_function(i));
1983 :
1984 125145 : value = (*temp).dot(get_rb_evaluation().get_basis_function(j));
1985 125145 : get_rb_evaluation().RB_Aq_vector[q_a](j,i) = value;
1986 : }
1987 : }
1988 : }
1989 : }
1990 4956 : }
1991 :
1992 :
1993 4673 : void RBConstruction::update_residual_terms(bool compute_inner_products)
1994 : {
1995 396 : LOG_SCOPE("update_residual_terms()", "RBConstruction");
1996 :
1997 132 : libMesh::out << "Updating RB residual terms" << std::endl;
1998 :
1999 4673 : unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
2000 :
2001 132 : if (store_untransformed_basis)
2002 : {
2003 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
2004 : }
2005 :
2006 18692 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
2007 : {
2008 28038 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
2009 : {
2010 : // Initialize the vector in which we'll store the representor
2011 14019 : if (!get_rb_evaluation().Aq_representor[q_a][i])
2012 : {
2013 27246 : get_rb_evaluation().Aq_representor[q_a][i] = NumericVector<Number>::build(this->comm());
2014 14019 : get_rb_evaluation().Aq_representor[q_a][i]->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2015 : }
2016 :
2017 396 : libmesh_assert(get_rb_evaluation().Aq_representor[q_a][i]->size() == this->n_dofs() &&
2018 : get_rb_evaluation().Aq_representor[q_a][i]->local_size() == this->n_local_dofs() );
2019 :
2020 14019 : rhs->zero();
2021 14019 : if (!store_untransformed_basis)
2022 : {
2023 14019 : get_Aq(q_a)->vector_mult(*rhs, get_rb_evaluation().get_basis_function(i));
2024 : }
2025 : else
2026 : {
2027 0 : get_Aq(q_a)->vector_mult(*rhs, *_untransformed_basis_functions[i]);
2028 : }
2029 14019 : rhs->scale(-1.);
2030 :
2031 14019 : if (!is_quiet())
2032 : {
2033 0 : libMesh::out << "Starting solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
2034 0 : << Utility::get_timestamp() << std::endl;
2035 : }
2036 :
2037 14415 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2038 :
2039 14019 : if (assert_convergence)
2040 14019 : check_convergence(*inner_product_solver);
2041 :
2042 14019 : if (!is_quiet())
2043 : {
2044 0 : libMesh::out << "Finished solve [q_a][i]=[" << q_a <<"]["<< i << "] in RBConstruction::update_residual_terms() at "
2045 0 : << Utility::get_timestamp() << std::endl;
2046 0 : libMesh::out << this->n_linear_iterations() << " iterations, final residual "
2047 0 : << this->final_linear_residual() << std::endl;
2048 : }
2049 :
2050 : // Store the representor
2051 14019 : *get_rb_evaluation().Aq_representor[q_a][i] = *solution;
2052 : }
2053 : }
2054 :
2055 : // Now compute and store the inner products (if requested)
2056 4673 : if (compute_inner_products)
2057 : {
2058 :
2059 19286 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2060 : {
2061 14613 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector,*Fq_representor[q_f]);
2062 :
2063 58452 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
2064 : {
2065 87678 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
2066 : {
2067 43839 : get_rb_evaluation().Fq_Aq_representor_innerprods[q_f][q_a][i] =
2068 43839 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a][i]);
2069 : }
2070 : }
2071 : }
2072 :
2073 132 : unsigned int q=0;
2074 18692 : for (unsigned int q_a1=0; q_a1<get_rb_theta_expansion().get_n_A_terms(); q_a1++)
2075 : {
2076 42057 : for (unsigned int q_a2=q_a1; q_a2<get_rb_theta_expansion().get_n_A_terms(); q_a2++)
2077 : {
2078 56076 : for (unsigned int i=(RB_size-delta_N); i<RB_size; i++)
2079 : {
2080 278364 : for (unsigned int j=0; j<RB_size; j++)
2081 : {
2082 250326 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][j]);
2083 250326 : get_rb_evaluation().Aq_Aq_representor_innerprods[q][i][j] =
2084 250326 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][i]);
2085 :
2086 250326 : if (i != j)
2087 : {
2088 222288 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *get_rb_evaluation().Aq_representor[q_a2][i]);
2089 222288 : get_rb_evaluation().Aq_Aq_representor_innerprods[q][j][i] =
2090 222288 : inner_product_storage_vector->dot(*get_rb_evaluation().Aq_representor[q_a1][j]);
2091 : }
2092 : }
2093 : }
2094 28038 : q++;
2095 : }
2096 : }
2097 : } // end if (compute_inner_products)
2098 4673 : }
2099 :
2100 576 : SparseMatrix<Number> & RBConstruction::get_matrix_for_output_dual_solves()
2101 : {
2102 576 : return *inner_product_matrix;
2103 : }
2104 :
2105 284 : void RBConstruction::compute_output_dual_innerprods()
2106 : {
2107 : // Skip calculations if we've already computed the output dual norms
2108 284 : if (!output_dual_innerprods_computed)
2109 : {
2110 : // Short circuit if we don't have any outputs
2111 284 : if (get_rb_theta_expansion().get_n_outputs() == 0)
2112 : {
2113 142 : output_dual_innerprods_computed = true;
2114 142 : return;
2115 : }
2116 :
2117 : // Only log if we get to here
2118 12 : LOG_SCOPE("compute_output_dual_innerprods()", "RBConstruction");
2119 :
2120 4 : libMesh::out << "Compute output dual inner products" << std::endl;
2121 :
2122 : // Find out the largest value of Q_l
2123 4 : unsigned int max_Q_l = 0;
2124 710 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2125 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;
2126 :
2127 142 : std::vector<std::unique_ptr<NumericVector<Number>>> L_q_representor(max_Q_l);
2128 284 : for (unsigned int q=0; q<max_Q_l; q++)
2129 : {
2130 142 : L_q_representor[q] = NumericVector<Number>::build(this->comm());
2131 146 : L_q_representor[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2132 : }
2133 :
2134 710 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2135 : {
2136 1136 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
2137 : {
2138 568 : rhs->zero();
2139 568 : rhs->add(1., *get_output_vector(n,q_l));
2140 :
2141 568 : if (!is_quiet())
2142 0 : libMesh::out << "Starting solve n=" << n << ", q_l=" << q_l
2143 0 : << " in RBConstruction::compute_output_dual_innerprods() at "
2144 0 : << Utility::get_timestamp() << std::endl;
2145 :
2146 : // Use the main linear solver here instead of the inner_product solver, since
2147 : // get_matrix_for_output_dual_solves() may not return the inner product matrix.
2148 568 : solve_for_matrix_and_rhs(*get_linear_solver(), get_matrix_for_output_dual_solves(), *rhs);
2149 :
2150 : // We possibly perform multiple solves here with the same matrix, hence
2151 : // set reuse_preconditioner(true) (and set it back to false again below
2152 : // at the end of this function).
2153 568 : linear_solver->reuse_preconditioner(true);
2154 :
2155 568 : if (assert_convergence)
2156 568 : check_convergence(*get_linear_solver());
2157 :
2158 568 : if (!is_quiet())
2159 : {
2160 0 : libMesh::out << "Finished solve n=" << n << ", q_l=" << q_l
2161 0 : << " in RBConstruction::compute_output_dual_innerprods() at "
2162 0 : << Utility::get_timestamp() << std::endl;
2163 :
2164 0 : libMesh::out << this->n_linear_iterations()
2165 0 : << " iterations, final residual "
2166 0 : << this->final_linear_residual() << std::endl;
2167 : }
2168 :
2169 600 : *L_q_representor[q_l] = *solution;
2170 : }
2171 :
2172 16 : unsigned int q=0;
2173 1136 : for (unsigned int q_l1=0; q_l1<get_rb_theta_expansion().get_n_output_terms(n); q_l1++)
2174 : {
2175 568 : get_matrix_for_output_dual_solves().vector_mult(*inner_product_storage_vector, *L_q_representor[q_l1]);
2176 :
2177 1136 : for (unsigned int q_l2=q_l1; q_l2<get_rb_theta_expansion().get_n_output_terms(n); q_l2++)
2178 : {
2179 600 : output_dual_innerprods[n][q] = L_q_representor[q_l2]->dot(*inner_product_storage_vector);
2180 64 : libMesh::out << "output_dual_innerprods[" << n << "][" << q << "] = " << output_dual_innerprods[n][q] << std::endl;
2181 :
2182 568 : q++;
2183 : }
2184 : }
2185 : }
2186 :
2187 : // We may not need to use linear_solver again (e.g. this would happen if we use
2188 : // extra_linear_solver for the truth_solves). As a result, let's clear linear_solver
2189 : // to release any memory it may be taking up. If we do need it again, it will
2190 : // be initialized when necessary.
2191 142 : linear_solver->clear();
2192 142 : linear_solver->reuse_preconditioner(false);
2193 :
2194 142 : output_dual_innerprods_computed = true;
2195 134 : }
2196 :
2197 142 : get_rb_evaluation().output_dual_innerprods = output_dual_innerprods;
2198 : }
2199 :
2200 284 : void RBConstruction::compute_Fq_representor_innerprods(bool compute_inner_products)
2201 : {
2202 :
2203 : // Skip calculations if we've already computed the Fq_representors
2204 284 : if (!Fq_representor_innerprods_computed)
2205 : {
2206 : // Only log if we get to here
2207 16 : LOG_SCOPE("compute_Fq_representor_innerprods()", "RBConstruction");
2208 :
2209 1278 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2210 : {
2211 1022 : if (!Fq_representor[q_f])
2212 : {
2213 1932 : Fq_representor[q_f] = NumericVector<Number>::build(this->comm());
2214 1022 : Fq_representor[q_f]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2215 : }
2216 :
2217 28 : libmesh_assert(Fq_representor[q_f]->size() == this->n_dofs() &&
2218 : Fq_representor[q_f]->local_size() == this->n_local_dofs() );
2219 :
2220 994 : rhs->zero();
2221 994 : rhs->add(1., *get_Fq(q_f));
2222 :
2223 994 : if (!is_quiet())
2224 0 : libMesh::out << "Starting solve q_f=" << q_f
2225 0 : << " in RBConstruction::update_residual_terms() at "
2226 0 : << Utility::get_timestamp() << std::endl;
2227 :
2228 1022 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2229 :
2230 994 : if (assert_convergence)
2231 994 : check_convergence(*inner_product_solver);
2232 :
2233 994 : if (!is_quiet())
2234 : {
2235 0 : libMesh::out << "Finished solve q_f=" << q_f
2236 0 : << " in RBConstruction::update_residual_terms() at "
2237 0 : << Utility::get_timestamp() << std::endl;
2238 :
2239 0 : libMesh::out << this->n_linear_iterations()
2240 0 : << " iterations, final residual "
2241 0 : << this->final_linear_residual() << std::endl;
2242 : }
2243 :
2244 1050 : *Fq_representor[q_f] = *solution;
2245 : }
2246 :
2247 284 : if (compute_inner_products)
2248 : {
2249 8 : unsigned int q=0;
2250 :
2251 1278 : for (unsigned int q_f1=0; q_f1<get_rb_theta_expansion().get_n_F_terms(); q_f1++)
2252 : {
2253 994 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *Fq_representor[q_f1]);
2254 :
2255 4118 : for (unsigned int q_f2=q_f1; q_f2<get_rb_theta_expansion().get_n_F_terms(); q_f2++)
2256 : {
2257 3212 : Fq_representor_innerprods[q] = inner_product_storage_vector->dot(*Fq_representor[q_f2]);
2258 :
2259 3124 : q++;
2260 : }
2261 : }
2262 : } // end if (compute_inner_products)
2263 :
2264 284 : Fq_representor_innerprods_computed = true;
2265 : }
2266 :
2267 284 : get_rb_evaluation().Fq_representor_innerprods = Fq_representor_innerprods;
2268 284 : }
2269 :
2270 214 : void RBConstruction::load_rb_solution()
2271 : {
2272 18 : LOG_SCOPE("load_rb_solution()", "RBConstruction");
2273 :
2274 214 : solution->zero();
2275 :
2276 214 : libmesh_error_msg_if(get_rb_evaluation().RB_solution.size() > get_rb_evaluation().get_n_basis_functions(),
2277 : "ERROR: System contains " << get_rb_evaluation().get_n_basis_functions() << " basis functions."
2278 : << " RB_solution vector contains " << get_rb_evaluation().RB_solution.size() << " entries."
2279 : << " RB_solution in RBConstruction::load_rb_solution is too long!");
2280 :
2281 3776 : for (auto i : make_range(get_rb_evaluation().RB_solution.size()))
2282 3556 : solution->add(get_rb_evaluation().RB_solution(i), get_rb_evaluation().get_basis_function(i));
2283 :
2284 214 : update();
2285 214 : }
2286 :
2287 : // The slow (but simple, non-error prone) way to compute the residual dual norm
2288 : // Useful for error checking
2289 0 : Real RBConstruction::compute_residual_dual_norm_slow(const unsigned int N)
2290 : {
2291 0 : LOG_SCOPE("compute_residual_dual_norm_slow()", "RBConstruction");
2292 :
2293 : // Put the residual in rhs in order to compute the norm of the Riesz representor
2294 : // Note that this only works in serial since otherwise each processor will
2295 : // have a different parameter value during the Greedy training.
2296 :
2297 0 : std::unique_ptr<NumericVector<Number>> RB_sol = NumericVector<Number>::build(comm());
2298 0 : RB_sol->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2299 :
2300 0 : std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(comm());
2301 0 : temp->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2302 :
2303 0 : if (store_untransformed_basis)
2304 : {
2305 0 : libmesh_assert_equal_to(_untransformed_basis_functions.size(), get_rb_evaluation().get_n_basis_functions());
2306 : }
2307 :
2308 0 : for (unsigned int i=0; i<N; i++)
2309 : {
2310 0 : if (!store_untransformed_basis)
2311 : {
2312 0 : RB_sol->add(get_rb_evaluation().RB_solution(i), get_rb_evaluation().get_basis_function(i));
2313 : }
2314 : else
2315 : {
2316 0 : RB_sol->add(get_rb_evaluation().RB_solution(i), *_untransformed_basis_functions[i]);
2317 : }
2318 : }
2319 :
2320 0 : this->truth_assembly();
2321 0 : matrix->vector_mult(*temp, *RB_sol);
2322 0 : rhs->add(-1., *temp);
2323 :
2324 : // Then solve to get the Reisz representor
2325 0 : matrix->zero();
2326 0 : matrix->add(1., *inner_product_matrix);
2327 :
2328 0 : solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);
2329 0 : get_non_dirichlet_inner_product_matrix_if_avail()->vector_mult(*inner_product_storage_vector, *solution);
2330 0 : Number slow_residual_norm_sq = solution->dot(*inner_product_storage_vector);
2331 :
2332 0 : return std::sqrt( libmesh_real(slow_residual_norm_sq) );
2333 0 : }
2334 :
2335 1059389 : SparseMatrix<Number> * RBConstruction::get_inner_product_matrix()
2336 : {
2337 1059389 : libmesh_error_msg_if(!store_dirichlet_operators,
2338 : "Error: Must have store_dirichlet_operators==true to access inner_product_matrix.");
2339 1059389 : return inner_product_matrix.get();
2340 : }
2341 :
2342 28400 : const SparseMatrix<Number> * RBConstruction::get_inner_product_matrix() const
2343 : {
2344 28400 : libmesh_error_msg_if(!store_dirichlet_operators,
2345 : "Error: Must have store_dirichlet_operators==true to access inner_product_matrix.");
2346 28400 : return inner_product_matrix.get();
2347 : }
2348 :
2349 0 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix()
2350 : {
2351 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2352 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_inner_product_matrix.");
2353 :
2354 0 : return non_dirichlet_inner_product_matrix.get();
2355 : }
2356 :
2357 0 : const SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix() const
2358 : {
2359 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2360 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_inner_product_matrix.");
2361 :
2362 0 : return non_dirichlet_inner_product_matrix.get();
2363 : }
2364 :
2365 1059388 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix_if_avail()
2366 : {
2367 1059388 : if (store_non_dirichlet_operators)
2368 : {
2369 0 : return get_non_dirichlet_inner_product_matrix();
2370 : }
2371 :
2372 1059388 : return get_inner_product_matrix();
2373 : }
2374 :
2375 28400 : const SparseMatrix<Number> * RBConstruction::get_non_dirichlet_inner_product_matrix_if_avail() const
2376 : {
2377 28400 : if (store_non_dirichlet_operators)
2378 : {
2379 0 : return get_non_dirichlet_inner_product_matrix();
2380 : }
2381 :
2382 28400 : return get_inner_product_matrix();
2383 : }
2384 :
2385 1130697 : SparseMatrix<Number> * RBConstruction::get_Aq(unsigned int q)
2386 : {
2387 1130697 : libmesh_error_msg_if(!store_dirichlet_operators,
2388 : "Error: Must have store_dirichlet_operators==true to access non_dirichlet_Aq.");
2389 :
2390 1130697 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
2391 : "Error: We must have q < Q_a in get_Aq.");
2392 :
2393 1162917 : return Aq_vector[q].get();
2394 : }
2395 :
2396 0 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_Aq(unsigned int q)
2397 : {
2398 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2399 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Aq.");
2400 :
2401 0 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
2402 : "Error: We must have q < Q_a in get_Aq.");
2403 :
2404 0 : return non_dirichlet_Aq_vector[q].get();
2405 : }
2406 :
2407 265158 : SparseMatrix<Number> * RBConstruction::get_non_dirichlet_Aq_if_avail(unsigned int q)
2408 : {
2409 265158 : if (store_non_dirichlet_operators)
2410 : {
2411 0 : return get_non_dirichlet_Aq(q);
2412 : }
2413 :
2414 265158 : return get_Aq(q);
2415 : }
2416 :
2417 171800 : NumericVector<Number> * RBConstruction::get_Fq(unsigned int q)
2418 : {
2419 171800 : libmesh_error_msg_if(!store_dirichlet_operators,
2420 : "Error: Must have store_dirichlet_operators==true to access non_dirichlet_Fq.");
2421 :
2422 171800 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_F_terms(),
2423 : "Error: We must have q < Q_f in get_Fq.");
2424 :
2425 176696 : return Fq_vector[q].get();
2426 : }
2427 :
2428 0 : NumericVector<Number> * RBConstruction::get_non_dirichlet_Fq(unsigned int q)
2429 : {
2430 0 : libmesh_error_msg_if(!store_non_dirichlet_operators,
2431 : "Error: Must have store_non_dirichlet_operators==true to access non_dirichlet_Fq.");
2432 :
2433 0 : libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_F_terms(),
2434 : "Error: We must have q < Q_f in get_Fq.");
2435 :
2436 0 : return non_dirichlet_Fq_vector[q].get();
2437 : }
2438 :
2439 15606 : NumericVector<Number> * RBConstruction::get_non_dirichlet_Fq_if_avail(unsigned int q)
2440 : {
2441 15606 : if (store_non_dirichlet_operators)
2442 : {
2443 0 : return get_non_dirichlet_Fq(q);
2444 : }
2445 :
2446 15606 : return get_Fq(q);
2447 : }
2448 :
2449 584312 : NumericVector<Number> * RBConstruction::get_output_vector(unsigned int n, unsigned int q_l)
2450 : {
2451 584312 : libmesh_error_msg_if((n >= get_rb_theta_expansion().get_n_outputs()) ||
2452 : (q_l >= get_rb_theta_expansion().get_n_output_terms(n)),
2453 : "Error: We must have n < n_outputs and "
2454 : "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_output_vector.");
2455 :
2456 617688 : return outputs_vector[n][q_l].get();
2457 : }
2458 :
2459 0 : NumericVector<Number> * RBConstruction::get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l)
2460 : {
2461 0 : libmesh_error_msg_if((n >= get_rb_theta_expansion().get_n_outputs()) ||
2462 : (q_l >= get_rb_theta_expansion().get_n_output_terms(n)),
2463 : "Error: We must have n < n_outputs and "
2464 : "q_l < get_rb_theta_expansion().get_n_output_terms(n) in get_non_dirichlet_output_vector.");
2465 :
2466 0 : return non_dirichlet_outputs_vector[n][q_l].get();
2467 : }
2468 :
2469 0 : void RBConstruction::get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices)
2470 : {
2471 0 : all_matrices.clear();
2472 :
2473 0 : if (store_dirichlet_operators)
2474 0 : all_matrices["inner_product"] = get_inner_product_matrix();
2475 :
2476 0 : if (store_non_dirichlet_operators)
2477 0 : all_matrices["non_dirichlet_inner_product"] = get_non_dirichlet_inner_product_matrix();
2478 :
2479 0 : for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
2480 : {
2481 0 : std::stringstream matrix_name;
2482 0 : matrix_name << "A" << q_a;
2483 :
2484 0 : if (store_dirichlet_operators)
2485 0 : all_matrices[matrix_name.str()] = get_Aq(q_a);
2486 :
2487 0 : if (store_non_dirichlet_operators)
2488 : {
2489 0 : matrix_name << "_non_dirichlet";
2490 0 : all_matrices[matrix_name.str()] = get_non_dirichlet_Aq(q_a);
2491 : }
2492 0 : }
2493 0 : }
2494 :
2495 0 : void RBConstruction::get_all_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors)
2496 : {
2497 0 : all_vectors.clear();
2498 :
2499 0 : get_output_vectors(all_vectors);
2500 :
2501 0 : for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
2502 : {
2503 0 : std::stringstream F_vector_name;
2504 0 : F_vector_name << "F" << q_f;
2505 :
2506 0 : if (store_dirichlet_operators)
2507 0 : all_vectors[F_vector_name.str()] = get_Fq(q_f);
2508 :
2509 0 : if (store_non_dirichlet_operators)
2510 : {
2511 0 : F_vector_name << "_non_dirichlet";
2512 0 : all_vectors[F_vector_name.str()] = get_non_dirichlet_Fq(q_f);
2513 : }
2514 0 : }
2515 0 : }
2516 :
2517 0 : void RBConstruction::get_output_vectors(std::map<std::string, NumericVector<Number> *> & output_vectors)
2518 : {
2519 0 : output_vectors.clear();
2520 :
2521 0 : for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
2522 0 : for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
2523 : {
2524 0 : std::stringstream output_name;
2525 0 : if (store_dirichlet_operators)
2526 : {
2527 0 : output_name << "output_" << n << "_"<< q_l;
2528 0 : output_vectors[output_name.str()] = get_output_vector(n,q_l);
2529 : }
2530 :
2531 0 : if (store_non_dirichlet_operators)
2532 : {
2533 0 : output_name << "_non_dirichlet";
2534 0 : output_vectors[output_name.str()] = get_non_dirichlet_output_vector(n,q_l);
2535 : }
2536 0 : }
2537 0 : }
2538 :
2539 : #ifdef LIBMESH_ENABLE_DIRICHLET
2540 :
2541 568 : std::unique_ptr<DirichletBoundary> RBConstruction::build_zero_dirichlet_boundary_object()
2542 : {
2543 32 : ZeroFunction<> zf;
2544 :
2545 32 : std::set<boundary_id_type> dirichlet_ids;
2546 16 : std::vector<unsigned int> variables;
2547 :
2548 : // The DirichletBoundary constructor clones zf, so it's OK that zf is only in local scope
2549 1136 : return std::make_unique<DirichletBoundary>(dirichlet_ids, variables, &zf);
2550 : }
2551 :
2552 : #endif
2553 :
2554 0 : void RBConstruction::write_riesz_representors_to_files(const std::string & riesz_representors_dir,
2555 : const bool write_binary_residual_representors)
2556 : {
2557 0 : LOG_SCOPE("write_riesz_representors_to_files()", "RBConstruction");
2558 :
2559 : // Write out Riesz representors. These are useful to have when restarting,
2560 : // so you don't have to recompute them all over again.
2561 :
2562 : // First we write out the Fq representors, these are independent of an RBEvaluation object.
2563 0 : libMesh::out << "Writing out the Fq_representors..." << std::endl;
2564 :
2565 0 : std::ostringstream file_name;
2566 0 : const std::string riesz_representor_suffix = (write_binary_residual_representors ? ".xdr" : ".dat");
2567 : struct stat stat_info;
2568 :
2569 : // Residual representors written out to their own separate directory
2570 0 : if ( this->processor_id() == 0)
2571 0 : if ( Utility::mkdir(riesz_representors_dir.c_str()) != 0)
2572 0 : libMesh::out << "Skipping creating residual_representors directory: " << strerror(errno) << std::endl;
2573 :
2574 0 : for (auto i : index_range(Fq_representor))
2575 : {
2576 0 : if (Fq_representor[i])
2577 : {
2578 0 : file_name.str(""); // reset filename
2579 0 : file_name << riesz_representors_dir << "/Fq_representor" << i << riesz_representor_suffix;
2580 :
2581 : // Check to see if file exists, if so, don't overwrite it, we assume it was
2582 : // there from a previous call to this function. Note: if stat returns zero
2583 : // it means it successfully got the file attributes (and therefore the file
2584 : // exists). Because of the following factors:
2585 : // 1.) write_serialized_data takes longer for proc 0 than it does for others,
2586 : // so processors can get out of sync.
2587 : // 2.) The constructor for Xdr opens a (0 length) file on *all* processors,
2588 : // there are typically hundreds of 0-length files created during this loop,
2589 : // and that screws up checking for the existence of files. One way to stay
2590 : // in sync is to but a barrier at each iteration of the loop -- not sure how
2591 : // bad this will affect performance, but it can't be much worse than serialized
2592 : // I/O already is :)
2593 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2594 :
2595 0 : if ( (stat_result != 0) || // file definitely doesn't already exist
2596 0 : (stat_info.st_size == 0)) // file exists, but has zero length (can happen if another proc already opened it!)
2597 : {
2598 : // No need to copy!
2599 : // *solution = *(Fq_representor[i]);
2600 : // std::swap doesn't work on pointers
2601 : //std::swap(solution.get(), Fq_representor[i]);
2602 0 : Fq_representor[i]->swap(*solution);
2603 :
2604 0 : Xdr fqr_data(file_name.str(),
2605 0 : write_binary_residual_representors ? ENCODE : WRITE);
2606 :
2607 0 : write_serialized_data(fqr_data, false);
2608 :
2609 : // Synchronize before moving on
2610 0 : this->comm().barrier();
2611 :
2612 : // Swap back.
2613 0 : Fq_representor[i]->swap(*solution);
2614 :
2615 : // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2616 : // for the system call, be sure to do it only on one processor, etc.
2617 0 : }
2618 : }
2619 : }
2620 :
2621 :
2622 : // Next, write out the Aq representors associated with rb_eval.
2623 0 : libMesh::out << "Writing out the Aq_representors..." << std::endl;
2624 :
2625 0 : const unsigned int jstop = get_rb_evaluation().get_n_basis_functions();
2626 0 : const unsigned int jstart = jstop-get_delta_N();
2627 :
2628 0 : RBEvaluation & rbe = get_rb_evaluation();
2629 0 : for (auto i : index_range(rbe.Aq_representor))
2630 0 : for (unsigned int j=jstart; j<jstop; ++j)
2631 : {
2632 0 : libMesh::out << "Writing out Aq_representor[" << i << "][" << j << "]..." << std::endl;
2633 0 : libmesh_assert(get_rb_evaluation().Aq_representor[i][j]);
2634 :
2635 0 : file_name.str(""); // reset filename
2636 : file_name << riesz_representors_dir
2637 0 : << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2638 :
2639 : {
2640 : // No need to copy! Use swap instead.
2641 : // *solution = *(Aq_representor[i][j]);
2642 0 : rbe.Aq_representor[i][j]->swap(*solution);
2643 :
2644 0 : Xdr aqr_data(file_name.str(),
2645 0 : write_binary_residual_representors ? ENCODE : WRITE);
2646 :
2647 0 : write_serialized_data(aqr_data, false);
2648 :
2649 : // Synchronize before moving on
2650 0 : this->comm().barrier();
2651 :
2652 : // Swap back.
2653 0 : rbe.Aq_representor[i][j]->swap(*solution);
2654 :
2655 : // TODO: bzip the resulting file? See $LIBMESH_DIR/src/mesh/unstructured_mesh.C
2656 : // for the system call, be sure to do it only on one processor, etc.
2657 0 : }
2658 : }
2659 0 : }
2660 :
2661 :
2662 :
2663 0 : void RBConstruction::read_riesz_representors_from_files(const std::string & riesz_representors_dir,
2664 : const bool read_binary_residual_representors)
2665 : {
2666 0 : LOG_SCOPE("read_riesz_representors_from_files()", "RBConstruction");
2667 :
2668 0 : libMesh::out << "Reading in the Fq_representors..." << std::endl;
2669 :
2670 0 : const std::string riesz_representor_suffix = (read_binary_residual_representors ? ".xdr" : ".dat");
2671 0 : std::ostringstream file_name;
2672 : struct stat stat_info;
2673 :
2674 : // Read in the Fq_representors. There should be Q_f of these. FIXME:
2675 : // should we be worried about leaks here?
2676 0 : for (const auto & rep : Fq_representor)
2677 0 : libmesh_error_msg_if(rep, "Error, must delete existing Fq_representor before reading in from file.");
2678 :
2679 0 : for (auto i : index_range(Fq_representor))
2680 : {
2681 0 : file_name.str(""); // reset filename
2682 : file_name << riesz_representors_dir
2683 0 : << "/Fq_representor" << i << riesz_representor_suffix;
2684 :
2685 : // On processor zero check to be sure the file exists
2686 0 : if (this->processor_id() == 0)
2687 : {
2688 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2689 :
2690 0 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
2691 : }
2692 :
2693 0 : Xdr fqr_data(file_name.str(),
2694 0 : read_binary_residual_representors ? DECODE : READ);
2695 :
2696 0 : read_serialized_data(fqr_data, false);
2697 :
2698 0 : Fq_representor[i] = NumericVector<Number>::build(this->comm());
2699 0 : Fq_representor[i]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
2700 :
2701 : // No need to copy, just swap
2702 : // *Fq_representor[i] = *solution;
2703 0 : Fq_representor[i]->swap(*solution);
2704 0 : }
2705 :
2706 : // Alert the update_residual_terms() function that we don't need to recompute
2707 : // the Fq_representors as we have already read them in from file!
2708 0 : Fq_representor_innerprods_computed = true;
2709 :
2710 :
2711 0 : libMesh::out << "Reading in the Aq_representors..." << std::endl;
2712 :
2713 : // Read in the Aq representors. The class makes room for [Q_a][Nmax] of these. We are going to
2714 : // read in [Q_a][get_rb_evaluation().get_n_basis_functions()]. FIXME:
2715 : // should we be worried about leaks in the locations where we're about to fill entries?
2716 0 : RBEvaluation & rbe = get_rb_evaluation();
2717 0 : for (const auto & row : rbe.Aq_representor)
2718 0 : for (const auto & rep : row)
2719 0 : libmesh_error_msg_if(rep, "Error, must delete existing Aq_representor before reading in from file.");
2720 :
2721 : // Now ready to read them in from file!
2722 0 : for (auto i : index_range(rbe.Aq_representor))
2723 0 : for (unsigned int j=0; j<rbe.get_n_basis_functions(); ++j)
2724 : {
2725 0 : file_name.str(""); // reset filename
2726 : file_name << riesz_representors_dir
2727 0 : << "/Aq_representor" << i << "_" << j << riesz_representor_suffix;
2728 :
2729 : // On processor zero check to be sure the file exists
2730 0 : if (this->processor_id() == 0)
2731 : {
2732 0 : int stat_result = stat(file_name.str().c_str(), &stat_info);
2733 :
2734 0 : libmesh_error_msg_if(stat_result != 0, "File does not exist: " << file_name.str());
2735 : }
2736 :
2737 0 : Xdr aqr_data(file_name.str(), read_binary_residual_representors ? DECODE : READ);
2738 :
2739 0 : read_serialized_data(aqr_data, false);
2740 :
2741 0 : rbe.Aq_representor[i][j] = NumericVector<Number>::build(this->comm());
2742 0 : rbe.Aq_representor[i][j]->init (n_dofs(), n_local_dofs(),
2743 0 : false, PARALLEL);
2744 :
2745 : // No need to copy, just swap
2746 0 : rbe.Aq_representor[i][j]->swap(*solution);
2747 0 : }
2748 0 : }
2749 :
2750 160467 : void RBConstruction::check_convergence(LinearSolver<Number> & input_solver)
2751 : {
2752 : libMesh::LinearConvergenceReason conv_flag;
2753 :
2754 160467 : conv_flag = input_solver.get_converged_reason();
2755 :
2756 160467 : libmesh_error_msg_if(conv_flag < 0, "Convergence error. Error id: " << conv_flag);
2757 160467 : }
2758 :
2759 0 : bool RBConstruction::get_convergence_assertion_flag() const
2760 : {
2761 0 : return assert_convergence;
2762 : }
2763 :
2764 0 : void RBConstruction::set_convergence_assertion_flag(bool flag)
2765 : {
2766 0 : assert_convergence = flag;
2767 0 : }
2768 :
2769 813684 : bool RBConstruction::get_preevaluate_thetas_flag() const
2770 : {
2771 813684 : return _preevaluate_thetas_flag;
2772 : }
2773 :
2774 0 : void RBConstruction::set_preevaluate_thetas_flag(bool flag)
2775 : {
2776 0 : _preevaluate_thetas_flag = flag;
2777 0 : }
2778 :
2779 0 : unsigned int RBConstruction::get_current_training_parameter_index() const
2780 : {
2781 0 : return _current_training_parameter_index;
2782 : }
2783 :
2784 0 : void RBConstruction::set_current_training_parameter_index(unsigned int index)
2785 : {
2786 0 : _current_training_parameter_index = index;
2787 0 : }
2788 :
2789 : const std::vector<Number> &
2790 0 : RBConstruction::get_evaluated_thetas(unsigned int training_parameter_index) const
2791 : {
2792 0 : const numeric_index_type first_index = get_first_local_training_index();
2793 0 : libmesh_assert(training_parameter_index >= first_index);
2794 :
2795 0 : const numeric_index_type local_index = training_parameter_index - first_index;
2796 0 : libmesh_assert(local_index < _evaluated_thetas.size());
2797 :
2798 0 : return _evaluated_thetas[local_index];
2799 : }
2800 :
2801 0 : void RBConstruction::preevaluate_thetas()
2802 : {
2803 0 : LOG_SCOPE("preevaluate_thetas()", "RBConstruction");
2804 :
2805 0 : _evaluated_thetas.resize(get_local_n_training_samples());
2806 :
2807 : // Early return if we've already preevaluated thetas.
2808 0 : if (_preevaluate_thetas_completed)
2809 0 : return;
2810 :
2811 0 : if ( get_local_n_training_samples() == 0 )
2812 0 : return;
2813 :
2814 0 : auto & rb_theta_expansion = get_rb_evaluation().get_rb_theta_expansion();
2815 0 : const unsigned int n_A_terms = rb_theta_expansion.get_n_A_terms();
2816 0 : const unsigned int n_F_terms = rb_theta_expansion.get_n_F_terms();
2817 0 : const unsigned int n_outputs = rb_theta_expansion.get_total_n_output_terms();
2818 :
2819 : // Collect all training parameters
2820 : // TODO: Here instead of using a vector of RBParameters objects,
2821 : // we could use a single RBParameters object with multiple samples.
2822 : // This would save memory over the current approach, but that may
2823 : // not be a big deal in practice unless the number of training samples
2824 : // is very large for some reason.
2825 0 : std::vector<RBParameters> mus(get_local_n_training_samples());
2826 0 : const numeric_index_type first_index = get_first_local_training_index();
2827 0 : for (unsigned int i=0; i<get_local_n_training_samples(); i++)
2828 : {
2829 : // Load training parameter i, this is only loaded
2830 : // locally since the RB solves are local.
2831 0 : set_params_from_training_set( first_index+i );
2832 0 : mus[i] = get_parameters();
2833 0 : _evaluated_thetas[i].resize(n_A_terms + n_F_terms + n_outputs);
2834 : }
2835 :
2836 : // Evaluate thetas for all training parameters simultaneously
2837 0 : for (unsigned int q_a=0; q_a<n_A_terms; q_a++)
2838 : {
2839 0 : const auto A_vals = rb_theta_expansion.eval_A_theta(q_a, mus);
2840 0 : for (auto i : make_range(get_local_n_training_samples()))
2841 0 : _evaluated_thetas[i][q_a] = A_vals[i];
2842 : }
2843 :
2844 0 : for (unsigned int q_f=0; q_f<n_F_terms; q_f++)
2845 : {
2846 0 : const auto F_vals = rb_theta_expansion.eval_F_theta(q_f, mus);
2847 0 : for (auto i : make_range(get_local_n_training_samples()))
2848 0 : _evaluated_thetas[i][n_A_terms + q_f] = F_vals[i];
2849 : }
2850 :
2851 : {
2852 0 : unsigned int output_counter = 0;
2853 0 : for (unsigned int n=0; n<rb_theta_expansion.get_n_outputs(); n++)
2854 0 : for (unsigned int q_l=0; q_l<rb_theta_expansion.get_n_output_terms(n); q_l++)
2855 : {
2856 : // Evaluate the current output functional term for all
2857 : // training parameters simultaneously.
2858 0 : const auto output_vals = rb_theta_expansion.eval_output_theta(n, q_l, mus);
2859 :
2860 : // TODO: the size of _evaluated_thetas is currently assumed to be
2861 : // the same as get_local_n_training_samples(), but this won't be
2862 : // the case if we use RBParameters objects that have multiple samples.
2863 : // So just make sure that's the case for now.
2864 0 : libmesh_error_msg_if(output_vals.size() != get_local_n_training_samples(),
2865 : "We currently only support single-sample RBParameters "
2866 : "objects during the training stage.");
2867 :
2868 0 : for (auto i : index_range(output_vals))
2869 0 : _evaluated_thetas[i][n_A_terms + n_F_terms + output_counter] = output_vals[i];
2870 :
2871 : // Go to next output term
2872 0 : output_counter++;
2873 : }
2874 : }
2875 :
2876 0 : _preevaluate_thetas_completed = true;
2877 0 : }
2878 :
2879 0 : void RBConstruction::reset_preevaluate_thetas_completed()
2880 : {
2881 0 : _preevaluate_thetas_completed = false;
2882 0 : }
2883 :
2884 : } // namespace libMesh
|
