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 )
85 class ParameterDefinition
88 static const uint16_type ParameterSpaceDimension = 2;
89 typedef ParameterSpace<ParameterSpaceDimension> parameterspace_type;
92 class FunctionSpaceDefinition
95 static const uint16_type Order = 1;
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;
102 typedef Simplex<CONVECTION_DIM> entity_type;
103 typedef Mesh<entity_type> mesh_type;
106 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
107 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
108 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
110 #if defined( FEELPP_USE_LM )
111 typedef Lagrange<0, Scalar> basis_l_type;
112 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
114 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
117 typedef FunctionSpace<mesh_type, basis_type> space_type;
119 typedef bases< Lagrange<Order_p, Scalar> > single_basis_type;
120 typedef FunctionSpace<mesh_type, single_basis_type> mono_space_type;
126 template <
typename ParameterDefinition,
typename FunctionSpaceDefinition >
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;
134 typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fun_type;
135 typedef EIMFunctionBase<mono_space_type, space_type , parameterspace_type> fund_type;
147 class ConvectionCrb :
public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition , EimDefinition< ParameterDefinition, FunctionSpaceDefinition> >,
148 public boost::enable_shared_from_this< ConvectionCrb >
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;
156 static const uint16_type Order = 1;
157 static const uint16_type ParameterSpaceDimension = 2;
158 static const bool is_time_dependent =
false;
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;
167 typedef Simplex<CONVECTION_DIM> entity_type;
168 typedef Mesh<entity_type> mesh_type;
169 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
171 typedef Backend<double> backend_type;
172 typedef boost::shared_ptr<backend_type> backend_ptrtype;
175 typedef backend_type::sparse_matrix_type sparse_matrix_type;
176 typedef backend_type::vector_type vector_type;
178 typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
179 typedef backend_type::vector_ptrtype vector_ptrtype;
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;
187 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
188 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
189 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
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;
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;
206 #if defined( FEELPP_USE_LM )
207 typedef Lagrange<0, Scalar> basis_l_type;
208 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
210 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
216 typedef FunctionSpace<mesh_type, basis_type> space_type;
220 typedef ReducedBasisSpace<ConvectionCrb, mesh_type, basis_type, value_type> rbfunctionspace_type;
221 typedef boost::shared_ptr< rbfunctionspace_type > rbfunctionspace_ptrtype;
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;
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;
238 typedef space_type functionspace_type;
239 typedef space_ptrtype functionspace_ptrtype;
241 typedef boost::shared_ptr<element_type> element_ptrtype;
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;
249 typedef Exporter<mesh_type> export_type;
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;
256 typedef Eigen::MatrixXd matrixN_type;
263 Feel::gmsh_ptrtype createMesh();
269 std::string modelName()
271 std::ostringstream ostr;
272 ostr <<
"naturalconvection" ;
283 parameter_type refParameter()
288 affine_decomposition_type computeAffineDecomposition()
290 return boost::make_tuple( M_Aqm, M_Fqm );
310 int Ql(
int l )
const;
313 int mMaxF(
int output_index,
int q );
319 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
320 computeBetaQm( parameter_type
const& mu,
double time=0 ) ;
322 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
323 computeBetaQm(element_type
const& T, parameter_type
const& mu,
double time=0 )
325 return computeBetaQm( mu , time );
328 void update( parameter_type
const& mu );
330 void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
337 void solve( parameter_type
const& mu, element_ptrtype& T );
342 element_type solve( parameter_type
const& mu );
347 void l2solve( vector_ptrtype& u, vector_ptrtype
const& f );
352 void exportResults( element_type& u );
353 void exportResults( element_ptrtype& U,
int t );
354 void exportResults( element_type& U,
double t );
360 double scalarProduct( vector_ptrtype
const& X, vector_ptrtype
const& Y );
365 double scalarProduct( vector_type
const& x, vector_type
const& y );
375 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
381 value_type output(
int output_index, parameter_type
const& mu , element_type& unknown,
bool need_to_solve=
false);
383 sparse_matrix_ptrtype newMatrix()
const
385 return M_backend->newMatrix( Xh, Xh );
388 vector_ptrtype newVector()
const
390 return M_backend->newVector( Xh );
394 space_ptrtype functionSpace()
408 void setMeshSize(
double s )
414 po::options_description
const& optionsDescription()
const
425 po::variables_map
const&
vm()
const
432 sparse_matrix_ptrtype innerProduct()
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 );
446 po::options_description M_desc;
448 po::variables_map M_vm;
449 backend_ptrtype M_backend;
453 rbfunctionspace_ptrtype RbXh;
454 boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
456 oplin_ptrtype M_oplin;
459 sparse_matrix_ptrtype M_L;
460 sparse_matrix_ptrtype M_D;
463 sparse_matrix_ptrtype D,M;
466 boost::shared_ptr<export_type> exporter;
469 std::map<std::string,std::pair<boost::timer,double> > timers;
471 std::vector <double> Grashofs;
472 double M_current_Grashofs;
473 double M_current_Prandtl;
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;
486 parameterspace_ptrtype M_Dmu;
487 beta_vector_type M_betaAqm;
488 std::vector<beta_vector_type> M_betaFqm;
490 element_type M_unknown;
493 sparse_matrix_ptrtype M_A_tril;
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