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_EIM_CONSTRUCTION_H
21 : #define LIBMESH_RB_EIM_CONSTRUCTION_H
22 :
23 : // rbOOmit includes
24 : #include "libmesh/rb_construction.h"
25 : #include "libmesh/rb_assembly_expansion.h"
26 : #include "libmesh/rb_eim_assembly.h"
27 : #include "libmesh/rb_eim_evaluation.h"
28 :
29 : // libMesh includes
30 : #include "libmesh/mesh_function.h"
31 : #include "libmesh/coupling_matrix.h"
32 :
33 : // C++ includes
34 : #include <unordered_map>
35 : #include <map>
36 : #include <string>
37 : #include <memory>
38 : #include <vector>
39 :
40 : namespace libMesh
41 : {
42 :
43 : /**
44 : * This struct is used to encapsulate the arguments required
45 : * to specify an EIM point that we may add to our list of
46 : * interpolation points.
47 : */
48 : struct EimPointData
49 : {
50 : dof_id_type elem_id;
51 : dof_id_type node_id;
52 : unsigned int side_index;
53 : unsigned int comp_index;
54 : unsigned int qp_index;
55 : };
56 :
57 : /**
58 : * This class is part of the rbOOmit framework.
59 : *
60 : * RBEIMConstruction implements the Construction stage of the
61 : * Empirical Interpolation Method (EIM). This can be used to
62 : * create an approximation to parametrized functions. In the context
63 : * of the reduced basis (RB) method, the EIM approximation is typically
64 : * used to create an affine approximation to non-affine operators,
65 : * so that the standard RB method can be applied in that case.
66 : */
67 0 : class RBEIMConstruction : public RBConstructionBase<System>
68 : {
69 : public:
70 :
71 : enum BEST_FIT_TYPE { PROJECTION_BEST_FIT, EIM_BEST_FIT, POD_BEST_FIT };
72 :
73 : /**
74 : * Constructor. Optionally initializes required
75 : * data structures.
76 : */
77 : RBEIMConstruction (EquationSystems & es,
78 : const std::string & name,
79 : const unsigned int number);
80 :
81 : /**
82 : * Special functions.
83 : * - This class has the same restrictions/defaults as its base class.
84 : * - Destructor is defaulted out-of-line
85 : */
86 : RBEIMConstruction (RBEIMConstruction &&) = default;
87 : RBEIMConstruction (const RBEIMConstruction &) = delete;
88 : RBEIMConstruction & operator= (const RBEIMConstruction &) = delete;
89 : RBEIMConstruction & operator= (RBEIMConstruction &&) = delete;
90 : virtual ~RBEIMConstruction ();
91 :
92 : /**
93 : * Type of the data structure used to map from (elem id) -> [n_vars][n_qp] data.
94 : */
95 : typedef RBEIMEvaluation::QpDataMap QpDataMap;
96 :
97 : /**
98 : * Type of the data structure used to map from (elem id,side_index) -> [n_vars][n_qp] data.
99 : */
100 : typedef RBEIMEvaluation::SideQpDataMap SideQpDataMap;
101 :
102 : /**
103 : * Type of the data structure used to map from node id -> [n_vars] data.
104 : */
105 : typedef RBEIMEvaluation::NodeDataMap NodeDataMap;
106 :
107 : /**
108 : * Clear this object.
109 : */
110 : virtual void clear() override;
111 :
112 : /**
113 : * Set the RBEIMEvaluation object.
114 : */
115 : void set_rb_eim_evaluation(RBEIMEvaluation & rb_eim_eval_in);
116 :
117 : /**
118 : * Get a reference to the RBEvaluation object.
119 : */
120 : RBEIMEvaluation & get_rb_eim_evaluation();
121 :
122 : /**
123 : * Get a const reference to the RBEvaluation object.
124 : */
125 : const RBEIMEvaluation & get_rb_eim_evaluation() const;
126 :
127 : /**
128 : * Perform initialization of this object to prepare for running
129 : * train_eim_approximation().
130 : */
131 : void initialize_eim_construction();
132 :
133 : /**
134 : * Read parameters in from file and set up this system
135 : * accordingly.
136 : */
137 : virtual void process_parameters_file (const std::string & parameters_filename);
138 :
139 : /**
140 : * Set the state of this RBConstruction object based on the arguments
141 : * to this function.
142 : */
143 : void set_rb_construction_parameters(unsigned int n_training_samples_in,
144 : bool deterministic_training_in,
145 : int training_parameters_random_seed_in,
146 : bool quiet_mode_in,
147 : unsigned int Nmax_in,
148 : Real rel_training_tolerance_in,
149 : Real abs_training_tolerance_in,
150 : const RBParameters & mu_min_in,
151 : const RBParameters & mu_max_in,
152 : const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
153 : const std::map<std::string,bool> & log_scaling,
154 : std::map<std::string, std::vector<RBParameter>> * training_sample_list=nullptr);
155 :
156 : /**
157 : * Specify which type of "best fit" we use to guide the EIM
158 : * greedy algorithm.
159 : */
160 : virtual void set_best_fit_type_flag (const std::string & best_fit_type_string);
161 :
162 : /**
163 : * Print out info that describes the current setup of this RBConstruction.
164 : */
165 : virtual void print_info();
166 :
167 : /**
168 : * Rescale solution snapshots so that they all have unity norm. This is relevant
169 : * if training samples have differing magnitudes and we want to approximate them
170 : * all with equal accuracy.
171 : */
172 : void apply_normalization_to_solution_snapshots();
173 :
174 : /**
175 : * Generate the EIM approximation for the specified parametrized function
176 : * using either POD or the Greedy Algorithm. Return the final tolerance.
177 : */
178 : virtual Real train_eim_approximation();
179 :
180 : /**
181 : * Generate the EIM approximation for the specified parametrized function
182 : * using the Greedy Algorithm. Return the final tolerance.
183 : * Method is virtual so that behavior can be specialized further in
184 : * subclasses, if needed.
185 : */
186 : virtual Real train_eim_approximation_with_greedy();
187 :
188 : /**
189 : * Generate the EIM approximation for the specified parametrized function
190 : * using Proper Orthogonal Decomposition (POD). Return the final tolerance.
191 : * Method is virtual so that behavior can be specialized further in
192 : * subclasses, if needed.
193 : */
194 : virtual Real train_eim_approximation_with_POD();
195 :
196 : /**
197 : * Build a vector of ElemAssembly objects that accesses the basis
198 : * functions stored in this RBEIMConstruction object. This is useful
199 : * for performing the Offline stage of the Reduced Basis method where
200 : * we want to use assembly functions based on this EIM approximation.
201 : */
202 : virtual void initialize_eim_assembly_objects();
203 :
204 : /**
205 : * \returns The vector of assembly objects that point to this RBEIMConstruction.
206 : */
207 : std::vector<std::unique_ptr<ElemAssembly>> & get_eim_assembly_objects();
208 :
209 : /**
210 : * Build an element assembly object that will access basis function
211 : * \p bf_index.
212 : * This is pure virtual, override in subclasses to specify the appropriate
213 : * ElemAssembly object.
214 : */
215 : virtual std::unique_ptr<ElemAssembly> build_eim_assembly(unsigned int bf_index) = 0;
216 :
217 : /**
218 : * Pre-request FE data needed for calculations.
219 : */
220 : virtual void init_context(FEMContext &);
221 :
222 : /**
223 : * Get/set the relative tolerance for the basis training.
224 : */
225 : void set_rel_training_tolerance(Real new_training_tolerance);
226 : Real get_rel_training_tolerance();
227 :
228 : /**
229 : * Get/set the absolute tolerance for the basis training.
230 : */
231 : void set_abs_training_tolerance(Real new_training_tolerance);
232 : Real get_abs_training_tolerance();
233 :
234 : /**
235 : * Get/set Nmax, the maximum number of RB
236 : * functions we are willing to compute.
237 : */
238 : unsigned int get_Nmax() const;
239 : virtual void set_Nmax(unsigned int Nmax);
240 :
241 : /**
242 : * Call this method to set _set_Nmax_from_n_snapshots=true and
243 : * _Nmax_from_n_snapshots_increment=increment. This means that
244 : * we will overrule Nmax to be n_snapshots + increment, where
245 : * increment can be positive or negative (typically it should be
246 : * negative to limit Nmax to be less that the number of snapshots).
247 : */
248 : void enable_set_Nmax_from_n_snapshots(int increment);
249 :
250 : /**
251 : * Call this method to set _set_Nmax_from_n_snapshots=false and
252 : * reset _Nmax_from_n_snapshots_increment to 0.
253 : */
254 : void disable_set_Nmax_from_n_snapshots();
255 :
256 : /**
257 : * Get the maximum value (across all processors) from
258 : * the parametrized functions in the training set.
259 : */
260 : Real get_max_abs_value_in_training_set() const;
261 :
262 : /**
263 : * Get the EIM solution vector at all parametrized functions in the training
264 : * set. In some cases we want to store this data for future use. For example
265 : * this is useful in the case that the parametrized function is defined
266 : * based on a look-up table rather than an analytical function, since
267 : * if we store the EIM solution data, we can do Online solves without
268 : * initializing the look-up table data.
269 : */
270 : void store_eim_solutions_for_training_set();
271 :
272 : /**
273 : * Get a const reference to the specified parametrized function from
274 : * the training set.
275 : */
276 : const QpDataMap & get_parametrized_function_from_training_set(unsigned int training_index) const;
277 : const SideQpDataMap & get_side_parametrized_function_from_training_set(unsigned int training_index) const;
278 : const NodeDataMap & get_node_parametrized_function_from_training_set(unsigned int training_index) const;
279 :
280 : /**
281 : * Get the interior and side quadrature weights.
282 : */
283 : const std::unordered_map<dof_id_type, std::vector<Real> > & get_local_quad_point_JxW();
284 : const std::map<std::pair<dof_id_type,unsigned int>, std::vector<Real> > & get_local_side_quad_point_JxW();
285 :
286 : /**
287 : * Get the number of parametrized functions used for training.
288 : */
289 : unsigned int get_n_parametrized_functions_for_training() const;
290 :
291 : /**
292 : * Zero the _eim_projection_matrix and resize it to be get_Nmax() x get_Nmax().
293 : */
294 : void reinit_eim_projection_matrix();
295 :
296 : /**
297 : * Enum that indicates which type of "best fit" algorithm
298 : * we should use.
299 : * a) projection: Find the best fit in the inner product
300 : * b) eim: Use empirical interpolation to find a "best fit"
301 : */
302 : BEST_FIT_TYPE best_fit_type_flag;
303 :
304 : protected:
305 :
306 : /**
307 : * Implementation of enrich_eim_approximation() for the case of element sides.
308 : *
309 : * If \p add_basis_function is true, then we add an extra basis function to the
310 : * EIM basis. If it is false, then we only store the data associated with the
311 : * interpolation point that we identify, which can be relevant when setting up
312 : * data for the error indicator, for example.
313 : *
314 : * If \p eim_point_data is not nullptr, then we add the extra point that is
315 : * specified rather than looking for the "optimal point" in \p interior_pf.
316 : *
317 : * @return true if \p side_pf is linearly dependent to the existing basis,
318 : * and in this case we skip adding the basis function since we do not want
319 : * to add linearly dependent data to the basis.
320 : */
321 : bool enrich_eim_approximation_on_sides(const SideQpDataMap & side_pf,
322 : bool add_basis_function,
323 : EimPointData * eim_point_data);
324 :
325 : /**
326 : * Implementation of enrich_eim_approximation() for the case of element nodes.
327 : */
328 : bool enrich_eim_approximation_on_nodes(const NodeDataMap & node_pf,
329 : bool add_basis_function,
330 : EimPointData * eim_point_data);
331 :
332 : /**
333 : * Implementation of enrich_eim_approximation() for the case of element interiors.
334 : */
335 : bool enrich_eim_approximation_on_interiors(const QpDataMap & interior_pf,
336 : bool add_basis_function,
337 : EimPointData * eim_point_data);
338 :
339 : /**
340 : * Update the matrices used in training the EIM approximation.
341 : *
342 : * If \p set_eim_error_indicator is true then we add data corresponding
343 : * to the EIM error indicator.
344 : */
345 : void update_eim_matrices(bool set_eim_error_indicator);
346 :
347 : private:
348 :
349 : /**
350 : * Find the training sample that has the largest EIM approximation error
351 : * based on the current EIM approximation. Return the maximum error, and
352 : * the training sample index at which it occurred.
353 : */
354 : std::pair<Real, unsigned int> compute_max_eim_error();
355 :
356 : /**
357 : * Compute and store the parametrized function for each
358 : * parameter in the training set at all the stored qp locations.
359 : */
360 : void initialize_parametrized_functions_in_training_set();
361 :
362 : /**
363 : * Initialize the data associated with each quad point (location, JxW, etc.)
364 : * so that we can use this in evaluation of the parametrized functions.
365 : */
366 : void initialize_qp_data();
367 :
368 : /**
369 : * Evaluate the inner product of vec1 and vec2 which specify values at
370 : * quadrature points. The inner product includes the JxW contributions
371 : * stored in _local_quad_point_JxW, so that this is equivalent to
372 : * computing w^t M v, where M is the mass matrix.
373 : *
374 : * If \p apply_comp_scaling then we will incorporate the scaling from
375 : * _component_scaling_in_training_set in the inner product.
376 : */
377 : Number inner_product(const QpDataMap & v, const QpDataMap & w, bool apply_comp_scaling);
378 :
379 : /**
380 : * Same as inner_product() except for side data.
381 : */
382 : Number side_inner_product(const SideQpDataMap & v, const SideQpDataMap & w, bool apply_comp_scaling);
383 :
384 : /**
385 : * Same as inner_product() except for node data.
386 : */
387 : Number node_inner_product(const NodeDataMap & v, const NodeDataMap & w, bool apply_comp_scaling);
388 :
389 : /**
390 : * Get the maximum absolute value from a vector stored in the format that we use
391 : * for basis functions.
392 : */
393 : template <class DataMap>
394 0 : Real get_max_abs_value(const DataMap & v) const
395 : {
396 0 : Real max_value = 0.;
397 :
398 0 : for (const auto & pr : v)
399 : {
400 0 : const auto & v_comp_and_qp = pr.second;
401 :
402 0 : for (const auto & comp : index_range(v_comp_and_qp))
403 : {
404 0 : Real comp_scaling = 1.;
405 0 : if (get_rb_eim_evaluation().scale_components_in_enrichment().count(comp))
406 : {
407 : // Make sure that _component_scaling_in_training_set is initialized
408 0 : libmesh_error_msg_if(comp >= _component_scaling_in_training_set.size(),
409 : "Invalid vector index");
410 0 : comp_scaling = _component_scaling_in_training_set[comp];
411 : }
412 :
413 0 : const std::vector<Number> & v_qp = v_comp_and_qp[comp];
414 0 : for (Number value : v_qp)
415 0 : max_value = std::max(max_value, std::abs(value * comp_scaling));
416 : }
417 : }
418 :
419 0 : comm().max(max_value);
420 0 : return max_value;
421 : }
422 :
423 : /**
424 : * Get the maximum absolute value from a vector stored in the format that we use
425 : * for basis functions. This case handles NodeDataMap.
426 : */
427 : Real get_node_max_abs_value(const NodeDataMap & v) const;
428 :
429 : /**
430 : * Add a new basis function to the EIM approximation.
431 : */
432 : void enrich_eim_approximation(unsigned int training_index,
433 : bool add_basis_function,
434 : EimPointData * eim_point_data);
435 :
436 : /**
437 : * Scale all values in \p pf by \p scaling_factor
438 : */
439 : template<class DataMap>
440 0 : static void scale_parametrized_function(DataMap & local_pf,
441 : Number scaling_factor)
442 : {
443 0 : for (auto & pr : local_pf)
444 : {
445 0 : auto & comp_and_qp = pr.second;
446 :
447 0 : for (unsigned int comp : index_range(comp_and_qp))
448 : {
449 0 : std::vector<Number> & qp_values = comp_and_qp[comp];
450 :
451 0 : for (unsigned int qp : index_range(qp_values))
452 : {
453 0 : qp_values[qp] *= scaling_factor;
454 : }
455 : }
456 : }
457 0 : }
458 :
459 : /**
460 : * Scale all values in \p pf by \p scaling_factor
461 : * The templated function above handles the elem and side cases,
462 : * and this separate case handles the node case.
463 : */
464 : static void scale_node_parametrized_function(NodeDataMap & local_pf,
465 : Number scaling_factor);
466 :
467 : /**
468 : * Static helper function that is used by get_random_point().
469 : */
470 : static unsigned int get_random_int_0_to_n(unsigned int n);
471 :
472 : /**
473 : * Helper function that identifies a random EIM point from \p v.
474 : */
475 : EimPointData get_random_point(const QpDataMap & v);
476 : EimPointData get_random_point(const SideQpDataMap & v);
477 : EimPointData get_random_point(const NodeDataMap & v);
478 :
479 : /**
480 : * Get a random point using the 0^th training sample as input to
481 : * get_random_point().
482 : */
483 : EimPointData get_random_point_from_training_sample();
484 :
485 : /**
486 : * Maximum number of EIM basis functions we are willing to use.
487 : */
488 : unsigned int _Nmax;
489 :
490 : /**
491 : * If _set_Nmax_from_n_snapshots=true, then we overrule Nmax to be
492 : * Nmax += _Nmax_from_n_snapshots_increment. Note that the "increment
493 : * can be positive or negative. Typically we would want to set the
494 : * increment to be negative or 0 to limit Nmax based on the number
495 : * of available snapshots, but in some rare cases it could make sense
496 : * to set it to a positive value, e.g. if we are appending to a basis
497 : * that has already been generated via a previous training.
498 : */
499 : bool _set_Nmax_from_n_snapshots;
500 : int _Nmax_from_n_snapshots_increment;
501 :
502 : /**
503 : * Relative and absolute tolerances for training the EIM approximation.
504 : */
505 : Real _rel_training_tolerance;
506 : Real _abs_training_tolerance;
507 :
508 : /**
509 : * The matrix we use in order to perform L2 projections of
510 : * parametrized functions as part of EIM training.
511 : */
512 : DenseMatrix<Number> _eim_projection_matrix;
513 :
514 : /**
515 : * The RBEIMEvaluation object that we use to perform the EIM training.
516 : */
517 : RBEIMEvaluation * _rb_eim_eval;
518 :
519 : /**
520 : * The vector of assembly objects that are created to point to
521 : * this RBEIMConstruction.
522 : */
523 : std::vector<std::unique_ptr<ElemAssembly>> _rb_eim_assembly_objects;
524 :
525 : /**
526 : * The parametrized functions that are used for training. We pre-compute and
527 : * store all of these functions, rather than recompute them at each iteration
528 : * of the training.
529 : *
530 : * We store values at quadrature points on elements that are local to this processor.
531 : * The indexing is as follows:
532 : * basis function index --> element ID --> variable --> quadrature point --> value
533 : * We use a map to index the element ID, since the IDs on this processor in
534 : * generally will not start at zero.
535 : */
536 : std::vector<QpDataMap> _local_parametrized_functions_for_training;
537 :
538 : /**
539 : * Same as _local_parametrized_functions_for_training except for side data.
540 : * The indexing is as follows:
541 : * basis function index --> (element ID,side index) --> variable --> quadrature point --> value
542 : */
543 : std::vector<SideQpDataMap> _local_side_parametrized_functions_for_training;
544 :
545 : /**
546 : * Same as _local_parametrized_functions_for_training except for node data.
547 : * The indexing is as follows:
548 : * basis function index --> node ID --> variable --> value
549 : */
550 : std::vector<NodeDataMap> _local_node_parametrized_functions_for_training;
551 :
552 : /**
553 : * Maximum value in _local_parametrized_functions_for_training across all processors.
554 : * This can be used for normalization purposes, for example.
555 : */
556 : Real _max_abs_value_in_training_set;
557 :
558 : /**
559 : * The training sample index at which we found _max_abs_value_in_training_set.
560 : */
561 : unsigned int _max_abs_value_in_training_set_index;
562 :
563 : /**
564 : * Keep track of a scaling factor for each component of the parametrized functions in
565 : * the training set which "scales up" each component to have a similar magnitude as
566 : * the largest component encountered in the training set. This can give more uniform
567 : * scaling across all components and is helpful in cases where components have widely
568 : * varying magnitudes.
569 : */
570 : std::vector<Real> _component_scaling_in_training_set;
571 :
572 : /**
573 : * The quadrature point locations, quadrature point weights (JxW), and subdomain IDs
574 : * on every element local to this processor.
575 : *
576 : * The indexing is as follows:
577 : * element ID --> quadrature point --> xyz
578 : * element ID --> quadrature point --> JxW
579 : * element ID --> subdomain_id
580 : * We use a map to index the element ID, since the IDs on this processor in
581 : * generally will not start at zero.
582 : */
583 : std::unordered_map<dof_id_type, std::vector<Point> > _local_quad_point_locations;
584 : std::unordered_map<dof_id_type, std::vector<Real> > _local_quad_point_JxW;
585 : std::unordered_map<dof_id_type, subdomain_id_type > _local_quad_point_subdomain_ids;
586 :
587 : /**
588 : * EIM approximations often arise when applying a geometric mapping to a Reduced Basis
589 : * formulation. In this context, we often need to approximate derivates of the mapping
590 : * function via EIM. In order to enable this, we also optionally store perturbations
591 : * about each point in _local_quad_point_locations to enable finite difference approximation
592 : * to the mapping function derivatives.
593 : */
594 : std::unordered_map<dof_id_type, std::vector<std::vector<Point>> > _local_quad_point_locations_perturbations;
595 :
596 : /**
597 : * Same as above except for side data.
598 : */
599 : std::map<std::pair<dof_id_type,unsigned int>, std::vector<Point> > _local_side_quad_point_locations;
600 : std::map<std::pair<dof_id_type,unsigned int>, std::vector<Real> > _local_side_quad_point_JxW;
601 : std::map<std::pair<dof_id_type,unsigned int>, subdomain_id_type > _local_side_quad_point_subdomain_ids;
602 : std::map<std::pair<dof_id_type,unsigned int>, boundary_id_type > _local_side_quad_point_boundary_ids;
603 : std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Point>> > _local_side_quad_point_locations_perturbations;
604 :
605 : /**
606 : * Same as above except for node data.
607 : */
608 : std::unordered_map<dof_id_type, Point > _local_node_locations;
609 : std::unordered_map<dof_id_type, boundary_id_type > _local_node_boundary_ids;
610 :
611 : /**
612 : * For side data, we also store "side type" info. This is used to distinguish between
613 : * data that is stored on a "shellface" vs. a "standard side". The convention we use
614 : * here is:
615 : * 0 --> standard side
616 : * 1 --> shellface
617 : */
618 : std::map<std::pair<dof_id_type,unsigned int>, unsigned int > _local_side_quad_point_side_types;
619 :
620 : };
621 :
622 : } // namespace libMesh
623 :
624 : #endif // LIBMESH_RB_EIM_CONSTRUCTION_H
|
