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"
#include "libmesh/rb_evaluation.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/elem_assembly.h"
#include "libmesh/dirichlet_boundaries.h"

// C++ includes

namespace libMesh
{

class RBThetaExpansion;
class RBAssemblyExpansion;

// ------------------------------------------------------------
// RBConstruction class definition

class RBConstruction : public RBConstructionBase<LinearImplicitSystem>
{
public:

  RBConstruction (EquationSystems& es,
                  const std::string& name,
                  const unsigned int number);

  virtual ~RBConstruction ();

  typedef RBConstruction sys_type;

  void set_rb_evaluation(RBEvaluation& rb_eval_in);

  RBEvaluation& get_rb_evaluation();

  bool is_rb_eval_initialized() const;

  RBThetaExpansion& get_rb_theta_expansion();

  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 ();

  virtual std::string system_type () const;

  virtual Real truth_solve(int plot_solution);

  virtual Real train_reduced_basis(const std::string& directory_name = "offline_data",
                                   const bool resize_rb_eval_data=true);

  virtual Real compute_max_error_bound();

  const RBParameters& get_greedy_parameter(unsigned int i);

  void set_training_tolerance(Real new_training_tolerance)
  {this->training_tolerance = new_training_tolerance; }
  Real get_training_tolerance() { return training_tolerance; }

  unsigned int get_Nmax() const    { return Nmax; }
  virtual void set_Nmax(unsigned int Nmax);

  void set_quiet_mode(bool quiet_mode_in)
  { this->quiet_mode = quiet_mode_in; }

  bool is_quiet() const
  { return this->quiet_mode; }

  virtual void load_basis_function(unsigned int i);

  virtual void load_rb_solution();

  SparseMatrix<Number>* get_inner_product_matrix();

  SparseMatrix<Number>* get_non_dirichlet_inner_product_matrix();

  SparseMatrix<Number>* get_Aq(unsigned int q);

  SparseMatrix<Number>* get_non_dirichlet_Aq(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_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_constraint_matrix(SparseMatrix<Number>* input_matrix);

  void assemble_and_add_constraint_matrix(SparseMatrix<Number>* input_matrix);

  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,
                                      std::string alternative_solver_in,
                                      bool reuse_preconditioner_in,
                                      bool use_relative_bound_in_greedy_in,
                                      bool write_data_during_training_in,
                                      unsigned int training_parameters_random_seed_in,
                                      bool quiet_mode_in,
                                      unsigned int Nmax_in,
                                      Real training_tolerance_in,
                                      RBParameters mu_min_in,
                                      RBParameters mu_max_in,
                                      std::map< std::string, std::vector<Real> > discrete_parameter_values_in,
                                      std::map<std::string,bool> log_scaling);

  virtual void print_info();

  void print_basis_function_orthogonality();

  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_constraint_assembly(ElemAssembly& constraint_assembly_in);

  ElemAssembly& get_constraint_assembly();

  void zero_constrained_dofs_on_vector(NumericVector<Number>& vector);

  static UniquePtr<DirichletBoundary> build_zero_dirichlet_boundary_object();

  void set_convergence_assertion_flag(bool flag);

  //----------- PUBLIC DATA MEMBERS -----------//

  std::vector<Real> training_error_bounds;

  UniquePtr< LinearSolver<Number> > inner_product_solver;

  LinearSolver<Number> *original_linear_solver;

  bool use_inner_product_solver;

  UniquePtr< SparseMatrix<Number> > inner_product_matrix;

  UniquePtr< SparseMatrix<Number> > non_dirichlet_inner_product_matrix;

  UniquePtr< SparseMatrix<Number> > constraint_matrix;

  std::vector< Number > truth_outputs;

  std::vector< std::vector< Number > > output_dual_innerprods;

  std::vector< NumericVector<Number>* > Fq_representor;

  std::vector<Number> Fq_representor_innerprods;

  bool constrained_problem;

  bool reuse_preconditioner;

  bool use_relative_bound_in_greedy;

  bool exit_on_repeated_greedy_parameters;

  bool write_data_during_training;

  bool impose_internal_dirichlet_BCs;

  bool impose_internal_fluxes;

  bool compute_RB_inner_product;

  bool store_non_dirichlet_operators;

  bool enforce_constraints_exactly;

  bool use_empty_rb_solve_in_greedy;

  bool Fq_representor_innerprods_computed;

protected:

  virtual void allocate_data_structures();

  virtual void truth_assembly();

  virtual UniquePtr<DGFEMContext> build_context();

  virtual void assemble_matrix_for_output_dual_solves();

  virtual bool greedy_termination_test(Real training_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 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& ) {}

  bool get_convergence_assertion_flag() const;

  void check_convergence();

  //----------- PROTECTED DATA MEMBERS -----------//

  unsigned int Nmax;

  unsigned int delta_N;

  bool quiet_mode;

  bool output_dual_innerprods_computed;

  bool assert_convergence;

private:

  //----------- PRIVATE DATA MEMBERS -----------//

  RBEvaluation* rb_eval;

  RBAssemblyExpansion* rb_assembly_expansion;

  ElemAssembly* inner_product_assembly;

  ElemAssembly* constraint_assembly;

  std::vector< SparseMatrix<Number>* > Aq_vector;

  std::vector< NumericVector<Number>* > Fq_vector;

  std::vector< std::vector< NumericVector<Number>* > > outputs_vector;

  std::vector< SparseMatrix<Number>* > non_dirichlet_Aq_vector;
  std::vector< NumericVector<Number>* > non_dirichlet_Fq_vector;
  std::vector< std::vector< NumericVector<Number>* > > non_dirichlet_outputs_vector;

  Real training_tolerance;

};

} // namespace libMesh

#endif // LIBMESH_RB_CONSTRUCTION_H