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