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_EVALUATION_H
21 : #define LIBMESH_RB_EIM_EVALUATION_H
22 :
23 : // libMesh includes
24 : #include "libmesh/point.h"
25 : #include "libmesh/rb_theta_expansion.h"
26 : #include "libmesh/rb_parametrized.h"
27 : #include "libmesh/parallel_object.h"
28 : #include "libmesh/dense_matrix.h"
29 : #include "libmesh/dense_vector.h"
30 : #include "libmesh/rb_parametrized_function.h"
31 : #include "libmesh/fe_type.h"
32 :
33 : // C++ includes
34 : #include <memory>
35 : #include <map>
36 : #include <vector>
37 : #include <string>
38 :
39 : namespace libMesh
40 : {
41 :
42 : class RBParameters;
43 : class RBParametrizedFunction;
44 : class RBTheta;
45 : class System;
46 : class EquationSystems;
47 : class Elem;
48 :
49 : /**
50 : * This struct encapsulates data that specifies how we will
51 : * perform plotting for EIM variable groups.
52 : */
53 0 : struct EIMVarGroupPlottingInfo
54 : {
55 :
56 : /**
57 : * Default constructor. Initialize values that would otherwise
58 : * be uninitialized. We do not initialize the FEType or string
59 : * members since they already have default constructors.
60 : */
61 : EIMVarGroupPlottingInfo()
62 : :
63 : first_eim_var_index(0),
64 : n_eim_vars(0),
65 : enforce_min_value(false),
66 : enforce_max_value(false),
67 : min_value(0.),
68 : max_value(0.)
69 : {
70 : }
71 :
72 : /**
73 : * The index for the first EIM variable in this variable group.
74 : */
75 : unsigned int first_eim_var_index;
76 :
77 : /**
78 : * The number of EIM variables in the group. The variables
79 : * are assumed to be numbered contiguously.
80 : */
81 : unsigned int n_eim_vars;
82 :
83 : /**
84 : * The name of the System we use for plotting this EIM variable group.
85 : */
86 : std::string eim_sys_name;
87 :
88 : /**
89 : * These booleans indicate if we should clamp the resulting output
90 : * to be above a min value or below a max value. This can be relevant
91 : * if we want to satisfy some physical constraints on the outputs, for
92 : * example, since these constraints may not be exactly satisfied by
93 : * the EIM output.
94 : */
95 : bool enforce_min_value;
96 : bool enforce_max_value;
97 :
98 : /**
99 : * The min (resp. max) value that we enforce if enforce_min_value
100 : * (resp. enforce_max_value) is true.
101 : */
102 : Real min_value;
103 : Real max_value;
104 : };
105 :
106 : /**
107 : * This class enables evaluation of an Empirical Interpolation Method (EIM)
108 : * approximation. RBEvaluation plays an analogous role in the context of
109 : * the regular reduced basis method.
110 : */
111 0 : class RBEIMEvaluation : public RBParametrized,
112 : public ParallelObject
113 : {
114 : public:
115 :
116 : /**
117 : * Constructor.
118 : */
119 : RBEIMEvaluation(const Parallel::Communicator & comm);
120 :
121 : /**
122 : * Special functions.
123 : * - This class contains unique_ptrs, so it can't be default copy
124 : constructed/assigned.
125 : * - The destructor is defaulted out of line.
126 : */
127 : RBEIMEvaluation (RBEIMEvaluation &&) = default;
128 : RBEIMEvaluation (const RBEIMEvaluation &) = delete;
129 : RBEIMEvaluation & operator= (const RBEIMEvaluation &) = delete;
130 : RBEIMEvaluation & operator= (RBEIMEvaluation &&) = default;
131 : virtual ~RBEIMEvaluation ();
132 :
133 : /**
134 : * Type of the data structure used to map from (elem id) -> [n_vars][n_qp] data.
135 : */
136 : typedef std::map<dof_id_type, std::vector<std::vector<Number>>> QpDataMap;
137 :
138 : /**
139 : * Type of the data structure used to map from (elem id, side index) -> [n_vars][n_qp] data.
140 : */
141 : typedef std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Number>>> SideQpDataMap;
142 :
143 : /**
144 : * Type of the data structure used to map from (node id) -> [n_vars] data.
145 : */
146 : typedef std::map<dof_id_type, std::vector<Number>> NodeDataMap;
147 :
148 : /**
149 : * Clear this object.
150 : */
151 : virtual void clear() override;
152 :
153 : /**
154 : * Resize the data structures for storing data associated
155 : * with this object.
156 : */
157 : void resize_data_structures(const unsigned int Nmax);
158 :
159 : /**
160 : * Set the parametrized function that we will approximate
161 : * using the Empirical Interpolation Method. This object
162 : * will take ownership of the unique pointer.
163 : */
164 : void set_parametrized_function(std::unique_ptr<RBParametrizedFunction> pf);
165 :
166 : /**
167 : * Get a reference to the parametrized function.
168 : */
169 : RBParametrizedFunction & get_parametrized_function();
170 :
171 : /**
172 : * Get a const reference to the parametrized function.
173 : */
174 : const RBParametrizedFunction & get_parametrized_function() const;
175 :
176 : /**
177 : * Calculate the EIM approximation for the given
178 : * right-hand side vector \p EIM_rhs. Store the
179 : * solution coefficients in the member _eim_solution.
180 : */
181 : DenseVector<Number> rb_eim_solve(DenseVector<Number> & EIM_rhs);
182 :
183 : /**
184 : * Perform rb_eim_solves at each mu in \p mus and store the results
185 : * in _rb_eim_solutions.
186 : */
187 : void rb_eim_solves(const std::vector<RBParameters> & mus, unsigned int N);
188 :
189 : /**
190 : * Initialize _interpolation_points_spatial_indices. Once this data is
191 : * initialized, we can store it in the training data, and read it back
192 : * in during the Online stage to be used in solves.
193 : */
194 : void initialize_interpolation_points_spatial_indices();
195 :
196 : /**
197 : * The Online counterpart of initialize_interpolation_points_spatial_indices().
198 : * This is used to initialize the spatial indices data in _parametrized_function
199 : * so that the _parametrized_function can be used in the Online stage without
200 : * reconstructing the spatial indices on every element or node in the mesh.
201 : */
202 : void initialize_param_fn_spatial_indices();
203 :
204 : /**
205 : * Return the current number of EIM basis functions.
206 : */
207 : unsigned int get_n_basis_functions() const;
208 :
209 : /**
210 : * Return the number of interpolation points. If we're not using the EIM error
211 : * indicator, then this matches get_n_basis_functions(), but if we are using
212 : * the EIM error indicator then we should have one extra interpolation point.
213 : */
214 : unsigned int get_n_interpolation_points() const;
215 :
216 : /**
217 : * Return the number of unique elements containing interpolation points
218 : */
219 : unsigned int get_n_elems() const;
220 :
221 : /**
222 : * Return the number of properties stored in the rb_property_map
223 : */
224 : unsigned int get_n_properties() const;
225 :
226 : /**
227 : * Set the number of basis functions. Useful when reading in
228 : * stored data.
229 : */
230 : void set_n_basis_functions(unsigned int n_bfs);
231 :
232 : /**
233 : * Subtract coeffs[i]*basis_function[i] from \p v.
234 : */
235 : void decrement_vector(QpDataMap & v,
236 : const DenseVector<Number> & coeffs);
237 :
238 : /**
239 : * Same as decrement_vector() except for Side data.
240 : */
241 : void side_decrement_vector(SideQpDataMap & v,
242 : const DenseVector<Number> & coeffs);
243 :
244 : /**
245 : * Same as decrement_vector() except for node data.
246 : */
247 : void node_decrement_vector(NodeDataMap & v,
248 : const DenseVector<Number> & coeffs);
249 :
250 : /**
251 : * Build a vector of RBTheta objects that accesses the components
252 : * of the RB_solution member variable of this RBEvaluation.
253 : * Store these objects in the member vector rb_theta_objects.
254 : */
255 : void initialize_eim_theta_objects();
256 :
257 : /**
258 : * \returns The vector of theta objects that point to this RBEIMEvaluation.
259 : */
260 : std::vector<std::unique_ptr<RBTheta>> & get_eim_theta_objects();
261 :
262 : /**
263 : * Build a theta object corresponding to EIM index \p index.
264 : * The default implementation builds an RBEIMTheta object, possibly
265 : * override in subclasses if we need more specialized behavior.
266 : */
267 : virtual std::unique_ptr<RBTheta> build_eim_theta(unsigned int index);
268 :
269 : /**
270 : * Fill up values by evaluating the parametrized function \p pf for all quadrature
271 : * points on element \p elem_id and component \p comp.
272 : */
273 : static void get_parametrized_function_values_at_qps(
274 : const QpDataMap & pf,
275 : dof_id_type elem_id,
276 : unsigned int comp,
277 : std::vector<Number> & values);
278 :
279 : /**
280 : * Same as get_parametrized_function_values_at_qps() except for side data.
281 : */
282 : static void get_parametrized_function_side_values_at_qps(
283 : const SideQpDataMap & pf,
284 : dof_id_type elem_id,
285 : unsigned int side_index,
286 : unsigned int comp,
287 : std::vector<Number> & values);
288 :
289 : /**
290 : * Same as get_parametrized_function_values_at_qps() except for node data.
291 : * Note that this does not do any parallel communication, so it is only
292 : * applicable to looking up local values.
293 : */
294 : static Number get_parametrized_function_node_local_value(
295 : const NodeDataMap & pf,
296 : dof_id_type node_id,
297 : unsigned int comp);
298 :
299 : /**
300 : * Same as above, except that we just return the value at the qp^th
301 : * quadrature point.
302 : */
303 : static Number get_parametrized_function_value(
304 : const Parallel::Communicator & comm,
305 : const QpDataMap & pf,
306 : dof_id_type elem_id,
307 : unsigned int comp,
308 : unsigned int qp);
309 :
310 : /**
311 : * Same as get_parametrized_function_value() except for side data.
312 : */
313 : static Number get_parametrized_function_side_value(
314 : const Parallel::Communicator & comm,
315 : const SideQpDataMap & pf,
316 : dof_id_type elem_id,
317 : unsigned int side_index,
318 : unsigned int comp,
319 : unsigned int qp);
320 :
321 : /**
322 : * Same as get_parametrized_function_value() except for node data.
323 : * Unlike get_parametrized_function_node_local_value(), this does
324 : * parallel communication, and therefore if can be used to look up
325 : * values regardless of whether or not \p node_id is local.
326 : */
327 : static Number get_parametrized_function_node_value(
328 : const Parallel::Communicator & comm,
329 : const NodeDataMap & pf,
330 : dof_id_type node_id,
331 : unsigned int comp);
332 :
333 : /**
334 : * Fill up \p values with the basis function values for basis function
335 : * \p basis_function_index and variable \p var, at all quadrature points
336 : * on element \p elem_id. Each processor stores data for only the
337 : * elements local to that processor, so if elem_id is not on this processor
338 : * then \p values will be empty.
339 : */
340 : void get_eim_basis_function_values_at_qps(unsigned int basis_function_index,
341 : dof_id_type elem_id,
342 : unsigned int var,
343 : std::vector<Number> & values) const;
344 :
345 : /**
346 : * Same as get_eim_basis_function_values_at_qps() except for side data.
347 : */
348 : void get_eim_basis_function_side_values_at_qps(unsigned int basis_function_index,
349 : dof_id_type elem_id,
350 : unsigned int side_index,
351 : unsigned int var,
352 : std::vector<Number> & values) const;
353 :
354 : /**
355 : * Same as get_eim_basis_function_values_at_qps() except for node data.
356 : * Note that this does not do any parallel communication, it just looks
357 : * up the value from _local_node_eim_basis_functions.
358 : */
359 : Number get_eim_basis_function_node_local_value(unsigned int basis_function_index,
360 : dof_id_type node_id,
361 : unsigned int var) const;
362 :
363 : /**
364 : * Same as above, except that we just return the value at the qp^th
365 : * quadrature point.
366 : */
367 : Number get_eim_basis_function_value(unsigned int basis_function_index,
368 : dof_id_type elem_id,
369 : unsigned int comp,
370 : unsigned int qp) const;
371 :
372 : /**
373 : * Same as get_eim_basis_function_value() except for side data.
374 : */
375 : Number get_eim_basis_function_side_value(unsigned int basis_function_index,
376 : dof_id_type elem_id,
377 : unsigned int side_index,
378 : unsigned int comp,
379 : unsigned int qp) const;
380 :
381 : /**
382 : * Same as get_eim_basis_function_value() except for node data.
383 : * Note that unlike get_eim_basis_function_node_local_value(),
384 : * this does do parallel communication so that it can be called
385 : * on any processor regardless of whether \p node_id is local or not.
386 : */
387 : Number get_eim_basis_function_node_value(unsigned int basis_function_index,
388 : dof_id_type node_id,
389 : unsigned int var) const;
390 :
391 : /**
392 : * Get a reference to the i^th basis function.
393 : */
394 : const QpDataMap & get_basis_function(unsigned int i) const;
395 :
396 : /**
397 : * Get a reference to the i^th side basis function.
398 : */
399 : const SideQpDataMap & get_side_basis_function(unsigned int i) const;
400 :
401 : /**
402 : * Get a reference to the i^th node basis function.
403 : */
404 : const NodeDataMap & get_node_basis_function(unsigned int i) const;
405 :
406 : /**
407 : * Set _rb_eim_solutions. Normally we update _rb_eim_solutions by performing
408 : * and EIM solve, but in some cases we want to set the EIM solution coefficients
409 : * elsewhere, so this setter enables us to do that.
410 : */
411 : void set_rb_eim_solutions(const std::vector<DenseVector<Number>> & rb_eim_solutions);
412 :
413 : /**
414 : * Return the EIM solution coefficients from the most recent call to rb_eim_solves().
415 : */
416 : const std::vector<DenseVector<Number>> & get_rb_eim_solutions() const;
417 :
418 : /**
419 : * Return entry \p index for each solution in _rb_eim_solutions.
420 : */
421 : std::vector<Number> get_rb_eim_solutions_entries(unsigned int index) const;
422 :
423 : /**
424 : * Return a const reference to the EIM solutions for the parameters in the training set.
425 : */
426 : const std::vector<DenseVector<Number>> & get_eim_solutions_for_training_set() const;
427 :
428 : /**
429 : * Return a writeable reference to the EIM solutions for the parameters in the training set.
430 : */
431 : std::vector<DenseVector<Number>> & get_eim_solutions_for_training_set();
432 :
433 : /**
434 : * Return the EIM error indicator values from the most recent call to rb_eim_solves().
435 : * The first entry of each pair is the normalized error indicator, and the second
436 : * entry is the normalization factor.
437 : */
438 : const std::vector<std::pair<Real,Real>> & get_rb_eim_error_indicators() const;
439 :
440 : /**
441 : * Set the data associated with EIM interpolation points.
442 : */
443 : void add_interpolation_points_xyz(Point p);
444 : void add_interpolation_points_comp(unsigned int comp);
445 : void add_interpolation_points_subdomain_id(subdomain_id_type sbd_id);
446 : void add_interpolation_points_boundary_id(boundary_id_type b_id);
447 : void add_interpolation_points_xyz_perturbations(const std::vector<Point> & perturbs);
448 : void add_interpolation_points_elem_id(dof_id_type elem_id);
449 : void add_interpolation_points_side_index(unsigned int side_index);
450 : void add_interpolation_points_node_id(dof_id_type node_id);
451 : void add_interpolation_points_qp(unsigned int qp);
452 : void add_interpolation_points_elem_type(ElemType elem_type);
453 : void add_interpolation_points_phi_i_qp(const std::vector<Real> & phi_i_qp);
454 : void add_interpolation_points_JxW_all_qp(const std::vector<Real> & JxW_all_qp);
455 : void add_interpolation_points_phi_i_all_qp(const std::vector<std::vector<Real>> & phi_i_all_qp);
456 : void add_interpolation_points_qrule_order(Order qrule_order);
457 : void add_elem_center_dxyzdxi(const Point & dxyzdxi);
458 : void add_elem_center_dxyzdeta(const Point & dxyzdxi);
459 : void add_interpolation_points_spatial_indices(const std::vector<unsigned int> & spatial_indices);
460 : void add_elem_id_local_index_map_entry(dof_id_type elem_id, unsigned int local_index);
461 : void add_rb_property_map_entry(std::string & property_name, std::set<dof_id_type> & entity_ids);
462 :
463 : /**
464 : * Get the data associated with EIM interpolation points.
465 : */
466 : Point get_interpolation_points_xyz(unsigned int index) const;
467 : unsigned int get_interpolation_points_comp(unsigned int index) const;
468 : subdomain_id_type get_interpolation_points_subdomain_id(unsigned int index) const;
469 : boundary_id_type get_interpolation_points_boundary_id(unsigned int index) const;
470 : const std::vector<Point> & get_interpolation_points_xyz_perturbations(unsigned int index) const;
471 : dof_id_type get_interpolation_points_elem_id(unsigned int index) const;
472 : unsigned int get_interpolation_points_side_index(unsigned int index) const;
473 : dof_id_type get_interpolation_points_node_id(unsigned int index) const;
474 : unsigned int get_interpolation_points_qp(unsigned int index) const;
475 : ElemType get_interpolation_points_elem_type(unsigned int index) const;
476 : const std::vector<Real> & get_interpolation_points_phi_i_qp(unsigned int index) const;
477 : const std::vector<Real> & get_interpolation_points_JxW_all_qp(unsigned int index) const;
478 : const std::vector<std::vector<Real>> & get_interpolation_points_phi_i_all_qp(unsigned int index) const;
479 : Order get_interpolation_points_qrule_order(unsigned int index) const;
480 : const Point & get_elem_center_dxyzdxi(unsigned int index) const;
481 : const Point & get_elem_center_dxyzdeta(unsigned int index) const;
482 : const std::vector<unsigned int> & get_interpolation_points_spatial_indices(unsigned int index) const;
483 : const std::map<dof_id_type, unsigned int> & get_elem_id_to_local_index_map() const;
484 : const std::unordered_map<std::string, std::set<dof_id_type>> & get_rb_property_map() const;
485 :
486 : /**
487 : * _interpolation_points_spatial_indices is optional data, so we need to be able to
488 : * check how many _interpolation_points_spatial_indices values have actually been set
489 : * since it may not match the number of interpolation points.
490 : */
491 : unsigned int get_n_interpolation_points_spatial_indices() const;
492 :
493 : /**
494 : * Set entry of the EIM interpolation matrix.
495 : */
496 : void set_interpolation_matrix_entry(unsigned int i, unsigned int j, Number value);
497 :
498 : /**
499 : * Get the EIM interpolation matrix.
500 : */
501 : const DenseMatrix<Number> & get_interpolation_matrix() const;
502 :
503 : /**
504 : * Add \p bf to our EIM basis.
505 : */
506 : void add_basis_function(
507 : const QpDataMap & bf);
508 :
509 : /**
510 : * Add interpolation data associated with a new basis function.
511 : */
512 : void add_interpolation_data(
513 : Point p,
514 : unsigned int comp,
515 : dof_id_type elem_id,
516 : subdomain_id_type subdomain_id,
517 : unsigned int qp,
518 : const std::vector<Point> & perturbs,
519 : const std::vector<Real> & phi_i_qp,
520 : ElemType elem_type,
521 : const std::vector<Real> & JxW_all_qp,
522 : const std::vector<std::vector<Real>> & phi_i_all_qp,
523 : Order qrule_order,
524 : const Point & dxyz_dxi_elem_center,
525 : const Point & dxyz_deta_elem_center);
526 :
527 : /**
528 : * Add \p side_bf to our EIM basis.
529 : */
530 : void add_side_basis_function(
531 : const SideQpDataMap & side_bf);
532 :
533 : /**
534 : * Add interpolation data associated with a new basis function.
535 : */
536 : void add_side_interpolation_data(
537 : Point p,
538 : unsigned int comp,
539 : dof_id_type elem_id,
540 : unsigned int side_index,
541 : subdomain_id_type subdomain_id,
542 : boundary_id_type boundary_id,
543 : unsigned int qp,
544 : const std::vector<Point> & perturbs,
545 : const std::vector<Real> & phi_i_qp);
546 :
547 : /**
548 : * Add \p node_bf to our EIM basis.
549 : */
550 : void add_node_basis_function(
551 : const NodeDataMap & node_bf);
552 :
553 : /**
554 : * Add interpolation data associated with a new basis function.
555 : */
556 : void add_node_interpolation_data(
557 : Point p,
558 : unsigned int comp,
559 : dof_id_type node_id,
560 : boundary_id_type boundary_id);
561 :
562 : /**
563 : * Set _preserve_rb_eim_solutions.
564 : */
565 : void set_preserve_rb_eim_solutions(bool preserve_rb_eim_solutions);
566 :
567 : /**
568 : * Get _preserve_rb_eim_solutions.
569 : */
570 : bool get_preserve_rb_eim_solutions() const;
571 :
572 : /**
573 : * Write out all the basis functions to file.
574 : * \p sys is used for file IO
575 : * \p directory_name specifies which directory to write files to
576 : * \p read_binary_basis_functions indicates whether to write
577 : * binary or ASCII data
578 : *
579 : * Note: this is not currently a virtual function and is not related
580 : * to the RBEvaluation function of the same name.
581 : */
582 : void write_out_basis_functions(const std::string & directory_name = "offline_data",
583 : bool write_binary_basis_functions = true);
584 :
585 : /**
586 : * Read in all the basis functions from file.
587 : *
588 : * \param sys The Mesh in this System determines the parallel distribution of the basis functions.
589 : * \param directory_name Specifies which directory to write files to.
590 : * \param read_binary_basis_functions Indicates whether to expect binary or ASCII data.
591 : *
592 : * Note: this is not a virtual function and is not related to the
593 : * RBEvaluation function of the same name.
594 : */
595 : void read_in_basis_functions(const System & sys,
596 : const std::string & directory_name = "offline_data",
597 : bool read_binary_basis_functions = true);
598 :
599 : /**
600 : * Project the EIM basis function data stored in \p bf_data onto
601 : * sys.solution. The intent of this function is to work with the
602 : * data format provided by get_interior_basis_functions_as_vecs().
603 : * That format can be easily serialized, if needed, and hence can
604 : * be used to provide \p bf_data after reading in data from disk,
605 : * for example.
606 : *
607 : * \p extra_options can be used to pass extra information to this
608 : * method, e.g. options related to how to perform the projection.
609 : *
610 : * This is a no-op by default, implement in sub-classes if needed.
611 : */
612 : virtual void project_qp_data_vector_onto_system(System & sys,
613 : const std::vector<Number> & bf_data,
614 : const EIMVarGroupPlottingInfo & eim_vargroup,
615 : const std::map<std::string,std::string> & extra_options);
616 :
617 : /**
618 : * Get _eim_vars_to_project_and_write.
619 : */
620 : const std::vector<EIMVarGroupPlottingInfo> & get_eim_vars_to_project_and_write() const;
621 :
622 : /**
623 : * Get _scale_components_in_enrichment.
624 : */
625 : const std::set<unsigned int> & scale_components_in_enrichment() const;
626 :
627 : /**
628 : * Virtual function to indicate if we use the EIM error indicator in this case.
629 : * This indicates if we will generate the data during the Offline training for
630 : * the EIM error indicator. In EIM solves, the error indicator will only
631 : * be used if set_eim_error_indicator_active() is set to true, since we want
632 : * to be able to enable or disable the error indicator depending on the type
633 : * of solve we are doing (e.g. EIM solves during training do not need the
634 : * error indicator).
635 : */
636 : virtual bool use_eim_error_indicator() const;
637 :
638 : /**
639 : * Activate/decative the error indicator in EIM solves. We need this option since
640 : * in some cases (e.g. during EIM training) we do not want to activate the EIM
641 : * error indicator, whereas in "online solves" we do want to activate it.
642 : */
643 : void set_eim_error_indicator_active(bool is_active);
644 :
645 : /**
646 : * Get/set _extra_points_interpolation_matrix.
647 : */
648 : const DenseVector<Number> & get_error_indicator_interpolation_row() const;
649 : void set_error_indicator_interpolation_row(const DenseVector<Number> & error_indicator_row);
650 :
651 : /**
652 : * Evaluates the EIM error indicator based on \p error_indicator_rhs, \p eim_solution,
653 : * and _error_indicator_interpolation_row. We also pass in \p eim_rhs since this is
654 : * used to normalize the error indicator.
655 : *
656 : * We return a pair that specifies the relative error indicator, and the normalization
657 : * that was used to compute the relative error indicator. We can then recover the
658 : * absolute error indicator via rel. indicator x normalization.
659 : */
660 : std::pair<Real,Real> get_eim_error_indicator(
661 : Number error_indicator_rhs,
662 : const DenseVector<Number> & eim_solution,
663 : const DenseVector<Number> & eim_rhs);
664 :
665 : /**
666 : * Get the VectorizedEvalInput data.
667 : */
668 : const VectorizedEvalInput & get_vec_eval_input() const;
669 :
670 : /**
671 : * Initialize the rb_property_map of RBEIMEvaluation (this) from the rb_property_map stored in
672 : * RBParametrizedFunction with empty entries but identical keys.
673 : */
674 : void initialize_rb_property_map();
675 :
676 : /**
677 : * Get all interior basis functions in the form of std::vectors.
678 : * This can provide a convenient format for processing the basis
679 : * function data, or writing it to disk.
680 : *
681 : * Indexing is as follows:
682 : * basis function index --> variable index --> data at qps per elems
683 : *
684 : * In parallel we gather basis functions to processor 0 and hence
685 : * we only return non-empty data on processor 0.
686 : */
687 : std::vector<std::vector<std::vector<Number>>> get_interior_basis_functions_as_vecs();
688 :
689 : /**
690 : * Get data that defines the sizes of interior EIM basis functions.
691 : */
692 : std::map<std::string,std::size_t>
693 : get_interior_basis_function_sizes(std::vector<unsigned int> & n_qp_per_elem);
694 :
695 : /**
696 : * Here we store an enum that defines the type of EIM error indicator
697 : * normalization that we use in get_eim_error_indicator(). The enum
698 : * is public so that it can be set in user code.
699 : *
700 : * RESIDUAL_SUM: Use the sum of the terms in the EIM residual to determine
701 : * the error indicator normalization. This ensures that the error indicator
702 : * value will be at most 1.0, which may be a desirable property of the
703 : * indicator.
704 : *
705 : * RHS: Use only the right-hand side value for the EIM residual to
706 : * determine the error indicator normalization.
707 : *
708 : * MAX_RHS: Use the maximum value in the EIM RHS vector to determine
709 : * the error indicator normalization (default). This is helpful when
710 : * the values at some EIM points are much larger than others, since in
711 : * this scenario we typically want to normalize the error indicator
712 : * based on the largest values in order to avoid overestimating the
713 : * error.
714 : */
715 : enum EimErrorIndicatorNormalization { RESIDUAL_SUM, RESIDUAL_RHS, MAX_RHS };
716 :
717 : EimErrorIndicatorNormalization eim_error_indicator_normalization;
718 :
719 : /**
720 : * If this boolean is true then we clamp EIM error indicator values to
721 : * be at most 1. This is often desirable since we typically think of
722 : * the error indicator as a percentage, and a value of 1 means 100%
723 : * error. We set this to true by default.
724 : */
725 : bool limit_eim_error_indicator_to_one;
726 :
727 : protected:
728 :
729 : /**
730 : * This vector specifies which EIM variables we want to write to disk and/or
731 : * project to nodes for plotting purposes. By default this is an empty
732 : * set, but can be updated in subclasses to specify the EIM variables that
733 : * are relevant for visualization.
734 : *
735 : * We identify groups of variables with one or more variables in a group.
736 : * The purpose of using a group is often we plot multiple components of
737 : * a tensor-valued or vector-valued quantity, so it makes sense to refer
738 : * to the entire group of variables together in those cases.
739 : */
740 : std::vector<EIMVarGroupPlottingInfo> _eim_vars_to_project_and_write;
741 :
742 : /**
743 : * This set that specifies which EIM variables will be scaled during EIM
744 : * enrichment so that their maximum value matches the maximum value across
745 : * all variables. This is helpful in cases where some components are much
746 : * smaller in magnitude than others, since in those cases if we do not apply
747 : * component scaling to the small components then the accuracy of the EIM
748 : * approximation for those components will not be controlled well by the
749 : * EIM enrichment process.
750 : */
751 : std::set<unsigned int> _scale_components_in_enrichment;
752 :
753 : private:
754 :
755 : /**
756 : * Method that writes out element interior EIM basis functions. This may be called by
757 : * write_out_basis_functions().
758 : */
759 : void write_out_interior_basis_functions(const std::string & directory_name,
760 : bool write_binary_basis_functions);
761 :
762 : /**
763 : * Method that writes out element side EIM basis functions. This may be called by
764 : * write_out_basis_functions().
765 : */
766 : void write_out_side_basis_functions(const std::string & directory_name,
767 : bool write_binary_basis_functions);
768 :
769 : /**
770 : * Method that writes out element node EIM basis functions. This may be called by
771 : * write_out_basis_functions().
772 : */
773 : void write_out_node_basis_functions(const std::string & directory_name,
774 : bool write_binary_basis_functions);
775 :
776 : /**
777 : * Method that reads in element interior EIM basis functions. This may be called by
778 : * read_in_basis_functions().
779 : */
780 : void read_in_interior_basis_functions(const System & sys,
781 : const std::string & directory_name,
782 : bool read_binary_basis_functions);
783 :
784 : /**
785 : * Method that reads in element side EIM basis functions. This may be called by
786 : * read_in_basis_functions().
787 : */
788 : void read_in_side_basis_functions(const System & sys,
789 : const std::string & directory_name,
790 : bool read_binary_basis_functions);
791 :
792 : /**
793 : * Method that reads in element node EIM basis functions. This may be called by
794 : * read_in_basis_functions().
795 : */
796 : void read_in_node_basis_functions(const System & sys,
797 : const std::string & directory_name,
798 : bool read_binary_basis_functions);
799 :
800 : /**
801 : * Helper function called by write_out_interior_basis_functions() to
802 : * get basis function \p bf_index stored as a std::vector per variable.
803 : */
804 : std::vector<std::vector<Number>> get_interior_basis_function_as_vec_helper(
805 : unsigned int n_vars,
806 : unsigned int n_qp_data,
807 : unsigned int bf_index);
808 :
809 : /**
810 : * The EIM solution coefficients from the most recent call to rb_eim_solves().
811 : */
812 : std::vector<DenseVector<Number>> _rb_eim_solutions;
813 :
814 : /**
815 : * If we're using the EIM error indicator, then we store the error indicator
816 : * values corresponding to _rb_eim_solutions here.
817 : */
818 : std::vector<std::pair<Real,Real>> _rb_eim_error_indicators;
819 :
820 : /**
821 : * Storage for EIM solutions from the training set. This is typically used in
822 : * the case that we have is_lookup_table==true in our RBParametrizedFunction,
823 : * since in that case we need to store all the EIM solutions on the training
824 : * set so that we do not always need to refer to the lookup table itself
825 : * (since in some cases, like in the Online stage, the lookup table is not
826 : * available).
827 : */
828 : std::vector<DenseVector<Number>> _eim_solutions_for_training_set;
829 :
830 : /**
831 : * The parameters and the number of basis functions that were used in the
832 : * most recent call to rb_eim_solves(). We store this so that we can
833 : * check if we can skip calling rb_eim_solves() again if the inputs
834 : * haven't changed.
835 : */
836 : std::vector<RBParameters> _rb_eim_solves_mus;
837 : unsigned int _rb_eim_solves_N;
838 :
839 : /**
840 : * Dense matrix that stores the lower triangular
841 : * interpolation matrix that can be used
842 : */
843 : DenseMatrix<Number> _interpolation_matrix;
844 :
845 : /**
846 : * We store the EIM interpolation point data in this object.
847 : */
848 : VectorizedEvalInput _vec_eval_input;
849 :
850 : /**
851 : * In the case of a "vector-valued" EIM, this vector determines which
852 : * component of the parameterized function we sample at each EIM point.
853 : */
854 : std::vector<unsigned int> _interpolation_points_comp;
855 :
856 : /**
857 : * Here we store the spatial indices that were initialized by
858 : * initialize_spatial_indices_at_interp_pts(). These are relevant
859 : * in the case that _parametrized_function is defined by indexing
860 : * into separate data based on the mesh-based data.
861 : */
862 : std::vector<std::vector<unsigned int>> _interpolation_points_spatial_indices;
863 :
864 : /**
865 : * Store the parametrized function that will be approximated
866 : * by this EIM system. Note that the parametrized function
867 : * may have more than one component, and each component is
868 : * approximated by a separate variable in the EIM system.
869 : */
870 : std::unique_ptr<RBParametrizedFunction> _parametrized_function;
871 :
872 : /**
873 : * The vector of RBTheta objects that are created to point to
874 : * this RBEIMEvaluation.
875 : */
876 : std::vector<std::unique_ptr<RBTheta>> _rb_eim_theta_objects;
877 :
878 : /**
879 : * The EIM basis functions. We store values at quadrature points
880 : * on elements that are local to this processor. The indexing
881 : * is as follows:
882 : * basis function index --> element ID --> variable --> quadrature point --> value
883 : * We use a map to index the element ID, since the IDs on this processor in
884 : * general will not start at zero.
885 : */
886 : std::vector<QpDataMap> _local_eim_basis_functions;
887 :
888 : /**
889 : * The EIM basis functions on element sides. We store values at quadrature points
890 : * on elements that are local to this processor. The indexing
891 : * is as follows:
892 : * basis function index --> (element ID,side index) --> variable --> quadrature point --> value
893 : * We use a map to index the element ID, since the IDs on this processor in
894 : * general will not start at zero.
895 : */
896 : std::vector<SideQpDataMap> _local_side_eim_basis_functions;
897 :
898 : /**
899 : * The EIM basis functions on element nodes (e.g. on a nodeset). We store values at nodes
900 : * that are local to this processor. The indexing
901 : * is as follows:
902 : * basis function index --> node ID --> variable --> value
903 : * We use a map to index the node ID, since the IDs on this processor in
904 : * general will not start at zero.
905 : */
906 : std::vector<NodeDataMap> _local_node_eim_basis_functions;
907 :
908 : /**
909 : * Print the contents of _local_eim_basis_functions to libMesh::out.
910 : * Helper function mainly useful for debugging.
911 : */
912 : void print_local_eim_basis_functions() const;
913 :
914 : /**
915 : * Helper function that gathers the contents of
916 : * _local_eim_basis_functions to processor 0 in preparation for
917 : * printing to file.
918 : */
919 : void gather_bfs();
920 :
921 : /**
922 : * Same as gather_bfs() except for side data.
923 : */
924 : void side_gather_bfs();
925 :
926 : /**
927 : * Same as gather_bfs() except for node data.
928 : */
929 : void node_gather_bfs();
930 :
931 : /**
932 : * Helper function that distributes the entries of
933 : * _local_eim_basis_functions to their respective processors after
934 : * they are read in on processor 0.
935 : */
936 : void distribute_bfs(const System & sys);
937 :
938 : /**
939 : * Same as distribute_bfs() except for side data.
940 : */
941 : void side_distribute_bfs(const System & sys);
942 :
943 : /**
944 : * Same as distribute_bfs() except for node data.
945 : */
946 : void node_distribute_bfs(const System & sys);
947 :
948 : /**
949 : * Boolean to indicate if we skip updating _rb_eim_solutions in rb_eim_solves().
950 : * This is relevant for cases when we set up _rb_eim_solutions elsewhere and we
951 : * want to avoid changing it.
952 : */
953 : bool _preserve_rb_eim_solutions;
954 :
955 : /**
956 : * Indicate if the EIM error indicator is active in RB EIM solves. Note that
957 : * this is distinct from use_eim_error_indicator(), since use_eim_error_indicator()
958 : * indicates if this RBEIMEvaluation has an EIM error indicator defined,
959 : * whereas _is_eim_error_indicator_active is used to turn on or off the
960 : * error indicator. This primary purpose of _is_eim_error_indicator_active
961 : * is to turn the error indicator off during EIM training (when it is not relevant)
962 : * and to turn it on during "online solves".
963 : */
964 : bool _is_eim_error_indicator_active;
965 :
966 : /**
967 : * Here we store an extra row of the interpolation matrix which is used to
968 : * compute the EIM error indicator. This stores the EIM basis function
969 : * values at the extra point associated with the error indicator.
970 : */
971 : DenseVector<Number> _error_indicator_interpolation_row;
972 :
973 : };
974 :
975 : }
976 :
977 : #endif // LIBMESH_RB_EIM_EVALUATION_H
|
