30 #define __CRBModel_H 1
32 #include <boost/shared_ptr.hpp>
38 #include <feel/feeldiscr/bdf2.hpp>
42 #include <feel/feeldiscr/operatorlinearcomposite.hpp>
44 #include <feel/feeldiscr/fsfunctionallinearcomposite.hpp>
47 #include <feel/feelcore/pslogger.hpp>
52 enum class CRBModelMode
54 PFEM = 0, SCM = 1, CRB = 2, SCM_ONLINE=3, CRB_ONLINE=4
72 template<
typename ModelType>
73 class CRBModel :
public boost::enable_shared_from_this<CRBModel<ModelType> >
82 static const uint16_type ParameterSpaceDimension = ModelType::ParameterSpaceDimension;
83 static const bool is_time_dependent = ModelType::is_time_dependent;
93 typedef boost::shared_ptr<ModelType> model_ptrtype;
108 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
112 typedef typename model_type::element_ptrtype element_ptrtype;
114 typedef typename model_type::backend_type backend_type;
115 typedef boost::shared_ptr<backend_type> backend_ptrtype;
116 typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
117 typedef typename model_type::vector_ptrtype vector_ptrtype;
118 typedef typename model_type::vector_type vector_type;
120 typedef typename model_type::eigen_matrix_type eigen_matrix_type;
122 typedef typename model_type::parameterspace_type parameterspace_type;
123 typedef typename model_type::parameterspace_ptrtype parameterspace_ptrtype;
124 typedef typename model_type::parameter_type parameter_type;
125 typedef typename model_type::parameter_ptrtype parameter_ptrtype;
126 typedef typename model_type::sampling_type sampling_type;
127 typedef typename model_type::sampling_ptrtype sampling_ptrtype;
130 typedef typename std::vector< std::vector < element_ptrtype > > initial_guess_type;
132 typedef Eigen::VectorXd vectorN_type;
133 typedef std::vector< std::vector< double > > beta_vector_type;
135 typedef typename boost::tuple<sparse_matrix_ptrtype,
136 sparse_matrix_ptrtype,
137 std::vector<vector_ptrtype>
138 > offline_merge_type;
141 typedef typename boost::tuple<std::vector< std::vector<sparse_matrix_ptrtype> >,
142 std::vector< std::vector<sparse_matrix_ptrtype> >,
143 std::vector< std::vector< std::vector<vector_ptrtype> > >
144 > affine_decomposition_type;
147 typedef typename boost::tuple< beta_vector_type,
149 std::vector<beta_vector_type>
154 typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
156 typedef OperatorLinearComposite< space_type , space_type > operatorcomposite_type;
157 typedef boost::shared_ptr<operatorcomposite_type> operatorcomposite_ptrtype;
159 typedef FsFunctionalLinearComposite< space_type > functionalcomposite_type;
160 typedef boost::shared_ptr<functionalcomposite_type> functionalcomposite_ptrtype;
163 typedef boost::shared_ptr<preconditioner_type> preconditioner_ptrtype;
165 static const int nb_spaces = functionspace_type::nSpaces;
166 typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
167 typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
168 typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<4> >, fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3> >,
169 fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
170 >::type >::type >::type index_vector_type;
184 M_is_initialized( false ),
185 M_mode( CRBModelMode::PFEM ),
187 M_backend( backend_type::build( BACKEND_PETSC ) ),
193 CRBModel( po::variables_map
const&
vm, CRBModelMode mode = CRBModelMode::PFEM )
197 M_InitialGuessVector(),
200 M_is_initialized( false ),
217 M_InitialGuessVector(),
220 M_is_initialized( false ),
222 M_mode( CRBModelMode::PFEM ),
224 M_backend( backend_type::build( model->
vm ) ),
230 CRBModel( model_ptrtype & model , CRBModelMode mode )
234 M_InitialGuessVector(),
237 M_is_initialized( false ),
253 M_InitialGuessV( o.M_InitialGuessV ),
254 M_InitialGuessVector( o.M_InitialGuessVector ),
257 M_is_initialized( o.M_is_initialized ),
260 M_model( o.M_model ),
261 M_backend( o.M_backend ),
274 if ( M_is_initialized )
279 _pcfactormatsolverpackage=(MatSolverPackageType) M_backend->matSolverPackageEnumType(),
280 _worldcomm=M_backend->comm(),
281 _prefix=M_backend->prefix() ,
283 M_is_initialized=
true;
285 if( ! M_model->isInitialized() )
287 LOG( INFO ) <<
"CRBModel Model is not initialized";
288 M_model->initModel();
289 M_model->setInitialized(
true );
294 if ( M_mode != CRBModelMode::CRB_ONLINE &&
295 M_mode != CRBModelMode::SCM_ONLINE )
302 auto Xh = M_model->functionSpace();
323 M_backend = o.M_backend;
338 po::variables_map
vm()
const
349 return M_model->newMatrix();
358 return M_model->newVector();
382 return M_model->innerProductForMassMatrix();
390 return M_model->innerProductForMassMatrix();
399 sparse_matrix_ptrtype
const&
h1()
const
404 sparse_matrix_ptrtype
h1()
412 return M_model->functionSpace();
418 return M_model->Qa();
426 return Qm( mpl::bool_<model_type::is_time_dependent>() );
430 return M_model->Qm();
434 if( M_model->constructOperatorCompositeM() )
435 return functionspace_type::nSpaces;
440 int QInitialGuess()
const
442 return M_model->QInitialGuess();
447 return M_model->mMaxA( q );
453 return mMaxM(q, mpl::bool_<model_type::is_time_dependent>() );
455 int mMaxM(
int q, mpl::bool_<true> )
457 return M_model->mMaxM( q );
459 int mMaxM(
int q , mpl::bool_<false> )
464 int mMaxInitialGuess(
int q )
const
466 return M_model->mMaxInitialGuess( q );
469 int mMaxF(
int output_index,
int q )
471 return M_model->mMaxF( output_index, q );
477 return M_model->Nl();
483 return M_model->Ql( l );
489 return M_model->parameterSpace();
493 parameter_type refParameter()
495 return M_model->refParameter();
507 void setMeshSize(
double s )
509 M_model->setMeshSize( s );
519 beta_vector_type computeBetaInitialGuess( parameter_type
const& mu )
const
521 return M_model->computeBetaInitialGuess( mu );
529 return computeBetaQm( mu , mpl::bool_<model_type::is_time_dependent>(), time );
531 betaqm_type
computeBetaQm( parameter_type
const& mu , mpl::bool_<true>,
double time=0 )
533 return M_model->computeBetaQm( mu , time );
535 betaqm_type
computeBetaQm( parameter_type
const& mu , mpl::bool_<false>,
double time=0 )
537 beta_vector_type betaAqm;
538 beta_vector_type betaMqm;
539 std::vector<beta_vector_type> betaFqm;
542 std::vector<beta_vector_type> >
545 steady_beta = M_model->computeBetaQm( mu , time );
546 betaAqm = steady_beta.get<0>();
547 betaFqm = steady_beta.get<1>();
549 int nspace = functionspace_type::nSpaces;
551 if ( M_model->constructOperatorCompositeM() )
553 betaMqm.resize( nspace );
554 for(
int q=0; q<nspace; q++)
556 betaMqm[q].resize(1);
563 betaMqm[0].resize(1);
567 return boost::make_tuple( betaMqm, betaAqm, betaFqm );
573 return computeBetaQm( T , mu , mpl::bool_<model_type::is_time_dependent>(), time );
577 return M_model->computeBetaQm( T, mu, time );
581 beta_vector_type betaAqm, betaMqm ;
582 std::vector<beta_vector_type> betaFqm;
585 std::vector<beta_vector_type> >
588 steady_beta = M_model->computeBetaQm(T, mu , time );
589 betaAqm = steady_beta.get<0>();
590 betaFqm = steady_beta.get<1>();
592 int nspace = functionspace_type::nSpaces;
593 if ( M_model->constructOperatorCompositeM() )
595 betaMqm.resize( nspace );
596 for(
int q=0; q<nspace; q++)
598 betaMqm[q].resize(1);
605 betaMqm[0].resize(1);
610 return boost::make_tuple( betaMqm, betaAqm, betaFqm );
614 element_ptrtype assembleInitialGuess( parameter_type
const& mu );
619 offline_merge_type
update( parameter_type
const& mu,
double time=0 )
622 offline_merge_type offline_merge;
624 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
625 offline_merge = offlineMerge( all_beta , mu );
627 offline_merge = offlineMergeOnFly( all_beta, mu );
629 return offline_merge;
631 offline_merge_type
update( parameter_type
const& mu,
element_type const& T,
double time=0 )
635 offline_merge_type offline_merge;
637 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
638 offline_merge = offlineMerge( all_beta , mu );
640 offline_merge = offlineMergeOnFly( all_beta, mu );
642 return offline_merge;
646 element_type solveFemUsingAffineDecompositionFixedPoint( parameter_type
const& mu );
647 element_type solveFemDualUsingAffineDecompositionFixedPoint( parameter_type
const& mu );
648 element_type solveFemUsingOfflineEim( parameter_type
const& mu );
656 return M_model->initModel();
658 void setInitialized(
const bool & b)
660 return M_model->setInitialized( b );
664 return M_model->isInitialized();
674 return M_model->scalarContinuousEim();
682 return M_model->scalarDiscontinuousEim();
686 struct ComputeNormL2InCompositeCase
689 ComputeNormL2InCompositeCase(
element_type const composite_u1 ,
692 M_composite_u1 ( composite_u1 ),
693 M_composite_u2 ( composite_u2 )
696 template<
typename T >
698 operator()(
const T& t )
const
704 M_vec.conservativeResize( i+1 );
706 auto u1 = M_composite_u1.template element< T::value >();
707 auto u2 = M_composite_u2.template element< T::value >();
709 double norm = normL2(_range=
elements( mesh ),_expr=( idv(u1)-idv(u2) ) );
718 mutable vectorN_type M_vec;
726 static const bool is_composite = functionspace_type::is_composite;
727 return computeNormL2( u1, u2 , mpl::bool_< is_composite >() );
731 double norm=normL2(_range=
elements( u2.mesh() ),_expr=( idv(u1)-idv(u2) ) );
737 ComputeNormL2InCompositeCase compute_normL2_in_composite_case( u1, u2 );
738 index_vector_type index_vector;
739 fusion::for_each( index_vector, compute_normL2_in_composite_case );
740 return compute_normL2_in_composite_case.norm();
756 return M_model->operatorCompositeA();
759 std::vector< functionalcomposite_ptrtype > functionalCompositeF()
const
761 return M_model->functionalCompositeF();
764 operatorcomposite_ptrtype operatorCompositeM()
const
766 return operatorCompositeM( mpl::bool_<model_type::is_time_dependent>() );
768 operatorcomposite_ptrtype operatorCompositeM( mpl::bool_<true> )
const
770 return M_model->operatorCompositeM();
772 operatorcomposite_ptrtype operatorCompositeM( mpl::bool_<false> )
const
774 bool constructed_by_model = M_model->constructOperatorCompositeM();
775 if( constructed_by_model )
776 return M_model->operatorCompositeM();
778 return preAssembleMassMatrix();
797 boost::tie(
M_Mqm,
M_Aqm,
M_Fqm ) = M_model->computeAffineDecomposition();
799 if(
M_Aqm.size() == 0 )
801 auto compositeM = operatorCompositeM();
802 int q_max = this->
Qm();
803 M_Mqm.resize( q_max);
804 for(
int q=0; q<q_max; q++)
806 int m_max = this->mMaxM(q);
807 M_Mqm[q].resize(m_max);
808 for(
int m=0; m<m_max;m++)
810 auto operatorfree = compositeM->operatorlinear(q,m);
811 size_type pattern = operatorfree->pattern();
812 auto trial = operatorfree->domainSpace();
813 auto test=operatorfree->dualImageSpace();
814 M_Mqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
815 operatorfree->matPtr(
M_Mqm[q][m]);
821 M_Aqm.resize( q_max);
822 for(
int q=0; q<q_max; q++)
824 int m_max = this->mMaxA(q);
825 M_Aqm[q].resize(m_max);
826 for(
int m=0; m<m_max;m++)
828 auto operatorfree = compositeA->operatorlinear(q,m);
829 size_type pattern = operatorfree->pattern();
830 auto trial = operatorfree->domainSpace();
831 auto test=operatorfree->dualImageSpace();
832 M_Aqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
833 operatorfree->matPtr(
M_Aqm[q][m]);
837 auto vector_compositeF = functionalCompositeF();
838 int number_outputs = vector_compositeF.size();
839 M_Fqm.resize(number_outputs);
842 auto composite_f = vector_compositeF[
output];
845 for(
int q=0; q<q_max; q++)
847 int m_max = this->mMaxF(
output,q);
849 for(
int m=0; m<m_max;m++)
851 auto operatorfree = composite_f->functionallinear(q,m);
852 auto space = operatorfree->space();
853 M_Fqm[
output][q][m]= M_backend->newVector( space );
864 initial_guess_type initial_guess;
865 boost::tie(
M_Aqm,
M_Fqm ) = M_model->computeAffineDecomposition();
867 if (
M_Aqm.size() > 0 )
869 assembleMassMatrix();
873 auto compositeM = operatorCompositeM();
874 int q_max = this->
Qm();
875 M_Mqm.resize( q_max);
876 for(
int q=0; q<q_max; q++)
878 int m_max = this->mMaxM(q);
879 M_Mqm[q].resize(m_max);
880 for(
int m=0; m<m_max;m++)
882 auto operatorfree = compositeM->operatorlinear(q,m);
883 size_type pattern = operatorfree->pattern();
884 auto trial = operatorfree->domainSpace();
885 auto test=operatorfree->dualImageSpace();
886 M_Mqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
887 operatorfree->matPtr(
M_Mqm[q][m]);
893 M_Aqm.resize( q_max);
894 for(
int q=0; q<q_max; q++)
896 int m_max = this->mMaxA(q);
897 M_Aqm[q].resize(m_max);
898 for(
int m=0; m<m_max;m++)
900 auto operatorfree = compositeA->operatorlinear(q,m);
901 size_type pattern = operatorfree->pattern();
902 auto trial = operatorfree->domainSpace();
903 auto test=operatorfree->dualImageSpace();
904 M_Aqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
905 operatorfree->matPtr(
M_Aqm[q][m]);
909 auto vector_compositeF = functionalCompositeF();
910 int number_outputs = M_model->Nl();
911 M_Fqm.resize(number_outputs);
914 auto composite_f = vector_compositeF[
output];
917 for(
int q=0; q<q_max; q++)
919 int m_max = this->mMaxF(
output,q);
921 for(
int m=0; m<m_max;m++)
923 auto functionalfree = composite_f->functionallinear(q,m);
924 auto space = functionalfree->space();
925 M_Fqm[
output][q][m]= M_backend->newVector( space );
936 std::vector< std::vector<element_ptrtype> > computeInitialGuessAffineDecomposition( )
938 return M_model->computeInitialGuessAffineDecomposition( );
941 std::vector< std::vector<element_ptrtype> > computeInitialGuessVAffineDecomposition( )
943 initial_guess_type initial_guess_v;
944 initial_guess_v = M_model->computeInitialGuessAffineDecomposition();
945 this->assembleInitialGuessV( initial_guess_v);
946 for(
int q=0; q<M_InitialGuessVector.size(); q++)
948 for(
int m=0; m<M_InitialGuessVector[q].size(); m++)
949 *M_InitialGuessV[q][m] = *M_InitialGuessVector[q][m];
951 return M_InitialGuessV;
965 return M_B->energy( xi_j, xi_i );
978 return M_B->energy( xi_i, xi_i );
991 sparse_matrix_ptrtype
Aqm( uint16_type q, uint16_type m,
bool transpose =
false )
993 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
997 return M_Aqm[q][m]->transpose();
1004 auto opfree = composite->operatorlinear(q,m);
1006 auto trial = opfree->domainSpace();
1007 auto test = opfree->dualImageSpace();
1008 auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1009 opfree->matPtr( matrix );
1012 return matrix->transpose();
1027 const sparse_matrix_ptrtype
Mqm( uint16_type q, uint16_type m,
bool transpose =
false )
const
1029 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1033 return M_Mqm[q][m]->transpose();
1039 auto composite = operatorCompositeM();
1040 auto opfree = composite->operatorlinear(q,m);
1042 auto trial = opfree->domainSpace();
1043 auto test = opfree->dualImageSpace();
1044 auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1045 opfree->matPtr( matrix );
1048 return matrix->transpose();
1062 sparse_matrix_ptrtype
Mqm( uint16_type q, uint16_type m,
bool transpose =
false )
1064 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1068 return M_Mqm[q][m]->transpose();
1074 auto composite = operatorCompositeM();
1075 auto opfree = composite->operatorlinear(q,m);
1077 auto trial = opfree->domainSpace();
1078 auto test = opfree->dualImageSpace();
1079 auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1080 opfree->matPtr( matrix );
1083 return matrix->transpose();
1089 const vector_ptrtype InitialGuessVector( uint16_type q, uint16_type m )
const
1091 return M_InitialGuessVector[q][m];
1094 vector_ptrtype InitialGuessVector( uint16_type q, uint16_type m )
1096 return M_InitialGuessVector[q][m];
1112 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1115 return M_Aqm[q][m]->energy( xi_j, xi_i, transpose );
1120 auto opfree = composite->operatorlinear(q,m);
1122 auto trial = opfree->domainSpace();
1123 auto test = opfree->dualImageSpace();
1124 auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1125 opfree->matPtr( matrix );
1127 return matrix->energy( xi_j, xi_i, transpose );
1143 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1146 return M_Mqm[q][m]->energy( xi_j, xi_i, transpose );
1150 auto composite = operatorCompositeM();
1151 auto opfree = composite->operatorlinear(q,m);
1153 auto trial = opfree->domainSpace();
1154 auto test = opfree->dualImageSpace();
1155 auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1156 opfree->matPtr( matrix );
1158 return matrix->energy( xi_j, xi_i, transpose );
1165 beta_vector_type
const& betaInitialGuessQm( mpl::bool_<true> )
const
1167 return M_model->betaInitialGuessQm();
1176 vector_ptrtype
Fqm( uint16_type l, uint16_type q,
int m )
const
1178 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1180 return M_Fqm[l][q][m];
1184 auto vector_composite = functionalCompositeF();
1185 auto composite = vector_composite[l];
1186 auto functional = composite->functionallinear(q,m);
1188 auto vector = M_backend->newVector( space );
1194 element_ptrtype InitialGuessQm( uint16_type q,
int m )
const
1196 return M_InitialGuessV[q][m];
1215 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1221 auto vector_composite = functionalCompositeF();
1222 auto composite = vector_composite[l];
1223 auto functional = composite->functionallinear(q,m);
1225 auto vector = M_backend->newVector( space );
1245 value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_ptrtype
const& xi )
1249 if( option(_name=
"crb.stock-matrices").
template as<bool>() )
1255 auto vector_composite = functionalCompositeF();
1256 auto composite = vector_composite[l];
1257 auto functional = composite->functionallinear(q,m);
1259 auto vector = M_backend->newVector( space );
1273 return M_model->scalarProduct( X, Y );
1280 return M_model->scalarProduct( X, Y );
1289 return M_model->scalarProductForMassMatrix( X, Y );
1296 return M_model->scalarProductForMassMatrix( X, Y );
1309 return M_model->scalarProductForPod( X, Y );
1311 double scalarProductForPod( vector_type
const& X, vector_type
const& Y , mpl::bool_<false> )
1324 double scalarProductForPod( vector_ptrtype
const& X, vector_ptrtype
const& Y , mpl::bool_<true> )
1326 return M_model->scalarProductForPod( X, Y );
1328 double scalarProductForPod( vector_ptrtype
const& X, vector_ptrtype
const& Y , mpl::bool_<false> )
1340 return M_model->solve( mu );
1347 void solve( parameter_type
const& mu, element_ptrtype& u )
1349 return M_model->solve( mu, u );
1356 void solve( parameter_type
const& mu, element_ptrtype& u, vector_ptrtype
const& L,
bool transpose =
false )
1358 return M_model->solve( mu, u, L, transpose );
1365 void l2solve( vector_ptrtype& u, vector_ptrtype
const& f )
1367 return M_model->l2solve( u, f );
1383 boost::tuple<int, value_type>
1386 vector_ptrtype
const& f,
1387 bool tranpose =
false
1400 exportResults(
double time, std::vector<element_ptrtype>
const& v )
1408 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
1416 return M_model->output( output_index, mu , u , need_to_solve );
1419 int computeNumberOfSnapshots()
1421 return computeNumberOfSnapshots( mpl::bool_<model_type::is_time_dependent>() );
1423 int computeNumberOfSnapshots( mpl::bool_<true> )
1425 return M_model->computeNumberOfSnapshots();
1427 int computeNumberOfSnapshots( mpl::bool_<false> )
1432 vectorN_type computeStatistics ( Eigen::VectorXd vector , std::string name )
1434 return M_model->computeStatistics( vector , name );
1439 return timeStep( mpl::bool_<model_type::is_time_dependent>() );
1441 double timeStep( mpl::bool_<true> )
1445 if ( M_model->isSteady() ) timestep=1e30;
1447 else timestep = M_model->timeStep();
1451 double timeStep( mpl::bool_<false> )
1456 double timeInitial()
1458 return timeInitial( mpl::bool_<model_type::is_time_dependent>() );
1460 double timeInitial( mpl::bool_<true> )
1462 return M_model->timeInitial();
1464 double timeInitial( mpl::bool_<false> )
1471 return timeFinal( mpl::bool_<model_type::is_time_dependent>() );
1473 double timeFinal( mpl::bool_<true> )
1477 if ( M_model->isSteady() ) timefinal=1e30;
1479 else timefinal = M_model->timeFinal();
1483 double timeFinal( mpl::bool_<false> )
1490 return timeOrder( mpl::bool_<model_type::is_time_dependent>() );
1492 int timeOrder( mpl::bool_<true> )
1494 return M_model->timeOrder();
1496 int timeOrder( mpl::bool_<false> )
1504 return isSteady( mpl::bool_<model_type::is_time_dependent>() );
1506 bool isSteady( mpl::bool_<true> )
1508 return M_model->isSteady();
1510 bool isSteady( mpl::bool_<false> )
1516 void initializationField( element_ptrtype& initial_field,parameter_type
const& mu )
1518 return initializationField( initial_field,mu,mpl::bool_<model_type::is_time_dependent>() );
1520 void initializationField( element_ptrtype& initial_field,parameter_type
const& mu,mpl::bool_<true> )
1522 return M_model->initializationField( initial_field,mu );
1524 void initializationField( element_ptrtype& initial_field,parameter_type
const& mu,mpl::bool_<false> ) {};
1535 std::vector< std::vector<sparse_matrix_ptrtype> >
M_Aqm;
1537 mutable std::vector< std::vector<element_ptrtype> > M_InitialGuessV;
1538 mutable std::vector< std::vector<vector_ptrtype> > M_InitialGuessVector;
1541 mutable std::vector< std::vector<sparse_matrix_ptrtype> >
M_Mqm;
1544 std::vector< std::vector<std::vector<vector_ptrtype> > >
M_Fqm;
1548 bool M_is_initialized;
1551 po::variables_map M_vm;
1554 CRBModelMode M_mode;
1557 model_ptrtype M_model;
1559 backend_ptrtype M_backend;
1562 sparse_matrix_ptrtype M_B;
1563 sparse_matrix_ptrtype M_H1;
1565 beta_vector_type M_dummy_betaMqm;
1576 offline_merge_type offlineMerge( betaqm_type
const& all_beta, parameter_type
const& mu );
1577 offline_merge_type offlineMergeOnFly( betaqm_type
const& all_beta , parameter_type
const& mu );
1579 void assembleMassMatrix( );
1580 void assembleMassMatrix( mpl::bool_<true> );
1581 void assembleMassMatrix( mpl::bool_<false> );
1583 operatorcomposite_ptrtype preAssembleMassMatrix( )
const ;
1584 operatorcomposite_ptrtype preAssembleMassMatrix( mpl::bool_<true> )
const ;
1585 operatorcomposite_ptrtype preAssembleMassMatrix( mpl::bool_<false> )
const ;
1587 void assembleInitialGuessV( initial_guess_type & initial_guess );
1588 void assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<true> );
1589 void assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<false> );
1592 preconditioner_ptrtype M_preconditioner;
1597 template <
typename ModelType>
1598 struct PreAssembleMassMatrixInCompositeCase
1603 typedef typename ModelType::mesh_type mesh_type;
1606 typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1609 typedef typename ModelType::space_type space_type;
1612 typedef typename ModelType::functionspace_type functionspace_type;
1613 typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1616 typedef typename ModelType::element_type element_type;
1617 typedef typename ModelType::element_ptrtype element_ptrtype;
1620 typedef OperatorLinearComposite<space_type, space_type> operatorcomposite_type;
1621 typedef boost::shared_ptr<operatorcomposite_type> operatorcomposite_ptrtype;
1623 PreAssembleMassMatrixInCompositeCase( element_type
const u ,
1624 element_type
const v )
1626 M_composite_u ( u ),
1629 auto Xh = M_composite_u.functionSpace();
1630 op_mass = opLinearComposite( _imageSpace=Xh, _domainSpace=Xh );
1633 template<
typename T >
1635 operator()(
const T& t )
const
1638 using namespace Feel::vf;
1640 auto u = M_composite_u.template element< T::value >();
1641 auto v = M_composite_v.template element< T::value >();
1642 auto Xh = M_composite_u.functionSpace();
1643 mesh_ptrtype mesh = Xh->mesh();
1645 auto expr =
integrate( _range=
elements( mesh ) , _expr=trans( idt( u ) )*
id( v ) ) ;
1646 auto opfree = opLinearFree( _imageSpace=Xh, _domainSpace=Xh, _expr=expr );
1652 auto tuple = boost::make_tuple( q , m );
1653 op_mass->addElement( tuple , opfree );
1657 operatorcomposite_ptrtype opmass()
1662 element_type M_composite_u;
1663 element_type M_composite_v;
1664 operatorcomposite_ptrtype op_mass;
1669 template <
typename ModelType>
1670 struct AssembleMassMatrixInCompositeCase
1674 typedef typename ModelType::mesh_type mesh_type;
1677 typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1680 typedef typename ModelType::space_type space_type;
1683 typedef typename ModelType::functionspace_type functionspace_type;
1684 typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1687 typedef typename ModelType::element_type element_type;
1688 typedef typename ModelType::element_ptrtype element_ptrtype;
1691 AssembleMassMatrixInCompositeCase( element_type
const u ,
1692 element_type
const v ,
1693 boost::shared_ptr<CRBModel<ModelType> > crb_model)
1695 M_composite_u ( u ),
1696 M_composite_v ( v ),
1697 M_crb_model( crb_model )
1700 template<
typename T >
1702 operator()(
const T& t )
const
1705 using namespace Feel::vf;
1707 auto u = M_composite_u.template element< T::value >();
1708 auto v = M_composite_v.template element< T::value >();
1709 auto Xh = M_composite_u.functionSpace();
1710 mesh_ptrtype mesh = Xh->mesh();
1712 form2( _test=Xh, _trial=Xh, _matrix=M_crb_model->Mqm(0,0) ) +=
1717 element_type M_composite_u;
1718 element_type M_composite_v;
1719 mutable boost::shared_ptr<CRBModel<ModelType> > M_crb_model;
1722 template <
typename ModelType>
1723 struct AssembleInitialGuessVInCompositeCase
1727 typedef typename ModelType::mesh_type mesh_type;
1730 typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1733 typedef typename ModelType::space_type space_type;
1736 typedef typename ModelType::functionspace_type functionspace_type;
1737 typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1740 typedef typename ModelType::element_type element_type;
1741 typedef typename ModelType::element_ptrtype element_ptrtype;
1743 typedef typename std::vector< std::vector < element_ptrtype > > initial_guess_type;
1745 AssembleInitialGuessVInCompositeCase( element_type
const v ,
1746 initial_guess_type
const initial_guess ,
1747 boost::shared_ptr<CRBModel<ModelType> > crb_model)
1749 M_composite_v ( v ),
1750 M_composite_initial_guess ( initial_guess ),
1751 M_crb_model ( crb_model )
1754 template<
typename T >
1756 operator()(
const T& t )
const
1758 auto v = M_composite_v.template element< T::value >();
1759 auto Xh = M_composite_v.functionSpace();
1760 mesh_ptrtype mesh = Xh->mesh();
1761 int q_max = M_crb_model->QInitialGuess();
1762 for(
int q = 0; q < q_max; q++)
1764 int m_max = M_crb_model->mMaxInitialGuess(q);
1765 for(
int m = 0; m < m_max ; m++)
1767 auto initial_guess_qm = M_crb_model->InitialGuessVector(q,m);
1768 auto view = M_composite_initial_guess[q][m]->template element< T::value >();
1769 form1( _test=Xh, _vector=initial_guess_qm ) +=
1770 integrate ( _range=
elements( mesh ), _expr=trans( Feel::vf::idv( view ) )*Feel::vf::id( v ) );
1776 element_type M_composite_v;
1777 initial_guess_type M_composite_initial_guess;
1778 mutable boost::shared_ptr<CRBModel<ModelType> > M_crb_model;
1784 template<
typename TruthModelType>
1786 CRBModel<TruthModelType>::initB()
1790 M_B = M_model->innerProduct();
1792 LOG(INFO) <<
"[CRBModel::initB] initialize scalar product\n";
1793 M_B = M_backend->newMatrix( M_model->functionSpace(), M_model->functionSpace() );
1794 using namespace Feel::vf;
1795 typename functionspace_type::element_type u( M_model->functionSpace() );
1796 form2( M_model->functionSpace(), M_model->functionSpace(), M_B, _init=true ) =
1798 gradt( u )*trans( grad( u ) ) );
1802 auto M = M_backend->newMatrix( M_model->functionSpace(), M_model->functionSpace() );
1803 form2( M_model->functionSpace(), M_model->functionSpace(), M, _init=true ) =
1806 M_B->printMatlab(
"ipB.m" );
1807 M->printMatlab(
"ipM.m" );
1809 LOG(INFO) <<
"[CRBModel::initB] starting eigen solve\n";
1811 SolverEigen<double>::eigenmodes_type modesmin=
1817 _spectrum=SMALLEST_MAGNITUDE,
1819 _ncv=M_vm[
"solvereigen-ncv"].as<int>(),
1820 _nev=M_vm[
"solvereigen-nev"].as<int>(),
1821 _tolerance=M_vm[
"solvereigen-tol"].as<double>(),
1822 _maxit=M_vm[
"solvereigen-maxiter"].as<int>()
1826 if ( modesmin.empty() || modesmin.begin()->second.get<0>()<1e-6 )
1828 LOG(INFO) <<
"coercivity constant not computed, taking 1\n";
1833 eigmin = modesmin.begin()->second.get<0>();
1836 LOG(INFO) <<
"[CRBModel::initB] coercivity constant (tau) = " << eigmin <<
"\n";
1840 M_B->addMatrix( eigmin, M );
1849 template<
typename TruthModelType>
1850 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1851 CRBModel<TruthModelType>::preAssembleMassMatrix()
const
1853 static const bool is_composite = functionspace_type::is_composite;
1854 return preAssembleMassMatrix( mpl::bool_< is_composite >() );
1857 template<
typename TruthModelType>
1858 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1859 CRBModel<TruthModelType>::preAssembleMassMatrix( mpl::bool_<false> )
const
1862 auto Xh = M_model->functionSpace();
1863 auto mesh = Xh->mesh();
1866 auto op_mass = opLinearComposite( _domainSpace=Xh , _imageSpace=Xh );
1867 auto opfree = opLinearFree( _domainSpace=Xh , _imageSpace=Xh , _expr=expr );
1868 opfree->setName(
"mass operator (automatically created)");
1870 op_mass->addElement( boost::make_tuple(0,0) , opfree );
1874 template<
typename TruthModelType>
1875 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1876 CRBModel<TruthModelType>::preAssembleMassMatrix( mpl::bool_<true> )
const
1878 auto Xh = M_model->functionSpace();
1880 index_vector_type index_vector;
1881 PreAssembleMassMatrixInCompositeCase<TruthModelType> preassemble_mass_matrix_in_composite_case ( u , v );
1882 fusion::for_each( index_vector, preassemble_mass_matrix_in_composite_case );
1884 auto op_mass = preassemble_mass_matrix_in_composite_case.opmass();
1889 template<
typename TruthModelType>
1891 CRBModel<TruthModelType>::assembleMassMatrix(
void )
1893 static const bool is_composite = functionspace_type::is_composite;
1894 return assembleMassMatrix( mpl::bool_< is_composite >() );
1897 template<
typename TruthModelType>
1899 CRBModel<TruthModelType>::assembleMassMatrix( mpl::bool_<false> )
1901 using namespace Feel::vf;
1902 auto Xh = M_model->functionSpace();
1904 M_Mqm[0].resize( 1 );
1905 M_Mqm[0][0] = M_backend->newMatrix( _test=Xh , _trial=Xh );
1906 auto mesh = Xh->mesh();
1907 form2( _test=Xh, _trial=Xh, _matrix=M_Mqm[0][0] ) =
1909 M_Mqm[0][0]->close();
1912 template<
typename TruthModelType>
1914 CRBModel<TruthModelType>::assembleMassMatrix( mpl::bool_<true> )
1916 auto Xh = M_model->functionSpace();
1918 index_vector_type index_vector;
1919 int size = functionspace_type::nSpaces;
1922 M_Mqm[0][0]=M_backend->newMatrix( _test=Xh , _trial=Xh );
1924 AssembleMassMatrixInCompositeCase<TruthModelType> assemble_mass_matrix_in_composite_case ( u , v , this->shared_from_this());
1925 fusion::for_each( index_vector, assemble_mass_matrix_in_composite_case );
1927 M_Mqm[0][0]->close();
1932 template<
typename TruthModelType>
1934 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess )
1936 static const bool is_composite = functionspace_type::is_composite;
1937 return assembleInitialGuessV( initial_guess, mpl::bool_< is_composite >() );
1940 template<
typename TruthModelType>
1942 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<true> )
1944 auto Xh = M_model->functionSpace();
1945 auto mesh = Xh->mesh();
1947 int q_max= this->QInitialGuess();
1948 M_InitialGuessV.resize( q_max );
1949 M_InitialGuessVector.resize( q_max );
1950 for(
int q = 0; q < q_max; q++ )
1952 int m_max= this->mMaxInitialGuess(q);
1953 M_InitialGuessV[q].resize( m_max );
1954 M_InitialGuessVector[q].resize( m_max );
1955 for(
int m = 0; m < m_max; m++ )
1957 M_InitialGuessV[q][m] = Xh->elementPtr();
1958 M_InitialGuessVector[q][m] = M_model->newVector();
1962 index_vector_type index_vector;
1963 AssembleInitialGuessVInCompositeCase<TruthModelType> assemble_initial_guess_v_in_composite_case ( v , initial_guess , this->shared_from_this());
1964 fusion::for_each( index_vector, assemble_initial_guess_v_in_composite_case );
1966 for(
int q = 0; q < q_max; q++ )
1968 int m_max = this->mMaxInitialGuess(q) ;
1969 for(
int m = 0; m < m_max; m++ )
1970 M_InitialGuessVector[q][m]->close();
1975 template<
typename TruthModelType>
1977 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<false> )
1979 using namespace Feel::vf;
1980 auto Xh = M_model->functionSpace();
1981 auto mesh = Xh->mesh();
1983 int q_max= this->QInitialGuess();
1984 M_InitialGuessV.resize( q_max );
1985 M_InitialGuessVector.resize( q_max );
1986 for(
int q = 0; q < q_max; q++ )
1988 int m_max= this->mMaxInitialGuess(q);
1989 M_InitialGuessV[q].resize( m_max );
1990 M_InitialGuessVector[q].resize( m_max );
1991 for(
int m = 0; m < m_max; m++ )
1993 M_InitialGuessV[q][m] = Xh->elementPtr();
1994 M_InitialGuessVector[q][m] = M_model->newVector();
1995 form1( _test=Xh, _vector=M_InitialGuessVector[q][m]) =
1997 M_InitialGuessVector[q][m]->close();
2003 template<
typename TruthModelType>
2004 typename CRBModel<TruthModelType>::element_ptrtype
2005 CRBModel<TruthModelType>::assembleInitialGuess( parameter_type
const& mu )
2007 auto Xh = M_model->functionSpace();
2008 element_ptrtype initial_guess = Xh->elementPtr();
2009 initial_guess_type vector_initial_guess;
2010 beta_vector_type beta;
2011 vector_initial_guess = M_model->computeInitialGuessAffineDecomposition();
2012 beta = M_model->computeBetaInitialGuess( mu );
2014 int q_max = vector_initial_guess.size();
2017 int m_max=vector_initial_guess[q].size();
2018 for (
size_type m = 0; m < m_max ; ++m )
2020 element_type temp = Xh->element();
2021 temp = *vector_initial_guess[q][m];
2022 temp.scale( beta[q][m] );
2023 *initial_guess += temp;
2026 return initial_guess;
2030 template<
typename TruthModelType>
2031 typename CRBModel<TruthModelType>::offline_merge_type
2032 CRBModel<TruthModelType>::offlineMergeOnFly(betaqm_type
const& all_beta, parameter_type
const& mu )
2035 auto compositeA = this->operatorCompositeA();
2037 auto compositeM = this->operatorCompositeM();
2038 auto vector_compositeF = this->functionalCompositeF();
2041 auto beta_M = all_beta.template get<0>();
2042 auto beta_A = all_beta.template get<1>();
2044 auto beta_F = all_beta.template get<2>();
2047 compositeA->setScalars( beta_A );
2048 compositeM->setScalars( beta_M );
2051 auto A = M_model->newMatrix();
2052 auto M = M_model->newMatrix();
2053 compositeA->sumAllMatrices( A );
2056 compositeM->sumAllMatrices( M );
2058 std::vector<vector_ptrtype> F( Nl() );
2060 for(
int output=0; output<Nl(); output++)
2062 auto compositeF = vector_compositeF[output];
2063 compositeF->setScalars( beta_F[output] );
2064 F[output] = M_model->newVector();
2065 compositeF->sumAllVectors( F[output] );
2068 return boost::make_tuple( M, A, F );
2072 template<
typename TruthModelType>
2073 typename CRBModel<TruthModelType>::offline_merge_type
2074 CRBModel<TruthModelType>::offlineMerge( betaqm_type
const& all_beta , parameter_type
const& mu )
2078 sparse_matrix_ptrtype A( M_backend->newMatrix(
2079 _test=M_model->functionSpace(),
2080 _trial=M_model->functionSpace()
2083 sparse_matrix_ptrtype M( M_backend->newMatrix(
2084 _test=M_model->functionSpace(),
2085 _trial=M_model->functionSpace()
2090 auto A = this->newMatrix();
2091 auto M = this->newMatrix();
2095 std::vector<vector_ptrtype> F( Nl() );
2098 auto beta_M = all_beta.template get<0>();
2099 auto beta_A = all_beta.template get<1>();
2101 auto beta_F = all_beta.template get<2>();
2106 for(
size_type m = 0; m < mMaxA(q); ++m )
2107 A->addMatrix( beta_A[q][m], M_Aqm[q][m] );
2114 for(
size_type m = 0; m < mMaxM(q) ; ++m )
2115 M->addMatrix( beta_M[q][m] , M_Mqm[q][m] );
2121 F[l] = M_backend->newVector( M_model->functionSpace() );
2124 for (
size_type q = 0; q < Ql( l ); ++q )
2126 for (
size_type m = 0; m < mMaxF(l,q); ++m )
2127 F[l]->add( beta_F[l][q][m] , M_Fqm[l][q][m] );
2132 return boost::make_tuple( M, A, F );
2136 template<
typename TruthModelType>
2137 typename CRBModel<TruthModelType>::element_type
2138 CRBModel<TruthModelType>::solveFemUsingOfflineEim( parameter_type
const& mu )
2140 auto Xh= this->functionSpace();
2143 mybdf = bdf( _space=Xh, _vm=this->vm() , _name=
"mybdf" );
2144 sparse_matrix_ptrtype A;
2145 sparse_matrix_ptrtype M;
2146 std::vector<vector_ptrtype> F;
2147 element_ptrtype InitialGuess = Xh->elementPtr();
2148 vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2149 auto u = Xh->element();
2151 double time_initial;
2155 if ( this->isSteady() )
2163 time_initial=this->timeInitial();
2164 time_step=this->timeStep();
2165 time_final=this->timeFinal();
2166 this->initializationField( InitialGuess, mu );
2169 mybdf->setTimeInitial( time_initial );
2170 mybdf->setTimeStep( time_step );
2171 mybdf->setTimeFinal( time_final );
2174 auto vec_bdf_poly = M_backend->newVector( Xh );
2176 for( mybdf->start(*InitialGuess); !mybdf->isFinished(); mybdf->next() )
2178 bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2179 auto bdf_poly = mybdf->polyDeriv();
2180 *vec_bdf_poly = bdf_poly;
2181 boost::tie(M, A, F) = this->update( mu , mybdf->time() );
2185 A->addMatrix( bdf_coeff, M );
2186 Rhs->addVector( *vec_bdf_poly, *M );
2188 M_backend->solve( _matrix=A , _solution=u, _rhs=Rhs );
2189 mybdf->shiftRight(u);
2196 template<
typename TruthModelType>
2197 typename CRBModel<TruthModelType>::element_type
2198 CRBModel<TruthModelType>::solveFemUsingAffineDecompositionFixedPoint( parameter_type
const& mu )
2200 auto Xh= this->functionSpace();
2203 mybdf = bdf( _space=Xh, _vm=this->vm() , _name=
"mybdf" );
2204 sparse_matrix_ptrtype A;
2205 sparse_matrix_ptrtype M;
2206 std::vector<vector_ptrtype> F;
2207 element_ptrtype InitialGuess = Xh->elementPtr();
2208 auto u = Xh->element();
2209 auto uold = Xh->element();
2210 vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2212 double time_initial;
2216 if ( this->isSteady() )
2222 InitialGuess = this->assembleInitialGuess( mu ) ;
2226 time_initial=this->timeInitial();
2227 time_step=this->timeStep();
2228 time_final=this->timeFinal();
2229 this->initializationField( InitialGuess, mu );
2232 mybdf->setTimeInitial( time_initial );
2233 mybdf->setTimeStep( time_step );
2234 mybdf->setTimeFinal( time_final );
2241 auto vec_bdf_poly = M_backend->newVector( Xh );
2243 int max_fixedpoint_iterations = this->vm()[
"crb.max-fixedpoint-iterations"].template as<int>();
2244 double increment_fixedpoint_tol = this->vm()[
"crb.increment-fixedpoint-tol"].template as<double>();
2245 for( mybdf->start(*InitialGuess); !mybdf->isFinished(); mybdf->next() )
2248 bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2249 auto bdf_poly = mybdf->polyDeriv();
2250 *vec_bdf_poly = bdf_poly;
2252 boost::tie(M, A, F) = this->update( mu , u , mybdf->time() );
2257 A->addMatrix( bdf_coeff, M );
2258 Rhs->addVector( *vec_bdf_poly, *M );
2261 M_preconditioner->setMatrix( A );
2262 M_backend->solve( _matrix=A , _solution=u, _rhs=Rhs , _prec=M_preconditioner);
2263 norm = this->computeNormL2( uold , u );
2265 }
while( norm > increment_fixedpoint_tol && iter<max_fixedpoint_iterations );
2266 mybdf->shiftRight(u);
2271 template<
typename TruthModelType>
2272 typename CRBModel<TruthModelType>::element_type
2273 CRBModel<TruthModelType>::solveFemDualUsingAffineDecompositionFixedPoint( parameter_type
const& mu )
2275 int output_index = option(_name=
"crb.output-index").template as<int>();
2277 auto Xh= this->functionSpace();
2280 mybdf = bdf( _space=Xh, _vm=this->vm() , _name=
"mybdf" );
2281 sparse_matrix_ptrtype A,Adu;
2282 sparse_matrix_ptrtype M;
2283 std::vector<vector_ptrtype> F;
2284 auto udu = Xh->element();
2285 auto uold = Xh->element();
2286 vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2287 auto dual_initial_field = Xh->elementPtr();
2289 double time_initial;
2293 if ( this->isSteady() )
2302 time_initial=this->timeFinal()+this->timeStep();
2303 time_step=-this->timeStep();
2304 time_final=this->timeInitial()+this->timeStep();
2307 mybdf->setTimeInitial( time_initial );
2308 mybdf->setTimeStep( time_step );
2309 mybdf->setTimeFinal( time_final );
2315 auto vec_bdf_poly = M_backend->newVector( Xh );
2317 if ( this->isSteady() )
2321 boost::tie( M, A, F) = this->update( mu , mybdf->timeInitial() );
2322 *Rhs=*F[output_index];
2323 M_preconditioner->setMatrix( M );
2324 M_backend->solve( _matrix=M, _solution=dual_initial_field, _rhs=Rhs, _prec=M_preconditioner );
2325 udu=*dual_initial_field;
2329 int max_fixedpoint_iterations = this->vm()[
"crb.max-fixedpoint-iterations"].template as<int>();
2330 double increment_fixedpoint_tol = this->vm()[
"crb.increment-fixedpoint-tol"].template as<double>();
2331 for( mybdf->start(udu); !mybdf->isFinished(); mybdf->next() )
2334 bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2335 auto bdf_poly = mybdf->polyDeriv();
2336 *vec_bdf_poly = bdf_poly;
2338 boost::tie(M, A, F) = this->update( mu , udu , mybdf->time() );
2342 A->addMatrix( bdf_coeff, M );
2344 *vec_bdf_poly = bdf_poly;
2345 Rhs->addVector( *vec_bdf_poly, *M );
2349 *Rhs = *F[output_index];
2353 if( option(
"crb.use-symmetric-matrix").
template as<bool>() )
2356 A->transpose( Adu );
2359 M_preconditioner->setMatrix( Adu );
2360 M_backend->solve( _matrix=Adu , _solution=udu, _rhs=Rhs , _prec=M_preconditioner);
2361 norm = this->computeNormL2( uold , udu );
2363 }
while( norm > increment_fixedpoint_tol && iter<max_fixedpoint_iterations );
2364 mybdf->shiftRight(udu);
2370 template<
typename TruthModelType>
2374 M_model->run( X, N, Y, P );