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