rb_eim_evaluation.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_EVALUATION_H
#define LIBMESH_RB_EIM_EVALUATION_H
 
// libMesh includes
#include "libmesh/point.h"
#include "libmesh/rb_theta_expansion.h"
#include "libmesh/rb_parametrized.h"
#include "libmesh/parallel_object.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/rb_parametrized_function.h"
#include "libmesh/fe_type.h"
 
// C++ includes
#include <memory>
#include <map>
#include <vector>
#include <string>
 
namespace libMesh
{
 
class RBParameters;
class RBParametrizedFunction;
class RBTheta;
class System;
class EquationSystems;
class Elem;
 
struct EIMVarGroupPlottingInfo
{
 
  EIMVarGroupPlottingInfo()
  :
  first_eim_var_index(0),
  n_eim_vars(0),
  enforce_min_value(false),
  enforce_max_value(false),
  min_value(0.),
  max_value(0.)
  {
  }
 
  unsigned int first_eim_var_index;
 
  unsigned int n_eim_vars;
 
  std::string eim_sys_name;
 
  bool enforce_min_value;
  bool enforce_max_value;
 
  Real min_value;
  Real max_value;
};
 
class RBEIMEvaluation : public RBParametrized,
                        public ParallelObject
{
public:
 
  RBEIMEvaluation(const Parallel::Communicator & comm);
 
  RBEIMEvaluation (RBEIMEvaluation &&) = default;
  RBEIMEvaluation (const RBEIMEvaluation &) = delete;
  RBEIMEvaluation & operator= (const RBEIMEvaluation &) = delete;
  RBEIMEvaluation & operator= (RBEIMEvaluation &&) = default;
  virtual ~RBEIMEvaluation ();
 
  typedef std::map<dof_id_type, std::vector<std::vector<Number>>> QpDataMap;
 
  typedef std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Number>>> SideQpDataMap;
 
  typedef std::map<dof_id_type, std::vector<Number>> NodeDataMap;
 
  virtual void clear() override;
 
  void resize_data_structures(const unsigned int Nmax);
 
  void set_parametrized_function(std::unique_ptr<RBParametrizedFunction> pf);
 
  RBParametrizedFunction & get_parametrized_function();
 
  const RBParametrizedFunction & get_parametrized_function() const;
 
  DenseVector<Number> rb_eim_solve(DenseVector<Number> & EIM_rhs);
 
  void rb_eim_solves(const std::vector<RBParameters> & mus, unsigned int N);
 
  void initialize_interpolation_points_spatial_indices();
 
  void initialize_param_fn_spatial_indices();
 
  unsigned int get_n_basis_functions() const;
 
  unsigned int get_n_interpolation_points() const;
 
  unsigned int get_n_elems() const;
 
  unsigned int get_n_properties() const;
 
  void set_n_basis_functions(unsigned int n_bfs);
 
  void decrement_vector(QpDataMap & v,
                        const DenseVector<Number> & coeffs);
 
  void side_decrement_vector(SideQpDataMap & v,
                             const DenseVector<Number> & coeffs);
 
  void node_decrement_vector(NodeDataMap & v,
                             const DenseVector<Number> & coeffs);
 
  void initialize_eim_theta_objects();
 
  std::vector<std::unique_ptr<RBTheta>> & get_eim_theta_objects();
 
  virtual std::unique_ptr<RBTheta> build_eim_theta(unsigned int index);
 
  static void get_parametrized_function_values_at_qps(
    const QpDataMap & pf,
    dof_id_type elem_id,
    unsigned int comp,
    std::vector<Number> & values);
 
  static void get_parametrized_function_side_values_at_qps(
    const SideQpDataMap & pf,
    dof_id_type elem_id,
    unsigned int side_index,
    unsigned int comp,
    std::vector<Number> & values);
 
  static Number get_parametrized_function_node_local_value(
    const NodeDataMap & pf,
    dof_id_type node_id,
    unsigned int comp);
 
  static Number get_parametrized_function_value(
    const Parallel::Communicator & comm,
    const QpDataMap & pf,
    dof_id_type elem_id,
    unsigned int comp,
    unsigned int qp);
 
  static Number get_parametrized_function_side_value(
    const Parallel::Communicator & comm,
    const SideQpDataMap & pf,
    dof_id_type elem_id,
    unsigned int side_index,
    unsigned int comp,
    unsigned int qp);
 
  static Number get_parametrized_function_node_value(
    const Parallel::Communicator & comm,
    const NodeDataMap & pf,
    dof_id_type node_id,
    unsigned int comp);
 
  void get_eim_basis_function_values_at_qps(unsigned int basis_function_index,
                                            dof_id_type elem_id,
                                            unsigned int var,
                                            std::vector<Number> & values) const;
 
  void get_eim_basis_function_side_values_at_qps(unsigned int basis_function_index,
                                                 dof_id_type elem_id,
                                                 unsigned int side_index,
                                                 unsigned int var,
                                                 std::vector<Number> & values) const;
 
  Number get_eim_basis_function_node_local_value(unsigned int basis_function_index,
                                                 dof_id_type node_id,
                                                 unsigned int var) const;
 
  Number get_eim_basis_function_value(unsigned int basis_function_index,
                                      dof_id_type elem_id,
                                      unsigned int comp,
                                      unsigned int qp) const;
 
  Number get_eim_basis_function_side_value(unsigned int basis_function_index,
                                          dof_id_type elem_id,
                                          unsigned int side_index,
                                          unsigned int comp,
                                          unsigned int qp) const;
 
  Number get_eim_basis_function_node_value(unsigned int basis_function_index,
                                           dof_id_type node_id,
                                           unsigned int var) const;
 
  const QpDataMap & get_basis_function(unsigned int i) const;
 
  const SideQpDataMap & get_side_basis_function(unsigned int i) const;
 
  const NodeDataMap & get_node_basis_function(unsigned int i) const;
 
  void set_rb_eim_solutions(const std::vector<DenseVector<Number>> & rb_eim_solutions);
 
  const std::vector<DenseVector<Number>> & get_rb_eim_solutions() const;
 
  std::vector<Number> get_rb_eim_solutions_entries(unsigned int index) const;
 
  const std::vector<DenseVector<Number>> & get_eim_solutions_for_training_set() const;
 
  std::vector<DenseVector<Number>> & get_eim_solutions_for_training_set();
 
  const std::vector<std::pair<Real,Real>> & get_rb_eim_error_indicators() const;
 
  void add_interpolation_points_xyz(Point p);
  void add_interpolation_points_comp(unsigned int comp);
  void add_interpolation_points_subdomain_id(subdomain_id_type sbd_id);
  void add_interpolation_points_boundary_id(boundary_id_type b_id);
  void add_interpolation_points_xyz_perturbations(const std::vector<Point> & perturbs);
  void add_interpolation_points_elem_id(dof_id_type elem_id);
  void add_interpolation_points_side_index(unsigned int side_index);
  void add_interpolation_points_node_id(dof_id_type node_id);
  void add_interpolation_points_qp(unsigned int qp);
  void add_interpolation_points_elem_type(ElemType elem_type);
  void add_interpolation_points_phi_i_qp(const std::vector<Real> & phi_i_qp);
  void add_interpolation_points_JxW_all_qp(const std::vector<Real> & JxW_all_qp);
  void add_interpolation_points_phi_i_all_qp(const std::vector<std::vector<Real>> & phi_i_all_qp);
  void add_interpolation_points_qrule_order(Order qrule_order);
  void add_elem_center_dxyzdxi(const Point & dxyzdxi);
  void add_elem_center_dxyzdeta(const Point & dxyzdxi);
  void add_interpolation_points_spatial_indices(const std::vector<unsigned int> & spatial_indices);
  void add_elem_id_local_index_map_entry(dof_id_type elem_id, unsigned int local_index);
  void add_rb_property_map_entry(std::string & property_name, std::set<dof_id_type> & entity_ids);
 
  Point get_interpolation_points_xyz(unsigned int index) const;
  unsigned int get_interpolation_points_comp(unsigned int index) const;
  subdomain_id_type get_interpolation_points_subdomain_id(unsigned int index) const;
  boundary_id_type get_interpolation_points_boundary_id(unsigned int index) const;
  const std::vector<Point> & get_interpolation_points_xyz_perturbations(unsigned int index) const;
  dof_id_type get_interpolation_points_elem_id(unsigned int index) const;
  unsigned int get_interpolation_points_side_index(unsigned int index) const;
  dof_id_type get_interpolation_points_node_id(unsigned int index) const;
  unsigned int get_interpolation_points_qp(unsigned int index) const;
  ElemType get_interpolation_points_elem_type(unsigned int index) const;
  const std::vector<Real> & get_interpolation_points_phi_i_qp(unsigned int index) const;
  const std::vector<Real> & get_interpolation_points_JxW_all_qp(unsigned int index) const;
  const std::vector<std::vector<Real>> & get_interpolation_points_phi_i_all_qp(unsigned int index) const;
  Order get_interpolation_points_qrule_order(unsigned int index) const;
  const Point & get_elem_center_dxyzdxi(unsigned int index) const;
  const Point & get_elem_center_dxyzdeta(unsigned int index) const;
  const std::vector<unsigned int> & get_interpolation_points_spatial_indices(unsigned int index) const;
  const std::map<dof_id_type, unsigned int> & get_elem_id_to_local_index_map() const;
  const std::unordered_map<std::string, std::set<dof_id_type>> & get_rb_property_map() const;
 
  unsigned int get_n_interpolation_points_spatial_indices() const;
 
  void set_interpolation_matrix_entry(unsigned int i, unsigned int j, Number value);
 
  const DenseMatrix<Number> & get_interpolation_matrix() const;
 
  void add_basis_function(
    const QpDataMap & bf);
 
  void add_interpolation_data(
    Point p,
    unsigned int comp,
    dof_id_type elem_id,
    subdomain_id_type subdomain_id,
    unsigned int qp,
    const std::vector<Point> & perturbs,
    const std::vector<Real> & phi_i_qp,
    ElemType elem_type,
    const std::vector<Real> & JxW_all_qp,
    const std::vector<std::vector<Real>> & phi_i_all_qp,
    Order qrule_order,
    const Point & dxyz_dxi_elem_center,
    const Point & dxyz_deta_elem_center);
 
  void add_side_basis_function(
    const SideQpDataMap & side_bf);
 
  void add_side_interpolation_data(
    Point p,
    unsigned int comp,
    dof_id_type elem_id,
    unsigned int side_index,
    subdomain_id_type subdomain_id,
    boundary_id_type boundary_id,
    unsigned int qp,
    const std::vector<Point> & perturbs,
    const std::vector<Real> & phi_i_qp);
 
  void add_node_basis_function(
    const NodeDataMap & node_bf);
 
  void add_node_interpolation_data(
    Point p,
    unsigned int comp,
    dof_id_type node_id,
    boundary_id_type boundary_id);
 
  void set_preserve_rb_eim_solutions(bool preserve_rb_eim_solutions);
 
  bool get_preserve_rb_eim_solutions() const;
 
  void write_out_basis_functions(const std::string & directory_name = "offline_data",
                                 bool write_binary_basis_functions = true);
 
  void read_in_basis_functions(const System & sys,
                               const std::string & directory_name = "offline_data",
                               bool read_binary_basis_functions = true);
 
  virtual void project_qp_data_vector_onto_system(System & sys,
                                                  const std::vector<Number> & bf_data,
                                                  const EIMVarGroupPlottingInfo & eim_vargroup,
                                                  const std::map<std::string,std::string> & extra_options);
 
  const std::vector<EIMVarGroupPlottingInfo> & get_eim_vars_to_project_and_write() const;
 
  const std::set<unsigned int> & scale_components_in_enrichment() const;
 
  virtual bool use_eim_error_indicator() const;
 
  void set_eim_error_indicator_active(bool is_active);
 
  const DenseVector<Number> & get_error_indicator_interpolation_row() const;
  void set_error_indicator_interpolation_row(const DenseVector<Number> & error_indicator_row);
 
  std::pair<Real,Real> get_eim_error_indicator(
    Number error_indicator_rhs,
    const DenseVector<Number> & eim_solution,
    const DenseVector<Number> & eim_rhs);
 
  const VectorizedEvalInput & get_vec_eval_input() const;
 
  void initialize_rb_property_map();
 
  std::vector<std::vector<std::vector<Number>>> get_interior_basis_functions_as_vecs();
 
  std::map<std::string,std::size_t>
    get_interior_basis_function_sizes(std::vector<unsigned int> & n_qp_per_elem);
 
  enum EimErrorIndicatorNormalization { RESIDUAL_SUM, RESIDUAL_RHS, MAX_RHS };
 
  EimErrorIndicatorNormalization eim_error_indicator_normalization;
 
  bool limit_eim_error_indicator_to_one;
 
protected:
 
  std::vector<EIMVarGroupPlottingInfo> _eim_vars_to_project_and_write;
 
  std::set<unsigned int> _scale_components_in_enrichment;
 
private:
 
  void write_out_interior_basis_functions(const std::string & directory_name,
                                          bool write_binary_basis_functions);
 
  void write_out_side_basis_functions(const std::string & directory_name,
                                      bool write_binary_basis_functions);
 
  void write_out_node_basis_functions(const std::string & directory_name,
                                      bool write_binary_basis_functions);
 
  void read_in_interior_basis_functions(const System & sys,
                                        const std::string & directory_name,
                                        bool read_binary_basis_functions);
 
  void read_in_side_basis_functions(const System & sys,
                                    const std::string & directory_name,
                                    bool read_binary_basis_functions);
 
  void read_in_node_basis_functions(const System & sys,
                                    const std::string & directory_name,
                                    bool read_binary_basis_functions);
 
  std::vector<std::vector<Number>> get_interior_basis_function_as_vec_helper(
    unsigned int n_vars,
    unsigned int n_qp_data,
    unsigned int bf_index);
 
  std::vector<DenseVector<Number>> _rb_eim_solutions;
 
  std::vector<std::pair<Real,Real>> _rb_eim_error_indicators;
 
  std::vector<DenseVector<Number>> _eim_solutions_for_training_set;
 
  std::vector<RBParameters> _rb_eim_solves_mus;
  unsigned int _rb_eim_solves_N;
 
  DenseMatrix<Number> _interpolation_matrix;
 
  VectorizedEvalInput _vec_eval_input;
 
  std::vector<unsigned int> _interpolation_points_comp;
 
  std::vector<std::vector<unsigned int>> _interpolation_points_spatial_indices;
 
  std::unique_ptr<RBParametrizedFunction> _parametrized_function;
 
  std::vector<std::unique_ptr<RBTheta>> _rb_eim_theta_objects;
 
  std::vector<QpDataMap> _local_eim_basis_functions;
 
  std::vector<SideQpDataMap> _local_side_eim_basis_functions;
 
  std::vector<NodeDataMap> _local_node_eim_basis_functions;
 
  void print_local_eim_basis_functions() const;
 
  void gather_bfs();
 
  void side_gather_bfs();
 
  void node_gather_bfs();
 
  void distribute_bfs(const System & sys);
 
  void side_distribute_bfs(const System & sys);
 
  void node_distribute_bfs(const System & sys);
 
  bool _preserve_rb_eim_solutions;
 
  bool _is_eim_error_indicator_active;
 
  DenseVector<Number> _error_indicator_interpolation_row;
 
};
 
}
 
#endif // LIBMESH_RB_EIM_EVALUATION_H