rb_eim_construction.h

Go to the documentation of this file.

// rbOOmit: An implementation of the Certified Reduced Basis method.
// Copyright (C) 2009, 2010 David J. Knezevic
 
// This file is part of rbOOmit.
 
// rbOOmit is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// rbOOmit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
#ifndef LIBMESH_RB_EIM_CONSTRUCTION_H
#define LIBMESH_RB_EIM_CONSTRUCTION_H
 
// rbOOmit includes
#include "libmesh/rb_construction.h"
#include "libmesh/rb_assembly_expansion.h"
#include "libmesh/rb_eim_assembly.h"
#include "libmesh/rb_eim_evaluation.h"
 
// libMesh includes
#include "libmesh/mesh_function.h"
#include "libmesh/coupling_matrix.h"
 
// C++ includes
#include <unordered_map>
#include <map>
#include <string>
#include <memory>
#include <vector>
 
namespace libMesh
{
 
struct EimPointData
{
  dof_id_type elem_id;
  dof_id_type node_id;
  unsigned int side_index;
  unsigned int comp_index;
  unsigned int qp_index;
};
 
class RBEIMConstruction : public RBConstructionBase<System>
{
public:
 
  enum BEST_FIT_TYPE { PROJECTION_BEST_FIT, EIM_BEST_FIT, POD_BEST_FIT };
 
  RBEIMConstruction (EquationSystems & es,
                     const std::string & name,
                     const unsigned int number);
 
  RBEIMConstruction (RBEIMConstruction &&) = default;
  RBEIMConstruction (const RBEIMConstruction &) = delete;
  RBEIMConstruction & operator= (const RBEIMConstruction &) = delete;
  RBEIMConstruction & operator= (RBEIMConstruction &&) = delete;
  virtual ~RBEIMConstruction ();
 
  typedef RBEIMEvaluation::QpDataMap QpDataMap;
 
  typedef RBEIMEvaluation::SideQpDataMap SideQpDataMap;
 
  typedef RBEIMEvaluation::NodeDataMap NodeDataMap;
 
  virtual void clear() override;
 
  void set_rb_eim_evaluation(RBEIMEvaluation & rb_eim_eval_in);
 
  RBEIMEvaluation & get_rb_eim_evaluation();
 
  const RBEIMEvaluation & get_rb_eim_evaluation() const;
 
  void initialize_eim_construction();
 
  virtual void process_parameters_file (const std::string & parameters_filename);
 
  void set_rb_construction_parameters(unsigned int n_training_samples_in,
                                      bool deterministic_training_in,
                                      int training_parameters_random_seed_in,
                                      bool quiet_mode_in,
                                      unsigned int Nmax_in,
                                      Real rel_training_tolerance_in,
                                      Real abs_training_tolerance_in,
                                      const RBParameters & mu_min_in,
                                      const RBParameters & mu_max_in,
                                      const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
                                      const std::map<std::string,bool> & log_scaling,
                                      std::map<std::string, std::vector<RBParameter>> * training_sample_list=nullptr);
 
  virtual void set_best_fit_type_flag (const std::string & best_fit_type_string);
 
  virtual void print_info();
 
  void apply_normalization_to_solution_snapshots();
 
  virtual Real train_eim_approximation();
 
  virtual Real train_eim_approximation_with_greedy();
 
  virtual Real train_eim_approximation_with_POD();
 
  virtual void initialize_eim_assembly_objects();
 
  std::vector<std::unique_ptr<ElemAssembly>> & get_eim_assembly_objects();
 
  virtual std::unique_ptr<ElemAssembly> build_eim_assembly(unsigned int bf_index) = 0;
 
  virtual void init_context(FEMContext &);
 
  void set_rel_training_tolerance(Real new_training_tolerance);
  Real get_rel_training_tolerance();
 
  void set_abs_training_tolerance(Real new_training_tolerance);
  Real get_abs_training_tolerance();
 
  unsigned int get_Nmax() const;
  virtual void set_Nmax(unsigned int Nmax);
 
  void enable_set_Nmax_from_n_snapshots(int increment);
 
  void disable_set_Nmax_from_n_snapshots();
 
  Real get_max_abs_value_in_training_set() const;
 
  void store_eim_solutions_for_training_set();
 
  const QpDataMap & get_parametrized_function_from_training_set(unsigned int training_index) const;
  const SideQpDataMap & get_side_parametrized_function_from_training_set(unsigned int training_index) const;
  const NodeDataMap & get_node_parametrized_function_from_training_set(unsigned int training_index) const;
 
  const std::unordered_map<dof_id_type, std::vector<Real> > & get_local_quad_point_JxW();
  const std::map<std::pair<dof_id_type,unsigned int>, std::vector<Real> > & get_local_side_quad_point_JxW();
 
  unsigned int get_n_parametrized_functions_for_training() const;
 
  void reinit_eim_projection_matrix();
 
  BEST_FIT_TYPE best_fit_type_flag;
 
protected:
 
  bool enrich_eim_approximation_on_sides(const SideQpDataMap & side_pf,
                                         bool add_basis_function,
                                         EimPointData * eim_point_data);
 
  bool enrich_eim_approximation_on_nodes(const NodeDataMap & node_pf,
                                         bool add_basis_function,
                                         EimPointData * eim_point_data);
 
  bool enrich_eim_approximation_on_interiors(const QpDataMap & interior_pf,
                                             bool add_basis_function,
                                             EimPointData * eim_point_data);
 
  void update_eim_matrices(bool set_eim_error_indicator);
 
private:
 
  std::pair<Real, unsigned int> compute_max_eim_error();
 
  void initialize_parametrized_functions_in_training_set();
 
  void initialize_qp_data();
 
  Number inner_product(const QpDataMap & v, const QpDataMap & w, bool apply_comp_scaling);
 
  Number side_inner_product(const SideQpDataMap & v, const SideQpDataMap & w, bool apply_comp_scaling);
 
  Number node_inner_product(const NodeDataMap & v, const NodeDataMap & w, bool apply_comp_scaling);
 
  template <class DataMap>
  Real get_max_abs_value(const DataMap & v) const
  {
    Real max_value = 0.;
 
    for (const auto & pr : v)
      {
        const auto & v_comp_and_qp = pr.second;
 
        for (const auto & comp : index_range(v_comp_and_qp))
          {
            Real comp_scaling = 1.;
            if (get_rb_eim_evaluation().scale_components_in_enrichment().count(comp))
              {
                // Make sure that _component_scaling_in_training_set is initialized
                libmesh_error_msg_if(comp >= _component_scaling_in_training_set.size(),
                                    "Invalid vector index");
                comp_scaling = _component_scaling_in_training_set[comp];
              }
 
            const std::vector<Number> & v_qp = v_comp_and_qp[comp];
            for (Number value : v_qp)
              max_value = std::max(max_value, std::abs(value * comp_scaling));
          }
      }
 
    comm().max(max_value);
    return max_value;
  }
 
  Real get_node_max_abs_value(const NodeDataMap & v) const;
 
  void enrich_eim_approximation(unsigned int training_index,
                                bool add_basis_function,
                                EimPointData * eim_point_data);
 
  template<class DataMap>
  static void scale_parametrized_function(DataMap & local_pf,
                                          Number scaling_factor)
  {
    for (auto & pr : local_pf)
      {
        auto & comp_and_qp = pr.second;
 
        for (unsigned int comp : index_range(comp_and_qp))
          {
            std::vector<Number> & qp_values = comp_and_qp[comp];
 
            for (unsigned int qp : index_range(qp_values))
              {
                qp_values[qp] *= scaling_factor;
              }
          }
      }
  }
 
  static void scale_node_parametrized_function(NodeDataMap & local_pf,
                                               Number scaling_factor);
 
  static unsigned int get_random_int_0_to_n(unsigned int n);
 
  EimPointData get_random_point(const QpDataMap & v);
  EimPointData get_random_point(const SideQpDataMap & v);
  EimPointData get_random_point(const NodeDataMap & v);
 
  EimPointData get_random_point_from_training_sample();
 
  unsigned int _Nmax;
 
  bool _set_Nmax_from_n_snapshots;
  int _Nmax_from_n_snapshots_increment;
 
  Real _rel_training_tolerance;
  Real _abs_training_tolerance;
 
  DenseMatrix<Number> _eim_projection_matrix;
 
  RBEIMEvaluation * _rb_eim_eval;
 
  std::vector<std::unique_ptr<ElemAssembly>> _rb_eim_assembly_objects;
 
  std::vector<QpDataMap> _local_parametrized_functions_for_training;
 
  std::vector<SideQpDataMap> _local_side_parametrized_functions_for_training;
 
  std::vector<NodeDataMap> _local_node_parametrized_functions_for_training;
 
  Real _max_abs_value_in_training_set;
 
  unsigned int _max_abs_value_in_training_set_index;
 
  std::vector<Real> _component_scaling_in_training_set;
 
  std::unordered_map<dof_id_type, std::vector<Point> > _local_quad_point_locations;
  std::unordered_map<dof_id_type, std::vector<Real> > _local_quad_point_JxW;
  std::unordered_map<dof_id_type, subdomain_id_type > _local_quad_point_subdomain_ids;
 
  std::unordered_map<dof_id_type, std::vector<std::vector<Point>> > _local_quad_point_locations_perturbations;
 
  std::map<std::pair<dof_id_type,unsigned int>, std::vector<Point> > _local_side_quad_point_locations;
  std::map<std::pair<dof_id_type,unsigned int>, std::vector<Real> > _local_side_quad_point_JxW;
  std::map<std::pair<dof_id_type,unsigned int>, subdomain_id_type > _local_side_quad_point_subdomain_ids;
  std::map<std::pair<dof_id_type,unsigned int>, boundary_id_type > _local_side_quad_point_boundary_ids;
  std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Point>> > _local_side_quad_point_locations_perturbations;
 
  std::unordered_map<dof_id_type, Point > _local_node_locations;
  std::unordered_map<dof_id_type, boundary_id_type > _local_node_boundary_ids;
 
  std::map<std::pair<dof_id_type,unsigned int>, unsigned int > _local_side_quad_point_side_types;
 
};
 
} // namespace libMesh
 
#endif // LIBMESH_RB_EIM_CONSTRUCTION_H