Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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>
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
50 #include <feel/feeldiscr/operatorlagrangep1.hpp>
51 // exporter
52 #include <feel/feelfilters/exporter.hpp>
53 
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 
63 // use the Feel namespace
64 using namespace Feel;
65 using namespace Feel::vf;
66 
67 #if !defined( CONVECTION_DIM )
68 #define CONVECTION_DIM 2
69 #endif
70 #if !defined( CONVECTION_ORDER_U )
71 #define CONVECTION_ORDER_U 2
72 #endif
73 #if !defined( CONVECTION_ORDER_P )
74 #define CONVECTION_ORDER_P 1
75 #endif
76 #if !defined( CONVECTION_ORDER_T )
77 #define CONVECTION_ORDER_T 2
78 #endif
79 #if !defined( CRB_SOLVER )
80 #define CRB_SOLVER 0
81 #endif
82 
83 
84 class ParameterDefinition
85 {
86 public :
87  static const uint16_type ParameterSpaceDimension = 2;
88  typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
89 };
90 
91 class FunctionSpaceDefinition
92 {
93 public :
94  static const uint16_type Order = 1;
95 
96  static const int Order_s = CONVECTION_ORDER_U;
97  static const int Order_p = CONVECTION_ORDER_P;
98  static const int Order_t = CONVECTION_ORDER_T;
99 
100  // Definitions pour mesh
101  typedef Simplex<CONVECTION_DIM> entity_type;
102  typedef Mesh<entity_type> mesh_type;
103 
104  // space and associated elements definitions
105  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
106  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
107  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
108 
109 #if defined( FEELPP_USE_LM )
110  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
111  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
112 #else
113  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
114 #endif
115 
116  typedef FunctionSpace<mesh_type, basis_type> space_type;
117 };
118 
119 
128 //template< int Order_s, int Order_p, int Order_t >
129 class ConvectionCrb : public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition >
130 {
131 public:
132 
133  typedef ModelCrbBase<ParameterDefinition,FunctionSpaceDefinition> super_type;
134  typedef typename super_type::funs_type funs_type;
135  typedef typename super_type::funsd_type funsd_type;
136 
137  static const uint16_type Order = 1;
138  static const uint16_type ParameterSpaceDimension = 2;
139  static const bool is_time_dependent = false;
140 
141  //typedef Convection<Order_s, Order_p, Order_t> self_type;
142  static const int Order_s = CONVECTION_ORDER_U;
143  static const int Order_p = CONVECTION_ORDER_P;
144  static const int Order_t = CONVECTION_ORDER_T;
145  typedef ConvectionCrb self_type;
146 
147  // Definitions pour mesh
150  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
151 
153  typedef boost::shared_ptr<backend_type> backend_ptrtype;
154 
155  /*matrix*/
158 
159  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
160  typedef backend_type::vector_ptrtype vector_ptrtype;
161 
162  typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> eigen_matrix_type;
163  typedef eigen_matrix_type ematrix_type;
164  typedef boost::shared_ptr<eigen_matrix_type> eigen_matrix_ptrtype;
165 
166 
167  // space and associated elements definitions
168  typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type; // velocity space
169  typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type; // pressure space
170  typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type; // temperature space
171 
173  typedef boost::shared_ptr<U_space_type> U_space_ptrtype;
175  typedef boost::shared_ptr<T_space_type> T_space_ptrtype;
176 
177  /* parameter space */
179  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
181  typedef parameterspace_type::element_ptrtype parameter_ptrtype;
183  typedef parameterspace_type::sampling_ptrtype sampling_ptrtype;
184  typedef std::vector< std::vector< double > > beta_vector_type;
185 
186 #if defined( FEELPP_USE_LM )
187  typedef Lagrange<0, Scalar> basis_l_type; // multipliers for pressure space
188  typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
189 #else
190  typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
191 #endif
192 
194  typedef double value_type;
195 
197 
198  /* EIM */
199  //typedef EIMFunctionBase<U_space_type, space_type , parameterspace_type> fun_type;
200  //typedef boost::shared_ptr<fun_type> fun_ptrtype;
201  //typedef std::vector<fun_ptrtype> funs_type;
202 
203  typedef boost::shared_ptr<space_type> space_ptrtype;
204  typedef typename space_type::element_type element_type;
205  typedef typename element_type:: sub_element<0>::type element_0_type;
206  //typedef typename space_type::sub_functionspace<0>::type::element_type::element_type E0;
207  typedef typename element_type:: sub_element<1>::type element_1_type;
208  typedef typename element_type:: sub_element<2>::type element_2_type;
209 #if defined( FEELPP_USE_LM )
210  typedef typename element_type:: sub_element<3>::type element_3_type;
211 #endif
212 
214  typedef space_ptrtype functionspace_ptrtype;
215 
216  typedef boost::shared_ptr<element_type> element_ptrtype;
217 
219  typedef boost::shared_ptr<oplin_type> oplin_ptrtype;
220  typedef FsFunctionalLinear<space_type> funlin_type;
221  typedef boost::shared_ptr<funlin_type> funlin_ptrtype;
222 
223  // Definition pour les exportations
225 
226  typedef boost::tuple<
227  std::vector< std::vector<sparse_matrix_ptrtype> >,
228  std::vector< std::vector<std::vector<vector_ptrtype> > >
229  > affine_decomposition_type;
230 
231  typedef Eigen::MatrixXd matrixN_type;
232 
233  // Constructeur
234  ConvectionCrb( );
235  ConvectionCrb( po::variables_map const& vm );
236 
237  // generate the mesh
238  Feel::gmsh_ptrtype createMesh();
239 
240  // Functions usefull for crb resolution :
241 
242  void initModel();
243 
244  std::string modelName()
245  {
246  std::ostringstream ostr;
247  ostr << "naturalconvection" ;
248  return ostr.str();
249  }
250 
251 
253  parameterspace_ptrtype parameterSpace() const
254  {
255  return M_Dmu;
256  };
257 
258  parameter_type refParameter()
259  {
260  return M_Dmu->min();
261  }
262 
263  affine_decomposition_type computeAffineDecomposition()
264  {
265  return boost::make_tuple( M_Aqm, M_Fqm );
266  }
267 
268  // \return the number of terms in affine decomposition of left hand
269  // side bilinear form
270  int Qa() const;
271  int QaTri() const;
272 
279  int Nl() const;
280 
285  int Ql( int l ) const;
286 
287  int mMaxA( int q );
288  int mMaxF( int output_index, int q );
289 
294  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
295  computeBetaQm( parameter_type const& mu, double time=0 ) ;
296 
297  boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
298  computeBetaQm(element_type const& T, parameter_type const& mu, double time=0 )
299  {
300  return computeBetaQm( mu , time );
301  }
302 
303  void update( parameter_type const& mu );
304 
305  void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
306 
312  void solve( parameter_type const& mu, element_ptrtype& T );
313 
317  element_type solve( parameter_type const& mu );
318 
322  void l2solve( vector_ptrtype& u, vector_ptrtype const& f );
323 
327  void exportResults( element_type& u );
328  void exportResults( element_ptrtype& U, int t );
329  void exportResults( element_type& U, double t );
330 
335  double scalarProduct( vector_ptrtype const& X, vector_ptrtype const& Y );
336 
340  double scalarProduct( vector_type const& x, vector_type const& y );
341 
350  void run( const double * X, unsigned long N, double * Y, unsigned long P );
351 
356  value_type output( int output_index, parameter_type const& mu , element_type& unknown, bool need_to_solve=false);
357 
358  sparse_matrix_ptrtype newMatrix() const
359  {
360  return M_backend->newMatrix( Xh, Xh );
361  };
362 
363  vector_ptrtype newVector() const
364  {
365  return M_backend->newVector( Xh );
366  }
367 
368 
369  space_ptrtype functionSpace()
370  {
371  return Xh;
372  };
373 
374  void setMeshSize( double s )
375  {
376  meshSize = s;
377  };
378 
379 
380  po::options_description const& optionsDescription() const
381  {
382  return M_desc;
383  }
384 
391  po::variables_map const& vm() const
392  {
393  return M_vm;
394  }
395 
396 
397 
398  sparse_matrix_ptrtype innerProduct()
399  {
400  return M;
401  }
402 
403  void updateJacobianWithoutAffineDecomposition( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
404  void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J );
405  void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );
406  sparse_matrix_ptrtype computeTrilinearForm( const element_type& X );
407  sparse_matrix_ptrtype jacobian( const element_type& X );
408  vector_ptrtype residual( const element_type& X );
409 
410 private:
411 
412  po::options_description M_desc;
413 
414  po::variables_map M_vm;
415  backend_ptrtype M_backend;
416 
417  space_ptrtype Xh;
418  boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
419 
420  oplin_ptrtype M_oplin;
421  funlin_ptrtype M_lf;
422 
423  sparse_matrix_ptrtype M_L;
424  sparse_matrix_ptrtype M_D;
425  vector_ptrtype F;
426 
427  sparse_matrix_ptrtype D,M;
428 
429  // Exporters
430  boost::shared_ptr<export_type> exporter;
431 
432  // Timers
433  std::map<std::string,std::pair<boost::timer,double> > timers;
434 
435  std::vector <double> Grashofs;
436  double M_current_Grashofs;
437  double M_current_Prandtl;
438 
439  // Variables usefull for crb resolution :
440 
441  double meshSize;
442 
443  element_ptrtype pT;
444 
445  std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
446  std::vector< std::vector<sparse_matrix_ptrtype> > M_Mqm;
447  std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
448 
449 
450  parameterspace_ptrtype M_Dmu;
451  beta_vector_type M_betaAqm;
452  std::vector<beta_vector_type> M_betaFqm;
453 
454  element_type M_unknown;
455  parameter_type M_mu;
456 
457  sparse_matrix_ptrtype M_A_tril;
458  funs_type M_funs;
459 
460 
461 };
462 #endif /* __ConvectionCrb_H */

Generated on Fri Oct 25 2013 14:24:05 for Feel++ by doxygen 1.8.4