26 #ifndef __ConvectionCrb_H
27 #define __ConvectionCrb_H 1
35 #include <feel/options.hpp>
37 #include <feel/feelcore/application.hpp>
40 #include <feel/feelalg/backend.hpp>
43 #include <feel/feelpoly/im.hpp>
46 #include <feel/feeldiscr/functionspace.hpp>
49 #include <feel/feeldiscr/operatorlinear.hpp>
50 #include <feel/feeldiscr/operatorlagrangep1.hpp>
52 #include <feel/feelfilters/exporter.hpp>
54 #include <feel/feelcrb/parameterspace.hpp>
55 #include <feel/feelcrb/eim.hpp>
59 #include <Eigen/Dense>
61 #include <feel/feelcrb/modelcrbbase.hpp>
62 #include <feel/feeldiscr/reducedbasisspace.hpp>
66 using namespace Feel::vf;
68 #if !defined( CONVECTION_DIM )
69 #define CONVECTION_DIM 2
71 #if !defined( CONVECTION_ORDER_U )
72 #define CONVECTION_ORDER_U 2
74 #if !defined( CONVECTION_ORDER_P )
75 #define CONVECTION_ORDER_P 1
77 #if !defined( CONVECTION_ORDER_T )
78 #define CONVECTION_ORDER_T 2
80 #if !defined( CRB_SOLVER )
88 static const uint16_type ParameterSpaceDimension = 2;
89 typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
95 static const bool is_time_dependent =
false;
96 static const bool is_linear =
false;
98 static const uint16_type Order = 1;
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;
105 typedef Simplex<CONVECTION_DIM> entity_type;
106 typedef Mesh<entity_type> mesh_type;
109 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
110 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
111 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
113 #if defined( FEELPP_USE_LM )
114 typedef Lagrange<0, Scalar> basis_l_type;
115 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
117 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
120 typedef FunctionSpace<mesh_type, basis_type> space_type;
122 typedef bases< Lagrange<Order_p, Scalar> > single_basis_type;
123 typedef FunctionSpace<mesh_type, single_basis_type> mono_space_type;
129 template <
typename ParameterDefinition,
typename FunctionSpaceDefinition >
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;
137 typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fun_type;
138 typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fund_type;
150 class ConvectionCrb :
public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition , EimDefinition< ParameterDefinition, FunctionSpaceDefinition> >,
151 public boost::enable_shared_from_this< ConvectionCrb >
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;
159 static const uint16_type Order = 1;
160 static const uint16_type ParameterSpaceDimension = 2;
161 static const bool is_time_dependent =
false;
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;
170 typedef Simplex<CONVECTION_DIM> entity_type;
171 typedef Mesh<entity_type> mesh_type;
172 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
174 typedef Backend<double> backend_type;
175 typedef boost::shared_ptr<backend_type> backend_ptrtype;
178 typedef backend_type::sparse_matrix_type sparse_matrix_type;
179 typedef backend_type::vector_type vector_type;
181 typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
182 typedef backend_type::vector_ptrtype vector_ptrtype;
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;
190 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
191 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
192 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
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;
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;
209 #if defined( FEELPP_USE_LM )
210 typedef Lagrange<0, Scalar> basis_l_type;
211 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
213 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
219 typedef FunctionSpace<mesh_type, basis_type> space_type;
223 typedef ReducedBasisSpace<ConvectionCrb, mesh_type, basis_type, value_type> rbfunctionspace_type;
224 typedef boost::shared_ptr< rbfunctionspace_type > rbfunctionspace_ptrtype;
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;
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;
241 typedef space_type functionspace_type;
242 typedef space_ptrtype functionspace_ptrtype;
244 typedef boost::shared_ptr<element_type> element_ptrtype;
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;
252 typedef Exporter<mesh_type> export_type;
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;
259 typedef Eigen::MatrixXd matrixN_type;
266 Feel::gmsh_ptrtype createMesh();
272 std::string modelName()
274 std::ostringstream ostr;
275 ostr <<
"naturalconvection" ;
286 parameter_type refParameter()
291 affine_decomposition_type computeAffineDecomposition()
293 return boost::make_tuple( M_Aqm, M_Fqm );
313 int Ql(
int l )
const;
316 int mMaxF(
int output_index,
int q );
322 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
323 computeBetaQm( parameter_type
const& mu,
double time=0 ) ;
325 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
326 computeBetaQm(element_type
const& T, parameter_type
const& mu,
double time=0 )
328 return computeBetaQm( mu , time );
331 void update( parameter_type
const& mu );
333 void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
340 void solve( parameter_type
const& mu, element_ptrtype& T );
345 element_type solve( parameter_type
const& mu );
350 void l2solve( vector_ptrtype& u, vector_ptrtype
const& f );
355 void exportResults( element_type& u );
356 void exportResults( element_ptrtype& U,
int t );
357 void exportResults( element_type& U,
double t );
363 double scalarProduct( vector_ptrtype
const& X, vector_ptrtype
const& Y );
368 double scalarProduct( vector_type
const& x, vector_type
const& y );
378 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
384 value_type output(
int output_index, parameter_type
const& mu , element_type& unknown,
bool need_to_solve=
false);
386 sparse_matrix_ptrtype newMatrix()
const
388 return M_backend->newMatrix( Xh, Xh );
391 vector_ptrtype newVector()
const
393 return M_backend->newVector( Xh );
397 space_ptrtype functionSpace()
411 void setMeshSize(
double s )
417 po::options_description
const& optionsDescription()
const
428 po::variables_map
const&
vm()
const
435 sparse_matrix_ptrtype energyMatrix()
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 );
449 po::options_description M_desc;
451 po::variables_map M_vm;
452 backend_ptrtype M_backend;
456 rbfunctionspace_ptrtype RbXh;
457 boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
459 oplin_ptrtype M_oplin;
462 sparse_matrix_ptrtype M_L;
463 sparse_matrix_ptrtype M_D;
466 sparse_matrix_ptrtype D,M;
469 boost::shared_ptr<export_type> exporter;
472 std::map<std::string,std::pair<boost::timer,double> > timers;
474 std::vector <double> Grashofs;
475 double M_current_Grashofs;
476 double M_current_Prandtl;
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;
489 parameterspace_ptrtype M_Dmu;
490 beta_vector_type M_betaAqm;
491 std::vector<beta_vector_type> M_betaFqm;
493 element_type M_unknown;
496 sparse_matrix_ptrtype M_A_tril;
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
po::variables_map const & vm() const
Definition: convection_crb.hpp:428