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 Namespaces 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 
85 class ParameterDefinition
86 {
87 public :
88  static const uint16_type ParameterSpaceDimension = 2;
89  typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
90 };
91 
92 class FunctionSpaceDefinition
93 {
94 public :
95  static const uint16_type Order = 1;
96 
97  static const int Order_s = CONVECTION_ORDER_U;
98  static const int Order_p = CONVECTION_ORDER_P;
99  static const int Order_t = CONVECTION_ORDER_T;
100 
101  // Definitions pour mesh
102  typedef Simplex<CONVECTION_DIM> entity_type;
103  typedef Mesh<entity_type> mesh_type;
104 
105  // space and associated elements definitions
106  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
107  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
108  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
109 
110 #if defined( FEELPP_USE_LM )
111  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
112  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
113 #else
114  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
115 #endif
116 
117  typedef FunctionSpace<mesh_type, basis_type> space_type;
118 
119  typedef bases< Lagrange<Order_p, Scalar> > single_basis_type;
120  typedef FunctionSpace<mesh_type, single_basis_type> mono_space_type;
121 
122 };
123 
124 
125 //for compilation
126 template <typename ParameterDefinition, typename FunctionSpaceDefinition >
127 class EimDefinition
128 {
129 public :
130  typedef typename ParameterDefinition::parameterspace_type parameterspace_type;
131  typedef typename FunctionSpaceDefinition::mono_space_type mono_space_type;
132  typedef typename FunctionSpaceDefinition::space_type space_type;
133 
134  typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fun_type;
135  typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fund_type;
136 };
137 
146 //template< int Order_s, int Order_p, int Order_t >
147 class ConvectionCrb : public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition , EimDefinition< ParameterDefinition, FunctionSpaceDefinition> >,
148  public boost::enable_shared_from_this< ConvectionCrb >
149 {
150 public:
151 
152  typedef ModelCrbBase<ParameterDefinition,FunctionSpaceDefinition, EimDefinition<ParameterDefinition,FunctionSpaceDefinition> > super_type;
153  typedef typename super_type::funs_type funs_type;
154  typedef typename super_type::funsd_type funsd_type;
155 
156  static const uint16_type Order = 1;
157  static const uint16_type ParameterSpaceDimension = 2;
158  static const bool is_time_dependent = false;
159 
160  //typedef Convection<Order_s, Order_p, Order_t> self_type;
161  static const int Order_s = CONVECTION_ORDER_U;
162  static const int Order_p = CONVECTION_ORDER_P;
163  static const int Order_t = CONVECTION_ORDER_T;
164  typedef ConvectionCrb self_type;
165 
166  // Definitions pour mesh
167  typedef Simplex<CONVECTION_DIM> entity_type;
168  typedef Mesh<entity_type> mesh_type;
169  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
170 
171  typedef Backend<double> backend_type;
172  typedef boost::shared_ptr<backend_type> backend_ptrtype;
173 
174  /*matrix*/
175  typedef backend_type::sparse_matrix_type sparse_matrix_type;
176  typedef backend_type::vector_type vector_type;
177 
178  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
179  typedef backend_type::vector_ptrtype vector_ptrtype;
180 
181  typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> eigen_matrix_type;
182  typedef eigen_matrix_type ematrix_type;
183  typedef boost::shared_ptr<eigen_matrix_type> eigen_matrix_ptrtype;
184 
185 
186  // space and associated elements definitions
187  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
188  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
189  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
190 
191  typedef FunctionSpace<mesh_type, basis_u_type> U_space_type;
192  typedef boost::shared_ptr<U_space_type> U_space_ptrtype;
193  typedef FunctionSpace<mesh_type, basis_t_type> T_space_type;
194  typedef boost::shared_ptr<T_space_type> T_space_ptrtype;
195 
196 
197  /* parameter space */
198  typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
199  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
200  typedef parameterspace_type::element_type parameter_type;
201  typedef parameterspace_type::element_ptrtype parameter_ptrtype;
202  typedef parameterspace_type::sampling_type sampling_type;
203  typedef parameterspace_type::sampling_ptrtype sampling_ptrtype;
204  typedef std::vector< std::vector< double > > beta_vector_type;
205 
206 #if defined( FEELPP_USE_LM )
207  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
208  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
209 #else
210  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
211 #endif
212 
214  typedef double value_type;
215 
216  typedef FunctionSpace<mesh_type, basis_type> space_type;
217 
218 
219  /*reduced basis space*/
220  typedef ReducedBasisSpace<ConvectionCrb, mesh_type, basis_type, value_type> rbfunctionspace_type;
221  typedef boost::shared_ptr< rbfunctionspace_type > rbfunctionspace_ptrtype;
222 
223  /* EIM */
224  //typedef EIMFunctionBase<U_space_type, space_type , parameterspace_type> fun_type;
225  //typedef boost::shared_ptr<fun_type> fun_ptrtype;
226  //typedef std::vector<fun_ptrtype> funs_type;
227 
228  typedef boost::shared_ptr<space_type> space_ptrtype;
229  typedef typename space_type::element_type element_type;
230  typedef typename element_type:: sub_element<0>::type element_0_type;
231  //typedef typename space_type::sub_functionspace<0>::type::element_type::element_type E0;
232  typedef typename element_type:: sub_element<1>::type element_1_type;
233  typedef typename element_type:: sub_element<2>::type element_2_type;
234 #if defined( FEELPP_USE_LM )
235  typedef typename element_type:: sub_element<3>::type element_3_type;
236 #endif
237 
238  typedef space_type functionspace_type;
239  typedef space_ptrtype functionspace_ptrtype;
240 
241  typedef boost::shared_ptr<element_type> element_ptrtype;
242 
243  typedef OperatorLinear<space_type,space_type> oplin_type;
244  typedef boost::shared_ptr<oplin_type> oplin_ptrtype;
245  typedef FsFunctionalLinear<space_type> funlin_type;
246  typedef boost::shared_ptr<funlin_type> funlin_ptrtype;
247 
248  // Definition pour les exportations
249  typedef Exporter<mesh_type> export_type;
250 
251  typedef boost::tuple<
252  std::vector< std::vector<sparse_matrix_ptrtype> >,
253  std::vector< std::vector<std::vector<vector_ptrtype> > >
254  > affine_decomposition_type;
255 
256  typedef Eigen::MatrixXd matrixN_type;
257 
258  // Constructeur
259  ConvectionCrb( );
260  ConvectionCrb( po::variables_map const& vm );
261 
262  // generate the mesh
263  Feel::gmsh_ptrtype createMesh();
264 
265  // Functions usefull for crb resolution :
266 
267  void initModel();
268 
269  std::string modelName()
270  {
271  std::ostringstream ostr;
272  ostr << "naturalconvection" ;
273  return ostr.str();
274  }
275 
276 
278  parameterspace_ptrtype parameterSpace() const
279  {
280  return M_Dmu;
281  };
282 
283  parameter_type refParameter()
284  {
285  return M_Dmu->min();
286  }
287 
288  affine_decomposition_type computeAffineDecomposition()
289  {
290  return boost::make_tuple( M_Aqm, M_Fqm );
291  }
292 
293  // \return the number of terms in affine decomposition of left hand
294  // side bilinear form
295  int Qa() const;
296  int QaTri() const;
297 
304  int Nl() const;
305 
310  int Ql( int l ) const;
311 
312  int mMaxA( int q );
313  int mMaxF( int output_index, int q );
314 
319  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
320  computeBetaQm( parameter_type const& mu, double time=0 ) ;
321 
322  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
323  computeBetaQm(element_type const& T, parameter_type const& mu, double time=0 )
324  {
325  return computeBetaQm( mu , time );
326  }
327 
328  void update( parameter_type const& mu );
329 
330  void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
331 
337  void solve( parameter_type const& mu, element_ptrtype& T );
338 
342  element_type solve( parameter_type const& mu );
343 
347  void l2solve( vector_ptrtype& u, vector_ptrtype const& f );
348 
352  void exportResults( element_type& u );
353  void exportResults( element_ptrtype& U, int t );
354  void exportResults( element_type& U, double t );
355 
360  double scalarProduct( vector_ptrtype const& X, vector_ptrtype const& Y );
361 
365  double scalarProduct( vector_type const& x, vector_type const& y );
366 
375  void run( const double * X, unsigned long N, double * Y, unsigned long P );
376 
381  value_type output( int output_index, parameter_type const& mu , element_type& unknown, bool need_to_solve=false);
382 
383  sparse_matrix_ptrtype newMatrix() const
384  {
385  return M_backend->newMatrix( Xh, Xh );
386  };
387 
388  vector_ptrtype newVector() const
389  {
390  return M_backend->newVector( Xh );
391  }
392 
393 
394  space_ptrtype functionSpace()
395  {
396  return Xh;
397  }
398 
402  rbfunctionspace_ptrtype rBFunctionSpace()
403  {
404  return RbXh;
405  }
406 
407 
408  void setMeshSize( double s )
409  {
410  meshSize = s;
411  };
412 
413 
414  po::options_description const& optionsDescription() const
415  {
416  return M_desc;
417  }
418 
425  po::variables_map const& vm() const
426  {
427  return M_vm;
428  }
429 
430 
431 
432  sparse_matrix_ptrtype innerProduct()
433  {
434  return M;
435  }
436 
437  void updateJacobianWithoutAffineDecomposition( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
438  void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
439  void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );
440  sparse_matrix_ptrtype computeTrilinearForm( const element_type& X );
441  sparse_matrix_ptrtype jacobian( const element_type& X );
442  vector_ptrtype residual( const element_type& X );
443 
444 private:
445 
446  po::options_description M_desc;
447 
448  po::variables_map M_vm;
449  backend_ptrtype M_backend;
450 
451  space_ptrtype Xh;
452 
453  rbfunctionspace_ptrtype RbXh;
454  boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
455 
456  oplin_ptrtype M_oplin;
457  funlin_ptrtype M_lf;
458 
459  sparse_matrix_ptrtype M_L;
460  sparse_matrix_ptrtype M_D;
461  vector_ptrtype F;
462 
463  sparse_matrix_ptrtype D,M;
464 
465  // Exporters
466  boost::shared_ptr<export_type> exporter;
467 
468  // Timers
469  std::map<std::string,std::pair<boost::timer,double> > timers;
470 
471  std::vector <double> Grashofs;
472  double M_current_Grashofs;
473  double M_current_Prandtl;
474 
475  // Variables usefull for crb resolution :
476 
477  double meshSize;
478 
479  element_ptrtype pT;
480 
481  std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
482  std::vector< std::vector<sparse_matrix_ptrtype> > M_Mqm;
483  std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
484 
485 
486  parameterspace_ptrtype M_Dmu;
487  beta_vector_type M_betaAqm;
488  std::vector<beta_vector_type> M_betaFqm;
489 
490  element_type M_unknown;
491  parameter_type M_mu;
492 
493  sparse_matrix_ptrtype M_A_tril;
494  funs_type M_funs;
495 
496 
497 };
498 #endif /* __ConvectionCrb_H */
double value_type
numerical type is double
Definition: convection_crb.hpp:214
rbfunctionspace_ptrtype rBFunctionSpace()
Returns the reduced basis function space.
Definition: convection_crb.hpp:402
parameterspace_ptrtype parameterSpace() const
return the parameter space
Definition: convection_crb.hpp:278
Definition: convection_crb.hpp:147
po::variables_map const & vm() const
Definition: convection_crb.hpp:425