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 : #ifndef LIBMESH_RB_CONSTRUCTION_H
21 : #define LIBMESH_RB_CONSTRUCTION_H
22 :
23 : // rbOOmit includes
24 : #include "libmesh/rb_construction_base.h"
25 :
26 : // libMesh includes
27 : #include "libmesh/linear_implicit_system.h"
28 : #include "libmesh/dense_vector.h"
29 : #include "libmesh/dense_matrix.h"
30 : #include "libmesh/dg_fem_context.h"
31 : #include "libmesh/dirichlet_boundaries.h"
32 :
33 : // C++ includes
34 :
35 : namespace libMesh
36 : {
37 :
38 : class RBThetaExpansion;
39 : class RBAssemblyExpansion;
40 : class RBEvaluation;
41 : class ElemAssembly;
42 :
43 : /**
44 : * This class is part of the rbOOmit framework.
45 : *
46 : * RBConstruction implements the Construction stage
47 : * of the certified reduced basis method for
48 : * steady-state elliptic parametrized PDEs.
49 : *
50 : * \author David J. Knezevic
51 : * \date 2009
52 : */
53 96 : class RBConstruction : public RBConstructionBase<LinearImplicitSystem>
54 : {
55 : public:
56 :
57 : /**
58 : * Constructor. Optionally initializes required
59 : * data structures.
60 : */
61 : RBConstruction (EquationSystems & es,
62 : const std::string & name,
63 : const unsigned int number);
64 :
65 : /**
66 : * Special functions.
67 : * - This class has the same restrictions/defaults as its base class.
68 : * - Destructor is defaulted out-of-line
69 : */
70 : RBConstruction (RBConstruction &&) = default;
71 : RBConstruction (const RBConstruction &) = delete;
72 : RBConstruction & operator= (const RBConstruction &) = delete;
73 : RBConstruction & operator= (RBConstruction &&) = delete;
74 : virtual ~RBConstruction ();
75 :
76 : /**
77 : * The type of system.
78 : */
79 : typedef RBConstruction sys_type;
80 :
81 : /**
82 : * Assembles & solves the linear system A*x=b for the specified
83 : * matrix \p input_matrix and right-hand side \p rhs.
84 : */
85 : virtual void solve_for_matrix_and_rhs (LinearSolver<Number> & input_solver,
86 : SparseMatrix<Number> & input_matrix,
87 : NumericVector<Number> & input_rhs);
88 :
89 : /**
90 : * Set the RBEvaluation object.
91 : */
92 : void set_rb_evaluation(RBEvaluation & rb_eval_in);
93 :
94 : /**
95 : * Get a reference to the RBEvaluation object.
96 : */
97 : RBEvaluation & get_rb_evaluation();
98 : const RBEvaluation & get_rb_evaluation() const;
99 :
100 : /**
101 : * \returns \p true if rb_eval is initialized. False, otherwise.
102 : */
103 : bool is_rb_eval_initialized() const;
104 :
105 : /**
106 : * Get a reference to the RBThetaExpansion object that
107 : * that belongs to rb_eval.
108 : */
109 : RBThetaExpansion & get_rb_theta_expansion();
110 : const RBThetaExpansion & get_rb_theta_expansion() const;
111 :
112 : /**
113 : * Set the rb_assembly_expansion object.
114 : */
115 : void set_rb_assembly_expansion(RBAssemblyExpansion & rb_assembly_expansion_in);
116 :
117 : /**
118 : * \returns A reference to the rb_assembly_expansion object
119 : */
120 : RBAssemblyExpansion & get_rb_assembly_expansion();
121 :
122 : /**
123 : * \returns A reference to *this.
124 : */
125 : sys_type & system () { return *this; }
126 :
127 : /**
128 : * The type of the parent.
129 : */
130 : typedef RBConstructionBase<LinearImplicitSystem> Parent;
131 :
132 : /**
133 : * Clear all the data structures associated with
134 : * the system.
135 : */
136 : virtual void clear () override;
137 :
138 : /**
139 : * \returns A string indicating the type of the system.
140 : */
141 : virtual std::string system_type () const override;
142 :
143 : /**
144 : * Perform a "truth" solve, i.e. solve the finite element system at
145 : * at the parameters currently set in the system. This is used
146 : * extensively in training the reduced basis, since "truth snapshots"
147 : * are employed as basis functions.
148 : */
149 : virtual Real truth_solve(int plot_solution);
150 :
151 : /**
152 : * Train the reduced basis. This can use different approaches, e.g. Greedy
153 : * or POD, which are chosen using the RB_training_type member variable.
154 : *
155 : * In the case that we use Greedy training, this function returns the
156 : * final maximum a posteriori error bound on the training set.
157 : */
158 : virtual Real train_reduced_basis(const bool resize_rb_eval_data=true);
159 :
160 : /**
161 : * Train the reduced basis using the "Greedy algorithm."
162 : *
163 : * Each stage of the Greedy algorithm involves solving the reduced basis
164 : * over a large training set and selecting the parameter at which the
165 : * reduced basis error bound is largest, then performing a truth_solve
166 : * at that parameter and enriching the reduced basis with the corresponding
167 : * snapshot.
168 : *
169 : * \p resize_rb_eval_data is a boolean flag to indicate whether or not we
170 : * call rb_eval->resize_data_structures(Nmax). True by default, but we may
171 : * set it to false if, for example, we are continuing from a previous
172 : * training run and don't want to clobber the existing rb_eval data.
173 : *
174 : * \returns The final maximum a posteriori error bound on the training set.
175 : */
176 : Real train_reduced_basis_with_greedy(const bool resize_rb_eval_data);
177 :
178 : /**
179 : * This function computes one basis function for each rhs term. This is
180 : * useful in some cases since we can avoid doing a full greedy if we know
181 : * that we do not have any "left-hand side" parameters, for example.
182 : */
183 : void enrich_basis_from_rhs_terms(const bool resize_rb_eval_data=true);
184 :
185 : /**
186 : * Train the reduced basis using Proper Orthogonal Decomposition (POD).
187 : * This is an alternative to train_reduced_basis(), which uses the RB greedy
188 : * algorithm. In contrast to the RB greedy algorithm, POD requires us to
189 : * perform truth solves at all training samples, which can be computationally
190 : * intensive.
191 : *
192 : * The main advantage of using POD is that it does not rely on the RB error
193 : * indicator. The RB error indicator typically stagnates due to rounding
194 : * error at approximately square-root of machine precision, since it involves
195 : * taking the square-root of a sum of terms that cancel. This error indicator
196 : * stagnation puts a limit on the accuracy level that can be achieved with
197 : * the RB greedy algorithm, so for cases where we need higher accuracy, the
198 : * POD approach is a good alternative.
199 : */
200 : void train_reduced_basis_with_POD();
201 :
202 : /**
203 : * (i) Compute the a posteriori error bound for each set of parameters
204 : * in the training set, (ii) set current_parameters to the parameters that
205 : * maximize the error bound, and (iii) return the maximum error bound.
206 : */
207 : virtual Real compute_max_error_bound();
208 :
209 : /**
210 : * Return the parameters chosen during the i^th step of
211 : * the Greedy algorithm.
212 : */
213 : const RBParameters & get_greedy_parameter(unsigned int i);
214 :
215 : /**
216 : * Get/set the relative tolerance for the basis training.
217 : */
218 8 : void set_rel_training_tolerance(Real new_training_tolerance)
219 284 : {this->rel_training_tolerance = new_training_tolerance; }
220 276 : Real get_rel_training_tolerance() const { return rel_training_tolerance; }
221 :
222 : /**
223 : * Get/set the absolute tolerance for the basis training.
224 : */
225 8 : void set_abs_training_tolerance(Real new_training_tolerance)
226 284 : {this->abs_training_tolerance = new_training_tolerance; }
227 276 : Real get_abs_training_tolerance() const { return abs_training_tolerance; }
228 :
229 : /**
230 : * Get/set the boolean to indicate if we normalize the RB error in the greedy.
231 : */
232 8 : void set_normalize_rb_bound_in_greedy(bool normalize_rb_bound_in_greedy_in)
233 284 : {this->normalize_rb_bound_in_greedy = normalize_rb_bound_in_greedy_in; }
234 276 : bool get_normalize_rb_bound_in_greedy() const { return normalize_rb_bound_in_greedy; }
235 :
236 : /**
237 : * @return true if \p RB_training_type_in is a type of training that
238 : * requires a serial training set. For example, POD training generally
239 : * does require a serial training set.
240 : */
241 : virtual bool is_serial_training_type(const std::string & RB_training_type_in);
242 :
243 : /**
244 : * Get/set the string that determines the training type.
245 : */
246 : void set_RB_training_type(const std::string & RB_training_type_in);
247 : const std::string & get_RB_training_type() const;
248 :
249 : /**
250 : * Get/set Nmax, the maximum number of RB
251 : * functions we are willing to compute.
252 : */
253 10459 : unsigned int get_Nmax() const { return Nmax; }
254 : virtual void set_Nmax(unsigned int Nmax);
255 :
256 : /**
257 : * Load the i^th RB function into the RBConstruction
258 : * solution vector.
259 : */
260 : virtual void load_basis_function(unsigned int i);
261 :
262 : /**
263 : * Load the RB solution from the most recent solve with rb_eval
264 : * into this system's solution vector.
265 : */
266 : virtual void load_rb_solution();
267 :
268 : /**
269 : * The slow (but simple, non-error prone) way to compute the residual dual norm.
270 : * Useful for error checking.
271 : */
272 : Real compute_residual_dual_norm_slow(const unsigned int N);
273 :
274 : /**
275 : * Get a pointer to inner_product_matrix. Accessing via this
276 : * function, rather than directly through the class member allows
277 : * us to do error checking (e.g. inner_product_matrix is not
278 : * defined in low-memory mode).
279 : */
280 : SparseMatrix<Number> * get_inner_product_matrix();
281 : const SparseMatrix<Number> * get_inner_product_matrix() const;
282 :
283 : /**
284 : * Get the non-Dirichlet (or more generally no-constraints) version
285 : * of the inner-product matrix. This is useful for performing multiplications
286 : * on vectors that already have constraints enforced.
287 : */
288 : SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix();
289 : const SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix() const;
290 :
291 : /**
292 : * Get the non-Dirichlet inner-product matrix if it's available,
293 : * otherwise get the inner-product matrix with constraints.
294 : */
295 : SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix_if_avail();
296 : const SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix_if_avail() const;
297 :
298 : /**
299 : * Get a pointer to Aq.
300 : */
301 : SparseMatrix<Number> * get_Aq(unsigned int q);
302 :
303 : /**
304 : * Get a pointer to non_dirichlet_Aq.
305 : */
306 : SparseMatrix<Number> * get_non_dirichlet_Aq(unsigned int q);
307 :
308 : /**
309 : * Get a pointer to non_dirichlet_Aq if it's available, otherwise
310 : * get Aq.
311 : */
312 : SparseMatrix<Number> * get_non_dirichlet_Aq_if_avail(unsigned int q);
313 :
314 : /**
315 : * Allocate all the data structures necessary for the construction
316 : * stage of the RB method. This function also performs
317 : * matrix and vector assembly of the "truth" affine expansion.
318 : *
319 : * We can optionally skip the matrix or vector assembly steps by setting
320 : * skip_matrix_assembly = true, or skip_vector_assembly = true,
321 : * respectively.
322 : */
323 : virtual void initialize_rb_construction(bool skip_matrix_assembly=false,
324 : bool skip_vector_assembly=false);
325 :
326 : /**
327 : * Get a pointer to Fq.
328 : */
329 : NumericVector<Number> * get_Fq(unsigned int q);
330 :
331 : /**
332 : * Get a pointer to non-Dirichlet Fq.
333 : */
334 : NumericVector<Number> * get_non_dirichlet_Fq(unsigned int q);
335 :
336 : /**
337 : * Get a pointer to non_dirichlet_Fq if it's available, otherwise
338 : * get Fq.
339 : */
340 : NumericVector<Number> * get_non_dirichlet_Fq_if_avail(unsigned int q);
341 :
342 : /**
343 : * Get a pointer to the n^th output.
344 : */
345 : NumericVector<Number> * get_output_vector(unsigned int n, unsigned int q_l);
346 :
347 : /**
348 : * Get a pointer to non-Dirichlet output vector.
349 : */
350 : NumericVector<Number> * get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l);
351 :
352 : /**
353 : * Get a map that stores pointers to all of the matrices.
354 : */
355 : virtual void get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices);
356 :
357 : /**
358 : * Get a map that stores pointers to all of the vectors.
359 : */
360 : virtual void get_all_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors);
361 :
362 : /**
363 : * Get a map that stores pointers to all of the vectors.
364 : */
365 : virtual void get_output_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors);
366 :
367 : /**
368 : * Assemble the matrices and vectors for this system.
369 : * Optionally skip matrix or vector assembly (e.g. we may want to
370 : * read data in from disk instead).
371 : */
372 : virtual void assemble_affine_expansion(bool skip_matrix_assembly,
373 : bool skip_vector_assembly);
374 :
375 : /**
376 : * Assemble the inner product matrix and store it in input_matrix.
377 : */
378 : void assemble_inner_product_matrix(SparseMatrix<Number> * input_matrix, bool apply_dof_constraints=true);
379 :
380 : /**
381 : * Assemble the q^th affine matrix and store it in input_matrix.
382 : */
383 : void assemble_Aq_matrix(unsigned int q, SparseMatrix<Number> * input_matrix, bool apply_dof_constraints=true);
384 :
385 : /**
386 : * Assemble the q^th affine vector and store it in input_matrix.
387 : */
388 : void assemble_Fq_vector(unsigned int q, NumericVector<Number> * input_vector, bool apply_dof_constraints=true);
389 :
390 : /**
391 : * Add the scaled q^th affine matrix to input_matrix. If symmetrize==true, then
392 : * we symmetrize Aq before adding it.
393 : */
394 : void add_scaled_Aq(Number scalar, unsigned int q_a,
395 : SparseMatrix<Number> * input_matrix,
396 : bool symmetrize);
397 :
398 : /**
399 : * Write out all the Riesz representor data to files.
400 : * \p residual_representors_dir specifies which directory to write to.
401 : * \p write_binary_residual_representors specifies whether we write binary or ASCII data.
402 : */
403 : virtual void write_riesz_representors_to_files(const std::string & riesz_representors_dir,
404 : const bool write_binary_residual_representors);
405 :
406 : /**
407 : * Read in all the Riesz representor data from files.
408 : * \p directory_name specifies which directory to read from.
409 : * \p io_flag specifies whether we read in all data, or only
410 : * a basis (in)dependent subset.
411 : */
412 : virtual void read_riesz_representors_from_files(const std::string & riesz_representors_dir,
413 : const bool write_binary_residual_representors);
414 :
415 :
416 : /**
417 : * This function computes all of the residual representors, can be useful
418 : * when restarting a basis training computation.
419 : * If \p compute_inner_products is false, we just compute the residual Riesz
420 : * representors, whereas if true, we also compute all the corresponding inner
421 : * product terms.
422 : */
423 : virtual void recompute_all_residual_terms(const bool compute_inner_products=true);
424 :
425 : /**
426 : * Read in from the file specified by \p parameters_filename
427 : * and set the this system's member variables accordingly.
428 : */
429 : virtual void process_parameters_file(const std::string & parameters_filename);
430 :
431 : /**
432 : * Set the state of this RBConstruction object based on the arguments
433 : * to this function.
434 : */
435 : void set_rb_construction_parameters(unsigned int n_training_samples_in,
436 : bool deterministic_training_in,
437 : int training_parameters_random_seed_in,
438 : bool quiet_mode_in,
439 : unsigned int Nmax_in,
440 : Real rel_training_tolerance_in,
441 : Real abs_training_tolerance_in,
442 : bool normalize_rb_error_bound_in_greedy_in,
443 : const std::string & RB_training_type_in,
444 : const RBParameters & mu_min_in,
445 : const RBParameters & mu_max_in,
446 : const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
447 : const std::map<std::string,bool> & log_scaling,
448 : std::map<std::string, std::vector<RBParameter>> * training_sample_list=nullptr);
449 :
450 : /**
451 : * Print out info that describes the current setup of this RBConstruction.
452 : */
453 : virtual void print_info() const;
454 :
455 : /**
456 : * Print out a matrix that shows the orthogonality of the RB basis functions.
457 : * This is a helpful debugging tool, e.g. orthogonality can be degraded
458 : * due to finite precision arithmetic.
459 : */
460 : void print_basis_function_orthogonality() const;
461 :
462 : /**
463 : * Get delta_N, the number of basis functions we
464 : * add to the RB space per iteration of the greedy
465 : * algorithm. For steady-state systems, this should
466 : * be 1, but can be more than 1 for time-dependent
467 : * systems.
468 : */
469 2788 : unsigned int get_delta_N() const { return delta_N; }
470 :
471 : /**
472 : * Set the rb_assembly_expansion object.
473 : */
474 : void set_inner_product_assembly(ElemAssembly & inner_product_assembly_in);
475 :
476 : /**
477 : * \returns A reference to the inner product assembly object
478 : */
479 : ElemAssembly & get_inner_product_assembly();
480 :
481 : /**
482 : * Specify the coefficients of the A_q operators to be used in the
483 : * energy inner-product.
484 : */
485 : void set_energy_inner_product(const std::vector<Number> & energy_inner_product_coeffs_in);
486 :
487 : /**
488 : * It is sometimes useful to be able to zero vector entries
489 : * that correspond to constrained dofs.
490 : */
491 : void zero_constrained_dofs_on_vector(NumericVector<Number> & vector) const;
492 :
493 : /**
494 : * @return true if the most recent truth solve gave a zero solution.
495 : */
496 : virtual bool check_if_zero_truth_solve() const;
497 :
498 : #ifdef LIBMESH_ENABLE_DIRICHLET
499 : /**
500 : * It's helpful to be able to generate a DirichletBoundary that stores a ZeroFunction in order
501 : * to impose Dirichlet boundary conditions.
502 : */
503 : static std::unique_ptr<DirichletBoundary> build_zero_dirichlet_boundary_object();
504 : #endif
505 :
506 : /**
507 : * Setter for the flag determining if convergence should be
508 : * checked after each solve.
509 : */
510 : void set_convergence_assertion_flag(bool flag);
511 :
512 : /**
513 : * Get/set flag to pre-evaluate the theta functions
514 : */
515 : bool get_preevaluate_thetas_flag() const;
516 : void set_preevaluate_thetas_flag(bool flag);
517 :
518 :
519 : //----------- PUBLIC DATA MEMBERS -----------//
520 :
521 : /**
522 : * Vector storing the values of the error bound
523 : * for each parameter in the training set --- the
524 : * parameter giving the largest error bound is
525 : * chosen for the next snapshot in the Greedy basis
526 : * training.
527 : */
528 : std::vector<Real> training_error_bounds;
529 :
530 : /**
531 : * We store an extra linear solver object which we can optionally
532 : * use for solving all systems in which the system matrix is set
533 : * to inner_product_matrix.
534 : */
535 : std::unique_ptr<LinearSolver<Number>> inner_product_solver;
536 :
537 : /**
538 : * Also, we store a pointer to an extra linear solver. This can be
539 : * useful if we want to pass in the linear solver from somewhere
540 : * else. For example, if a solver is already primed elsewhere
541 : * then it can be more efficient to use that solver.
542 : */
543 : LinearSolver<Number> * extra_linear_solver;
544 :
545 : /**
546 : * The inner product matrix.
547 : */
548 : std::unique_ptr<SparseMatrix<Number>> inner_product_matrix;
549 :
550 : /**
551 : * Vector storing the truth output values from the most
552 : * recent truth solve.
553 : */
554 : std::vector<Number > truth_outputs;
555 :
556 : /**
557 : * The vector storing the dual norm inner product terms
558 : * for each output.
559 : */
560 : std::vector<std::vector<Number >> output_dual_innerprods;
561 :
562 : /**
563 : * Vector storing the residual representors associated with the
564 : * right-hand side.
565 : * These are basis independent and hence stored here, whereas
566 : * the Aq_representors are stored in RBEvaluation
567 : */
568 : std::vector<std::unique_ptr<NumericVector<Number>>> Fq_representor;
569 :
570 : /**
571 : * Vectors storing the residual representor inner products
572 : * to be used in computing the residuals online. We store
573 : * the Fq representor norms here because they are independent
574 : * of a reduced basis space. The basis dependent representors
575 : * are stored in RBEvaluation.
576 : */
577 : std::vector<Number> Fq_representor_innerprods;
578 :
579 : /**
580 : * Boolean flag to indicate if we skip residual calculations
581 : * in train_reduced_basis. This should only be used in
582 : * special cases, e.g. when we know a priori that we want
583 : * exactly one basis function and hence we do not need the
584 : * residual based error indicator.
585 : */
586 : bool skip_residual_in_train_reduced_basis;
587 :
588 : /**
589 : * Boolean flag to indicate whether we exit the greedy if
590 : * we select the same parameters twice in a row. In some
591 : * problems this indicates that the greedy has "saturated"
592 : * typically due to numerical rounding effects.
593 : */
594 : bool exit_on_repeated_greedy_parameters;
595 :
596 : /**
597 : * Boolean flag to indicate whether we impose "fluxes"
598 : * (i.e. element boundary contributions to the weak form)
599 : * on internal element boundaries in the assembly routines.
600 : */
601 : bool impose_internal_fluxes;
602 :
603 : /**
604 : * In some cases meshes are intentionally created with degenerate sides
605 : * as a way to represent, say, triangles using a hex-only mesh. In this
606 : * situation we should detect and skip any degenerate sides in order to
607 : * prevent zero or negative element Jacobian errors.
608 : */
609 : bool skip_degenerate_sides;
610 :
611 : /**
612 : * Boolean flag to indicate whether we compute the RB_inner_product_matrix.
613 : * This is false by default in RBConstruction since (in the default implementation)
614 : * the RB inner-product matrix will just be the identity. But we may need the
615 : * inner-product matrix subclasses.
616 : */
617 : bool compute_RB_inner_product;
618 :
619 : /**
620 : * Boolean flag to indicate whether we store a second copy of each
621 : * affine operator and vector which does not have Dirichlet bcs
622 : * enforced.
623 : */
624 : bool store_non_dirichlet_operators;
625 :
626 : /**
627 : * Boolean flag to indicate whether we store a second copy of the
628 : * basis without constraints or dof transformations applied to it.
629 : * This is necessary when we have dof transformations and need
630 : * to calculate the residual R(U) = C^T F - C^T A C U, since we
631 : * need to evaluate R(U) using the untransformed basis U rather
632 : * than C U to avoid "double applying" dof transformations in C.
633 : */
634 : bool store_untransformed_basis;
635 :
636 : /**
637 : * A boolean flag to indicate whether or not we initialize the
638 : * Greedy algorithm by performing rb_solves on the training set
639 : * with an "empty" (i.e. N=0) reduced basis space.
640 : */
641 : bool use_empty_rb_solve_in_greedy;
642 :
643 : /**
644 : * A boolean flag to indicate whether or not the Fq representor norms
645 : * have already been computed --- used to make sure that we don't
646 : * recompute them unnecessarily.
647 : */
648 : bool Fq_representor_innerprods_computed;
649 :
650 : protected:
651 :
652 : /**
653 : * Helper function that actually allocates all the data
654 : * structures required by this class.
655 : */
656 : virtual void allocate_data_structures();
657 :
658 : /**
659 : * Assemble the truth matrix and right-hand side
660 : * for current_parameters.
661 : */
662 : virtual void truth_assembly();
663 :
664 : /**
665 : * Builds a DGFEMContext object with enough information to do
666 : * evaluations on each element. We use DGFEMContext since it
667 : * allows for both DG and continuous Galerkin formulations.
668 : */
669 : virtual std::unique_ptr<DGFEMContext> build_context();
670 :
671 : /**
672 : * Return the matrix for the output residual dual
673 : * norm solves. By default we use the inner product matrix
674 : * for steady state problems.
675 : */
676 : virtual SparseMatrix<Number> & get_matrix_for_output_dual_solves();
677 :
678 : /**
679 : * Function that indicates when to terminate the Greedy
680 : * basis training. Override in subclasses to specialize.
681 : */
682 : virtual bool greedy_termination_test(Real abs_greedy_error,
683 : Real initial_greedy_error,
684 : int count);
685 :
686 : /**
687 : * Update the list of Greedily chosen parameters with
688 : * current_parameters.
689 : */
690 : void update_greedy_param_list();
691 :
692 : /**
693 : * This function loops over the mesh and applies the specified
694 : * interior and/or boundary assembly routines, then adds the
695 : * scaled result to input_matrix and/or input_vector.
696 : * If symmetrize==true then we assemble the symmetric part
697 : * of the matrix, 0.5*(A + A^T)
698 : */
699 : void add_scaled_matrix_and_vector(Number scalar,
700 : ElemAssembly * elem_assembly,
701 : SparseMatrix<Number> * input_matrix,
702 : NumericVector<Number> * input_vector,
703 : bool symmetrize=false,
704 : bool apply_dof_constraints=true);
705 :
706 : /**
707 : * This function is called from add_scaled_matrix_and_vector()
708 : * before each element matrix and vector are assembled into their
709 : * global counterparts. By default it is a no-op, but it could be
710 : * used to apply any user-defined transformations immediately prior
711 : * to assembly. We use DGFEMContext since it allows for both DG and
712 : * continuous Galerkin formulations.
713 : */
714 416000 : virtual void post_process_elem_matrix_and_vector
715 416000 : (DGFEMContext & /*context*/) {}
716 :
717 : /**
718 : * Similarly, provide an opportunity to post-process the truth
719 : * solution after the solve is complete. By default this is a no-op,
720 : * but it could be used to apply any required user-defined post
721 : * processing to the solution vector. Note: the truth solution is
722 : * stored in the "solution" member of this class, which is inherited
723 : * from the parent System class several levels up.
724 : */
725 3456 : virtual void post_process_truth_solution() {}
726 :
727 : /**
728 : * Set current_local_solution = vec so that we can access vec
729 : * from FEMContext during assembly. Override in subclasses if
730 : * different behavior is required.
731 : */
732 : virtual void set_context_solution_vec(NumericVector<Number> & vec);
733 :
734 : /**
735 : * Assemble and store all the inner-product
736 : * matrix, the constraint matrix (for constrained
737 : * problems) and the mass matrix (for time-dependent
738 : * problems).
739 : */
740 : virtual void assemble_misc_matrices();
741 :
742 : /**
743 : * Assemble and store all Q_a affine operators
744 : * as well as the inner-product matrix.
745 : */
746 : virtual void assemble_all_affine_operators();
747 :
748 : /**
749 : * Assemble and store the affine RHS vectors.
750 : */
751 : virtual void assemble_all_affine_vectors();
752 :
753 : /**
754 : * Assemble and store the output vectors.
755 : */
756 : virtual void assemble_all_output_vectors();
757 :
758 : /**
759 : * Compute and store the dual norm of each output functional.
760 : */
761 : virtual void compute_output_dual_innerprods();
762 :
763 : /**
764 : * Compute the terms that are combined `online'
765 : * to determine the dual norm of the residual. Here we
766 : * compute the terms associated with the right-hand side.
767 : * These terms are basis independent, hence we separate
768 : * them from the rest of the calculations that are done in
769 : * update_residual_terms.
770 : * By default,
771 : * inner product terms are also computed, but you can turn
772 : * this feature off e.g. if you are already reading in that
773 : * data from files.
774 : */
775 : virtual void compute_Fq_representor_innerprods(bool compute_inner_products=true);
776 :
777 : /**
778 : * Add a new basis function to the RB space. This is called
779 : * during train_reduced_basis.
780 : */
781 : virtual void enrich_RB_space();
782 :
783 : /**
784 : * Update the system after enriching the RB space; this calls
785 : * a series of functions to update the system properly.
786 : */
787 : virtual void update_system();
788 :
789 : /**
790 : * \returns The RB error bound for the current parameters.
791 : *
792 : * Used in the Greedy algorithm to select the next parameter.
793 : */
794 : virtual Real get_RB_error_bound();
795 :
796 : /**
797 : * Compute the reduced basis matrices for the current basis.
798 : */
799 : virtual void update_RB_system_matrices();
800 :
801 : /**
802 : * Compute the terms that are combined `online'
803 : * to determine the dual norm of the residual. By default,
804 : * inner product terms are also computed, but you can turn
805 : * this feature off e.g. if you are already reading in that
806 : * data from files.
807 : */
808 : virtual void update_residual_terms(bool compute_inner_products=true);
809 :
810 : /**
811 : * Initialize the FEMContext prior to performing
812 : * an element loop.
813 : * Reimplement this in derived classes in order to
814 : * call FE::get_*() as the particular physics requires.
815 : */
816 0 : virtual void init_context(FEMContext &)
817 : {
818 : // Failing to rederive init_context() means your FE objects don't
819 : // know what to compute.
820 : libmesh_deprecated();
821 0 : }
822 :
823 : /**
824 : * Getter for the flag determining if convergence should be
825 : * checked after each solve.
826 : */
827 : bool get_convergence_assertion_flag() const;
828 :
829 : /**
830 : * Check if the linear solver reports convergence.
831 : * Throw an error when that is not the case.
832 : */
833 : void check_convergence(LinearSolver<Number> & input_solver);
834 :
835 : /**
836 : * Get/set the current training parameter index
837 : */
838 : unsigned int get_current_training_parameter_index() const;
839 : void set_current_training_parameter_index(unsigned int index);
840 :
841 : /**
842 : * Return the evaluated theta functions at the given training parameter index.
843 : */
844 : const std::vector<Number> & get_evaluated_thetas(unsigned int training_parameter_index) const;
845 :
846 : /*
847 : * Pre-evaluate the theta functions on the entire (local) training parameter set.
848 : */
849 : virtual void preevaluate_thetas();
850 :
851 : /**
852 : * Reset the _preevaluate_thetas_completed flag to false. We can use this to force
853 : * us to recalculate preevaluate thetas, in cases where that is necessary.
854 : */
855 : void reset_preevaluate_thetas_completed();
856 :
857 : //----------- PROTECTED DATA MEMBERS -----------//
858 :
859 : /**
860 : * Maximum number of reduced basis
861 : * functions we are willing to use.
862 : */
863 : unsigned int Nmax;
864 :
865 : /**
866 : * The number of basis functions that we add at each greedy step.
867 : * This defaults to 1 in the steady case.
868 : */
869 : unsigned int delta_N;
870 :
871 : /**
872 : * A boolean flag to indicate whether or not the output dual norms
873 : * have already been computed --- used to make sure that we don't
874 : * recompute them unnecessarily.
875 : */
876 : bool output_dual_innerprods_computed;
877 :
878 : /**
879 : * A boolean flag to indicate whether to check for proper convergence
880 : * after each solve.
881 : */
882 : bool assert_convergence;
883 :
884 : private:
885 :
886 : //----------- PRIVATE DATA MEMBERS -----------//
887 :
888 : /**
889 : * The current RBEvaluation object we are using to
890 : * perform the Evaluation stage of the reduced basis
891 : * method.
892 : */
893 : RBEvaluation * rb_eval;
894 :
895 : /**
896 : * This member holds the (parameter independent) assembly functors
897 : * that define the "affine expansion" of the PDE that we are solving.
898 : */
899 : RBAssemblyExpansion * rb_assembly_expansion;
900 :
901 : /**
902 : * Pointer to inner product assembly.
903 : */
904 : ElemAssembly * inner_product_assembly;
905 :
906 : /**
907 : * Boolean to indicate whether we're using the energy inner-product.
908 : * If this is false then we use inner_product_assembly instead.
909 : */
910 : bool use_energy_inner_product;
911 :
912 : /**
913 : * We may optionally want to use the "energy inner-product" rather
914 : * than the inner-product assembly specified in inner_product_assembly.
915 : * In this case the inner-product will be defined by sum_q^Q k_q * A_q.
916 : * Here we provide the k_q values that will be used.
917 : * (Note that a true "energy-inner product" would obtain the k_q from
918 : * the theta_q's, but this is different for each parameter choice so
919 : * we just provide a fixed set of k_q's here to ensure that the
920 : * inner-product is parameter independent)
921 : */
922 : std::vector<Number> energy_inner_product_coeffs;
923 :
924 : /**
925 : * Vector storing the Q_a matrices from the affine expansion
926 : */
927 : std::vector<std::unique_ptr<SparseMatrix<Number>>> Aq_vector;
928 :
929 : /**
930 : * Vector storing the Q_f vectors in the affine decomposition
931 : * of the right-hand side.
932 : */
933 : std::vector<std::unique_ptr<NumericVector<Number>>> Fq_vector;
934 :
935 : /**
936 : * The libMesh vectors that define the output functionals.
937 : * Each row corresponds to the affine expansion of an output.
938 : */
939 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>>> outputs_vector;
940 :
941 : /**
942 : * We may also need a second set of matrices/vectors
943 : * that do not have the Dirichlet boundary conditions
944 : * enforced.
945 : */
946 : std::vector<std::unique_ptr<SparseMatrix<Number>>> non_dirichlet_Aq_vector;
947 : std::vector<std::unique_ptr<NumericVector<Number>>> non_dirichlet_Fq_vector;
948 : std::vector<std::vector<std::unique_ptr<NumericVector<Number>>>> non_dirichlet_outputs_vector;
949 : std::unique_ptr<SparseMatrix<Number>> non_dirichlet_inner_product_matrix;
950 :
951 : /**
952 : * Relative and absolute tolerances for training reduced basis
953 : * using the Greedy scheme.
954 : */
955 : Real rel_training_tolerance;
956 : Real abs_training_tolerance;
957 :
958 : /**
959 : * This boolean indicates if we normalize the RB error in the greedy using
960 : * RBEvaluation::get_error_bound_normalization().
961 : */
962 : bool normalize_rb_bound_in_greedy;
963 :
964 : /**
965 : * This string indicates the type of training that we will use.
966 : * Options are:
967 : * - Greedy: Reduced basis greedy algorithm
968 : * - POD: Proper Orthogonal Decomposition
969 : */
970 : std::string RB_training_type;
971 :
972 : /**
973 : * In cases where we have dof transformations such as a change of
974 : * coordinates at some nodes we need to store an extra set of basis
975 : * functions which have not had dof transformations applied to them.
976 : * These vectors are required in order to compute the residual in
977 : * the error indicator.
978 : */
979 : std::vector<std::unique_ptr<NumericVector<Number>>> _untransformed_basis_functions;
980 :
981 : /**
982 : * We also store a copy of the untransformed solution in order to
983 : * create _untransformed_basis_functions.
984 : */
985 : std::unique_ptr<NumericVector<Number>> _untransformed_solution;
986 :
987 : /**
988 : * Flag to indicate if we preevaluate the theta functions
989 : */
990 : bool _preevaluate_thetas_flag;
991 :
992 : /**
993 : * Flag to indicate if the preevaluate_thetas function has been called, since
994 : * this allows us to avoid calling preevaluate_thetas more than once, which is
995 : * typically unnecessary.
996 : */
997 : bool _preevaluate_thetas_completed;
998 :
999 : /**
1000 : * The current training parameter index during reduced basis training.
1001 : */
1002 : unsigned int _current_training_parameter_index;
1003 :
1004 : /**
1005 : * Storage of evaluated theta functions at a set of parameters. This
1006 : * can be used to store all of our theta functions at training samples
1007 : * instead of re-evaluating the same values repeatedly during training.
1008 : */
1009 : std::vector<std::vector<Number>> _evaluated_thetas;
1010 : };
1011 :
1012 : } // namespace libMesh
1013 :
1014 : #endif // LIBMESH_RB_CONSTRUCTION_H
|
