Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
convection_crb.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Samuel Quinodoz
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2009-02-25
8 
9  Copyright (C) 2007 Samuel Quinodoz
10  Copyright (C) 2009 Université Joseph Fourier (Grenoble I)
11 
12  This library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU Lesser General Public
14  License as published by the Free Software Foundation; either
15  version 3.0 of the License, or (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  Lesser General Public License for more details.
21 
22  You should have received a copy of the GNU Lesser General Public
23  License along with this library; if not, write to the Free Software
24  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 */
26 #ifndef __ConvectionCrb_H
27 #define __ConvectionCrb_H 1
28 
35 #include <feel/options.hpp>
36 //#include <feel/feelcore/applicationxml.hpp>
37 #include <feel/feelcore/application.hpp>
38 
39 // (non)linear algebra backend
40 #include <feel/feelalg/backend.hpp>
41 
42 // quadrature rules
43 #include <feel/feelpoly/im.hpp>
44 
45 // function space
46 #include <feel/feeldiscr/functionspace.hpp>
47 
48 // linear operators
49 #include <feel/feeldiscr/operatorlinear.hpp>
50 #include <feel/feeldiscr/operatorlagrangep1.hpp>
51 // exporter
52 #include <feel/feelfilters/exporter.hpp>
53 
54 #include <feel/feelcrb/parameterspace.hpp>
55 #include <feel/feelcrb/eim.hpp>
56 
57 #include <Eigen/Core>
58 #include <Eigen/LU>
59 #include <Eigen/Dense>
60 
61 #include <feel/feelcrb/modelcrbbase.hpp>
62 #include <feel/feeldiscr/reducedbasisspace.hpp>
63 
64 // use the Feel namespace
65 using namespace Feel;
66 using namespace Feel::vf;
67 
68 #if !defined( CONVECTION_DIM )
69 #define CONVECTION_DIM 2
70 #endif
71 #if !defined( CONVECTION_ORDER_U )
72 #define CONVECTION_ORDER_U 2
73 #endif
74 #if !defined( CONVECTION_ORDER_P )
75 #define CONVECTION_ORDER_P 1
76 #endif
77 #if !defined( CONVECTION_ORDER_T )
78 #define CONVECTION_ORDER_T 2
79 #endif
80 #if !defined( CRB_SOLVER )
81 #define CRB_SOLVER 0
82 #endif
83 
84 
86 {
87 public :
88  static const uint16_type ParameterSpaceDimension = 2;
89  typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
90 };
91 
93 {
94 public :
95  static const bool is_time_dependent = false;
96  static const bool is_linear = false;
97 
98  static const uint16_type Order = 1;
99 
100  static const int Order_s = CONVECTION_ORDER_U;
101  static const int Order_p = CONVECTION_ORDER_P;
102  static const int Order_t = CONVECTION_ORDER_T;
103 
104  // Definitions pour mesh
105  typedef Simplex<CONVECTION_DIM> entity_type;
106  typedef Mesh<entity_type> mesh_type;
107 
108  // space and associated elements definitions
109  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
110  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
111  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
112 
113 #if defined( FEELPP_USE_LM )
114  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
115  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
116 #else
117  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
118 #endif
119 
120  typedef FunctionSpace<mesh_type, basis_type> space_type;
121 
122  typedef bases< Lagrange<Order_p, Scalar> > single_basis_type;
123  typedef FunctionSpace<mesh_type, single_basis_type> mono_space_type;
124 
125 };
126 
127 
128 //for compilation
129 template <typename ParameterDefinition, typename FunctionSpaceDefinition >
131 {
132 public :
133  typedef typename ParameterDefinition::parameterspace_type parameterspace_type;
134  typedef typename FunctionSpaceDefinition::mono_space_type mono_space_type;
135  typedef typename FunctionSpaceDefinition::space_type space_type;
136 
137  typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fun_type;
138  typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fund_type;
139 };
140 
149 //template< int Order_s, int Order_p, int Order_t >
150 class ConvectionCrb : public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition , EimDefinition< ParameterDefinition, FunctionSpaceDefinition> >,
151  public boost::enable_shared_from_this< ConvectionCrb >
152 {
153 public:
154 
155  typedef ModelCrbBase<ParameterDefinition,FunctionSpaceDefinition, EimDefinition<ParameterDefinition,FunctionSpaceDefinition> > super_type;
156  typedef typename super_type::funs_type funs_type;
157  typedef typename super_type::funsd_type funsd_type;
158 
159  static const uint16_type Order = 1;
160  static const uint16_type ParameterSpaceDimension = 2;
161  static const bool is_time_dependent = false;
162 
163  //typedef Convection<Order_s, Order_p, Order_t> self_type;
164  static const int Order_s = CONVECTION_ORDER_U;
165  static const int Order_p = CONVECTION_ORDER_P;
166  static const int Order_t = CONVECTION_ORDER_T;
167  typedef ConvectionCrb self_type;
168 
169  // Definitions pour mesh
170  typedef Simplex<CONVECTION_DIM> entity_type;
171  typedef Mesh<entity_type> mesh_type;
172  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
173 
174  typedef Backend<double> backend_type;
175  typedef boost::shared_ptr<backend_type> backend_ptrtype;
176 
177  /*matrix*/
178  typedef backend_type::sparse_matrix_type sparse_matrix_type;
179  typedef backend_type::vector_type vector_type;
180 
181  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
182  typedef backend_type::vector_ptrtype vector_ptrtype;
183 
184  typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> eigen_matrix_type;
185  typedef eigen_matrix_type ematrix_type;
186  typedef boost::shared_ptr<eigen_matrix_type> eigen_matrix_ptrtype;
187 
188 
189  // space and associated elements definitions
190  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
191  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
192  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
193 
194  typedef FunctionSpace<mesh_type, basis_u_type> U_space_type;
195  typedef boost::shared_ptr<U_space_type> U_space_ptrtype;
196  typedef FunctionSpace<mesh_type, basis_t_type> T_space_type;
197  typedef boost::shared_ptr<T_space_type> T_space_ptrtype;
198 
199 
200  /* parameter space */
201  typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
202  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
203  typedef parameterspace_type::element_type parameter_type;
204  typedef parameterspace_type::element_ptrtype parameter_ptrtype;
205  typedef parameterspace_type::sampling_type sampling_type;
206  typedef parameterspace_type::sampling_ptrtype sampling_ptrtype;
207  typedef std::vector< std::vector< double > > beta_vector_type;
208 
209 #if defined( FEELPP_USE_LM )
210  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
211  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
212 #else
213  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
214 #endif
215 
217  typedef double value_type;
218 
219  typedef FunctionSpace<mesh_type, basis_type> space_type;
220 
221 
222  /*reduced basis space*/
223  typedef ReducedBasisSpace<ConvectionCrb, mesh_type, basis_type, value_type> rbfunctionspace_type;
224  typedef boost::shared_ptr< rbfunctionspace_type > rbfunctionspace_ptrtype;
225 
226  /* EIM */
227  //typedef EIMFunctionBase<U_space_type, space_type , parameterspace_type> fun_type;
228  //typedef boost::shared_ptr<fun_type> fun_ptrtype;
229  //typedef std::vector<fun_ptrtype> funs_type;
230 
231  typedef boost::shared_ptr<space_type> space_ptrtype;
232  typedef typename space_type::element_type element_type;
233  typedef typename element_type:: sub_element<0>::type element_0_type;
234  //typedef typename space_type::sub_functionspace<0>::type::element_type::element_type E0;
235  typedef typename element_type:: sub_element<1>::type element_1_type;
236  typedef typename element_type:: sub_element<2>::type element_2_type;
237 #if defined( FEELPP_USE_LM )
238  typedef typename element_type:: sub_element<3>::type element_3_type;
239 #endif
240 
241  typedef space_type functionspace_type;
242  typedef space_ptrtype functionspace_ptrtype;
243 
244  typedef boost::shared_ptr<element_type> element_ptrtype;
245 
246  typedef OperatorLinear<space_type,space_type> oplin_type;
247  typedef boost::shared_ptr<oplin_type> oplin_ptrtype;
248  typedef FsFunctionalLinear<space_type> funlin_type;
249  typedef boost::shared_ptr<funlin_type> funlin_ptrtype;
250 
251  // Definition pour les exportations
252  typedef Exporter<mesh_type> export_type;
253 
254  typedef boost::tuple<
255  std::vector< std::vector<sparse_matrix_ptrtype> >,
256  std::vector< std::vector<std::vector<vector_ptrtype> > >
257  > affine_decomposition_type;
258 
259  typedef Eigen::MatrixXd matrixN_type;
260 
261  // Constructeur
262  ConvectionCrb( );
263  ConvectionCrb( po::variables_map const& vm );
264 
265  // generate the mesh
266  Feel::gmsh_ptrtype createMesh();
267 
268  // Functions usefull for crb resolution :
269 
270  void initModel();
271 
272  std::string modelName()
273  {
274  std::ostringstream ostr;
275  ostr << "naturalconvection" ;
276  return ostr.str();
277  }
278 
279 
281  parameterspace_ptrtype parameterSpace() const
282  {
283  return M_Dmu;
284  };
285 
286  parameter_type refParameter()
287  {
288  return M_Dmu->min();
289  }
290 
291  affine_decomposition_type computeAffineDecomposition()
292  {
293  return boost::make_tuple( M_Aqm, M_Fqm );
294  }
295 
296  // \return the number of terms in affine decomposition of left hand
297  // side bilinear form
298  int Qa() const;
299  int QaTri() const;
300 
307  int Nl() const;
308 
313  int Ql( int l ) const;
314 
315  int mMaxA( int q );
316  int mMaxF( int output_index, int q );
317 
322  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
323  computeBetaQm( parameter_type const& mu, double time=0 ) ;
324 
325  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
326  computeBetaQm(element_type const& T, parameter_type const& mu, double time=0 )
327  {
328  return computeBetaQm( mu , time );
329  }
330 
331  void update( parameter_type const& mu );
332 
333  void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
334 
340  void solve( parameter_type const& mu, element_ptrtype& T );
341 
345  element_type solve( parameter_type const& mu );
346 
350  void l2solve( vector_ptrtype& u, vector_ptrtype const& f );
351 
355  void exportResults( element_type& u );
356  void exportResults( element_ptrtype& U, int t );
357  void exportResults( element_type& U, double t );
358 
363  double scalarProduct( vector_ptrtype const& X, vector_ptrtype const& Y );
364 
368  double scalarProduct( vector_type const& x, vector_type const& y );
369 
378  void run( const double * X, unsigned long N, double * Y, unsigned long P );
379 
384  value_type output( int output_index, parameter_type const& mu , element_type& unknown, bool need_to_solve=false);
385 
386  sparse_matrix_ptrtype newMatrix() const
387  {
388  return M_backend->newMatrix( Xh, Xh );
389  };
390 
391  vector_ptrtype newVector() const
392  {
393  return M_backend->newVector( Xh );
394  }
395 
396 
397  space_ptrtype functionSpace()
398  {
399  return Xh;
400  }
401 
405  rbfunctionspace_ptrtype rBFunctionSpace()
406  {
407  return RbXh;
408  }
409 
410 
411  void setMeshSize( double s )
412  {
413  meshSize = s;
414  };
415 
416 
417  po::options_description const& optionsDescription() const
418  {
419  return M_desc;
420  }
421 
428  po::variables_map const& vm() const
429  {
430  return M_vm;
431  }
432 
433 
434 
435  sparse_matrix_ptrtype energyMatrix()
436  {
437  return M;
438  }
439 
440  void updateJacobianWithoutAffineDecomposition( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
441  void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
442  void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );
443  sparse_matrix_ptrtype computeTrilinearForm( const element_type& X );
444  sparse_matrix_ptrtype jacobian( const element_type& X );
445  vector_ptrtype residual( const element_type& X );
446 
447 private:
448 
449  po::options_description M_desc;
450 
451  po::variables_map M_vm;
452  backend_ptrtype M_backend;
453 
454  space_ptrtype Xh;
455 
456  rbfunctionspace_ptrtype RbXh;
457  boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
458 
459  oplin_ptrtype M_oplin;
460  funlin_ptrtype M_lf;
461 
462  sparse_matrix_ptrtype M_L;
463  sparse_matrix_ptrtype M_D;
464  vector_ptrtype F;
465 
466  sparse_matrix_ptrtype D,M;
467 
468  // Exporters
469  boost::shared_ptr<export_type> exporter;
470 
471  // Timers
472  std::map<std::string,std::pair<boost::timer,double> > timers;
473 
474  std::vector <double> Grashofs;
475  double M_current_Grashofs;
476  double M_current_Prandtl;
477 
478  // Variables usefull for crb resolution :
479 
480  double meshSize;
481 
482  element_ptrtype pT;
483 
484  std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
485  std::vector< std::vector<sparse_matrix_ptrtype> > M_Mqm;
486  std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
487 
488 
489  parameterspace_ptrtype M_Dmu;
490  beta_vector_type M_betaAqm;
491  std::vector<beta_vector_type> M_betaFqm;
492 
493  element_type M_unknown;
494  parameter_type M_mu;
495 
496  sparse_matrix_ptrtype M_A_tril;
497  funs_type M_funs;
498 
499 
500 };
501 #endif /* __ConvectionCrb_H */
Definition: convection_crb.hpp:85
double value_type
numerical type is double
Definition: convection_crb.hpp:217
Definition: convection_crb.hpp:92
Definition: convection_crb.hpp:130
rbfunctionspace_ptrtype rBFunctionSpace()
Returns the reduced basis function space.
Definition: convection_crb.hpp:405
parameterspace_ptrtype parameterSpace() const
return the parameter space
Definition: convection_crb.hpp:281
Definition: convection_crb.hpp:150
Definition: qs_laplacian.doc:2
po::variables_map const & vm() const
Definition: convection_crb.hpp:428