Line data    Source code

       1             : // rbOOmit: An implementation of the Certified Reduced Basis method.
       2             : // Copyright (C) 2009, 2010 David J. Knezevic
       3             : 
       4             : // This file is part of rbOOmit.
       5             : 
       6             : // rbOOmit is free software; you can redistribute it and/or
       7             : // modify it under the terms of the GNU Lesser General Public
       8             : // License as published by the Free Software Foundation; either
       9             : // version 2.1 of the License, or (at your option) any later version.
      10             : 
      11             : // rbOOmit is distributed in the hope that it will be useful,
      12             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      14             : // Lesser General Public License for more details.
      15             : 
      16             : // You should have received a copy of the GNU Lesser General Public
      17             : // License along with this library; if not, write to the Free Software
      18             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      19             : 
      20             : #ifndef LIBMESH_TRANSIENT_RB_CONSTRUCTION_H
      21             : #define LIBMESH_TRANSIENT_RB_CONSTRUCTION_H
      22             : 
      23             : // rbOOmit includes
      24             : #include "libmesh/rb_construction.h"
      25             : #include "libmesh/transient_rb_evaluation.h"
      26             : #include "libmesh/rb_temporal_discretization.h"
      27             : 
      28             : // libMesh includes
      29             : #include "libmesh/transient_system.h"
      30             : 
      31             : // C++ includes
      32             : 
      33             : namespace libMesh
      34             : {
      35             : 
      36             : /**
      37             :  * This class is part of the rbOOmit framework.
      38             :  *
      39             :  * TransientRBConstruction extends RBConstruction to add
      40             :  * functionality relevant in the time-dependent case.
      41             :  *
      42             :  * We can handle time controls on the RHS as h(t)*f(x,\f$ \mu \f$).
      43             :  * See Martin Grepl's thesis for more details.
      44             :  *
      45             :  * \author David J. Knezevic
      46             :  * \date 2009
      47             :  */
      48          16 : class TransientRBConstruction : public TransientSystem<RBConstruction>, public RBTemporalDiscretization
      49             : {
      50             : public:
      51             : 
      52             :   /**
      53             :    * Constructor.  Optionally initializes required
      54             :    * data structures.
      55             :    */
      56             :   TransientRBConstruction (EquationSystems & es,
      57             :                            const std::string & name,
      58             :                            const unsigned int number);
      59             : 
      60             :   /**
      61             :    * Special functions.
      62             :    * - This class has the same restrictions/defaults as its base class.
      63             :    * - Destructor is defaulted out-of-line
      64             :    */
      65             :   TransientRBConstruction (TransientRBConstruction &&) = default;
      66             :   TransientRBConstruction (const TransientRBConstruction &) = delete;
      67             :   TransientRBConstruction & operator= (const TransientRBConstruction &) = delete;
      68             :   TransientRBConstruction & operator= (TransientRBConstruction &&) = delete;
      69             :   virtual ~TransientRBConstruction ();
      70             : 
      71             :   /**
      72             :    * The type of system.
      73             :    */
      74             :   typedef TransientRBConstruction sys_type;
      75             : 
      76             :   /**
      77             :    * The type of the parent.
      78             :    */
      79             :   typedef TransientSystem<RBConstruction> Parent;
      80             : 
      81             :   /**
      82             :    * Clear all the data structures associated with
      83             :    * the system.
      84             :    */
      85             :   virtual void clear () override;
      86             : 
      87             :   /**
      88             :    * Allocate all the data structures necessary for the construction
      89             :    * stage of the RB method. This function also performs
      90             :    * matrix and vector assembly of the "truth" affine expansion.
      91             :    *
      92             :    * Override to check that theta and assembly expansions are consistently
      93             :    * sized.
      94             :    */
      95             :   virtual void initialize_rb_construction(bool skip_matrix_assembly=false,
      96             :                                           bool skip_vector_assembly=false) override;
      97             : 
      98             :   /**
      99             :    * Perform a truth solve at the current parameter.
     100             :    */
     101             :   virtual Real truth_solve(int write_interval) override;
     102             : 
     103             :   /**
     104             :    * Train the reduced basis. Overridden so that we can set the
     105             :    * flag compute_truth_projection_error to true so that the calls
     106             :    * to truth_solve during the basis construction will compute the
     107             :    * projection error. Other calls to truth_solve generally do not
     108             :    * need to perform these projection calculations.
     109             :    */
     110             :   virtual Real train_reduced_basis(const bool resize_rb_eval_data=true) override;
     111             : 
     112             :   /**
     113             :    * Read in the parameters from file and set up the system
     114             :    * accordingly.
     115             :    */
     116             :   virtual void process_parameters_file (const std::string & parameters_filename) override;
     117             : 
     118             :   /**
     119             :    * Print out info that describes the current setup of this RBConstruction.
     120             :    */
     121             :   virtual void print_info() const override;
     122             : 
     123             :   /**
     124             :    * Function that indicates when to terminate the Greedy
     125             :    * basis training.
     126             :    */
     127             :   virtual bool greedy_termination_test(Real abs_greedy_error,
     128             :                                        Real initial_greedy_error,
     129             :                                        int count) override;
     130             : 
     131             :   /**
     132             :    * Assemble and store all the affine operators.
     133             :    * Override to assemble the mass matrix operators.
     134             :    */
     135             :   virtual void assemble_all_affine_operators() override;
     136             : 
     137             :   /**
     138             :    * Override to assemble the L2 matrix as well.
     139             :    */
     140             :   virtual void assemble_misc_matrices() override;
     141             : 
     142             :   /**
     143             :    * Assemble the L2 matrix.
     144             :    */
     145             :   void assemble_L2_matrix(SparseMatrix<Number> * input_matrix,
     146             :                           bool apply_dirichlet_bc=true);
     147             : 
     148             :   /**
     149             :    * Assemble the mass matrix at the current parameter
     150             :    * and store it in input_matrix.
     151             :    */
     152             :   void assemble_mass_matrix(SparseMatrix<Number> * input_matrix);
     153             : 
     154             :   /**
     155             :    * Add the scaled mass matrix (assembled for the current parameter)
     156             :    * to input_matrix.
     157             :    */
     158             :   void add_scaled_mass_matrix(Number scalar,
     159             :                               SparseMatrix<Number> * input_matrix);
     160             : 
     161             :   /**
     162             :    * Perform a matrix-vector multiplication with the current mass matrix
     163             :    * and store the result in dest.
     164             :    */
     165             :   void mass_matrix_scaled_matvec(Number scalar,
     166             :                                  NumericVector<Number> & dest,
     167             :                                  NumericVector<Number> & arg);
     168             : 
     169             :   /**
     170             :    * Set the L2 object.
     171             :    */
     172             :   void set_L2_assembly(ElemAssembly & L2_assembly_in);
     173             : 
     174             :   /**
     175             :    * \returns A reference to the L2 assembly object
     176             :    */
     177             :   ElemAssembly & get_L2_assembly();
     178             : 
     179             :   /**
     180             :    * Assemble the q^th affine term of the mass matrix and store it in input_matrix.
     181             :    */
     182             :   void assemble_Mq_matrix(unsigned int q,
     183             :                           SparseMatrix<Number> * input_matrix,
     184             :                           bool apply_dirichlet_bc=true);
     185             : 
     186             :   /**
     187             :    * Get a pointer to M_q.
     188             :    */
     189             :   SparseMatrix<Number> * get_M_q(unsigned int q);
     190             : 
     191             :   /**
     192             :    * Get a pointer to non_dirichlet_M_q.
     193             :    */
     194             :   SparseMatrix<Number> * get_non_dirichlet_M_q(unsigned int q);
     195             : 
     196             :   /**
     197             :    * Get a map that stores pointers to all of the matrices.
     198             :    */
     199             :   virtual void get_all_matrices(std::map<std::string, SparseMatrix<Number> *> & all_matrices) override;
     200             : 
     201             :   /**
     202             :    * Assemble the truth system in the transient linear case.
     203             :    */
     204             :   virtual void truth_assembly() override;
     205             : 
     206             :   /**
     207             :    * Get/set max_truth_solves, the maximum number of RB
     208             :    * truth solves we are willing to compute in the transient
     209             :    * case.
     210             :    *
     211             :    * \note In the steady state case, \p max_truth_solves is not needed
     212             :    * since it is equivalent to \p Nmax.
     213             :    */
     214        1360 :   int get_max_truth_solves() const                   { return max_truth_solves; }
     215          70 :   void set_max_truth_solves(int max_truth_solves_in) { this->max_truth_solves = max_truth_solves_in; }
     216             : 
     217             :   /**
     218             :    * Get/set POD_tol
     219             :    */
     220          68 :   Real get_POD_tol() const                { return POD_tol; }
     221          70 :   void set_POD_tol(const Real POD_tol_in) { this->POD_tol = POD_tol_in; }
     222             : 
     223             :   /**
     224             :    * Set delta_N, the number of basis functions we add to the
     225             :    * RB space from each POD
     226             :    */
     227          70 :   void set_delta_N(const unsigned int new_delta_N) { this->delta_N = new_delta_N; }
     228             : 
     229             :   /**
     230             :    * Load the RB solution from the current time-level
     231             :    * into the libMesh solution vector.
     232             :    */
     233             :   virtual void load_rb_solution() override;
     234             : 
     235             :   /**
     236             :    * Get the column of temporal_data corresponding to the current time level.
     237             :    * This gives access to the truth projection error data. If
     238             :    * the RB basis is empty, then this corresponds to the truth
     239             :    * solution data itself.
     240             :    */
     241             :   const NumericVector<Number> & get_error_temporal_data();
     242             : 
     243             :   /**
     244             :    * Compute the L2 projection of the initial condition
     245             :    * onto the RB space for 1 <= N <= RB_size and store
     246             :    * each projection in RB_initial_condition_matrix.
     247             :    */
     248             :   void update_RB_initial_condition_all_N();
     249             : 
     250             :   /**
     251             :    * Write out all the Riesz representor data to files. Override
     252             :    * to write out transient data too.
     253             :    */
     254             :   virtual void write_riesz_representors_to_files(const std::string & riesz_representors_dir,
     255             :                                                  const bool write_binary_residual_representors) override;
     256             : 
     257             :   /**
     258             :    * Write out all the Riesz representor data to files. Override
     259             :    * to read in transient data too.
     260             :    */
     261             :   virtual void read_riesz_representors_from_files(const std::string & riesz_representors_dir,
     262             :                                                   const bool write_binary_residual_representors) override;
     263             : 
     264             : 
     265             :   //----------- PUBLIC DATA MEMBERS -----------//
     266             : 
     267             :   /**
     268             :    * The L2 matrix.
     269             :    */
     270             :   std::unique_ptr<SparseMatrix<Number>> L2_matrix;
     271             : 
     272             :   /**
     273             :    * The L2 matrix without Dirichlet conditions enforced.
     274             :    * (This is only computed if store_non_dirichlet_operators == true.)
     275             :    */
     276             :   std::unique_ptr<SparseMatrix<Number>> non_dirichlet_L2_matrix;
     277             : 
     278             :   /**
     279             :    * Vector storing the Q_m matrices from the mass operator
     280             :    */
     281             :   std::vector<std::unique_ptr<SparseMatrix<Number>>> M_q_vector;
     282             : 
     283             :   /**
     284             :    * We sometimes also need a second set of M_q matrices
     285             :    * that do not have the Dirichlet boundary conditions
     286             :    * enforced.
     287             :    */
     288             :   std::vector<std::unique_ptr<SparseMatrix<Number>>> non_dirichlet_M_q_vector;
     289             : 
     290             :   /**
     291             :    * The truth outputs for all time-levels from the
     292             :    * most recent truth_solve.
     293             :    */
     294             :   std::vector<std::vector<Number>> truth_outputs_all_k;
     295             : 
     296             :   /**
     297             :    * Boolean flag to indicate whether we are using a non-zero initialization.
     298             :    * If we are, then an initialization function must be attached to the system.
     299             :    */
     300             :   bool nonzero_initialization;
     301             : 
     302             :   /**
     303             :    * Boolean flag that indicates whether we will compute the projection error
     304             :    * for the truth solution into the RB space (at every time level).
     305             :    * This typically only needs to true during a call to train_reduced_basis.
     306             :    */
     307             :   bool compute_truth_projection_error;
     308             : 
     309             :   /**
     310             :    * The filename of the file containing the initial
     311             :    * condition projected onto the truth mesh.
     312             :    */
     313             :   std::string init_filename;
     314             : 
     315             : protected:
     316             : 
     317             :   /**
     318             :    * Helper function that actually allocates all the data
     319             :    * structures required by this class.
     320             :    */
     321             :   virtual void allocate_data_structures() override;
     322             : 
     323             :   /**
     324             :    * Override assemble_affine_expansion to also initialize
     325             :    * RB_ic_proj_rhs_all_N, if necessary.
     326             :    */
     327             :   virtual void assemble_affine_expansion(bool skip_matrix_assembly,
     328             :                                          bool skip_vector_assembly) override;
     329             : 
     330             :   /**
     331             :    * This function imposes a truth initial condition,
     332             :    * defaults to zero initial condition if the flag
     333             :    * nonzero_initialization is true.
     334             :    */
     335             :   virtual void initialize_truth();
     336             : 
     337             :   /**
     338             :    * Override to return the L2 product matrix for output
     339             :    * dual norm solves for transient state problems.
     340             :    */
     341             :   virtual SparseMatrix<Number> & get_matrix_for_output_dual_solves() override;
     342             : 
     343             :   /**
     344             :    * Initialize RB space by adding the truth initial condition
     345             :    * as the first RB basis function.
     346             :    */
     347             :   void add_IC_to_RB_space();
     348             : 
     349             :   /**
     350             :    * Add a new basis functions to the RB space. In the transient
     351             :    * case we first perform a POD of the time-dependent "truth"
     352             :    * and then add a certain number of POD modes to the reduced basis.
     353             :    */
     354             :   virtual void enrich_RB_space() override;
     355             : 
     356             :   /**
     357             :    * Update the system after enriching the RB space.
     358             :    */
     359             :   virtual void update_system() override;
     360             : 
     361             :   /**
     362             :    * Compute the reduced basis matrices for the current basis.
     363             :    */
     364             :   virtual void update_RB_system_matrices() override;
     365             : 
     366             :   /**
     367             :    * Compute the terms that are combined `online'
     368             :    * to determine the dual norm of the residual.
     369             :    */
     370             :   virtual void update_residual_terms(bool compute_inner_products) override;
     371             : 
     372             :   /**
     373             :    * Set column k (i.e. the current time level) of temporal_data to the
     374             :    * difference between the current solution and the orthogonal
     375             :    * projection of the current solution onto the current RB space.
     376             :    */
     377             :   Number set_error_temporal_data();
     378             : 
     379             :   //----------- PROTECTED DATA MEMBERS -----------//
     380             : 
     381             :   /**
     382             :    * If positive, this tolerance determines the number of POD modes we
     383             :    * add to the space on a call to enrich_RB_space(). If negative, we add
     384             :    * delta_N POD modes.
     385             :    */
     386             :   Real POD_tol;
     387             : 
     388             :   /**
     389             :    * Maximum number of truth solves in the POD-Greedy. This can be
     390             :    * different from Nmax in the transient case since we may add
     391             :    * more than one basis function per truth solve.
     392             :    * If negative, it's ignored.
     393             :    */
     394             :   int max_truth_solves;
     395             : 
     396             :   /**
     397             :    * Function pointer for assembling the L2 matrix.
     398             :    */
     399             :   ElemAssembly * L2_assembly;
     400             : 
     401             :   /**
     402             :    * The vector that stores the right-hand side for the initial
     403             :    * condition projections.
     404             :    */
     405             :   DenseVector<Number> RB_ic_proj_rhs_all_N;
     406             : 
     407             : private:
     408             : 
     409             :   //----------- PRIVATE DATA MEMBERS -----------//
     410             : 
     411             :   /**
     412             :    * Dense matrix to store the data that we use for the temporal POD.
     413             :    */
     414             :   std::vector<std::unique_ptr<NumericVector<Number>>> temporal_data;
     415             : };
     416             : 
     417             : } // namespace libMesh
     418             : 
     419             : #endif // LIBMESH_TRANSIENT_RB_CONSTRUCTION_H