31 #ifndef __CRBTrilinear_H
32 #define __CRBTrilinear_H 1
34 #include <boost/multi_array.hpp>
35 #include <boost/tuple/tuple.hpp>
36 #include "boost/tuple/tuple_io.hpp"
37 #include <boost/format.hpp>
38 #include <boost/foreach.hpp>
39 #include <boost/bimap.hpp>
40 #include <boost/bimap/support/lambda.hpp>
41 #include <boost/archive/text_oarchive.hpp>
42 #include <boost/archive/text_iarchive.hpp>
43 #include <boost/math/special_functions/fpclassify.hpp>
47 #include <boost/serialization/vector.hpp>
48 #include <boost/serialization/list.hpp>
49 #include <boost/serialization/string.hpp>
50 #include <boost/serialization/version.hpp>
51 #include <boost/serialization/split_member.hpp>
57 #include <Eigen/Dense>
60 #include <feel/feelcore/feel.hpp>
62 #include <feel/feelcore/parameter.hpp>
65 #include <feel/feelcrb/crbscm.hpp>
67 #include <feel/feelfilters/exporter.hpp>
82 template<
typename TruthModelType>
99 typedef TruthModelType truth_model_type;
100 typedef truth_model_type model_type;
101 typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
104 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
106 typedef typename parameterspace_type::element_ptrtype parameter_ptrtype;
108 typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
110 typedef boost::bimap< int, boost::tuple<double,double,double> > convergence_type;
112 typedef double value_type;
114 typedef typename convergence_type::value_type convergence;
118 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
122 typedef typename model_type::element_ptrtype element_ptrtype;
124 typedef typename model_type::backend_type backend_type;
125 typedef boost::shared_ptr<backend_type> backend_ptrtype;
126 typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
127 typedef typename model_type::vector_ptrtype vector_ptrtype;
128 typedef typename model_type::beta_vector_type beta_vector_type;
130 typedef std::vector<element_type> wn_type;
131 typedef boost::tuple< std::vector<wn_type> , std::vector<std::string> > export_vector_wn_type;
133 typedef Eigen::VectorXd vectorN_type;
134 typedef Eigen::MatrixXd matrixN_type;
138 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
142 typedef boost::shared_ptr<space_type> space_ptrtype;
144 typedef Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> > map_dense_matrix_type;
145 typedef Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, 1> > map_dense_vector_type;
147 typedef boost::tuple< vectorN_type , vectorN_type > solutions_tuple;
151 typedef boost::shared_ptr<export_type> export_ptrtype;
158 typedef boost::shared_ptr<scm_type> scm_ptrtype;
170 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
175 M_error_type( CRB_NO_RESIDUAL ),
189 super( ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
191 ( boost::format(
"%1%-%2%-%3%" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
193 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
195 M_backend( backend_type::build( vm ) ),
196 M_output_index( vm[
"crb.output-index"].template as<int>() ),
197 M_tolerance( vm[
"crb.error-max"].template as<double>() ),
198 M_iter_max( vm[
"crb.dimension-max"].template as<int>() ),
199 M_error_type(
CRBErrorType( vm[
"crb.error-type"].template as<int>() ) ),
210 po::variables_map
const&
vm,
211 truth_model_ptrtype
const & model )
213 super( ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
215 ( boost::format(
"%1%-%2%-%3%" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
217 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
219 M_backend( backend_type::build( vm ) ),
220 M_output_index( vm[
"crb.output-index"].template as<int>() ),
221 M_tolerance( vm[
"crb.error-max"].template as<double>() ),
222 M_iter_max( vm[
"crb.dimension-max"].template as<int>() ),
223 M_error_type(
CRBErrorType( vm[
"crb.error-type"].template as<int>() ) ),
233 LOG(INFO) <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
241 M_output_index( o.M_output_index ),
242 M_tolerance( o.M_tolerance ),
243 M_iter_max( o.M_iter_max ),
244 M_error_type( o.M_error_type ),
245 M_maxerror( o.M_maxerror ),
249 M_WNmu_complement( o.M_WNmu_complement )
279 boost::tuple<double,double>
lb(
size_type N, parameter_type
const& mu, vectorN_type& uN )
const;
288 void computationalTimeStatistics( std::string );
300 void updateJacobian(
const map_dense_vector_type& X, map_dense_matrix_type& J , parameter_type
const& mu ,
int N )
const;
306 void updateResidual(
const map_dense_vector_type& X, map_dense_vector_type& R , parameter_type
const& mu ,
int N )
const;
309 void displayVector(
const map_dense_vector_type& V )
const ;
310 void displayVector(
const vectorN_type& V )
const ;
311 void displayMatrix(
const matrixN_type& M )
const ;
314 void orthonormalize(
size_type N, wn_type& wn,
int Nm = 1 );
315 void checkOrthonormality(
int N,
const wn_type& wn )
const;
322 bool print = this->
vm()[
"crb.print-error-during-rb-construction"].template as<bool>();
328 return CRB_NO_RESIDUAL;
337 sampling_ptrtype wnmu ( )
const
344 bool use = this->
vm()[
"crb.run-on-WNmu"].template as<bool>();
386 boost::tuple<double,double, solutions_tuple, double, double, double, boost::tuple<double,double,double> > run( parameter_type
const& mu,
double eps = 1e-6,
int N = -1 );
388 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P ){};
394 M_Dmu = M_model->parameterSpace();
401 LOG(INFO) <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
403 M_scm->setTruthModel( M_model );
423 bool show = this->
vm()[
"crb.show-mu-selection"].template as<bool>();
440 bool rebuild = this->
vm()[
"crb.rebuild-database"].template as<bool>();
454 WorldComm
const& worldComm()
const {
return Environment::worldComm() ; }
459 std::vector < std::vector < matrixN_type> > M_Aqm_tril_pr;
460 mutable matrixN_type M_bilinear_terms;
461 mutable vectorN_type M_linear_terms;
462 boost::shared_ptr<SolverNonLinear<double> > M_nlsolver;
467 truth_model_ptrtype M_model;
468 backend_ptrtype M_backend;
475 parameterspace_ptrtype M_Dmu;
477 sampling_ptrtype M_Xi;
479 sampling_ptrtype M_WNmu;
480 sampling_ptrtype M_WNmu_complement;
485 export_ptrtype exporter;
488 friend class boost::serialization::access;
492 template<
class Archive>
493 void save( Archive & ar,
const unsigned int version )
const;
495 template<
class Archive>
496 void load( Archive & ar,
const unsigned int version ) ;
498 BOOST_SERIALIZATION_SPLIT_MEMBER()
505 bool orthonormalize_primal;
506 bool orthonormalize_dual;
507 bool solve_dual_problem;
509 convergence_type M_rbconv;
512 std::vector < std::vector<matrixN_type> > M_Aqm_pr;
515 std::vector < std::vector<vectorN_type> > M_Fqm_pr;
517 std::vector < std::vector<vectorN_type> > M_Lqm_pr;
519 std::vector<
int> M_index;
522 parameter_type M_current_mu;
523 int M_no_residual_index;
529 template<typename TruthModelType>
534 int proc_number = this->worldComm().globalRank();
536 bool rebuild_database = this->
vm()[
"crb.rebuild-database"].template as<bool>() ;
537 orthonormalize_primal = this->
vm()[
"crb.orthonormalize-primal"].template as<bool>() ;
540 if( proc_number == 0 ) std::cout <<
"Offline CRB starts, this may take a while until Database is computed...\n";
541 LOG(INFO) <<
"[CRB::offline] Starting offline for output " << M_output_index <<
"\n";
542 LOG(INFO) <<
"[CRB::offline] initialize underlying finite element model\n";
544 LOG( INFO )<<
" -- model init done in " << ti.elapsed() <<
"s";
552 if ( rebuild_database || M_N == 0)
557 LOG(INFO) <<
"[CRB::offline] compute random sampling\n";
559 int sampling_size = this->
vm()[
"crb.sampling-size"].template as<int>();
560 std::string file_name = ( boost::format(
"M_Xi_%1%") % sampling_size ).str();
561 std::ifstream file ( file_name );
565 M_Xi->randomize( sampling_size );
567 M_Xi->writeOnFile(file_name);
572 M_Xi->readFromFile(file_name);
575 M_WNmu->setSuperSampling( M_Xi );
577 LOG( INFO )<<
"[CRB offline] M_error_type = "<<M_error_type<<std::endl;
579 LOG(INFO) <<
" -- sampling init done in " << ti.elapsed() <<
"s";
584 if( M_error_type == CRB_NO_RESIDUAL )
585 mu = M_Dmu->element();
589 boost::tie( mu, index ) = M_Xi->min();
593 int size = mu.size();
604 LOG(INFO) <<
"[CRB::offline] allocate reduced basis data structures\n";
605 M_Aqm_pr.resize( M_model->Qa() );
606 for(
int q=0; q<M_model->Qa(); q++)
608 M_Aqm_pr[q].resize( 1 );
611 M_Aqm_tril_pr.resize( M_model->QaTri() );
612 for(
int q=0; q<M_model->QaTri(); q++)
614 M_Fqm_pr.resize( M_model->Ql( 0 ) );
616 for(
int q=0; q<M_model->Ql( 0 ); q++)
618 M_Fqm_pr[q].resize( 1 );
621 M_Lqm_pr.resize( M_model->Ql( M_output_index ) );
622 for(
int q=0; q<M_model->Ql( M_output_index ); q++)
623 M_Lqm_pr[q].resize( 1 );
630 if( proc_number == 0 )
632 std::cout<<
"we are going to enrich the reduced basis"<<std::endl;
633 std::cout<<
"there are "<<M_N<<
" elements in the database"<<std::endl;
635 LOG(INFO) <<
"we are going to enrich the reduced basis"<<std::endl;
636 LOG(INFO) <<
"there are "<<M_N<<
" elements in the database"<<std::endl;
640 LOG(INFO) <<
"[CRBTrilinear::offline] compute affine decomposition\n";
641 std::vector< std::vector<sparse_matrix_ptrtype> > Aqm;
642 std::vector< std::vector<sparse_matrix_ptrtype> > Aqm_tril;
643 std::vector< std::vector<std::vector<vector_ptrtype> > > Fqm;
645 boost::tie( boost::tuples::ignore, Aqm, Fqm ) = M_model->computeAffineDecomposition();
647 element_ptrtype u(
new element_type( M_model->functionSpace() ) );
649 LOG(INFO) <<
"[CRBTrilinear::offline] starting offline adaptive loop\n";
651 bool reuse_prec = this->
vm()[
"crb.reuse-prec"].template as<bool>() ;
653 bool use_predefined_WNmu = this->
vm()[
"crb.use-predefined-WNmu"].template as<bool>() ;
654 int N_log_equi = this->
vm()[
"crb.use-logEquidistributed-WNmu"].template as<int>() ;
655 int N_equi = this->
vm()[
"crb.use-equidistributed-WNmu"].template as<int>() ;
657 if( N_log_equi > 0 || N_equi > 0 )
658 use_predefined_WNmu =
true;
660 if ( use_predefined_WNmu )
662 std::string file_name = ( boost::format(
"SamplingWNmu") ).str();
663 std::ifstream file ( file_name );
667 std::vector< parameter_type > V;
669 __mu = M_Dmu->element();
670 __mu(0)= 1 ; __mu(1)= 1 ; V.push_back( __mu );
671 __mu(0)= 111112 ; __mu(1)= 1 ; V.push_back( __mu );
672 __mu(0)= 222223 ; __mu(1)= 1 ; V.push_back( __mu );
673 __mu(0)= 333334 ; __mu(1)= 1 ; V.push_back( __mu );
674 __mu(0)= 444445 , __mu(1)= 1 ; V.push_back( __mu );
675 __mu(0)= 555556 ; __mu(1)= 1 ; V.push_back( __mu );
676 __mu(0)= 666667 ; __mu(1)= 1 ; V.push_back( __mu );
677 __mu(0)= 777778 ; __mu(1)= 1 ; V.push_back( __mu );
678 __mu(0)= 888889 ; __mu(1)= 1 ; V.push_back( __mu );
679 __mu(0)= 1e+06 ; __mu(1)= 1 ; V.push_back( __mu );
680 __mu(0)= 8123 ; __mu(1)= 1 ; V.push_back( __mu );
681 __mu(0)= 9123 ; __mu(1)= 1 ; V.push_back( __mu );
682 __mu(0)=1.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
683 __mu(0)=2.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
684 __mu(0)=4.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
685 __mu(0)=912 ; __mu(1)= 1 ; V.push_back( __mu );
686 __mu(0)=1.123e3 ; __mu(1)= 1 ; V.push_back( __mu );
687 __mu(0)=4.123e3 ; __mu(1)= 1 ; V.push_back( __mu );
688 __mu(0)=7.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
689 __mu(0)=2123 ; __mu(1)= 1 ; V.push_back( __mu );
690 __mu(0)=6.123e3 ; __mu(1)= 1 ; V.push_back( __mu );
691 __mu(0)=3.123e3 ; __mu(1)= 1 ; V.push_back( __mu );
692 __mu(0)=3.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
693 __mu(0)=5.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
694 __mu(0)=9.123e4 ; __mu(1)= 1 ; V.push_back( __mu );
695 __mu(0)=812 ; __mu(1)= 1 ; V.push_back( __mu );
696 __mu(0)=5.111e3 ; __mu(1)= 1 ; V.push_back( __mu );
697 __mu(0)=5.124e2 ; __mu(1)= 1 ; V.push_back( __mu );
698 M_WNmu->setElements( V );
699 M_iter_max = M_WNmu->size();
704 int sampling_size = M_WNmu->readFromFile(file_name);
705 M_iter_max = sampling_size;
707 mu = M_WNmu->at( M_N );
710 LOG( INFO )<<
"[CRB::offline] read WNmu ( sampling size : "<<M_iter_max<<
" )";
714 if( M_error_type == CRB_NO_RESIDUAL || use_predefined_WNmu )
721 LOG(INFO) <<
"[CRBTrilinear::offline] strategy "<< M_error_type <<
"\n";
723 while ( M_maxerror > M_tolerance && M_N < M_iter_max )
726 boost::timer timer, timer2;
727 LOG(INFO) <<
"========================================"<<
"\n";
728 if( proc_number == this->worldComm().masterRank() )
729 std::cout<<
"construction of "<<M_N<<
"/"<<M_iter_max<<
" basis "<<std::endl;
730 LOG(INFO) <<
"N=" << M_N <<
"/" << M_iter_max <<
"( nb proc : "<<worldComm().globalSize()<<
")";
733 u->setName( ( boost::format(
"fem-primal-N%1%-proc%2%" ) % (M_N) % proc_number ).str() );
740 LOG(INFO) <<
"[CRB::offline] solving primal" <<
"\n";
741 *u = M_model->solve( mu );
747 if( ! use_predefined_WNmu )
748 M_WNmu->push_back( mu, index );
750 M_WNmu_complement = M_WNmu->complement();
752 M_WN.push_back( *u );
753 int number_of_added_elements=1;
754 M_N+=number_of_added_elements;
756 if ( orthonormalize_primal )
758 orthonormalize( M_N, M_WN, number_of_added_elements );
759 orthonormalize( M_N, M_WN, number_of_added_elements );
760 orthonormalize( M_N, M_WN, number_of_added_elements );
763 LOG(INFO) <<
"[CRB::offline] compute Aq_pr, Aq_du, Aq_pr_du" <<
"\n";
764 for (
size_type q = 0; q < M_model->Qa(); ++q )
766 M_Aqm_pr[q][0].conservativeResize( M_N, M_N );
769 for (
size_type i = M_N-number_of_added_elements; i < M_N; i++ )
773 M_Aqm_pr[q][0]( i, j ) = Aqm[q][0]->energy( M_WN[i], M_WN[j] );
777 for (
size_type j=M_N-number_of_added_elements; j < M_N; j++ )
781 M_Aqm_pr[q][0]( i, j ) = Aqm[q][0]->energy( M_WN[i], M_WN[j] );
787 LOG(INFO) <<
"[CRBTrilinear::offline] compute Fq_pr" <<
"\n";
789 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
791 M_Fqm_pr[q][0].conservativeResize( M_N );
793 for (
size_type l = 1; l <= number_of_added_elements; ++l )
796 M_Fqm_pr[q][0]( index ) = M_model->Fqm( 0, q, 0, M_WN[index] );
800 LOG(INFO) <<
"[CRB::offline] compute Lq_pr" <<
"\n";
802 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
804 M_Lqm_pr[q][0].conservativeResize( M_N );
806 for (
size_type l = 1; l <= number_of_added_elements; ++l )
809 M_Lqm_pr[q][0]( index ) = M_model->Fqm( M_output_index, q, 0, M_WN[index] );
813 sparse_matrix_ptrtype trilinear_form;
814 for (
size_type q = 0; q < M_model->QaTri(); ++q )
816 M_Aqm_tril_pr[q].resize( M_N );
817 for (
int k=0 ; k<M_N; k++)
821 trilinear_form = M_model->computeTrilinearForm( M_WN[k] );
823 M_Aqm_tril_pr[q][k].conservativeResize( M_N, M_N );
824 for (
int i = 0; i < M_N; ++i )
826 for (
int j = 0; j < M_N; ++j )
828 M_Aqm_tril_pr[q][k]( i, j ) = trilinear_form->energy( M_WN[j], M_WN[i] );
836 if ( ! use_predefined_WNmu )
844 mu = M_Dmu->element();
846 BOOST_FOREACH(
auto _mu, *M_WNmu )
852 while( already_exist );
858 if( M_N < M_iter_max )
860 mu = M_WNmu->at( M_N );
865 M_rbconv.insert( convergence( M_N, boost::make_tuple(M_maxerror,delta_pr,delta_du) ) );
868 LOG(INFO) <<
"time: " << timer.elapsed() <<
"\n";
869 LOG(INFO) <<
"========================================"<<
"\n";
876 LOG( INFO )<<
"number of elements in the reduced basis : "<<M_N<<
" ( nb proc : "<<worldComm().globalSize()<<
")";
877 LOG(INFO) <<
" index choosen : ";
878 BOOST_FOREACH(
auto id, M_index )
881 bool visualize_basis = this->
vm()[
"crb.visualize-basis"].template as<bool>() ;
883 if ( visualize_basis )
885 std::vector<wn_type> wn;
886 std::vector<std::string> names;
887 wn.push_back( M_WN );
888 names.push_back(
"primal" );
891 if ( orthonormalize_primal )
892 LOG(INFO)<<
"[CRB::offline] Basis functions have been exported but warning elements have been orthonormalized";
895 LOG( INFO ) <<
"Offline CRB is done";
902 template<
typename TruthModelType>
906 std::vector<int> index_vector;
907 sampling_ptrtype S = M_WNmu->searchNearestNeighbors( mu, 1 , index_vector);
908 neighbor = S->at( 0 );
909 index = index_vector[0];
913 template<
typename TruthModelType>
914 boost::tuple<double,double>
918 if ( N > M_N ) N = M_N;
919 beta_vector_type betaAqm;
920 std::vector<beta_vector_type> betaFqm, betaLqm;
927 boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu );
938 vectorN_type R( (
int) N );
939 matrixN_type J( (
int) N , (
int) N );
944 uN.setZero( (
int) N );
948 if( this->
vm()[
"crb.cvg-study"].
template as<bool>() ==
true )
955 double *r_data = R.data();
956 double *j_data = J.data();
957 double *uN_data = uN.data();
959 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , 1> > map_R ( r_data, N );
960 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , 1> > map_uN ( uN_data, N );
961 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , Eigen::Dynamic> > map_J ( j_data, N , N );
966 bool use_continuity = this->
vm()[
"crb.use-continuity"].template as<bool>();
970 Nmax=std::max( 1.0,std::max( std::ceil( std::log( gr ) ),std::ceil( std::log( pr )-std::log( 1.e-2 ) ) ) );
972 for (
int i = 0; i < Nmax; ++i )
975 double current_Grashofs;
976 double current_Prandtl;
979 int denom = ( Nmax==1 )?1:Nmax-1;
980 current_Grashofs = math::exp( math::log( 1. )+i*( math::log( gr )-math::log( 1. ) )/denom );
981 current_Prandtl = math::exp( math::log( 1.e-2 )+i*( math::log( pr )-math::log( 1.e-2 ) )/denom );
982 LOG( INFO ) <<
"[CRBTrilinear::lb] i/N = " << i+1 <<
"/" << Nmax ;
983 LOG( INFO ) <<
"[CRBTrilinear::lb] intermediary Grashof = " << current_Grashofs;
984 LOG( INFO ) <<
"[CRBTrilinear::lb] and Prandtl = " << current_Prandtl ;
988 current_Grashofs = gr;
989 current_Prandtl = pr;
992 current_mu << current_Grashofs, current_Prandtl;
997 M_nlsolver->map_dense_jacobian = boost::bind( &
self_type::updateJacobian, boost::ref( *
this ), _1, _2 , current_mu , N );
998 M_nlsolver->map_dense_residual = boost::bind( &
self_type::updateResidual, boost::ref( *
this ), _1, _2 , current_mu , N );
999 M_nlsolver->setType( TRUST_REGION );
1000 M_nlsolver->solve( map_J , map_uN , map_R, 1e-12, 100);
1003 LOG(INFO) <<
"[CRBTrilinear::lb] solve with Newton done";
1005 double condition_number = 0;
1007 LOG( INFO ) <<
"[CRBTrilinear::lb] The condition number of jacobian done is not computed ";
1008 vectorN_type L ( (
int )N );
1011 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
1013 L += betaFqm[M_output_index][q][0]*M_Lqm_pr[q][0].head( N );
1015 output = L.dot( uN );
1016 LOG(INFO) <<
"[CRBTrilinear::lb] computation of the output done";
1020 return boost::make_tuple( output, condition_number );
1023 template<
typename TruthModelType>
1028 LOG(INFO) <<
"update linear terms \n";
1030 beta_vector_type betaAqm;
1031 std::vector<beta_vector_type> betaFqm, betaLqm;
1033 boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu );
1035 M_bilinear_terms.setZero( N , N );
1037 for (
size_type q = 0; q < M_model->Qa(); ++q )
1039 M_bilinear_terms += betaAqm[q][0]*M_Aqm_pr[q][0].block( 0, 0, N, N );
1042 M_linear_terms.setZero( N );
1044 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
1046 M_linear_terms += betaFqm[0][q][0]*M_Fqm_pr[q][0].head( N );
1050 template<
typename TruthModelType>
1054 LOG(INFO) <<
"updateJacobian \n";
1055 map_J = M_bilinear_terms;
1057 bool enable = this->
vm()[
"crb.enable-convection-terms"].template as<bool>();
1060 for (
size_type q = 0; q < M_model->QaTri(); ++q )
1062 for (
int k = 0 ; k < N; ++k)
1064 for (
int i = 0; i < N; ++i )
1066 map_J( i, k ) += ( M_Aqm_tril_pr[q][k].row( i ).head( N ) ).dot(map_X);
1067 map_J( i, k ) += ( M_Aqm_tril_pr[q][k].col( i ).head( N ) ).dot(map_X);
1073 if ( this->
vm()[
"crb.compute-error-on-reduced-residual-jacobian"].template as<bool>() )
1076 auto expansionX =
expansion( map_X , N , M_WN);
1077 auto J = M_model->jacobian( expansionX );
1078 matrixN_type model_reduced_jacobian( N , N );
1079 for(
int i=0; i<N; i++)
1081 for(
int j=0; j<N; j++)
1082 model_reduced_jacobian(i,j) = J->energy( M_WN[i], M_WN[j]);
1085 matrixN_type diff = map_J - model_reduced_jacobian;
1086 double max = diff.maxCoeff();
1087 std::cout<<std::setprecision(14)<<
"[CRB::updateJacobian] with X : "; this->displayVector(map_X); std::cout<<
" the max coeff of the difference jacobian matrix : "<<max<<std::endl;
1088 if( math::abs(max) > 1e-14 )
1090 std::cout<<
"here is the jacobian matrix containing difference between reduced jacobian from CRB and the model";
1091 this->displayMatrix( diff );
1093 std::cout<<std::endl;
1099 template<
typename TruthModelType>
1103 LOG(INFO) <<
" updateResidual \n";
1105 matrixN_type temp ( N , N );
1106 temp.setZero( N , N );
1108 map_R = M_linear_terms;
1109 map_R += M_bilinear_terms * map_X ;
1111 bool enable = this->
vm()[
"crb.enable-convection-terms"].template as<bool>();
1114 for (
size_type q = 0; q < M_model->QaTri(); ++q )
1116 for (
int k = 0 ; k < N; ++k)
1118 for (
int i = 0; i < N; ++i )
1120 temp( k, i ) = map_X.dot( M_Aqm_tril_pr[q][k].row( i ).head( N ) ) ;
1123 map_R += temp * map_X ;
1127 if ( this->
vm()[
"crb.compute-error-on-reduced-residual-jacobian"].
template as<bool>() )
1130 auto expansionX =
expansion( map_X , N , M_WN);
1131 auto R = M_model->residual( expansionX );
1132 vectorN_type model_reduced_residual( N );
1133 element_ptrtype eltR(
new element_type( M_model->functionSpace() ) );
1134 for(
int i=0; i<eltR->localSize();i++)
1135 eltR->operator()(i)=R->operator()(i);
1136 for(
int i=0; i<N; i++)
1137 model_reduced_residual(i) =
inner_product( *eltR , M_WN[i] );
1139 vectorN_type diff = map_R - model_reduced_residual;
1140 double max = diff.maxCoeff();
1141 std::cout<<std::setprecision(14)<<
"[CRB::updateResidual] with X :"; this->displayVector(map_X); std::cout<<
" Residual :";this->displayVector(map_R);
1142 std::cout<<
" the max error : "<<max<<std::endl;
1143 if( math::abs(max) > 1e-14 )
1145 std::cout<<
"here is the residual containing difference between reduced residual from CRB and the model";
1146 this->displayVector( diff );
1148 std::cout<<std::endl;
1152 template<
typename TruthModelType>
1157 std::cout<<std::setprecision(14)<<
" ( ";
1158 for(
int i=0; i<size-1;i++)
1159 std::cout<<V(i)<<
" , ";
1160 std::cout<<V(size-1)<<
" ) ";
1163 template<
typename TruthModelType>
1165 CRBTrilinear<TruthModelType>::displayVector(
const vectorN_type& V )
const
1168 std::cout<<std::setprecision(14)<<
" ( ";
1169 for(
int i=0; i<size-1;i++)
1170 std::cout<<V(i)<<
" , ";
1171 std::cout<<V(size-1)<<
" ) ";
1174 template<
typename TruthModelType>
1176 CRBTrilinear<TruthModelType>::displayMatrix(
const matrixN_type& M )
const
1178 std::cout<<std::setprecision(14)<<
"\n[ ";
1179 int cols = M.cols();
1180 int rows = M.rows();
1181 for(
int i=0; i<rows-1;i++)
1183 for(
int j=0; j<cols-1;j++)
1184 std::cout<<M(i,j)<<
" , ";
1185 std::cout<<M(i,cols-1)<<
" ]"<<std::endl;std::cout<<
"[ ";
1187 for(
int j=0; j<cols-1;j++)
1188 std::cout<<M(rows-1,j)<<
" , ";
1189 std::cout<<M(rows-1,cols-1)<<
" ]"<<std::endl;std::cout<<
" ";
1197 template<
typename TruthModelType>
1199 CRBTrilinear<TruthModelType>::orthonormalize(
size_type N, wn_type& wn,
int Nm )
1201 int proc_number = this->worldComm().globalRank();
1202 if( proc_number == 0 ) std::cout <<
" -- orthonormalization (Gram-Schmidt)\n";
1203 DVLOG(2) <<
"[CRB::orthonormalize] orthonormalize basis for N=" << N <<
"\n";
1204 DVLOG(2) <<
"[CRB::orthonormalize] orthonormalize basis for WN="
1205 << wn.size() <<
"\n";
1206 DVLOG(2) <<
"[CRB::orthonormalize] starting ...\n";
1210 for (
size_type j = std::max( i+1,N-Nm ); j < N; ++j )
1212 value_type __rij_pr = M_model->scalarProduct( wn[i], wn[ j ] );
1213 wn[j].add( -__rij_pr, wn[i] );
1220 value_type __rii_pr = math::sqrt( M_model->scalarProduct( wn[i], wn[i] ) );
1221 wn[i].scale( 1./__rii_pr );
1224 DVLOG(2) <<
"[CRB::orthonormalize] finished ...\n";
1225 DVLOG(2) <<
"[CRB::orthonormalize] copying back results in basis\n";
1227 if ( this->
vm()[
"crb.check.gs"].
template as<int>() )
1228 checkOrthonormality( N , wn );
1232 template <
typename TruthModelType>
1234 CRBTrilinear<TruthModelType>::checkOrthonormality (
int N,
const wn_type& wn )
const
1239 throw std::logic_error(
"[CRB::checkOrthonormality] ERROR : size of wn is zero" );
1242 if ( orthonormalize_primal*orthonormalize_dual==0 && this->worldComm().globalRank() == this->worldComm().masterRank() )
1244 std::cout<<
"Warning : calling checkOrthonormality is called but ";
1245 std::cout<<
" orthonormalize_dual = "<<orthonormalize_dual;
1246 std::cout<<
" and orthonormalize_primal = "<<orthonormalize_primal<<std::endl;
1251 I.setIdentity( N, N );
1253 for (
int i = 0; i < N; ++i )
1255 for (
int j = 0; j < N; ++j )
1257 A( i, j ) = M_model->scalarProduct( wn[i], wn[j] );
1262 DVLOG(2) <<
"orthonormalization: " << A.norm() <<
"\n";
1263 if( this->worldComm().globalRank() == this->worldComm().masterRank() )
1264 std::cout <<
" o check : " << A.norm() <<
" (should be 0)\n";
1268 template <
typename TruthModelType>
1274 std::vector<wn_type> vect_wn=export_vector_wn.template get<0>();
1275 std::vector<std::string> vect_names=export_vector_wn.template get<1>();
1277 if ( vect_wn.size()==0 )
1279 throw std::logic_error(
"[CRB::exportBasisFunctions] ERROR : there are no wn_type to export" );
1283 auto first_wn = vect_wn[0];
1284 auto first_element = first_wn[0];
1286 exporter->step( 0 )->setMesh( first_element.functionSpace()->mesh() );
1288 BOOST_FOREACH(
auto wn , vect_wn )
1293 throw std::logic_error(
"[CRB::exportBasisFunctions] ERROR : there are no element to export" );
1296 int element_number=0;
1299 BOOST_FOREACH(
auto element, wn )
1302 std::string basis_name = vect_names[basis_number];
1303 std::string number = ( boost::format(
"%1%_with_parameters" ) %element_number ).str() ;
1304 mu = M_WNmu->at( element_number );
1307 for (
int i=0; i<mu.size(); i++ )
1309 mu_str= mu_str + ( boost::format(
"_%1%" ) %mu[i] ).str() ;
1312 std::string
name = basis_name + number + mu_str;
1313 exporter->step( 0 )->add( name, element );
1324 template<
typename TruthModelType>
1329 LOG(INFO)<<
" List of parameter selectionned during the offline algorithm \n";
1330 for(
int k=0;k<M_WNmu->size();k++)
1332 LOG(INFO)<<
" mu "<<k<<
" = [ ";
1334 for(
int i=0; i<_mu.size()-1; i++ )
1335 LOG(INFO)<<_mu(i)<<
" , ";
1336 LOG(INFO)<<_mu( _mu.size()-1 )<<
" ] \n";
1340 template<
typename TruthModelType>
1341 boost::tuple<double,double, typename CRBTrilinear<TruthModelType>::solutions_tuple, double, double, double, boost::tuple<double,double,double> >
1346 int Nwn_max =
vm()[
"crb.dimension-max"].template as<int>();
1362 auto o =
lb( Nwn, mu, uN );
1363 double output = o.template get<0>();
1365 double condition_number = o.template get<1>();
1366 auto upper_bounds = boost::make_tuple(1 , 1, 1 );
1367 auto solutions = boost::make_tuple( uN , uN );
1368 return boost::make_tuple( output , Nwn , solutions, condition_number, 1, 1, upper_bounds);
1372 template<
typename TruthModelType>
1373 typename CRBTrilinear<TruthModelType>::element_type
1385 auto o =
lb( Nwn, mu, uN );
1386 int size = uN.size();
1387 return Feel::expansion( M_WN, uN , Nwn);
1391 template<
typename TruthModelType>
1405 FEELPP_ASSERT( Nwn <= M_WN.size() )( Nwn )( M_WN.size() ).error(
"invalid expansion size ( N and M_WN ) ");
1406 return Feel::expansion( WN, u, Nwn );
1410 template<
typename TruthModelType>
1415 double min=0,max=0,mean=0,standard_deviation=0;
1417 int n_eval = option(_name=
"crb.computational-time-neval").template as<int>();
1419 Eigen::Matrix<double, Eigen::Dynamic, 1> time_crb;
1420 time_crb.resize( n_eval );
1422 sampling_ptrtype Sampling(
new sampling_type( M_Dmu ) );
1423 Sampling->logEquidistribute( n_eval );
1425 bool cvg = option(_name=
"crb.cvg-study").template as<bool>();
1427 double tol = option(_name=
"crb.online-tolerance").template as<double>();
1434 int proc_number = Environment::worldComm().globalRank();
1435 int master = Environment::worldComm().masterRank();
1438 std::string file_name =
"cvg-timing-crb.dat";
1441 if( proc_number == master )
1443 conv.open(file_name, std::ios::app);
1444 conv <<
"NbBasis" <<
"\t" <<
"min" <<
"\t"<<
"max" <<
"\t"<<
"mean"<<
"\t"<<
"standard_deviation" <<
"\n";
1452 BOOST_FOREACH(
auto mu, *Sampling )
1454 boost::mpi::timer tcrb;
1455 auto o = this->run( mu, tol , N);
1456 time_crb( mu_number ) = tcrb.elapsed() ;
1460 auto stat = M_model->computeStatistics( time_crb , appname );
1465 standard_deviation=stat(3);
1467 if( proc_number == master )
1468 conv << N <<
"\t" << min <<
"\t" << max<<
"\t"<< mean<<
"\t"<< standard_deviation<<
"\n";
1473 template<
typename TruthModelType>
1474 template<
class Archive>
1476 CRBTrilinear<TruthModelType>::save( Archive & ar,
const unsigned int version )
const
1478 int proc_number = this->worldComm().globalRank();
1480 LOG(INFO) <<
"[CRBTrilinear::save] version : "<<version<<std::endl;
1482 auto mesh = mesh_type::New();
1483 auto is_mesh_loaded = mesh->load( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
1485 if ( ! is_mesh_loaded )
1487 auto first_element = M_WN[0];
1488 mesh = first_element.functionSpace()->mesh() ;
1489 mesh->save( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
1492 auto Xh = space_type::New( mesh );
1494 ar & boost::serialization::base_object<super>( *this );
1495 ar & BOOST_SERIALIZATION_NVP( M_output_index );
1496 ar & BOOST_SERIALIZATION_NVP( M_N );
1497 ar & BOOST_SERIALIZATION_NVP( M_rbconv );
1498 ar & BOOST_SERIALIZATION_NVP( M_error_type );
1499 ar & BOOST_SERIALIZATION_NVP( M_Xi );
1500 ar & BOOST_SERIALIZATION_NVP( M_WNmu );
1501 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr );
1502 ar & BOOST_SERIALIZATION_NVP( M_Aqm_tril_pr );
1503 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
1504 ar & BOOST_SERIALIZATION_NVP( M_Lqm_pr );
1506 ar & BOOST_SERIALIZATION_NVP( M_current_mu );
1507 ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
1509 for(
int i=0; i<M_N; i++)
1510 ar & BOOST_SERIALIZATION_NVP( M_WN[i] );
1512 ar & BOOST_SERIALIZATION_NVP( M_maxerror );
1515 template<
typename TruthModelType>
1516 template<
class Archive>
1518 CRBTrilinear<TruthModelType>::load( Archive & ar,
const unsigned int version )
1521 int proc_number = this->worldComm().globalRank();
1523 LOG(INFO) <<
"[CRBTrilinear::load] version"<< version <<std::endl;
1530 LOG(INFO) <<
"[load] model not initialized, loading fdb files...\n";
1531 mesh = mesh_type::New();
1533 bool is_mesh_loaded = mesh->load( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
1535 Xh = space_type::New( mesh );
1536 LOG(INFO) <<
"[load] loading fdb files done.\n";
1540 LOG(INFO) <<
"[load] get mesh/Xh from model...\n";
1541 mesh = M_model->functionSpace()->mesh();
1542 Xh = M_model->functionSpace();
1543 LOG(INFO) <<
"[load] get mesh/Xh from model done.\n";
1546 ar & boost::serialization::base_object<super>( *this );
1547 ar & BOOST_SERIALIZATION_NVP( M_output_index );
1548 ar & BOOST_SERIALIZATION_NVP( M_N );
1550 ar & BOOST_SERIALIZATION_NVP( M_rbconv );
1552 ar & BOOST_SERIALIZATION_NVP( M_error_type );
1553 ar & BOOST_SERIALIZATION_NVP( M_Xi );
1554 ar & BOOST_SERIALIZATION_NVP( M_WNmu );
1555 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr );
1556 ar & BOOST_SERIALIZATION_NVP( M_Aqm_tril_pr );
1557 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
1558 ar & BOOST_SERIALIZATION_NVP( M_Lqm_pr );
1561 ar & BOOST_SERIALIZATION_NVP( M_current_mu );
1562 ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
1568 for(
int i = 0 ; i < M_N ; i++ )
1570 temp.setName( (boost::format(
"fem-primal-%1%" ) % ( i ) ).str() );
1571 ar & BOOST_SERIALIZATION_NVP( temp );
1575 ar & BOOST_SERIALIZATION_NVP( M_maxerror );
1577 LOG(INFO) <<
"[CRBTrilinear::load] end of load function" << std::endl;
1581 template<
typename TruthModelType>
1589 boost::archive::text_oarchive oa( ofs );
1596 template<
typename TruthModelType>
1600 if ( this->rebuildDB() )
1608 if ( !fs::exists( db ) )
1612 fs::ifstream ifs( db );
1616 boost::archive::text_iarchive ia( ifs );
1620 this->setIsLoaded(
true );