rb_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_CONSTRUCTION_H
#define LIBMESH_RB_CONSTRUCTION_H
 
// rbOOmit includes
#include "libmesh/rb_construction_base.h"
 
// libMesh includes
#include "libmesh/linear_implicit_system.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dg_fem_context.h"
#include "libmesh/dirichlet_boundaries.h"
 
// C++ includes
 
namespace libMesh
{
 
class RBThetaExpansion;
class RBAssemblyExpansion;
class RBEvaluation;
class ElemAssembly;
 
class RBConstruction : public RBConstructionBase<LinearImplicitSystem>
{
public:
 
  RBConstruction (EquationSystems & es,
                  const std::string & name,
                  const unsigned int number);
 
  RBConstruction (RBConstruction &&) = default;
  RBConstruction (const RBConstruction &) = delete;
  RBConstruction & operator= (const RBConstruction &) = delete;
  RBConstruction & operator= (RBConstruction &&) = delete;
  virtual ~RBConstruction ();
 
  typedef RBConstruction sys_type;
 
  virtual void solve_for_matrix_and_rhs (LinearSolver<Number> & input_solver,
                                         SparseMatrix<Number> & input_matrix,
                                         NumericVector<Number> & input_rhs);
 
  void set_rb_evaluation(RBEvaluation & rb_eval_in);
 
  RBEvaluation & get_rb_evaluation();
  const RBEvaluation & get_rb_evaluation() const;
 
  bool is_rb_eval_initialized() const;
 
  RBThetaExpansion & get_rb_theta_expansion();
  const RBThetaExpansion & get_rb_theta_expansion() const;
 
  void set_rb_assembly_expansion(RBAssemblyExpansion & rb_assembly_expansion_in);
 
  RBAssemblyExpansion & get_rb_assembly_expansion();
 
  sys_type & system () { return *this; }
 
  typedef RBConstructionBase<LinearImplicitSystem> Parent;
 
  virtual void clear () override;
 
  virtual std::string system_type () const override;
 
  virtual Real truth_solve(int plot_solution);
 
  virtual Real train_reduced_basis(const bool resize_rb_eval_data=true);
 
  Real train_reduced_basis_with_greedy(const bool resize_rb_eval_data);
 
  void enrich_basis_from_rhs_terms(const bool resize_rb_eval_data=true);
 
  void train_reduced_basis_with_POD();
 
  virtual Real compute_max_error_bound();
 
  const RBParameters & get_greedy_parameter(unsigned int i);
 
  void set_rel_training_tolerance(Real new_training_tolerance)
  {this->rel_training_tolerance = new_training_tolerance; }
  Real get_rel_training_tolerance() const { return rel_training_tolerance; }
 
  void set_abs_training_tolerance(Real new_training_tolerance)
  {this->abs_training_tolerance = new_training_tolerance; }
  Real get_abs_training_tolerance() const { return abs_training_tolerance; }
 
  void set_normalize_rb_bound_in_greedy(bool normalize_rb_bound_in_greedy_in)
  {this->normalize_rb_bound_in_greedy = normalize_rb_bound_in_greedy_in; }
  bool get_normalize_rb_bound_in_greedy() const { return normalize_rb_bound_in_greedy; }
 
  virtual bool is_serial_training_type(const std::string & RB_training_type_in);
 
  void set_RB_training_type(const std::string & RB_training_type_in);
  const std::string & get_RB_training_type() const;
 
  unsigned int get_Nmax() const    { return Nmax; }
  virtual void set_Nmax(unsigned int Nmax);
 
  virtual void load_basis_function(unsigned int i);
 
  virtual void load_rb_solution();
 
  Real compute_residual_dual_norm_slow(const unsigned int N);
 
  SparseMatrix<Number> * get_inner_product_matrix();
  const SparseMatrix<Number> * get_inner_product_matrix() const;
 
  SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix();
  const SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix() const;
 
  SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix_if_avail();
  const SparseMatrix<Number> * get_non_dirichlet_inner_product_matrix_if_avail() const;
 
  SparseMatrix<Number> * get_Aq(unsigned int q);
 
  SparseMatrix<Number> * get_non_dirichlet_Aq(unsigned int q);
 
  SparseMatrix<Number> * get_non_dirichlet_Aq_if_avail(unsigned int q);
 
  virtual void initialize_rb_construction(bool skip_matrix_assembly=false,
                                          bool skip_vector_assembly=false);
 
  NumericVector<Number> * get_Fq(unsigned int q);
 
  NumericVector<Number> * get_non_dirichlet_Fq(unsigned int q);
 
  NumericVector<Number> * get_non_dirichlet_Fq_if_avail(unsigned int q);
 
  NumericVector<Number> * get_output_vector(unsigned int n, unsigned int q_l);
 
  NumericVector<Number> * get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l);
 
  virtual void get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices);
 
  virtual void get_all_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors);
 
  virtual void get_output_vectors(std::map<std::string, NumericVector<Number> *> & all_vectors);
 
  virtual void assemble_affine_expansion(bool skip_matrix_assembly,
                                         bool skip_vector_assembly);
 
  void assemble_inner_product_matrix(SparseMatrix<Number> * input_matrix, bool apply_dof_constraints=true);
 
  void assemble_Aq_matrix(unsigned int q, SparseMatrix<Number> * input_matrix, bool apply_dof_constraints=true);
 
  void assemble_Fq_vector(unsigned int q, NumericVector<Number> * input_vector, bool apply_dof_constraints=true);
 
  void add_scaled_Aq(Number scalar, unsigned int q_a,
                     SparseMatrix<Number> * input_matrix,
                     bool symmetrize);
 
  virtual void write_riesz_representors_to_files(const std::string & riesz_representors_dir,
                                                 const bool write_binary_residual_representors);
 
  virtual void read_riesz_representors_from_files(const std::string & riesz_representors_dir,
                                                  const bool write_binary_residual_representors);
 
 
  virtual void recompute_all_residual_terms(const bool compute_inner_products=true);
 
  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,
                                      bool normalize_rb_error_bound_in_greedy_in,
                                      const std::string & RB_training_type_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 print_info() const;
 
  void print_basis_function_orthogonality() const;
 
  unsigned int get_delta_N() const { return delta_N; }
 
  void set_inner_product_assembly(ElemAssembly & inner_product_assembly_in);
 
  ElemAssembly & get_inner_product_assembly();
 
  void set_energy_inner_product(const std::vector<Number> & energy_inner_product_coeffs_in);
 
  void zero_constrained_dofs_on_vector(NumericVector<Number> & vector) const;
 
  virtual bool check_if_zero_truth_solve() const;
 
#ifdef LIBMESH_ENABLE_DIRICHLET
  static std::unique_ptr<DirichletBoundary> build_zero_dirichlet_boundary_object();
#endif
 
  void set_convergence_assertion_flag(bool flag);
 
  bool get_preevaluate_thetas_flag() const;
  void set_preevaluate_thetas_flag(bool flag);
 
 
  //----------- PUBLIC DATA MEMBERS -----------//
 
  std::vector<Real> training_error_bounds;
 
  std::unique_ptr<LinearSolver<Number>> inner_product_solver;
 
  LinearSolver<Number> * extra_linear_solver;
 
  std::unique_ptr<SparseMatrix<Number>> inner_product_matrix;
 
  std::vector<Number > truth_outputs;
 
  std::vector<std::vector<Number >> output_dual_innerprods;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> Fq_representor;
 
  std::vector<Number> Fq_representor_innerprods;
 
  bool skip_residual_in_train_reduced_basis;
 
  bool exit_on_repeated_greedy_parameters;
 
  bool impose_internal_fluxes;
 
  bool skip_degenerate_sides;
 
  bool compute_RB_inner_product;
 
  bool store_dirichlet_operators;
 
  bool store_non_dirichlet_operators;
 
  bool store_untransformed_basis;
 
  bool use_empty_rb_solve_in_greedy;
 
  bool Fq_representor_innerprods_computed;
 
protected:
 
  virtual void allocate_data_structures();
 
  virtual void truth_assembly();
 
  virtual std::unique_ptr<DGFEMContext> build_context();
 
  virtual SparseMatrix<Number> & get_matrix_for_output_dual_solves();
 
  virtual bool greedy_termination_test(Real abs_greedy_error,
                                       Real initial_greedy_error,
                                       int count);
 
  void update_greedy_param_list();
 
  void add_scaled_matrix_and_vector(Number scalar,
                                    ElemAssembly * elem_assembly,
                                    SparseMatrix<Number> * input_matrix,
                                    NumericVector<Number> * input_vector,
                                    bool symmetrize=false,
                                    bool apply_dof_constraints=true);
 
  virtual void post_process_elem_matrix_and_vector
  (DGFEMContext & /*context*/) {}
 
  virtual void post_process_truth_solution() {}
 
  virtual void set_context_solution_vec(NumericVector<Number> & vec);
 
  virtual void assemble_misc_matrices();
 
  virtual void assemble_all_affine_operators();
 
  virtual void assemble_all_affine_vectors();
 
  virtual void assemble_all_output_vectors();
 
  virtual void compute_output_dual_innerprods();
 
  virtual void compute_Fq_representor_innerprods(bool compute_inner_products=true);
 
  virtual void enrich_RB_space();
 
  virtual void update_system();
 
  virtual Real get_RB_error_bound();
 
  virtual void update_RB_system_matrices();
 
  virtual void update_residual_terms(bool compute_inner_products=true);
 
  virtual void init_context(FEMContext &)
  {
    // Failing to rederive init_context() means your FE objects don't
    // know what to compute.
    libmesh_deprecated();
  }
 
  bool get_convergence_assertion_flag() const;
 
  void check_convergence(LinearSolver<Number> & input_solver);
 
  unsigned int get_current_training_parameter_index() const;
  void set_current_training_parameter_index(unsigned int index);
 
  const std::vector<Number> & get_evaluated_thetas(unsigned int training_parameter_index) const;
 
  /*
   * Pre-evaluate the theta functions on the entire (local) training parameter set.
   */
  virtual void preevaluate_thetas();
 
  void reset_preevaluate_thetas_completed();
 
  //----------- PROTECTED DATA MEMBERS -----------//
 
  unsigned int Nmax;
 
  unsigned int delta_N;
 
  bool output_dual_innerprods_computed;
 
  bool assert_convergence;
 
private:
 
  //----------- PRIVATE DATA MEMBERS -----------//
 
  RBEvaluation * rb_eval;
 
  RBAssemblyExpansion * rb_assembly_expansion;
 
  ElemAssembly * inner_product_assembly;
 
  bool use_energy_inner_product;
 
  std::vector<Number> energy_inner_product_coeffs;
 
  std::vector<std::unique_ptr<SparseMatrix<Number>>> Aq_vector;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> Fq_vector;
 
  std::vector<std::vector<std::unique_ptr<NumericVector<Number>>>> outputs_vector;
 
  std::vector<std::unique_ptr<SparseMatrix<Number>>> non_dirichlet_Aq_vector;
  std::vector<std::unique_ptr<NumericVector<Number>>> non_dirichlet_Fq_vector;
  std::vector<std::vector<std::unique_ptr<NumericVector<Number>>>> non_dirichlet_outputs_vector;
  std::unique_ptr<SparseMatrix<Number>> non_dirichlet_inner_product_matrix;
 
  Real rel_training_tolerance;
  Real abs_training_tolerance;
 
  bool normalize_rb_bound_in_greedy;
 
  std::string RB_training_type;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> _untransformed_basis_functions;
 
  std::unique_ptr<NumericVector<Number>> _untransformed_solution;
 
  bool _preevaluate_thetas_flag;
 
  bool _preevaluate_thetas_completed;
 
  unsigned int _current_training_parameter_index;
 
  std::vector<std::vector<Number>> _evaluated_thetas;
};
 
} // namespace libMesh
 
#endif // LIBMESH_RB_CONSTRUCTION_H