29 #ifndef _FEELPP_EIM_HPP
30 #define _FEELPP_EIM_HPP 1
38 #include <boost/ref.hpp>
39 #include <boost/next_prior.hpp>
40 #include <boost/type_traits.hpp>
41 #include <boost/tuple/tuple.hpp>
42 #if BOOST_VERSION >= 104700
43 #include <boost/math/special_functions/nonfinite_num_facets.hpp>
45 #include <boost/serialization/vector.hpp>
46 #include <boost/serialization/string.hpp>
47 #include <boost/serialization/export.hpp>
48 #include <boost/serialization/base_object.hpp>
62 class ModelCrbBaseBase {};
125 template<
typename ModelType>
132 static const uint16_type nDim = ModelType::nDim;
136 typedef ModelType model_type;
137 typedef typename ModelType::value_type value_type;
139 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_type;
140 typedef typename matrix_type::ColXpr column_type;
141 typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector_type;
142 typedef typename ModelType::node_type node_type;
144 typedef typename ModelType::functionspace_type functionspace_type;
145 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
146 typedef typename functionspace_type::element_type element_type;
147 typedef typename functionspace_type::element_ptrtype element_ptrtype;
149 typedef typename ModelType::solution_type solution_type;
151 typedef typename ModelType::parameterspace_type parameterspace_type;
152 typedef typename ModelType::parameter_type parameter_type;
153 typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
155 typedef boost::tuple<double,Eigen::Matrix<double,nDim,1> > space_residual_type;
156 typedef boost::tuple<double,parameter_type> parameter_residual_type;
169 M_is_written(
false ),
174 M_offline_done(
false ),
182 EIM( po::variables_map
const&
vm, model_type* model, sampling_ptrtype sampling,
double __tol = 1e-8 )
184 super(model->modelName(), model->name(), model->name(),
vm ),
187 M_is_written(
false ),
188 M_name( model->name() ),
189 M_trainset( sampling ),
192 M_WN( M_vm[
"eim.dimension-max"].
template as<int>() ),
193 M_offline_done(
false ),
201 if ( !
loadDB() || M_vm[
"eim.rebuild-database"].
template as<bool>() )
203 LOG(INFO) <<
"construct EIM approximation...\n";
204 M_offline_done =
false;
211 if( M_vm[
"eim.study-cvg"].
template as<bool>() )
218 M_offline_done =
false;
220 }
while( M_WN < M_vm[
"eim.dimension-max"].
template as<int>() );
229 M_is_read( __bbf.M_is_read ),
230 M_is_written( __bbf.M_is_written ),
231 M_name( __bbf.M_name ),
233 MM_max (__bbf.MM_max ),
235 M_offline_done( __bbf.M_offline_done ),
236 M_tol( __bbf.M_tol ),
240 M_index_max( __bbf.M_index_max ),
241 M_model( __bbf.M_model )
262 size_type nDOF()
const { FEELPP_ASSERT( M_model != 0 ).error(
"Invalid EIM model" );
return M_model->functionSpace()->nLocalDof(); }
267 std::vector<element_type>
const&
q()
const {
return M_q; }
275 FEELPP_ASSERT( __m >= 0 && __m < M_M )( __m )( M_M ).error(
"out of bounds access" );
280 size_type mMax()
const {
return MM_max; }
289 return M_offline_done;
301 void setTrainSet( sampling_ptrtype pset ) { M_trainset = pset; }
318 boost::archive::text_oarchive oa( ofs );
338 if ( !fs::exists( db ) )
342 fs::ifstream ifs( db );
346 boost::archive::text_iarchive ia( ifs );
350 this->setIsLoaded(
true );
376 vector_type
beta( parameter_type
const& mu )
const {
return beta( mu, this->mMax() ); }
377 vector_type beta( parameter_type
const& mu, solution_type
const& T )
const {
return beta( mu, T, this->mMax() ); }
378 vector_type
beta( parameter_type
const& mu,
size_type M )
const;
379 vector_type
beta( parameter_type
const& mu, solution_type
const& T,
size_type M )
const;
381 std::vector<double> studyConvergence( parameter_type
const & mu, solution_type & solution )
const;
383 void computationalTimeStatistics( std::string appname ) {
return M_model->computationalTimeStatistics(); }
384 element_type residual (
size_type M )
const;
386 parameter_residual_type computeBestFit( sampling_ptrtype trainset,
int __M );
388 element_type operator()( parameter_type
const& mu ,
int N)
const {
return expansion( M_q,
beta( mu ) , N); }
389 element_type operator()( parameter_type
const& mu, solution_type
const& T ,
int N )
const {
return expansion( M_q,
beta( mu, T ) , N ); }
395 friend class boost::serialization::access;
396 template<
class Archive>
397 void serialize(Archive & __ar,
const unsigned int __version )
399 LOG(INFO) <<
"serializing...\n";
401 __ar & BOOST_SERIALIZATION_NVP( M_name );
403 LOG(INFO) <<
"name saved/loaded...\n";
405 __ar & BOOST_SERIALIZATION_NVP( M_offline_done );
406 LOG(INFO) <<
"offline status...\n";
408 __ar & BOOST_SERIALIZATION_NVP( MM_max );
409 LOG(INFO) <<
"M saved/loaded\n";
413 __ar & BOOST_SERIALIZATION_NVP( M_index_max );
414 LOG(INFO) <<
"index saved/loaded\n";
417 __ar & BOOST_SERIALIZATION_NVP( M_t );
418 LOG(INFO) <<
"t saved/loaded\n";
421 if ( Archive::is_loading::value )
423 for(
int i = 0; i < MM_max; ++ i )
425 M_q.push_back( M_model->functionSpace()->element() );
427 for(
int i = 0; i < MM_max; ++ i )
428 __ar & BOOST_SERIALIZATION_NVP( M_q[i] );
430 LOG(INFO) <<
"q saved/loaded\n";
434 for(
int i = 0; i < MM_max; ++ i )
435 __ar & BOOST_SERIALIZATION_NVP( M_q[i] );
440 __ar & BOOST_SERIALIZATION_NVP( M_B );
441 LOG(INFO) <<
"B saved/loaded\n";
449 po::variables_map M_vm;
450 mutable bool M_is_read;
451 mutable bool M_is_written;
454 sampling_ptrtype M_trainset;
460 mutable bool M_offline_done;
464 std::vector<element_type> M_g;
466 std::vector<element_type> M_q;
470 std::vector<node_type> M_t;
472 std::vector<size_type> M_index_max;
518 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
521 template<
typename ModelType>
522 typename EIM<ModelType>::vector_type
526 vector_type __beta( __M );
527 for (
size_type __m = 0;__m < __M;++__m )
529 __beta[__m] = M_model->operator()( this->M_t[__m], mu );
531 this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(__beta);
535 template<
typename ModelType>
536 typename EIM<ModelType>::vector_type
540 vector_type __beta( __M );
541 for (
size_type __m = 0;__m < __M;++__m )
543 __beta[__m] = M_model->operator()( T, this->M_t[__m], mu );
545 this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(__beta);
550 template<
typename ModelType>
551 typename EIM<ModelType>::element_type
552 EIM<ModelType>::residual(
size_type __M )
const
554 LOG(INFO) <<
"compute residual for m=" << __M <<
"...\n";
555 vector_type rhs( __M );
558 for (
size_type __m = 0;__m < __M;++__m )
560 LOG(INFO) <<
"t[" << __m <<
"]=" << M_t[__m] <<
"\n";
561 rhs[__m]= M_g[__M]( M_t[__m] )(0,0,0);
563 this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(rhs);
564 LOG(INFO) <<
"solve B sol = rhs with rhs = " << rhs <<
"\n";
567 LOG(INFO) <<
"compute residual..." <<
"\n";
569 auto z = expansion( M_q, rhs, __M );
570 LOG(INFO) <<
"return residual..." <<
"\n";
571 return vf::project( _space=M_model->functionSpace(),
572 _expr=idv(M_g[__M])-idv( z ) );
575 template<
typename ModelType>
580 for (
size_type __i = 0;__i < __M-1; ++__i )
583 __Z[__M-1].add(- __s,__Z[__i]);
589 template<
typename ModelType>
590 typename EIM<ModelType>::parameter_residual_type
593 LOG(INFO) <<
"compute best fit for m=" << __M
594 <<
" and trainset of size " << trainset->size() <<
"...\n";
596 parameter_type mu = M_model->parameterSpace()->element();
601 LOG(INFO) <<
"Compute best fit M=" << __M <<
"\n";
602 BOOST_FOREACH( mu, *trainset )
604 LOG(INFO) <<
"compute best fit check mu...\n";
606 LOG_EVERY_N(INFO, 1 ) <<
" (every 10 mu) compute fit at mu="<< mu <<
"\n" ;
609 auto Z = M_model->operator()( mu );
612 for (
size_type __m = 0;__m < __M;++__m )
614 rhs[__m]= Z( M_t[__m] )(0,0,0);
616 this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(rhs);
617 auto res =
vf::project( _space=M_model->functionSpace(),
618 _expr=idv(Z)-idv( expansion( M_q, rhs, __M ) ) );
620 std::string norm_used = option(_name=
"eim.norm-used-for-residual").template as<std::string>();
621 bool check_name_norm =
false;
622 LOG( INFO ) <<
"[computeBestFit] norm used : "<<norm_used;
623 if( norm_used ==
"Linfty" )
625 check_name_norm=
true;
626 auto resmax = normLinf( _range=
elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(res) );
627 LOG_ASSERT( index < trainset->size() ) <<
"Invalid index " << index <<
" should be less than trainset size = " << trainset->size() <<
"\n";
628 maxerr( index++ ) = resmax.template get<0>();
630 if( norm_used ==
"L2" )
632 check_name_norm=
true;
633 double norm = math::sqrt(
integrate( _range=
elements( M_model->mesh() ) ,_expr=idv(res)*idv(res) ).evaluate()( 0,0 ) );
634 LOG_ASSERT( index < trainset->size() ) <<
"Invalid index " << index <<
" should be less than trainset size = " << trainset->size() <<
"\n";
635 maxerr( index++ ) = norm;
637 if( norm_used ==
"LinftyVec" )
639 check_name_norm=
true;
640 double norm = res.linftyNorm();
641 LOG_ASSERT( index < trainset->size() ) <<
"Invalid index " << index <<
" should be less than trainset size = " << trainset->size() <<
"\n";
642 maxerr( index++ ) = norm ;
644 CHECK( check_name_norm ) <<
"[EIM] The name of the norm "<<norm_used<<
" is not known\n";
647 auto err = maxerr.array().abs().maxCoeff( &index2 );
648 LOG_EVERY_N(INFO, 1 ) <<
" (every 10 mu) maxerr=" << err <<
" at index = " << index2 <<
" at mu = " << trainset->at(index2) <<
"\n";
650 LOG_ASSERT( index == trainset->size() ) <<
"Invalid index " << index <<
" should be equal to trainset size = " << trainset->size() <<
"\n";
651 auto err = maxerr.array().abs().maxCoeff( &index );
652 LOG(INFO)<<
"err=" << err <<
" reached at index " << index <<
" and mu=" << trainset->at(index) <<
"\n";
653 return boost::make_tuple( err, trainset->at(index) );
655 template<
typename ModelType>
657 EIM<ModelType>::offline( )
661 if ( this->isOfflineDone() )
663 LOG(INFO) <<
"[offline] starting offline stage ...\n";
665 LOG(INFO) <<
"[offline] create mu_1...\n";
667 auto mu = M_model->parameterSpace()->min();
670 M_trainset = M_model->parameterSpace()->sampling();
671 if ( M_trainset->empty() )
673 int sampling_size = M_vm[
"eim.sampling-size"].template as<int>();
674 std::string file_name = ( boost::format(
"eim_trainset_%1%") % sampling_size ).str();
675 std::ifstream file ( file_name );
678 M_trainset->randomize( sampling_size );
679 M_trainset->writeOnFile(file_name);
684 M_trainset->readFromFile(file_name);
689 auto res = M_model->functionSpace()->element();
693 LOG(INFO) <<
"compute finite element solution at mu_1...\n";
694 M_g.push_back( M_model->operator()( mu ) );
696 LOG(INFO) <<
"compute T^" << 0 <<
"...\n";
698 auto zmax = normLinf( _range=
elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(M_g[0]) );
700 M_t.push_back( zmax.template get<1>() );
701 LOG(INFO) <<
"norm Linf = " << zmax.template get<0>() <<
" at " << zmax.template get<1>() <<
"\n";
704 LOG(INFO) <<
"compute and insert q_0...\n";
708 q.scale( 1./ M_g[0]( M_t[0] )( 0, 0, 0 ) );
711 LOG(INFO) <<
"compute entry (0,0) of interpolation matrix...\n";
712 this->M_B.resize( 1, 1 );
713 this->M_B( 0, 0 ) = 1;
714 CHECK( math::abs( M_q[0]( M_t[0] )( 0, 0, 0 ) - 1 ) < 1e-10 )
715 <<
"q[0](t[0] != 1 " <<
"q[0] = " << M_q[0]( M_t[0] )( 0, 0, 0 )
716 <<
" t[0] = "<< M_t[0] <<
"\n";
731 LOG(INFO) <<
"start greedy algorithm...\n";
733 for( ; M_M < M_vm[
"eim.dimension-max"].template as<int>(); ++M_M )
736 LOG(INFO) <<
"M=" << M_M <<
"...\n";
738 LOG(INFO) <<
"compute best fit error...\n";
740 auto bestfit = computeBestFit( M_trainset, this->M_M-1 );
742 auto g_bestfit = M_model->operator()( bestfit.template get<1>() );
743 auto gmax = normLinf( _range=
elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(g_bestfit) );
745 DVLOG(2) <<
"best fit max error = " << bestfit.template get<0>() <<
" relative error = " << bestfit.template get<0>()/gmax.template get<0>() <<
" at mu = "
746 << bestfit.template get<1>() <<
" tolerance=" << M_vm[
"eim.error-max"].template as<double>() <<
"\n";
749 if ( (bestfit.template get<0>()/gmax.template get<0>()) < M_vm[
"eim.error-max"].template as<double>() && ! M_vm[
"eim.use-dimension-max-functions"].
template as<bool>() )
756 LOG(INFO) <<
"[offline] S(" << this->M_M-1 <<
") = " << bestfit.template get<1>() <<
"\n";
759 M_g.push_back( g_bestfit );
764 LOG(INFO) <<
"[offline] compute residual M="<< M_M <<
"..." <<
"\n";
765 res = this->residual(M_M-1);
767 LOG(INFO) <<
"[offline] compute arg sup |residual|..." <<
"\n";
768 auto resmax = normLinf( _range=
elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(res) );
769 LOG(INFO) <<
"[offline] store coordinates where max absolute value is attained " << resmax.template get<1>() <<
"..." <<
"\n";
771 M_t.push_back( resmax.template get<1>() );
773 LOG(INFO) <<
"[offline] scale new basis function by " << 1./resmax.template get<0>() <<
"..." <<
"\n";
774 res.scale( 1./res(resmax.template get<1>())(0,0,0) );
775 LOG(INFO) <<
"store new basis function..." <<
"\n";
776 M_q.push_back( res );
778 std::for_each( M_t.begin(), M_t.end(), []( node_type
const& t ) { DVLOG(2) <<
"t=" << t <<
"\n"; } );
780 M_B.conservativeResize( M_M, M_M );
781 for(
int __i = 0; __i < M_M; ++__i )
783 for(
int __j = 0; __j < M_M; ++__j )
785 this->M_B( __i, __j ) = M_q[__j]( M_t[__i] )(0,0,0);
789 DVLOG(2) <<
"[offline] Interpolation matrix: M_B = " << this->M_B <<
"\n";
791 for(
int __i = 0; __i < M_M; ++__i )
793 for(
int __j = __i+1; __j < M_M; ++__j )
795 LOG_ASSERT( math::abs( M_q[__j]( M_t[__i] )(0,0,0)) < 1e-10 )
796 <<
"q[j](t_i) when j > i should be 0, = "
797 << math::abs( M_q[__j]( M_t[__i] )(0,0,0)) <<
"\n";
801 LOG_ASSERT( math::abs( M_q[__i]( M_t[__i] )( 0, 0, 0 ) - 1 ) < 1e-10 )
802 <<
"q[ " << __i <<
"](t[" << __i <<
"] != 1 "
803 <<
"t[" << __i <<
"] = "<< M_t[__i] <<
"\n";
808 VLOG(2) <<
"================================================================================\n";
810 if ( resmax.template get<0>() < M_vm[
"eim.error-max"].template as<double>() && ! M_vm[
"eim.use-dimension-max-functions"].template as<bool>() )
817 LOG(INFO) <<
"[offline] M_max = " << M_M-1 <<
"...\n";
818 this->MM_max = this->M_M-1;
820 this->M_offline_done =
true;
821 LOG(INFO) <<
"[offline] saving DB...\n";
823 LOG(INFO) <<
"[offline] done with offline stage ...\n";
826 template<
typename ModelType>
828 EIM<ModelType>::studyConvergence( parameter_type
const & mu , solution_type & solution)
const
830 LOG(INFO) <<
" Convergence study \n";
831 int proc_number = Environment::worldComm().globalRank();
833 std::vector<double> l2ErrorVec(this->mMax(), 0.0);
836 for (
int i=0; i<mu.size(); i++ )
837 mu_str= mu_str + ( boost::format(
"_%1%" ) %mu(i) ).str() ;
839 std::string file_name =
"cvg-eim-"+M_model->name()+
"-"+mu_str+
".dat";
840 if( std::ifstream( file_name ) )
841 std::remove( file_name.c_str() );
844 if( proc_number == Environment::worldComm().masterRank() )
846 conv.open(file_name, std::ios::app);
847 conv <<
"#Nb_basis" <<
"\t" <<
"L2_error" <<
"\n";
850 bool use_expression = option(_name=
"eim.compute-error-with-truth-expression").template as<bool>();
851 int max = this->mMax();
852 for(
int N=1; N<=max; N++)
855 double exprl2norm = 0 , diffl2norm = 0 ;
859 exprl2norm =M_model->expressionL2Norm( solution , mu );
860 auto eim_approximation = this->operator()(mu , solution, N);
861 diffl2norm = M_model->diffL2Norm( solution , mu , eim_approximation );
865 exprl2norm =M_model->projExpressionL2Norm( solution , mu );
866 auto eim_approximation = this->operator()(mu , solution, N);
867 diffl2norm = M_model->projDiffL2Norm( solution , mu , eim_approximation );
870 double l2_error = diffl2norm / exprl2norm ;
872 double interpolation_error = M_model->interpolationError( solution , mu );
874 int size = mu.size();
875 LOG(INFO)<<
" mu = [ ";
876 for (
int i=0; i<size-1; i++ ) LOG(INFO)<< mu[i] <<
" , ";
877 LOG(INFO)<< mu[size-1]<<
" ]\n";
881 auto expression = M_model->operator()(mu);
883 auto eim_approximation = this->operator()(mu , N);
885 LOG(INFO) <<
"EIM name : " << M_model->name() <<
"\n";
886 double norm_l2_expression = expression.l2Norm();
887 LOG(INFO) <<
"norm_l2 expression = " << norm_l2_expression <<
"\n";
888 double norm_l2_approximation = eim_approximation.l2Norm();
889 LOG(INFO) <<
"norm_l2_approximation = " << norm_l2_approximation <<
"\n";
890 auto l2_error = math::abs( norm_l2_expression - norm_l2_approximation ) / norm_l2_expression;
891 LOG(INFO) <<
"norm l2 error = " << l2_error <<
"\n";
893 l2ErrorVec[N-1] = l2_error;
895 if( proc_number == Environment::worldComm().masterRank() )
896 conv << N <<
"\t" << l2_error <<
"\t" << interpolation_error <<
"\n";
904 template<
typename SpaceType,
typename ModelSpaceType,
typename ParameterSpaceType>
905 class EIMFunctionBase
909 typedef SpaceType functionspace_type;
910 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
911 typedef typename functionspace_type::element_type element_type;
912 typedef typename functionspace_type::element_ptrtype element_ptrtype;
913 typedef typename functionspace_type::mesh_type mesh_type;
914 typedef typename functionspace_type::mesh_ptrtype mesh_ptrtype;
915 typedef typename functionspace_type::value_type value_type;
917 typedef ModelSpaceType model_functionspace_type;
918 typedef typename model_functionspace_type::element_type solution_type;
920 static const uint16_type nDim = mesh_type::nDim;
922 typedef ParameterSpaceType parameterspace_type;
923 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
924 typedef typename parameterspace_type::element_type parameter_type;
925 typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
926 typedef Eigen::Matrix<double, nDim, 1> node_type;
928 typedef EIM<EIMFunctionBase<SpaceType,model_functionspace_type, ParameterSpaceType> > eim_type;
929 typedef typename eim_type::vector_type vector_type;
931 typedef boost::shared_ptr<eim_type> eim_ptrtype;
936 EIMFunctionBase( po::variables_map
const& vm,
937 functionspace_ptrtype fspace,
938 parameterspace_ptrtype pspace,
939 sampling_ptrtype sampling,
940 std::string
const& modelname,
941 std::string
const& name )
946 M_trainset( sampling ),
947 M_modelname( modelname ),
950 LOG(INFO)<<
"EimFunctionBase constructor\n";
952 virtual ~EIMFunctionBase()
954 std::string
const& name()
const {
return M_name; }
955 void setName( std::string
const& name ) { M_name = name; }
956 std::string modelName()
const {
return M_modelname; }
957 void setModelName( std::string
const& name ) { M_modelname = name; }
959 functionspace_ptrtype functionSpace()
const {
return M_fspace; }
960 functionspace_ptrtype functionSpace() {
return M_fspace; }
961 parameterspace_ptrtype parameterSpace()
const {
return M_pspace; }
962 parameterspace_ptrtype parameterSpace() {
return M_pspace; }
963 sampling_ptrtype trainSet() {
return M_trainset; }
964 sampling_ptrtype trainSet()
const {
return M_trainset; }
965 virtual void setTrainSet( sampling_ptrtype tset ) { M_trainset = tset; }
967 void addOnlineTime(
const double time )
969 int size = M_online_time.size();
970 M_online_time.conservativeResize( size+1 );
971 M_online_time( size ) = time;
973 Eigen::VectorXd onlineTime()
const {
return M_online_time; }
975 mesh_ptrtype mesh()
const {
return M_fspace->mesh(); }
976 mesh_ptrtype mesh() {
return M_fspace->mesh(); }
978 virtual element_type operator()( parameter_type
const& ) = 0;
979 virtual element_type operator()( solution_type
const& T, parameter_type
const& ) = 0;
980 virtual element_type interpolant( parameter_type
const& ) = 0;
981 value_type operator()( node_type
const& x, parameter_type
const& mu )
983 VLOG(2) <<
"calling EIMFunctionBase::operator()( x=" << x <<
", mu=" << mu <<
")\n";
984 element_type v = this->operator()( mu );
985 value_type res = v(x)(0,0,0);
986 VLOG(2) <<
"EIMFunctionBase::operator() v(x)=" << res <<
"\n";
993 value_type operator()( solution_type
const& T, node_type
const& x, parameter_type
const& mu )
995 VLOG(2) <<
"calling EIMFunctionBase::operator()( x=" << x <<
", mu=" << mu <<
")\n";
996 element_type v = this->operator()( T, mu );
997 value_type res = v(x)(0,0,0);
998 VLOG(2) <<
"EIMFunctionBase::operator() v(x)=" << res <<
"\n";
1002 virtual element_type
const& q(
int m )
const = 0;
1003 virtual vector_type beta( parameter_type
const& mu )
const = 0;
1004 virtual vector_type beta( parameter_type
const& mu, solution_type
const& T )
const = 0;
1007 virtual std::vector<double> studyConvergence( parameter_type
const & mu , solution_type & solution )
const = 0;
1008 virtual void computationalTimeStatistics( std::string appname ) = 0;
1009 virtual double expressionL2Norm( solution_type
const& T , parameter_type
const& mu )
const = 0;
1010 virtual double diffL2Norm( solution_type
const& T , parameter_type
const& mu , element_type
const& eim_expansion )
const = 0;
1011 virtual double projExpressionL2Norm( solution_type
const& T , parameter_type
const& mu )
const = 0;
1012 virtual double projDiffL2Norm( solution_type
const& T , parameter_type
const& mu , element_type
const& eim_expansion )
const = 0;
1013 virtual double interpolationError( solution_type
const& T , parameter_type
const& mu )
const = 0;
1015 po::variables_map M_vm;
1016 functionspace_ptrtype M_fspace;
1017 parameterspace_ptrtype M_pspace;
1018 sampling_ptrtype M_trainset;
1019 std::string M_modelname;
1021 Eigen::VectorXd M_online_time;
1026 template<
typename ModelType,
typename SpaceType,
typename ExprType>
1028 :
public EIMFunctionBase<SpaceType, typename ModelType::functionspace_type, typename ModelType::parameterspace_type>
1030 typedef EIMFunctionBase<SpaceType, typename ModelType::functionspace_type, typename ModelType::parameterspace_type> super;
1032 typedef ModelType model_type;
1034 typedef boost::shared_ptr<model_type> model_ptrtype;
1036 typedef SpaceType functionspace_type;
1037 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
1038 typedef typename super::element_type element_type;
1039 typedef typename super::element_ptrtype element_ptrtype;
1041 typedef typename ModelType::element_type solution_type;
1043 typedef typename super::parameterspace_type parameterspace_type;
1044 typedef typename super::parameter_type parameter_type;
1045 typedef typename super::sampling_ptrtype sampling_ptrtype;
1047 typedef ExprType expr_type;
1048 typedef boost::shared_ptr<expr_type> expr_ptrtype;
1050 typedef typename super::eim_type eim_type;
1051 typedef typename super::eim_ptrtype eim_ptrtype;
1053 typedef typename super::vector_type vector_type;
1055 typedef CRBModel<ModelType> crbmodel_type;
1056 typedef boost::shared_ptr<crbmodel_type> crbmodel_ptrtype;
1057 typedef CRB<crbmodel_type> crb_type;
1058 typedef boost::shared_ptr<crb_type> crb_ptrtype;
1060 EIMFunction( po::variables_map
const& vm,
1061 model_ptrtype model,
1062 functionspace_ptrtype space,
1066 sampling_ptrtype sampling,
1067 std::string
const& name )
1069 super( vm, space, model->parameterSpace(), sampling, model->modelName(), name ),
1075 M_eim( new eim_type( vm, this, sampling ) )
1078 element_type operator()( parameter_type
const& mu )
1081 #if !defined(NDEBUG)
1084 M_u = M_model->solve( mu );
1086 return vf::project( _space=this->functionSpace(), _expr=M_expr );
1088 element_type operator()( solution_type
const& T, parameter_type
const& mu )
1091 #if !defined(NDEBUG)
1097 return vf::project( _space=this->functionSpace(), _expr=M_expr );
1102 double expressionL2Norm( solution_type
const& T , parameter_type
const& mu )
const
1105 #if !defined(NDEBUG)
1109 auto mesh = this->functionSpace()->mesh();
1110 return math::sqrt(
integrate( _range=
elements( mesh ), _expr=M_expr*M_expr ).evaluate()( 0,0 ) );
1115 double diffL2Norm( solution_type
const& T , parameter_type
const& mu , element_type
const & eim_expansion )
const
1118 #if !defined(NDEBUG)
1122 auto mesh = this->functionSpace()->mesh();
1123 auto difference = M_expr - idv(eim_expansion);
1124 return math::sqrt(
integrate( _range=
elements( mesh ), _expr=difference*difference ).evaluate()( 0,0 ) );
1130 double projExpressionL2Norm( solution_type
const& T , parameter_type
const& mu )
const
1133 #if !defined(NDEBUG)
1137 auto mesh = this->functionSpace()->mesh();
1138 auto pi_g =
vf::project( _space=this->functionSpace(), _expr=M_expr );
1139 return math::sqrt(
integrate( _range=
elements( mesh ), _expr=idv(pi_g)*idv(pi_g) ).evaluate()( 0,0 ) );
1143 double projDiffL2Norm( solution_type
const& T , parameter_type
const& mu , element_type
const& eim_expansion )
const
1146 #if !defined(NDEBUG)
1150 auto mesh = this->functionSpace()->mesh();
1151 auto pi_g =
vf::project( _space=this->functionSpace(), _expr=M_expr );
1152 auto diff = pi_g - eim_expansion;
1153 return math::sqrt(
integrate( _range=
elements( mesh ), _expr=idv(diff)*idv(diff) ).evaluate()( 0,0 ) );
1156 double interpolationError( solution_type
const& T , parameter_type
const& mu )
const
1159 #if !defined(NDEBUG)
1163 auto mesh = this->functionSpace()->mesh();
1164 auto pi_g =
vf::project( _space=this->functionSpace(), _expr=M_expr );
1165 auto difference = M_expr - idv(pi_g);
1166 return math::sqrt(
integrate( _range=
elements( mesh ), _expr=difference*difference ).evaluate()( 0,0 ) );
1169 void computationalTimeStatistics( std::string appname )
1171 computationalTimeStatistics( appname,
typename boost::is_base_of<ModelCrbBaseBase,model_type>::type() );
1173 void computationalTimeStatistics( std::string appname, boost::mpl::bool_<false> )
1175 void computationalTimeStatistics( std::string appname, boost::mpl::bool_<true> )
1178 auto crbmodel = crbmodel_ptrtype(
new crbmodel_type( M_model , CRBModelMode::CRB ) );
1180 M_crb = crb_ptrtype(
new crb_type( appname,
1184 if ( !M_crb->isDBLoaded() || M_crb->rebuildDB() )
1186 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1187 LOG( INFO ) <<
"No CRB DB available, do crb offline computations...";
1188 M_crb->setOfflineStep(
true );
1192 int n_eval = option(_name=
"eim.computational-time-neval").template as<int>();
1194 Eigen::Matrix<double, Eigen::Dynamic, 1> time_crb;
1195 time_crb.resize( n_eval );
1197 typename crb_type::sampling_ptrtype Sampling(
new typename crb_type::sampling_type( M_model->parameterSpace() ) );
1198 Sampling->logEquidistribute( n_eval );
1201 int N = option(_name=
"crb.dimension").template as<int>();
1203 auto WN = M_crb->wn();
1205 BOOST_FOREACH(
auto mu, *Sampling )
1208 boost::mpi::timer tcrb;
1209 auto o = M_crb->run( mu, option(_name=
"crb.online-tolerance").
template as<double>() , N);
1210 auto solutions=o.template get<2>();
1211 auto uN = solutions.template get<0>();
1214 auto u_crb = M_crb->expansion( uN[size-1] , N , WN );
1216 boost::mpi::timer teim;
1217 this->beta( mu , u_crb );
1218 this->addOnlineTime( teim.elapsed() );
1219 time_crb( mu_number ) = tcrb.elapsed() ;
1223 M_model->computeStatistics( super::onlineTime() , super::name() );
1224 M_model->computeStatistics( time_crb , super::name()+
" - global crb timing" );
1228 void setTrainSet( sampling_ptrtype tset ) { M_eim->setTrainSet( tset ); }
1229 element_type interpolant( parameter_type
const& mu ) {
return M_eim->operator()( mu , M_eim->mMax() ); }
1231 element_type
const& q(
int m )
const {
return M_eim->q( m ); }
1233 vector_type beta( parameter_type
const& mu )
const {
return M_eim->beta( mu ); }
1234 vector_type beta( parameter_type
const& mu, solution_type
const& T )
const {
return M_eim->beta( mu, T ); }
1236 std::vector<double> studyConvergence( parameter_type
const & mu , solution_type & solution )
const {
return M_eim->studyConvergence( mu , solution ) ; };
1238 size_type mMax()
const {
return M_eim->mMax(); }
1241 po::variables_map M_vm;
1242 model_ptrtype M_model;
1245 parameter_type& M_mu;
1253 template<
typename Args>
1254 struct compute_eim_return
1256 typedef typename boost::remove_reference<typename boost::remove_pointer<typename parameter::binding<Args, tag::model>::type>::type>::type::element_type model1_type;
1257 typedef typename boost::remove_const<typename boost::remove_pointer<model1_type>::type>::type model_type;
1258 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::expr>::type>::type expr_type;
1259 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::space>::type>::type::element_type space_type;
1260 typedef EIMFunction<model_type, space_type, expr_type> type;
1261 typedef boost::shared_ptr<EIMFunction<model_type, space_type, expr_type> > ptrtype;
1265 BOOST_PARAMETER_FUNCTION(
1266 (
typename Feel::detail::compute_eim_return<Args>::ptrtype ),
1270 ( in_out(model), * )
1271 ( in_out(element), * )
1272 ( in_out(parameter), * )
1278 ( options, *, Environment::vm())
1281 ( sampling, *, model->parameterSpace()->sampling() )
1282 ( verbose, (
int), 0 )
1286 Feel::detail::ignore_unused_variable_warning( args );
1287 typedef typename Feel::detail::compute_eim_return<Args>::type eim_type;
1288 typedef typename Feel::detail::compute_eim_return<Args>::ptrtype eim_ptrtype;
1289 return eim_ptrtype(
new eim_type( options, model, space, element, parameter, expr, sampling, name ) );
1292 template<
typename ModelType>
1293 struct EimFunctionNoSolve
1295 typedef typename ModelType::functionspace_type functionspace_type;
1296 typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1297 typedef typename functionspace_type::element_type element_type;
1298 typedef typename ModelType::parameterspace_type parameterspace_type;
1299 typedef typename ModelType::parameterspace_ptrtype parameterspace_ptrtype;
1300 typedef typename parameterspace_type::element_type parameter_type;
1301 typedef typename parameterspace_type::sampling_type sampling_type;
1302 typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
1303 typedef typename functionspace_type::value_type value_type;
1305 typedef boost::shared_ptr<ModelType> model_ptrtype;
1307 static const int nb_spaces = functionspace_type::nSpaces;
1308 typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
1309 typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
1310 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> >,
1311 fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
1312 >::type >::type >::type index_vector_type;
1314 EimFunctionNoSolve( model_ptrtype model )
1317 M_elt( M_model->functionSpace()->element() )
1320 element_type solve( parameter_type
const& mu )
1322 DVLOG(2) <<
"no solve required\n";
1323 static const bool is_composite = functionspace_type::is_composite;
1324 return solve( mu , mpl::bool_< is_composite >() );
1326 element_type solve( parameter_type
const& mu , mpl::bool_<false> )
1328 value_type x = boost::lexical_cast<value_type>(
"inf");
1329 M_elt =
vf::project( _space=M_model->functionSpace(), _expr=cst(x) );
1332 element_type solve( parameter_type
const& mu , mpl::bool_<true> )
1334 ProjectInfCompositeCase project_inf_composite_case( M_elt );
1335 index_vector_type index_vector;
1336 fusion::for_each( index_vector, project_inf_composite_case );
1337 return project_inf_composite_case.element();
1340 std::string modelName()
const {
return M_model->modelName(); }
1341 functionspace_ptrtype functionSpace() {
return M_model->functionSpace(); }
1342 parameterspace_ptrtype parameterSpace() {
return M_model->parameterSpace(); }
1344 struct ProjectInfCompositeCase
1346 ProjectInfCompositeCase( element_type & composite_element)
1348 M_element( composite_element )
1351 template<
typename T >
1353 operator()(
const T& t )
const
1355 auto view = M_element.template element< T::value >();
1356 auto space = view.functionSpace();
1357 view =
vf::project( _space=space, _expr=cst( boost::lexical_cast<value_type>(
"inf") ) );
1360 element_type element()
1365 element_type M_element;
1369 model_ptrtype M_model;
1373 template<
typename ModelType>
1374 boost::shared_ptr<EimFunctionNoSolve<ModelType>>
1375 eim_no_solve( boost::shared_ptr<ModelType> model )
1377 return boost::shared_ptr<EimFunctionNoSolve<ModelType>>(
new EimFunctionNoSolve<ModelType>( model ) );
1381 po::options_description eimOptions( std::string
const& prefix =
"");