Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
crb_trilinear.hpp
Go to the documentation of this file.
1 /* -*- mode: c++ -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2009-11-24
7 
8  Copyright (C) 2009-2012 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
31 #ifndef __CRBTrilinear_H
32 #define __CRBTrilinear_H 1
33 
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>
44 #include <fstream>
45 
46 
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>
52 
53 #include <vector>
54 
55 #include <Eigen/Core>
56 #include <Eigen/LU>
57 #include <Eigen/Dense>
58 
60 #include <feel/feelcore/feel.hpp>
62 #include <feel/feelcore/parameter.hpp>
64 #include <feel/feelcrb/crbdb.hpp>
65 #include <feel/feelcrb/crbscm.hpp>
67 #include <feel/feelfilters/exporter.hpp>
68 
69 namespace Feel
70 {
71 
82 template<typename TruthModelType>
83 class CRBTrilinear : public CRBDB
84 {
85  typedef CRBDB super;
86 public:
87 
88 
92 
93 
95 
99  typedef TruthModelType truth_model_type;
100  typedef truth_model_type model_type;
101  typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
102 
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;
109 
110  typedef boost::bimap< int, boost::tuple<double,double,double> > convergence_type;
111 
112  typedef double value_type;
113 
114  typedef typename convergence_type::value_type convergence;
115 
117  typedef typename model_type::functionspace_type functionspace_type;
118  typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
119 
121  typedef typename model_type::element_type element_type;
122  typedef typename model_type::element_ptrtype element_ptrtype;
123 
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;
129 
130  typedef std::vector<element_type> wn_type;
131  typedef boost::tuple< std::vector<wn_type> , std::vector<std::string> > export_vector_wn_type;
132 
133  typedef Eigen::VectorXd vectorN_type;
134  typedef Eigen::MatrixXd matrixN_type;
135 
137  typedef typename model_type::mesh_type mesh_type;
138  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
139 
141  typedef typename model_type::space_type space_type;
142  typedef boost::shared_ptr<space_type> space_ptrtype;
143 
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;
146 
147  typedef boost::tuple< vectorN_type , vectorN_type > solutions_tuple;
148 
149  // ! export
151  typedef boost::shared_ptr<export_type> export_ptrtype;
152 
153  typedef CRBTrilinear self_type;
154 
155 
158  typedef boost::shared_ptr<scm_type> scm_ptrtype;
159 
161 
165 
168  :
169  super(),
170  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
171  M_model(),
172  M_output_index( 0 ),
173  M_tolerance( 1e-2 ),
174  M_iter_max( 3 ),
175  M_error_type( CRB_NO_RESIDUAL ),
176  M_Dmu( new parameterspace_type ),
177  M_Xi( new sampling_type( M_Dmu ) ),
178  M_WNmu( new sampling_type( M_Dmu, 1, M_Xi ) ),
179  M_WNmu_complement(),
180  M_scm( new scm_type( name, vm ) ),
181  exporter( Exporter<mesh_type>::New( "ensight" ) )
182  {
183 
184  }
185 
187  CRBTrilinear( std::string name, po::variables_map const& vm )
188  :
189  super( ( boost::format( "%1%" ) % vm["crb.error-type"].template as<int>() ).str(),
190  name,
191  ( boost::format( "%1%-%2%-%3%" ) % name % vm["crb.output-index"].template as<int>() % vm["crb.error-type"].template as<int>() ).str(),
192  vm ),
193  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
194  M_model(),
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>() ) ),
200  M_Dmu( new parameterspace_type ),
201  M_Xi( new sampling_type( M_Dmu ) ),
202  M_WNmu( new sampling_type( M_Dmu, 1, M_Xi ) ),
203  M_WNmu_complement(),
204  exporter( Exporter<mesh_type>::New( vm, "BasisFunction" ) )
205  {
206  }
207 
209  CRBTrilinear( std::string name,
210  po::variables_map const& vm,
211  truth_model_ptrtype const & model )
212  :
213  super( ( boost::format( "%1%" ) % vm["crb.error-type"].template as<int>() ).str(),
214  name,
215  ( boost::format( "%1%-%2%-%3%" ) % name % vm["crb.output-index"].template as<int>() % vm["crb.error-type"].template as<int>() ).str(),
216  vm ),
217  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
218  M_model(),
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>() ) ),
224  M_Dmu( new parameterspace_type ),
225  M_Xi( new sampling_type( M_Dmu ) ),
226  M_WNmu( new sampling_type( M_Dmu, 1, M_Xi ) ),
227  M_WNmu_complement(),
228  M_scm( new scm_type( name, vm ) ),
229  exporter( Exporter<mesh_type>::New( vm, "BasisFunction" ) )
230  {
231  this->setTruthModel( model );
232  if ( this->loadDB() )
233  LOG(INFO) << "Database " << this->lookForDB() << " available and loaded\n";
234  }
235 
236 
239  :
240  super( o ),
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 ),
246  M_Dmu( o.M_Dmu ),
247  M_Xi( o.M_Xi ),
248  M_WNmu( o.M_WNmu ),
249  M_WNmu_complement( o.M_WNmu_complement )
250  {}
251 
254  {}
255 
257 
261 
266  void findNearestNeighborInWNmu( parameter_type const& mu, parameter_type & neighbor, int & index ) const;
267 
279  boost::tuple<double,double> lb( size_type N, parameter_type const& mu, vectorN_type& uN ) const;
280 
286  convergence_type offline();
287 
288  void computationalTimeStatistics( std::string );
289 
294  void updateLinearTerms( parameter_type const& mu , int N ) const;
295 
300  void updateJacobian( const map_dense_vector_type& X, map_dense_matrix_type& J , parameter_type const& mu , int N ) const;
301 
306  void updateResidual( const map_dense_vector_type& X, map_dense_vector_type& R , parameter_type const& mu , int N ) const;
307 
308 
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 ;
312  // -------- from crb
313 
314  void orthonormalize( size_type N, wn_type& wn, int Nm = 1 );
315  void checkOrthonormality( int N, const wn_type& wn ) const;
316 
321  {
322  bool print = this->vm()["crb.print-error-during-rb-construction"].template as<bool>();
323  return print;
324  }
325 
326  CRBErrorType errorType() const
327  {
328  return CRB_NO_RESIDUAL;
329  }
330 
332  scm_ptrtype scm() const
333  {
334  return M_scm;
335  }
336 
337  sampling_ptrtype wnmu ( ) const
338  {
339  return M_WNmu;
340  }
341 
342  bool useWNmu()
343  {
344  bool use = this->vm()["crb.run-on-WNmu"].template as<bool>();
345  return use;
346  }
347 
348 
349  //only to compile
350  wn_type wndu() const
351  {
352  return M_WN;
353  }
354 
355 
356 
363 
370  element_type expansion( parameter_type const& mu , int N=-1);
371 
377  element_type expansion( vectorN_type const& u , int const N, wn_type const & WN ) const;
378 
379 
384  void exportBasisFunctions( const export_vector_wn_type& wn )const ;
385 
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 );
387 
388  void run( const double * X, unsigned long N, double * Y, unsigned long P ){};
389 
391  void setTruthModel( truth_model_ptrtype const& model )
392  {
393  M_model = model;
394  M_Dmu = M_model->parameterSpace();
395  M_Xi = sampling_ptrtype( new sampling_type( M_Dmu ) );
396 
397  if ( ! loadDB() )
398  M_WNmu = sampling_ptrtype( new sampling_type( M_Dmu ) );
399  else
400  {
401  LOG(INFO) << "Database " << this->lookForDB() << " available and loaded\n";
402  }
403  M_scm->setTruthModel( M_model );
404  }
405 
407  void setOfflineStep( bool b )
408  {
409  M_offline_step = b;
410  }
411 
412  wn_type wn() const
413  {
414  return M_WN;
415  }
416 
417 
422  {
423  bool show = this->vm()["crb.show-mu-selection"].template as<bool>();
424  return show;
425  }
426 
430  void printMuSelection( void );
431 
433  int dimension() const
434  {
435  return M_N;
436  }
437 
438  bool rebuildDB()
439  {
440  bool rebuild = this->vm()["crb.rebuild-database"].template as<bool>();
441  return rebuild;
442  }
443 
447  void saveDB();
448 
452  bool loadDB();
453 
454  WorldComm const& worldComm() const { return Environment::worldComm() ; }
456 
457 
458 private:
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;
463 
464 
465  // ------ from crb
466 
467  truth_model_ptrtype M_model;
468  backend_ptrtype M_backend;
469  int M_output_index;
470  double M_tolerance;
471  size_type M_iter_max;
472  CRBErrorType M_error_type;
473  double M_maxerror;
474  // parameter space
475  parameterspace_ptrtype M_Dmu;
476  // fine sampling of the parameter space
477  sampling_ptrtype M_Xi;
478  // sampling of parameter space to build WN
479  sampling_ptrtype M_WNmu;
480  sampling_ptrtype M_WNmu_complement;
481 
482  scm_ptrtype M_scm;
483 
484  //export
485  export_ptrtype exporter;
486 
487 
488  friend class boost::serialization::access;
489  // When the class Archive corresponds to an output archive, the
490  // & operator is defined similar to <<. Likewise, when the class Archive
491  // is a type of input archive the & operator is defined similar to >>.
492  template<class Archive>
493  void save( Archive & ar, const unsigned int version ) const;
494 
495  template<class Archive>
496  void load( Archive & ar, const unsigned int version ) ;
497 
498  BOOST_SERIALIZATION_SPLIT_MEMBER()
499 
500  // reduced basis space
501  wn_type M_WN;
502 
503  size_type M_N;
504 
505  bool orthonormalize_primal;
506  bool orthonormalize_dual;
507  bool solve_dual_problem;
508 
509  convergence_type M_rbconv;
510 
511  // left hand side ( bilinear )
512  std::vector < std::vector<matrixN_type> > M_Aqm_pr;
513 
514  // right hand side
515  std::vector < std::vector<vectorN_type> > M_Fqm_pr;
516  // output
517  std::vector < std::vector<vectorN_type> > M_Lqm_pr;
518 
519  std::vector<int> M_index;
520  int M_mode_number;
521 
522  parameter_type M_current_mu;
523  int M_no_residual_index;
524 
525  bool M_offline_step;
526 
527 };
528 
529 template<typename TruthModelType>
530 typename CRBTrilinear<TruthModelType>::convergence_type
531 CRBTrilinear<TruthModelType>::offline()
532 {
533 
534  int proc_number = this->worldComm().globalRank();
535 
536  bool rebuild_database = this->vm()["crb.rebuild-database"].template as<bool>() ;
537  orthonormalize_primal = this->vm()["crb.orthonormalize-primal"].template as<bool>() ;
538 
539  boost::timer ti;
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";
543  //M_model->initModel();
544  LOG( INFO )<< " -- model init done in " << ti.elapsed() << "s";
545 
546  parameter_type mu( M_Dmu );
547 
548  double delta_pr;
549  double delta_du;
550  size_type index;
551  //if M_N == 0 then there is not an already existing database
552  if ( rebuild_database || M_N == 0)
553  {
554 
555  ti.restart();
556 
557  LOG(INFO) << "[CRB::offline] compute random sampling\n";
558 
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 );
562  if( ! file )
563  {
564  // random sampling
565  M_Xi->randomize( sampling_size );
566  //M_Xi->equidistribute( this->vm()["crb.sampling-size"].template as<int>() );
567  M_Xi->writeOnFile(file_name);
568  }
569  else
570  {
571  M_Xi->clear();
572  M_Xi->readFromFile(file_name);
573  }
574 
575  M_WNmu->setSuperSampling( M_Xi );
576 
577  LOG( INFO )<<"[CRB offline] M_error_type = "<<M_error_type<<std::endl;
578 
579  LOG(INFO) << " -- sampling init done in " << ti.elapsed() << "s";
580  ti.restart();
581 
582  // empty sets
583  M_WNmu->clear();
584  if( M_error_type == CRB_NO_RESIDUAL )
585  mu = M_Dmu->element();
586  else
587  {
588  // start with M_C = { arg min mu, mu \in Xi }
589  boost::tie( mu, index ) = M_Xi->min();
590  }
591 
592 
593  int size = mu.size();
594  //std::cout << " -- WN size : " << M_WNmu->size() << "\n";
595 
596  // dimension of reduced basis space
597  M_N = 0;
598 
599  M_maxerror = 1e10;
600  delta_pr = 0;
601  delta_du = 0;
602  //boost::tie( M_maxerror, mu, index ) = maxErrorBounds( N );
603 
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++)
607  {
608  M_Aqm_pr[q].resize( 1 );
609  }
610 
611  M_Aqm_tril_pr.resize( M_model->QaTri() );
612  for(int q=0; q<M_model->QaTri(); q++)
613 
614  M_Fqm_pr.resize( M_model->Ql( 0 ) );
615 
616  for(int q=0; q<M_model->Ql( 0 ); q++)
617  {
618  M_Fqm_pr[q].resize( 1 );
619  }
620 
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 );
624 
625  }//end of if( rebuild_database )
626 #if 1
627  else
628  {
629  mu = M_current_mu;
630  if( proc_number == 0 )
631  {
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;
634  }
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;
637  }//end of else associated to if ( rebuild_databse )
638 #endif
639 
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;
644 
645  boost::tie( boost::tuples::ignore, Aqm, Fqm ) = M_model->computeAffineDecomposition();
646 
647  element_ptrtype u( new element_type( M_model->functionSpace() ) );
648 
649  LOG(INFO) << "[CRBTrilinear::offline] starting offline adaptive loop\n";
650 
651  bool reuse_prec = this->vm()["crb.reuse-prec"].template as<bool>() ;
652 
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>() ;
656 
657  if( N_log_equi > 0 || N_equi > 0 )
658  use_predefined_WNmu = true;
659 
660  if ( use_predefined_WNmu )
661  {
662  std::string file_name = ( boost::format("SamplingWNmu") ).str();
663  std::ifstream file ( file_name );
664  if( ! file )
665  {
666  M_WNmu->clear();
667  std::vector< parameter_type > V;
668  parameter_type __mu;
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();
700  }
701  else
702  {
703  M_WNmu->clear();
704  int sampling_size = M_WNmu->readFromFile(file_name);
705  M_iter_max = sampling_size;
706  }
707  mu = M_WNmu->at( M_N ); // first element
708  //std::cout<<" [use_predefined_WNmu] mu = \n"<<mu<<std::endl;
709 
710  LOG( INFO )<<"[CRB::offline] read WNmu ( sampling size : "<<M_iter_max<<" )";
711 
712  }
713 
714  if( M_error_type == CRB_NO_RESIDUAL || use_predefined_WNmu )
715  {
716  //in this case it makes no sens to check the estimated error
717  M_maxerror = 1e10;
718  }
719 
720 
721  LOG(INFO) << "[CRBTrilinear::offline] strategy "<< M_error_type <<"\n";
722 
723  while ( M_maxerror > M_tolerance && M_N < M_iter_max )
724  {
725 
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()<<")";
731 
732  // for a given parameter \p mu assemble the left and right hand side
733  u->setName( ( boost::format( "fem-primal-N%1%-proc%2%" ) % (M_N) % proc_number ).str() );
734 
735  mu.check();
736  u->zero();
737 
738  timer2.restart();
739 
740  LOG(INFO) << "[CRB::offline] solving primal" << "\n";
741  *u = M_model->solve( mu );
742 
743  //if( proc_number == this->worldComm().masterRank() ) std::cout << " -- primal problem solved in " << timer2.elapsed() << "s\n";
744  timer2.restart();
745 
746 
747  if( ! use_predefined_WNmu )
748  M_WNmu->push_back( mu, index );
749 
750  M_WNmu_complement = M_WNmu->complement();
751 
752  M_WN.push_back( *u );
753  int number_of_added_elements=1;
754  M_N+=number_of_added_elements;
755 
756  if ( orthonormalize_primal )
757  {
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 );
761  }
762 
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 )
765  {
766  M_Aqm_pr[q][0].conservativeResize( M_N, M_N );
767 
768  // only compute the last line and last column of reduced matrices
769  for ( size_type i = M_N-number_of_added_elements; i < M_N; i++ )
770  {
771  for ( size_type j = 0; j < M_N; ++j )
772  {
773  M_Aqm_pr[q][0]( i, j ) = Aqm[q][0]->energy( M_WN[i], M_WN[j] );
774  }
775  }
776 
777  for ( size_type j=M_N-number_of_added_elements; j < M_N; j++ )
778  {
779  for ( size_type i = 0; i < M_N; ++i )
780  {
781  M_Aqm_pr[q][0]( i, j ) = Aqm[q][0]->energy( M_WN[i], M_WN[j] );
782  }
783  }
784  }//loop over q
785 
786 
787  LOG(INFO) << "[CRBTrilinear::offline] compute Fq_pr" << "\n";
788 
789  for ( size_type q = 0; q < M_model->Ql( 0 ); ++q )
790  {
791  M_Fqm_pr[q][0].conservativeResize( M_N );
792 
793  for ( size_type l = 1; l <= number_of_added_elements; ++l )
794  {
795  int index = M_N-l;
796  M_Fqm_pr[q][0]( index ) = M_model->Fqm( 0, q, 0, M_WN[index] );
797  }
798  }//loop over q
799 
800  LOG(INFO) << "[CRB::offline] compute Lq_pr" << "\n";
801 
802  for ( size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
803  {
804  M_Lqm_pr[q][0].conservativeResize( M_N );
805 
806  for ( size_type l = 1; l <= number_of_added_elements; ++l )
807  {
808  int index = M_N-l;
809  M_Lqm_pr[q][0]( index ) = M_model->Fqm( M_output_index, q, 0, M_WN[index] );
810  }
811  }//loop over q
812 
813  sparse_matrix_ptrtype trilinear_form;
814  for (size_type q = 0; q < M_model->QaTri(); ++q )
815  {
816  M_Aqm_tril_pr[q].resize( M_N );
817  for (int k=0 ; k<M_N; k++)
818  {
819  //bring back the matrix associated to the trilinear form for a given basis function
820  //we do this here to use only one matrix
821  trilinear_form = M_model->computeTrilinearForm( M_WN[k] );
822 
823  M_Aqm_tril_pr[q][k].conservativeResize( M_N, M_N );
824  for ( int i = 0; i < M_N; ++i )
825  {
826  for ( int j = 0; j < M_N; ++j )
827  {
828  M_Aqm_tril_pr[q][k]( i, j ) = trilinear_form->energy( M_WN[j], M_WN[i] );
829  }//j
830  }//i
831  }//k
832  }// q
833 
834  timer2.restart();
835 
836  if ( ! use_predefined_WNmu )
837  {
838  bool already_exist;
839  do
840  {
841  //initialization
842  already_exist=false;
843  //pick randomly an element
844  mu = M_Dmu->element();
845  //make sure that the new mu is not already is M_WNmu
846  BOOST_FOREACH( auto _mu, *M_WNmu )
847  {
848  if( mu == _mu )
849  already_exist=true;
850  }
851  }
852  while( already_exist );
853  M_current_mu = mu;
854  }
855  else
856  {
857  //remmber that in this case M_iter_max = sampling size
858  if( M_N < M_iter_max )
859  {
860  mu = M_WNmu->at( M_N );
861  M_current_mu = mu;
862  }
863  }
864 
865  M_rbconv.insert( convergence( M_N, boost::make_tuple(M_maxerror,delta_pr,delta_du) ) );
866 
867  timer2.restart();
868  LOG(INFO) << "time: " << timer.elapsed() << "\n";
869  LOG(INFO) <<"========================================"<<"\n";
870 
871  //save DB after adding an element
872  this->saveDB();
873  }
874 
875 
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 )
879  LOG(INFO)<<id<<" ";
880  LOG(INFO)<<"\n";
881  bool visualize_basis = this->vm()["crb.visualize-basis"].template as<bool>() ;
882 
883  if ( visualize_basis )
884  {
885  std::vector<wn_type> wn;
886  std::vector<std::string> names;
887  wn.push_back( M_WN );
888  names.push_back( "primal" );
889  exportBasisFunctions( boost::make_tuple( wn ,names ) );
890 
891  if ( orthonormalize_primal )
892  LOG(INFO)<<"[CRB::offline] Basis functions have been exported but warning elements have been orthonormalized";
893  }
894 
895  LOG( INFO ) << "Offline CRB is done";
896 
897  return M_rbconv;
898 
899 }
900 
901 
902 template<typename TruthModelType>
903 void
905 {
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];
910  //std::cout<<"[CRBTrilinear::findNearestNeighborInWNmu] for Gr = "<<mu(0)<<" th nearest neighbor in WNmu is "<<neighbor(0)<<" at index "<<index<<std::endl;
911 }
912 
913 template<typename TruthModelType>
914 boost::tuple<double,double>
915 CRBTrilinear<TruthModelType>::lb( size_type N, parameter_type const& mu, vectorN_type & uN ) const
916 {
917 
918  if ( N > M_N ) N = M_N;
919  beta_vector_type betaAqm;
920  std::vector<beta_vector_type> betaFqm, betaLqm;
921 
922  //-- end of initialization step
923 
924  //vector containing outputs from time=time_step until time=time_for_output
925  double output;
926 
927  boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu );
928 
929  /*
930  ------> Here add Newton Method and continuous Gr,Pr
931  */
932  using namespace vf;
933  //Feel::ParameterSpace<2>::Element current_mu( mu );
934  parameter_type current_mu ;
935 
936  backend_ptrtype backend_primal_problem = backend_type::build( BACKEND_PETSC );
937 
938  vectorN_type R( (int) N );
939  matrixN_type J( (int) N , (int) N );
940 
941  //initialization of uN
942  //we look for the nearest neighbor of mu in the sampling WNmu
943  //let i the index of this neighbor in WNmu, we will set zeros in uN except at the i^th component where we will set 1
944  uN.setZero( (int) N );
945  parameter_type neighbor( M_Dmu );
946  int index;
947  findNearestNeighborInWNmu( mu, neighbor, index );
948  if( this->vm()["crb.cvg-study"].template as<bool>() == true )
949  {
950  //in this case, index may be smaller than uN.size
951  //so we do nothing
952  }
953  else
954  uN( index ) = 1;
955  double *r_data = R.data();
956  double *j_data = J.data();
957  double *uN_data = uN.data();
958 
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 );
962 
963  double gr = mu( 0 );
964  double pr = mu( 1 );
965 
966  bool use_continuity = this->vm()["crb.use-continuity"].template as<bool>();
967  int Nmax=1;
968 
969  if( use_continuity )
970  Nmax=std::max( 1.0,std::max( std::ceil( std::log( gr ) ),std::ceil( std::log( pr )-std::log( 1.e-2 ) ) ) );
971 
972  for ( int i = 0; i < Nmax; ++i )
973  {
974 
975  double current_Grashofs;
976  double current_Prandtl;
977  if( use_continuity )
978  {
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 ;
985  }
986  else
987  {
988  current_Grashofs = gr;
989  current_Prandtl = pr;
990  }
991 
992  current_mu << current_Grashofs, current_Prandtl;
993 
994  this->updateLinearTerms( current_mu , N );
995 
996  //M_nlsolver->setRelativeResidualTol( 1e-12 );
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);
1001  }
1002 
1003  LOG(INFO) << "[CRBTrilinear::lb] solve with Newton done";
1004 
1005  double condition_number = 0;
1006 
1007  LOG( INFO ) <<"[CRBTrilinear::lb] The condition number of jacobian done is not computed ";
1008  vectorN_type L ( ( int )N );
1009  L.setZero( N );
1010 
1011  for ( size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
1012  {
1013  L += betaFqm[M_output_index][q][0]*M_Lqm_pr[q][0].head( N );
1014  }
1015  output = L.dot( uN );
1016  LOG(INFO) << "[CRBTrilinear::lb] computation of the output done";
1017 
1018  //std::cout<<"[CRBTrilinear uN] : \n"<<uN<<std::endl;
1019 
1020  return boost::make_tuple( output, condition_number );
1021 
1022 }
1023 template<typename TruthModelType>
1024 void
1026 {
1027 
1028  LOG(INFO) << "update linear terms \n";
1029 
1030  beta_vector_type betaAqm;
1031  std::vector<beta_vector_type> betaFqm, betaLqm;
1032 
1033  boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu );
1034 
1035  M_bilinear_terms.setZero( N , N );
1036 
1037  for ( size_type q = 0; q < M_model->Qa(); ++q )
1038  {
1039  M_bilinear_terms += betaAqm[q][0]*M_Aqm_pr[q][0].block( 0, 0, N, N );
1040  }
1041 
1042  M_linear_terms.setZero( N );
1043 
1044  for ( size_type q = 0; q < M_model->Ql( 0 ); ++q )
1045  {
1046  M_linear_terms += betaFqm[0][q][0]*M_Fqm_pr[q][0].head( N );
1047  }
1048 
1049 }
1050 template<typename TruthModelType>
1051 void
1052 CRBTrilinear<TruthModelType>::updateJacobian( const map_dense_vector_type& map_X, map_dense_matrix_type& map_J , const parameter_type & mu , int N) const
1053 {
1054  LOG(INFO) << "updateJacobian \n";
1055  map_J = M_bilinear_terms;
1056 
1057  bool enable = this->vm()["crb.enable-convection-terms"].template as<bool>();
1058  if( enable )
1059  {
1060  for ( size_type q = 0; q < M_model->QaTri(); ++q )
1061  {
1062  for (int k = 0 ; k < N; ++k)
1063  {
1064  for ( int i = 0; i < N; ++i )
1065  {
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);
1068  }
1069  }
1070  }
1071  }
1072 
1073  if ( this->vm()["crb.compute-error-on-reduced-residual-jacobian"].template as<bool>() )
1074  {
1075  //bring the jacobian matrix from the model and then project it into the reduced basis
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++)
1080  {
1081  for(int j=0; j<N; j++)
1082  model_reduced_jacobian(i,j) = J->energy( M_WN[i], M_WN[j]);
1083  }
1084  //compute difference
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 )
1089  {
1090  std::cout<<"here is the jacobian matrix containing difference between reduced jacobian from CRB and the model";
1091  this->displayMatrix( diff );
1092  }
1093  std::cout<<std::endl;
1094  }
1095 }
1096 
1097 
1098 
1099 template<typename TruthModelType>
1100 void
1101 CRBTrilinear<TruthModelType>::updateResidual( const map_dense_vector_type& map_X, map_dense_vector_type& map_R , const parameter_type & mu, int N ) const
1102 {
1103  LOG(INFO) << " updateResidual \n";
1104 
1105  matrixN_type temp ( N , N );
1106  temp.setZero( N , N );
1107 
1108  map_R = M_linear_terms;
1109  map_R += M_bilinear_terms * map_X ;
1110 
1111  bool enable = this->vm()["crb.enable-convection-terms"].template as<bool>();
1112  if( enable )
1113  {
1114  for ( size_type q = 0; q < M_model->QaTri(); ++q )
1115  {
1116  for (int k = 0 ; k < N; ++k)
1117  {
1118  for ( int i = 0; i < N; ++i )
1119  {
1120  temp( k, i ) = map_X.dot( M_Aqm_tril_pr[q][k].row( i ).head( N ) ) ;
1121  }
1122  }
1123  map_R += temp * map_X ;
1124  }
1125  }
1126 
1127  if ( this->vm()["crb.compute-error-on-reduced-residual-jacobian"].template as<bool>() )
1128  {
1129  //bring the residual matrix from the model and then project it into the reduced basis
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] );
1138  //compute difference
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 )
1144  {
1145  std::cout<<"here is the residual containing difference between reduced residual from CRB and the model";
1146  this->displayVector( diff );
1147  }
1148  std::cout<<std::endl;
1149  }
1150 }
1151 
1152 template<typename TruthModelType>
1153 void
1154 CRBTrilinear<TruthModelType>::displayVector( const map_dense_vector_type& V ) const
1155 {
1156  int size=V.size();
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)<<" ) ";
1161 }
1162 
1163 template<typename TruthModelType>
1164 void
1165 CRBTrilinear<TruthModelType>::displayVector( const vectorN_type& V ) const
1166 {
1167  int size=V.size();
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)<<" ) ";
1172 }
1173 
1174 template<typename TruthModelType>
1175 void
1176 CRBTrilinear<TruthModelType>::displayMatrix( const matrixN_type& M ) const
1177 {
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++)
1182  {
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<<"[ ";
1186  }
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<<" ";
1190 
1191 
1192 }
1193 
1194 
1195 //---------------- from crb
1196 
1197 template<typename TruthModelType>
1198 void
1199 CRBTrilinear<TruthModelType>::orthonormalize( size_type N, wn_type& wn, int Nm )
1200 {
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";
1207 
1208  for ( size_type i = 0; i < N; ++i )
1209  {
1210  for ( size_type j = std::max( i+1,N-Nm ); j < N; ++j )
1211  {
1212  value_type __rij_pr = M_model->scalarProduct( wn[i], wn[ j ] );
1213  wn[j].add( -__rij_pr, wn[i] );
1214  }
1215  }
1216 
1217  // normalize
1218  for ( size_type i =N-Nm; i < N; ++i )
1219  {
1220  value_type __rii_pr = math::sqrt( M_model->scalarProduct( wn[i], wn[i] ) );
1221  wn[i].scale( 1./__rii_pr );
1222  }
1223 
1224  DVLOG(2) << "[CRB::orthonormalize] finished ...\n";
1225  DVLOG(2) << "[CRB::orthonormalize] copying back results in basis\n";
1226 
1227  if ( this->vm()["crb.check.gs"].template as<int>() )
1228  checkOrthonormality( N , wn );
1229 
1230 }
1231 
1232 template <typename TruthModelType>
1233 void
1234 CRBTrilinear<TruthModelType>::checkOrthonormality ( int N, const wn_type& wn ) const
1235 {
1236 
1237  if ( wn.size()==0 )
1238  {
1239  throw std::logic_error( "[CRB::checkOrthonormality] ERROR : size of wn is zero" );
1240  }
1241 
1242  if ( orthonormalize_primal*orthonormalize_dual==0 && this->worldComm().globalRank() == this->worldComm().masterRank() )
1243  {
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;
1247  }
1248 
1249  matrixN_type A, I;
1250  A.setZero( N, N );
1251  I.setIdentity( N, N );
1252 
1253  for ( int i = 0; i < N; ++i )
1254  {
1255  for ( int j = 0; j < N; ++j )
1256  {
1257  A( i, j ) = M_model->scalarProduct( wn[i], wn[j] );
1258  }
1259  }
1260 
1261  A -= I;
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";
1265  //FEELPP_ASSERT( A.norm() < 1e-14 )( A.norm() ).error( "orthonormalization failed.");
1266 }
1267 
1268 template <typename TruthModelType>
1269 void
1270 CRBTrilinear<TruthModelType>::exportBasisFunctions( const export_vector_wn_type& export_vector_wn )const
1271 {
1272 
1273 
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>();
1276 
1277  if ( vect_wn.size()==0 )
1278  {
1279  throw std::logic_error( "[CRB::exportBasisFunctions] ERROR : there are no wn_type to export" );
1280  }
1281 
1282 
1283  auto first_wn = vect_wn[0];
1284  auto first_element = first_wn[0];
1285 
1286  exporter->step( 0 )->setMesh( first_element.functionSpace()->mesh() );
1287  int basis_number=0;
1288  BOOST_FOREACH( auto wn , vect_wn )
1289  {
1290 
1291  if ( wn.size()==0 )
1292  {
1293  throw std::logic_error( "[CRB::exportBasisFunctions] ERROR : there are no element to export" );
1294  }
1295 
1296  int element_number=0;
1297  parameter_type mu;
1298 
1299  BOOST_FOREACH( auto element, wn )
1300  {
1301 
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 );
1305  std::string mu_str;
1306 
1307  for ( int i=0; i<mu.size(); i++ )
1308  {
1309  mu_str= mu_str + ( boost::format( "_%1%" ) %mu[i] ).str() ;
1310  }
1311 
1312  std::string name = basis_name + number + mu_str;
1313  exporter->step( 0 )->add( name, element );
1314  element_number++;
1315  }
1316  basis_number++;
1317  }
1318 
1319  exporter->save();
1320 
1321 }
1322 
1323 
1324 template<typename TruthModelType>
1325 void
1327 {
1328 
1329  LOG(INFO)<<" List of parameter selectionned during the offline algorithm \n";
1330  for(int k=0;k<M_WNmu->size();k++)
1331  {
1332  LOG(INFO)<<" mu "<<k<<" = [ ";
1333  parameter_type const& _mu = M_WNmu->at( k );
1334  for( int i=0; i<_mu.size()-1; i++ )
1335  LOG(INFO)<<_mu(i)<<" , ";
1336  LOG(INFO)<<_mu( _mu.size()-1 )<<" ] \n";
1337  }
1338 }
1339 
1340 template<typename TruthModelType>
1341 boost::tuple<double,double, typename CRBTrilinear<TruthModelType>::solutions_tuple, double, double, double, boost::tuple<double,double,double> >
1342 CRBTrilinear<TruthModelType>::run( parameter_type const& mu, double eps , int N)
1343 {
1344 
1345  //int Nwn = M_N;
1346  int Nwn_max = vm()["crb.dimension-max"].template as<int>();
1347 
1348  vectorN_type uN;
1349 
1350  int Nwn;
1351  if( N > 0 )
1352  {
1353  //in this case we want to fix Nwn
1354  Nwn = N;
1355  }
1356  else
1357  {
1358  //Nwn = Nwn_max;
1359  //M_N may be different of dimension-max
1360  Nwn = M_N;
1361  }
1362  auto o = lb( Nwn, mu, uN );
1363  double output = o.template get<0>();
1364 
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);
1369 }
1370 
1371 
1372 template<typename TruthModelType>
1373 typename CRBTrilinear<TruthModelType>::element_type
1375 {
1376  int Nwn;
1377 
1378  if( N > 0 )
1379  Nwn = N;
1380  else
1381  Nwn = M_N;
1382 
1383  vectorN_type uN;
1384 
1385  auto o = lb( Nwn, mu, uN );
1386  int size = uN.size();
1387  return Feel::expansion( M_WN, uN , Nwn);
1388 }
1389 
1390 
1391 template<typename TruthModelType>
1393 CRBTrilinear<TruthModelType>::expansion( vectorN_type const& u , int const N, wn_type const & WN) const
1394 {
1395  int Nwn;
1396 
1397  if( N > 0 )
1398  Nwn = N;
1399  else
1400  Nwn = M_N;
1401 
1402  //FEELPP_ASSERT( N == u.size() )( N )( u.size() ).error( "invalid expansion size");
1403  //FEELPP_ASSERT( Nwn == u.size() )( Nwn )( u.size() ).error( "invalid expansion size");
1404  //int size = uN.size();
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 );
1407 }
1408 
1409 
1410 template<typename TruthModelType>
1411 void
1413 {
1414 
1415  double min=0,max=0,mean=0,standard_deviation=0;
1416 
1417  int n_eval = option(_name="crb.computational-time-neval").template as<int>();
1418 
1419  Eigen::Matrix<double, Eigen::Dynamic, 1> time_crb;
1420  time_crb.resize( n_eval );
1421 
1422  sampling_ptrtype Sampling( new sampling_type( M_Dmu ) );
1423  Sampling->logEquidistribute( n_eval );
1424 
1425  bool cvg = option(_name="crb.cvg-study").template as<bool>();
1426  int dimension = this->dimension();
1427  double tol = option(_name="crb.online-tolerance").template as<double>();
1428 
1429  int N=dimension;//by default we perform only one time statistics
1430 
1431  if( cvg ) //if we want to compute time statistics for every crb basis then we start at 1
1432  N=1;
1433 
1434  int proc_number = Environment::worldComm().globalRank();
1435  int master = Environment::worldComm().masterRank();
1436 
1437  //write on a file
1438  std::string file_name = "cvg-timing-crb.dat";
1439 
1440  std::ofstream conv;
1441  if( proc_number == master )
1442  {
1443  conv.open(file_name, std::ios::app);
1444  conv << "NbBasis" << "\t" << "min" <<"\t"<< "max" <<"\t"<< "mean"<<"\t"<<"standard_deviation" << "\n";
1445  }
1446 
1447  //loop over basis functions (if cvg option)
1448  for(; N<=dimension; N++)
1449  {
1450 
1451  int mu_number = 0;
1452  BOOST_FOREACH( auto mu, *Sampling )
1453  {
1454  boost::mpi::timer tcrb;
1455  auto o = this->run( mu, tol , N);
1456  time_crb( mu_number ) = tcrb.elapsed() ;
1457  mu_number++;
1458  }
1459 
1460  auto stat = M_model->computeStatistics( time_crb , appname );
1461 
1462  min=stat(0);
1463  max=stat(1);
1464  mean=stat(2);
1465  standard_deviation=stat(3);
1466 
1467  if( proc_number == master )
1468  conv << N << "\t" << min << "\t" << max<< "\t"<< mean<< "\t"<< standard_deviation<<"\n";
1469  }//loop over basis functions
1470  conv.close();
1471 }
1472 
1473 template<typename TruthModelType>
1474 template<class Archive>
1475 void
1476 CRBTrilinear<TruthModelType>::save( Archive & ar, const unsigned int version ) const
1477 {
1478  int proc_number = this->worldComm().globalRank();
1479 
1480  LOG(INFO) <<"[CRBTrilinear::save] version : "<<version<<std::endl;
1481 
1482  auto mesh = mesh_type::New();
1483  auto is_mesh_loaded = mesh->load( _name="mymesh",_path=this->dbLocalPath(),_type="binary" );
1484 
1485  if ( ! is_mesh_loaded )
1486  {
1487  auto first_element = M_WN[0];
1488  mesh = first_element.functionSpace()->mesh() ;
1489  mesh->save( _name="mymesh",_path=this->dbLocalPath(),_type="binary" );
1490  }
1491 
1492  auto Xh = space_type::New( mesh );
1493 
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 );
1505 
1506  ar & BOOST_SERIALIZATION_NVP( M_current_mu );
1507  ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
1508 
1509  for(int i=0; i<M_N; i++)
1510  ar & BOOST_SERIALIZATION_NVP( M_WN[i] );
1511 
1512  ar & BOOST_SERIALIZATION_NVP( M_maxerror );
1513 }
1514 
1515 template<typename TruthModelType>
1516 template<class Archive>
1517 void
1518 CRBTrilinear<TruthModelType>::load( Archive & ar, const unsigned int version )
1519 {
1520 
1521  int proc_number = this->worldComm().globalRank();
1522 
1523  LOG(INFO) <<"[CRBTrilinear::load] version"<< version <<std::endl;
1524 
1525  mesh_ptrtype mesh;
1526  space_ptrtype Xh;
1527 
1528  if ( !M_model )
1529  {
1530  LOG(INFO) << "[load] model not initialized, loading fdb files...\n";
1531  mesh = mesh_type::New();
1532 
1533  bool is_mesh_loaded = mesh->load( _name="mymesh",_path=this->dbLocalPath(),_type="binary" );
1534 
1535  Xh = space_type::New( mesh );
1536  LOG(INFO) << "[load] loading fdb files done.\n";
1537  }
1538  else
1539  {
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";
1544  }
1545 
1546  ar & boost::serialization::base_object<super>( *this );
1547  ar & BOOST_SERIALIZATION_NVP( M_output_index );
1548  ar & BOOST_SERIALIZATION_NVP( M_N );
1549 
1550  ar & BOOST_SERIALIZATION_NVP( M_rbconv );
1551 
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 );
1559 
1560 
1561  ar & BOOST_SERIALIZATION_NVP( M_current_mu );
1562  ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
1563 
1564  element_type temp = Xh->element();
1565 
1566  M_WN.resize( M_N );
1567 
1568  for( int i = 0 ; i < M_N ; i++ )
1569  {
1570  temp.setName( (boost::format( "fem-primal-%1%" ) % ( i ) ).str() );
1571  ar & BOOST_SERIALIZATION_NVP( temp );
1572  M_WN[i] = temp;
1573  }
1574 
1575  ar & BOOST_SERIALIZATION_NVP( M_maxerror );
1576 
1577  LOG(INFO) << "[CRBTrilinear::load] end of load function" << std::endl;
1578 }
1579 
1580 
1581 template<typename TruthModelType>
1582 void
1584 {
1585  fs::ofstream ofs( this->dbLocalPath() / this->dbFilename() );
1586 
1587  if ( ofs )
1588  {
1589  boost::archive::text_oarchive oa( ofs );
1590  // write class instance to archive
1591  oa << *this;
1592  // archive and stream closed when destructors are called
1593  }
1594 }
1595 
1596 template<typename TruthModelType>
1597 bool
1599 {
1600  if ( this->rebuildDB() )
1601  return false;
1602 
1603  fs::path db = this->lookForDB();
1604 
1605  if ( db.empty() )
1606  return false;
1607 
1608  if ( !fs::exists( db ) )
1609  return false;
1610 
1611  //std::cout << "Loading " << db << "...\n";
1612  fs::ifstream ifs( db );
1613 
1614  if ( ifs )
1615  {
1616  boost::archive::text_iarchive ia( ifs );
1617  // write class instance to archive
1618  ia >> *this;
1619  //std::cout << "Loading " << db << " done...\n";
1620  this->setIsLoaded( true );
1621  // archive and stream closed when destructors are called
1622  return true;
1623  }
1624 
1625  return false;
1626 }
1627 
1628 
1629 } // Feel
1630 
1631 #endif /* __CRBTrilinear_H */

Generated on Sun Oct 20 2013 08:24:55 for Feel++ by doxygen 1.8.4