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
opusapp.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2011-06-18
7 
8  Copyright (C) 2011 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 */
29 #ifndef __OpusApp_H
30 #define __OpusApp_H 1
31 
32 #include <feel/feel.hpp>
33 
34 #include <boost/assign/std/vector.hpp>
35 #include <feel/feelcrb/crb.hpp>
36 #include <feel/feelcrb/eim.hpp>
38 #include <boost/serialization/version.hpp>
39 #include <boost/range/join.hpp>
40 #include <boost/regex.hpp>
41 #include <boost/assign/list_of.hpp>
42 
43 #include <Eigen/Core>
44 #include <Eigen/LU>
45 #include <Eigen/Dense>
46 
47 
48 namespace Feel
49 {
50 po::options_description opusapp_options( std::string const& prefix );
51 std::string _o( std::string const& prefix, std::string const& opt )
52 {
53  std::string o = prefix;
54 
55  if ( !o.empty() )
56  o += ".";
57 
58  return o + opt;
59 }
60 
61 
62 enum class SamplingMode
63 {
64  RANDOM = 0, EQUIDISTRIBUTED = 1, LOGEQUIDISTRIBUTED = 2
65 };
66 
67 #define prec 4
68 #define Pdim 7
69 #define fill ' '
70 #define dmanip std::scientific << std::setprecision( prec )
71 #define hdrmanip(N) std::setw(N) << std::setfill(fill) << std::right
72 #define tabmanip(N) std::setw(N) << std::setfill(fill) << std::right << dmanip
73 
74 
81  template<typename ModelType,
82  template < typename ReducedMethod > class RM=CRB,
83  template < typename ModelInterface > class Model=CRBModel>
84 class OpusApp : public Application
85 {
86  typedef Application super;
87 public:
88 
89  typedef double value_type;
90 
92  typedef ModelType model_type;
93  typedef boost::shared_ptr<ModelType> model_ptrtype;
94 
96  typedef typename model_type::functionspace_type functionspace_type;
97  typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
98 
99  typedef typename model_type::element_type element_type;
100 
101  typedef Eigen::VectorXd vectorN_type;
102 
103 #if 0
104  //old
105  typedef CRBModel<ModelType> crbmodel_type;
106  typedef boost::shared_ptr<crbmodel_type> crbmodel_ptrtype;
107  typedef CRB<crbmodel_type> crb_type;
108  typedef boost::shared_ptr<crb_type> crb_ptrtype;
109 #endif
110 
111  typedef Model<ModelType> crbmodel_type;
112  typedef boost::shared_ptr<crbmodel_type> crbmodel_ptrtype;
113  typedef RM<crbmodel_type> crb_type;
114  typedef boost::shared_ptr<crb_type> crb_ptrtype;
115 
117 
118  typedef typename ModelType::parameter_type parameter_type;
119  typedef std::vector< parameter_type > vector_parameter_type;
120 
121  typedef typename crb_type::sampling_ptrtype sampling_ptrtype;
122 
123  OpusApp()
124  :
125  super(),
126  M_mode( ( CRBModelMode )option(_name=_o( this->about().appName(),"run.mode" )).template as<int>() )
127  {
128  this->init();
129  }
130 
131  OpusApp( AboutData const& ad, po::options_description const& od )
132  :
133  super( ad, opusapp_options( ad.appName() ).add( od ).add( crbOptions() ).add( feel_options() ).add( eimOptions() ) ),
134  M_mode( ( CRBModelMode )option(_name=_o( this->about().appName(),"run.mode" )).template as<int>() )
135  {
136  this->init();
137  }
138 
139  OpusApp( AboutData const& ad, po::options_description const& od, CRBModelMode mode )
140  :
141  super( ad, opusapp_options( ad.appName() ).add( od ).add( crbOptions() ).add( feel_options() ).add( eimOptions() ) ),
142  M_mode( mode )
143  {
144  this->init();
145  }
146 
147  OpusApp( int argc, char** argv, AboutData const& ad, po::options_description const& od )
148  :
149  super( argc, argv, ad, opusapp_options( ad.appName() ).add( od ).add( crbOptions() ).add( feel_options() ).add( eimOptions() ) ),
150  M_mode( ( CRBModelMode )option(_name=_o( this->about().appName(),"run.mode" )).template as<int>() )
151  {
152  this->init();
153  }
154  OpusApp( int argc, char** argv, AboutData const& ad, po::options_description const& od, CRBModelMode mode )
155  :
156  super( argc, argv, ad, opusapp_options( ad.appName() ).add( od ).add( crbOptions() ).add( feel_options() ).add( eimOptions() ) ),
157  M_mode( mode )
158  {
159  this->init();
160  }
161  void init()
162  {
163  try
164  {
165  M_current_path = fs::current_path();
166 
167  std::srand( static_cast<unsigned>( std::time( 0 ) ) );
168  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
169  std::cout << this->about().appName() << std::endl;
170  LOG(INFO) << "[OpusApp] constructor " << this->about().appName() << "\n";
171 
172  // Check if user have given a name for result files repo
173  // Note : this name is also use for database storage
174  std::string results_repo_name;
175  if( this->vm().count("crb.results-repo-name") )
176  results_repo_name = option(_name="crb.results-repo-name").template as<std::string>();
177  else
178  results_repo_name = "default_repo";
179 
180  LOG(INFO) << "Name for results repo : " << results_repo_name << "\n";
181  this->changeRepository( boost::format( "%1%/%2%/" )
182  % this->about().appName()
183  % results_repo_name
184  );
185 
186  LOG(INFO) << "[OpusApp] ch repo" << "\n";
187  this->setLogs();
188  LOG(INFO) << "[OpusApp] set Logs" << "\n";
189  LOG(INFO) << "[OpusApp] mode:" << ( int )M_mode << "\n";
190 
191  model = crbmodel_ptrtype( new crbmodel_type( this->vm(),M_mode ) );
192  LOG(INFO) << "[OpusApp] get model done" << "\n";
193 
194  crb = crb_ptrtype( new crb_type( this->about().appName(),
195  this->vm() ,
196  model ) );
197  LOG(INFO) << "[OpusApp] get crb done" << "\n";
198 
199  //VLOG(1) << "[OpusApp] get crb done" << "\n";
200  //crb->setTruthModel( model );
201  //VLOG(1) << "[OpusApp] constructor done" << "\n";
202  }
203 
204  catch ( boost::bad_any_cast const& e )
205  {
206  std::cout << "[OpusApp] a bad any cast occured, probably a nonexistant or invalid command line/ config options\n";
207  std::cout << "[OpusApp] exception reason: " << e.what() << "\n";
208  }
209 
210  }
211 
212  void setMode( std::string const& mode )
213  {
214  if ( mode == "pfem" ) M_mode = CRBModelMode::PFEM;
215 
216  if ( mode == "crb" ) M_mode = CRBModelMode::CRB;
217 
218  if ( mode == "scm" ) M_mode = CRBModelMode::SCM;
219 
220  if ( mode == "scm_online" ) M_mode = CRBModelMode::SCM_ONLINE;
221 
222  if ( mode == "crb_online" ) M_mode = CRBModelMode::CRB_ONLINE;
223  }
224  void setMode( CRBModelMode mode )
225  {
226  M_mode = mode;
227  }
228 
229  void loadDB()
230  {
231  bool crb_use_predefined = option(_name="crb.use-predefined-WNmu").template as<bool>();
232  std::string file_name;
233  int NlogEquidistributed = option(_name="crb.use-logEquidistributed-WNmu").template as<int>();
234  int Nequidistributed = option(_name="crb.use-equidistributed-WNmu").template as<int>();
235  int NlogEquidistributedScm = option(_name="crb.scm.use-logEquidistributed-C").template as<int>();
236  int NequidistributedScm = option(_name="crb.scm.use-equidistributed-C").template as<int>();
237  typename crb_type::sampling_ptrtype Sampling( new typename crb_type::sampling_type( model->parameterSpace() ) );
238  if( NlogEquidistributed+Nequidistributed > 0 )
239  {
240  file_name = ( boost::format("SamplingWNmu") ).str();
241  if( NlogEquidistributed > 0 )
242  Sampling->logEquidistribute( NlogEquidistributed );
243  if( Nequidistributed > 0 )
244  Sampling->equidistribute( Nequidistributed );
245  Sampling->writeOnFile(file_name);
246  }
247  if( NlogEquidistributedScm+NequidistributedScm > 0 )
248  {
249  file_name = ( boost::format("SamplingC") ).str();
250  if( NlogEquidistributedScm > 0 )
251  Sampling->logEquidistribute( NlogEquidistributedScm );
252  if( NequidistributedScm > 0 )
253  Sampling->equidistribute( NequidistributedScm );
254  Sampling->writeOnFile(file_name);
255  }
256 
257  if ( M_mode == CRBModelMode::PFEM )
258  return;
259 
260  if ( !crb->scm()->isDBLoaded() || crb->scm()->rebuildDB() )
261  {
262  if ( M_mode == CRBModelMode::SCM )
263  {
264  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
265  std::cout << "No SCM DB available, do scm offline computations first...\n";
266  if( crb->scm()->doScmForMassMatrix() )
267  crb->scm()->setScmForMassMatrix( true );
268 
269  crb->scm()->offline();
270  }
271  }
272 
273  if ( !crb->isDBLoaded() || crb->rebuildDB() )
274  {
275  if ( M_mode == CRBModelMode::CRB )
276  //|| M_mode == CRBModelMode::SCM )
277  {
278  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
279  std::cout << "No CRB DB available, do crb offline computations...\n";
280  crb->setOfflineStep( true );
281  crb->offline();
282  }
283 
284  else if ( M_mode != CRBModelMode::SCM )
285  throw std::logic_error( "CRB/SCM Database could not be loaded" );
286  }
287 
288  if( crb->isDBLoaded() )
289  {
290  bool do_offline = false;
291  int current_dimension = crb->dimension();
292  int dimension_max = option(_name="crb.dimension-max").template as<int>();
293  int sampling_size = 0;
294  if( crb_use_predefined )
295  sampling_size = Sampling->readFromFile(file_name);
296 
297  if( sampling_size > current_dimension )
298  do_offline = true;
299 
300  if( current_dimension < dimension_max && !crb_use_predefined )
301  do_offline=true;
302 
303  if( ! do_offline )
304  {
305  crb->loadSCMDB();
306  }
307 
308  if( do_offline )
309  {
310  crb->setOfflineStep( true );
311  crb->offline();
312  }
313  }
314  }
315 
316 
317  FEELPP_DONT_INLINE
318  void run()
319  {
320 
321  bool export_solution = option(_name=_o( this->about().appName(),"export-solution" )).template as<bool>();
322 
323  int proc_number = Environment::worldComm().globalRank();
324 
325  if ( this->vm().count( "help" ) )
326  {
327  std::cout << this->optionsDescription() << "\n";
328  return;
329  }
330 
331  this->loadDB();
332  int run_sampling_size = option(_name=_o( this->about().appName(),"run.sampling.size" )).template as<int>();
333  SamplingMode run_sampling_type = ( SamplingMode )option(_name=_o( this->about().appName(),"run.sampling.mode" )).template as<int>();
334  int output_index = option(_name="crb.output-index").template as<int>();
335  //int output_index = option(_name=_o(this->about().appName(),"output.index")).template as<int>();
336 
337  typename crb_type::sampling_ptrtype Sampling( new typename crb_type::sampling_type( model->parameterSpace() ) );
338 
339  int n_eval_computational_time = option(_name="eim.computational-time-neval").template as<int>();
340  bool compute_fem = option(_name="crb.compute-fem-during-online").template as<bool>();
341  bool compute_stat = option(_name="crb.compute-stat").template as<bool>();
342 
343  if( n_eval_computational_time > 0 )
344  {
345  compute_fem = false;
346  auto eim_sc_vector = model->scalarContinuousEim();
347  auto eim_sd_vector = model->scalarDiscontinuousEim();
348  int size1 = eim_sc_vector.size();
349  int size2 = eim_sd_vector.size();
350  if( size1 + size2 == 0 )
351  throw std::logic_error( "[OpusApp] no eim object detected" );
352 
353  std::string appname = this->about().appName();
354  for(int i=0; i<size1; i++)
355  eim_sc_vector[i]->computationalTimeStatistics(appname);
356  for(int i=0; i<size2; i++)
357  eim_sd_vector[i]->computationalTimeStatistics(appname);
358 
359  run_sampling_size = 0;
360  }
361  n_eval_computational_time = option(_name="crb.computational-time-neval").template as<int>();
362  if( n_eval_computational_time > 0 )
363  {
364  if( ! option(_name="crb.cvg-study").template as<bool>() )
365  {
366  compute_fem = false;
367  run_sampling_size = 0;
368  }
369  std::string appname = this->about().appName();
370  crb->computationalTimeStatistics( appname );
371  }
372 
373 
374  switch ( run_sampling_type )
375  {
376  default:
377  case SamplingMode::RANDOM:
378  Sampling->randomize( run_sampling_size );
379  break;
380 
381  case SamplingMode::EQUIDISTRIBUTED:
382  Sampling->equidistribute( run_sampling_size );
383  break;
384 
385  case SamplingMode::LOGEQUIDISTRIBUTED:
386  Sampling->logEquidistribute( run_sampling_size );
387  break;
388  }
389 
390 
391  std::map<CRBModelMode,std::vector<std::string> > hdrs;
392  using namespace boost::assign;
393  std::vector<std::string> pfemhdrs = boost::assign::list_of( "FEM Output" )( "FEM Time" );
394  std::vector<std::string> crbhdrs = boost::assign::list_of( "FEM Output" )( "FEM Time" )( "RB Output" )( "Error Bounds" )( "CRB Time" )( "output error" )( "Conditionning" )( "l2_error" )( "h1_error" );
395  std::vector<std::string> scmhdrs = boost::assign::list_of( "Lb" )( "Lb Time" )( "Ub" )( "Ub Time" )( "FEM" )( "FEM Time" )( "output error" );
396  std::vector<std::string> crbonlinehdrs = boost::assign::list_of( "RB Output" )( "Error Bounds" )( "CRB Time" );
397  std::vector<std::string> scmonlinehdrs = boost::assign::list_of( "Lb" )( "Lb Time" )( "Ub" )( "Ub Time" )( "Rel.(FEM-Lb)" );
398  hdrs[CRBModelMode::PFEM] = pfemhdrs;
399  hdrs[CRBModelMode::CRB] = crbhdrs;
400  hdrs[CRBModelMode::SCM] = scmhdrs;
401  hdrs[CRBModelMode::CRB_ONLINE] = crbonlinehdrs;
402  hdrs[CRBModelMode::SCM_ONLINE] = scmonlinehdrs;
403  std::ostringstream ostr;
404 
405  //if( boost::is_same< crbmodel_type , crbmodelbilinear_type >::value )
406  {
407  if( crb->printErrorDuringOfflineStep() )
408  crb->printErrorsDuringRbConstruction();
409  if ( crb->showMuSelection() )
410  crb->printMuSelection();
411  }
412 
413  auto exporter = Exporter<typename crbmodel_type::mesh_type>::New( "ensight" );
414  if( export_solution )
415  exporter->step( 0 )->setMesh( model->functionSpace()->mesh() );
416 
417  printParameterHdr( ostr, model->parameterSpace()->dimension(), hdrs[M_mode] );
418 
419  int crb_error_type = option(_name="crb.error-type").template as<int>();
420 
421  int dim=0;
422  if( M_mode==CRBModelMode::CRB )
423  {
424  dim=crb->dimension();
425  if( crb->useWNmu() )
426  Sampling = crb->wnmu();
427 
428  if( option(_name="crb.run-on-scm-parameters").template as<bool>() )
429  {
430  Sampling = crb->scm()->c();
431  if( crb_error_type!=1 )
432  throw std::logic_error( "[OpusApp] The SCM has not been launched, you can't use the option crb.run-on-scm-parameters. Run the SCM ( option crb.error-type=1 ) or comment this option line." );
433  }
434  }
435  if( M_mode==CRBModelMode::SCM )
436  {
437  dim=crb->scm()->KMax();
438  if( option(_name="crb.scm.run-on-C").template as<bool>() )
439  Sampling = crb->scm()->c();
440  }
441 
442  std::ofstream file_summary_of_simulations( ( boost::format( "summary_of_simulations_%d" ) %dim ).str().c_str() ,std::ios::out | std::ios::app );
443 
444  int curpar = 0;
445 
446  /* Example of use of the setElements (but can use write in the file SamplingForTest)
447  vector_parameter_type V;
448  parameter_type UserMu( model->parameterSpace() );
449  double j=0.1;
450  //for(int i=1; i<101; i++) { UserMu(0)=j; UserMu(1)=1; UserMu(2)=1.5; UserMu(3)=2; UserMu(4)=3; UserMu(5)=4; UserMu(6)=4.5; UserMu(7)=5; UserMu(8)=6; V.push_back(UserMu ); j+=0.1;}
451  //for(int i=10; i<100; i+=10) { UserMu(0)=i; UserMu(1)=1; V.push_back(UserMu );}
452  //for(int i=1e2; i<1e3; i+=1e2) { UserMu(0)=i; UserMu(1)=1; V.push_back(UserMu );}
453  //for(int i=1e3; i<1e4; i+=1e3) { UserMu(0)=i; UserMu(1)=1; V.push_back(UserMu );}
454  //for(int i=1e4; i<1e5; i+=1e4) { UserMu(0)=i; UserMu(1)=1; V.push_back(UserMu );}
455  //for(int i=1e5; i<1e6; i+=1e5) { UserMu(0)=i; UserMu(1)=1; V.push_back(UserMu );}
456  //UserMu(0)=1e6; UserMu(1)=1; V.push_back(UserMu );
457  //Sampling->setElements( V );
458  */
459 
460  // Script write current mu in cfg => need to write in SamplingForTest
461  if( option(_name="crb.script-mode").template as<bool>() )
462  {
463  // Sampling will be the parameter given by OT
464  if( proc_number == Environment::worldComm().masterRank() )
465  buildSamplingFromCfg();
466 
467  Environment::worldComm().barrier();
468 
469  compute_fem = false;
470  compute_stat = false;
471  }
472 
478  if( option(_name="crb.use-predefined-test-sampling").template as<bool>() || option(_name="crb.script-mode").template as<bool>() )
479  {
480  std::string file_name = ( boost::format("SamplingForTest") ).str();
481  std::ifstream file ( file_name );
482  if( file )
483  {
484  Sampling->clear();
485  Sampling->readFromFile( file_name ) ;
486  }
487  else
488  {
489  VLOG(2) << "proc number : " << proc_number << "can't find file \n";
490  throw std::logic_error( "[OpusApp] file SamplingForTest was not found" );
491  }
492  }
493 
494  //Statistics
495  vectorN_type l2_error_vector;
496  vectorN_type h1_error_vector;
497  vectorN_type relative_error_vector;
498  vectorN_type time_fem_vector;
499  vectorN_type time_crb_vector;
500  vectorN_type relative_estimated_error_vector;
501 
502  vectorN_type scm_relative_error;
503 
504  bool solve_dual_problem = option(_name="crb.solve-dual-problem").template as<bool>();
505 
506  if (option(_name="crb.cvg-study").template as<bool>() && compute_fem )
507  {
508  int Nmax = crb->dimension();
509  vector_sampling_for_primal_efficiency_under_1.resize(Nmax);
510  for(int N=0; N<Nmax; N++)
511  {
512  sampling_ptrtype sampling_primal ( new typename crb_type::sampling_type( model->parameterSpace() ) );
513  vector_sampling_for_primal_efficiency_under_1[N]=sampling_primal;
514  }
515  if( solve_dual_problem )
516  {
517  vector_sampling_for_dual_efficiency_under_1.resize(Nmax);
518  for(int N=0; N<Nmax; N++)
519  {
520  sampling_ptrtype sampling_dual ( new typename crb_type::sampling_type( model->parameterSpace() ) );
521  vector_sampling_for_dual_efficiency_under_1[N]=sampling_dual;
522  }
523  }
524  }
525 
526  if( M_mode==CRBModelMode::CRB )
527  {
528 
529  l2_error_vector.resize( Sampling->size() );
530  h1_error_vector.resize( Sampling->size() );
531  relative_error_vector.resize( Sampling->size() );
532  time_fem_vector.resize( Sampling->size() );
533  time_crb_vector.resize( Sampling->size() );
534 
535  if( crb->errorType()!=2 )
536  relative_estimated_error_vector.resize( Sampling->size() );
537 
538  crb->setOfflineStep( false );
539 
540  if (option(_name="eim.cvg-study").template as<bool>() )
541  {
542  this->initializeConvergenceEimMap( Sampling->size() );
543  compute_fem=false;
544  }
545  if (option(_name="crb.cvg-study").template as<bool>() )
546  this->initializeConvergenceCrbMap( Sampling->size() );
547 
548  }
549  if( M_mode==CRBModelMode::SCM )
550  {
551  if (option(_name="crb.scm.cvg-study").template as<bool>() )
552  this->initializeConvergenceScmMap( Sampling->size() );
553 
554  scm_relative_error.resize( Sampling->size() );
555  }
556 
557  int crb_dimension = option(_name="crb.dimension").template as<int>();
558  int crb_dimension_max = option(_name="crb.dimension-max").template as<int>();
559  double crb_online_tolerance = option(_name="crb.online-tolerance").template as<double>();
560  bool crb_compute_variance = option(_name="crb.compute-variance").template as<bool>();
561 
562  double output_fem = -1;
563 
564 
565  //in the case we don't do the offline step, we need the affine decomposition
566  model->computeAffineDecomposition();
567 
568  //compute beta coeff for reference parameters
569  auto ref_mu = model->refParameter();
570  double dt = model->timeStep();
571  double ti = model->timeInitial();
572  double tf = model->timeFinal();
573  int K = ( tf - ti )/dt;
574  std::vector< std::vector< std::vector< double > > > ref_betaAqm;
575  for(int time_index=0; time_index<K; time_index++)
576  {
577  double time = time_index*dt;
578  ref_betaAqm.push_back( model->computeBetaQm( ref_mu , time ).template get<1>() );
579  }
580  auto ref_betaMqm = model->computeBetaQm( ref_mu , tf ).template get<0>() ;
581 
582 
583  BOOST_FOREACH( auto mu, *Sampling )
584  {
585  int size = mu.size();
586 
587  element_type u_crb; // expansion of reduced solution
588  element_type u_crb_dual; // expansion of reduced solution ( dual )
589 #if !NDEBUG
590  if( proc_number == Environment::worldComm().masterRank() )
591  {
592  std::cout << "(" << curpar << "/" << Sampling->size() << ") mu = [ ";
593  for ( int i=0; i<size-1; i++ ) std::cout<< mu[i] <<" , ";
594  std::cout<< mu[size-1]<<" ]\n ";
595  }
596 #endif
597  curpar++;
598 
599  std::ostringstream mu_str;
600  //if too many parameters, it will crash
601  int sizemax=7;
602  if( size < sizemax )
603  sizemax=size;
604  for ( int i=0; i<sizemax-1; i++ ) mu_str << std::scientific << std::setprecision( 5 ) << mu[i] <<",";
605  mu_str << std::scientific << std::setprecision( 5 ) << mu[size-1];
606 
607 #if !NDEBUG
608  LOG(INFO) << "mu=" << mu << "\n";
609  mu.check();
610 #endif
611  if( option(_name="crb.script-mode").template as<bool>() )
612  {
613  unsigned long N = mu.size() + 5;
614  unsigned long P = 2;
615  std::vector<double> X( N );
616  std::vector<double> Y( P );
617  for(int i=0; i<mu.size(); i++)
618  X[i] = mu[i];
619 
620  int N_dim = crb_dimension;
621  int N_dimMax = crb_dimension_max;
622  int Nwn;
623  if( N_dim > 0 )
624  Nwn = N_dim;
625  else
626  Nwn = N_dimMax;
627  X[N-5] = output_index;
628  X[N-4] = Nwn;
629  X[N-3] = crb_online_tolerance;
630  X[N-2] = crb_error_type;
631  //X[N-1] = option(_name="crb.compute-variance").template as<int>();
632  X[N-1] = 0;
633  bool compute_variance = crb_compute_variance;
634  if ( compute_variance )
635  X[N-1] = 1;
636 
637  this->run( X.data(), X.size(), Y.data(), Y.size() );
638  //std::cout << "output = " << Y[0] << std::endl;
639 
640  std::ofstream res(option(_name="result-file").template as<std::string>() );
641  res << "output="<< Y[0] << "\n";
642  }
643  else
644  {
645  switch ( M_mode )
646  {
647  case CRBModelMode::PFEM:
648  {
649  LOG(INFO) << "PFEM mode" << std::endl;
650  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
651  std::cout << "PFEM mode" << std::endl;
652  boost::mpi::timer ti;
653 
654  model->computeAffineDecomposition();
655  auto u_fem = model->solveFemUsingAffineDecompositionFixedPoint( mu );
656  std::ostringstream u_fem_str;
657  u_fem_str << "u_fem(" << mu_str.str() << ")";
658  u_fem.setName( u_fem_str.str() );
659 
660  LOG(INFO) << "compute output\n";
661  if( export_solution )
662  exporter->step(0)->add( u_fem.name(), u_fem );
663  //model->solve( mu );
664  std::vector<double> o = boost::assign::list_of( model->output( output_index,mu , u_fem, true) )( ti.elapsed() );
665  if(proc_number == Environment::worldComm().masterRank() ) std::cout << "output=" << o[0] << "\n";
666  printEntry( ostr, mu, o );
667 
668  std::ofstream res(option(_name="result-file").template as<std::string>() );
669  res << "output="<< o[0] << "\n";
670 
671  }
672  break;
673 
674  case CRBModelMode::CRB:
675  {
676  LOG(INFO) << "CRB mode\n";
677  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
678  std::cout << "CRB mode\n";
679 
680 
681  boost::mpi::timer ti;
682 
683  ti.restart();
684  LOG(INFO) << "solve crb\n";
685  google::FlushLogFiles(google::GLOG_INFO);
686 
687  //dimension of the RB (not necessarily the max)
688  int N = option(_name="crb.dimension").template as<int>();
689 
690  auto o = crb->run( mu, option(_name="crb.online-tolerance").template as<double>() , N);
691  double time_crb = ti.elapsed();
692 
693  auto WN = crb->wn();
694  auto WNdu = crb->wndu();
695  //auto u_crb = crb->expansion( mu , N );
696  auto solutions=o.template get<2>();
697  auto uN = solutions.template get<0>();
698  auto uNdu = solutions.template get<1>();
699 
700  int size = uN.size();
701 
702  //if( model->isSteady()) // Re-use uN given by lb in crb->run
703  u_crb = crb->expansion( uN[size-1] , N , WN ); // Re-use uN given by lb in crb->run
704  if( solve_dual_problem )
705  u_crb_dual = crb->expansion( uNdu[0] , N , WNdu );
706  //else
707  // u_crb = crb->expansion( mu , N , WN );
708 
709  std::ostringstream u_crb_str;
710  u_crb_str << "u_crb(" << mu_str.str() << ")";
711  u_crb.setName( u_crb_str.str() );
712  LOG(INFO) << "export u_crb \n";
713  if( export_solution )
714  exporter->step(0)->add( u_crb.name(), u_crb );
715 
716  double relative_error = -1;
717  double relative_estimated_error = -1;
718  double condition_number = o.template get<3>();
719  double l2_error = -1;
720  double h1_error = -1;
721  double l2_dual_error = -1;
722  double h1_dual_error = -1;
723  double time_fem = -1;
724 
725  element_type u_fem ;
726  element_type u_dual_fem ;
727 
728  auto all_upper_bounds = o.template get<6>();
729  double output_estimated_error = all_upper_bounds.template get<0>();
730  double solution_estimated_error = all_upper_bounds.template get<1>();
731  double dual_solution_estimated_error = all_upper_bounds.template get<2>();
732 
733  if ( compute_fem )
734  {
735  bool use_newton = option(_name="crb.use-newton").template as<bool>();
736 
737  ti.restart();
738  LOG(INFO) << "solve u_fem\n";
739 
740  //auto u_fem = model->solveRB( mu );
741  //auto u_fem = model->solveFemUsingOfflineEim( mu );
742 
743  if( boost::is_same< crbmodel_type , crbmodelbilinear_type >::value && ! use_newton )
744  //use affine decomposition
745  u_fem = model->solveFemUsingAffineDecompositionFixedPoint( mu );
746  else
747  u_fem = model->solve( mu );
748 
749  std::ostringstream u_fem_str;
750  u_fem_str << "u_fem(" << mu_str.str() << ")";
751  u_fem.setName( u_fem_str.str() );
752 
753  if( export_solution )
754  {
755  LOG(INFO) << "export u_fem \n";
756  exporter->step(0)->add( u_fem.name(), u_fem );
757  }
758  std::vector<double> ofem = boost::assign::list_of( model->output( output_index,mu, u_fem ) )( ti.elapsed() );
759 
760  double ocrb = o.template get<0>();
761  relative_error = std::abs( ofem[0]- ocrb) /ofem[0];
762  relative_estimated_error = output_estimated_error / ofem[0];
763 
764  //compute || u_fem - u_crb||_L2
765  LOG(INFO) << "compute error \n";
766  auto u_error = model->functionSpace()->element();
767  auto u_dual_error = model->functionSpace()->element();
768  std::ostringstream u_error_str;
769  u_error = (( u_fem - u_crb ).pow(2)).sqrt() ;
770  u_error_str << "u_error(" << mu_str.str() << ")";
771  u_error.setName( u_error_str.str() );
772  if( export_solution )
773  exporter->step(0)->add( u_error.name(), u_error );
774  LOG(INFO) << "L2(fem)=" << l2Norm( u_fem ) << "\n";
775  LOG(INFO) << "H1(fem)=" << h1Norm( u_fem ) << "\n";
776  l2_error = l2Norm( u_error )/l2Norm( u_fem );
777  h1_error = h1Norm( u_error )/h1Norm( u_fem );
778 
779  output_fem = ofem[0];
780  time_fem = ofem[1];
781 
782  if( boost::is_same< crbmodel_type , crbmodelbilinear_type >::value && ! use_newton )
783  {
784  if( solve_dual_problem )
785  {
786  u_dual_fem = model->solveFemDualUsingAffineDecompositionFixedPoint( mu );
787 
788  u_dual_error = model->functionSpace()->element();
789  u_dual_error = (( u_dual_fem - u_crb_dual ).pow(2)).sqrt() ;
790  l2_dual_error = l2Norm( u_dual_error )/l2Norm( u_dual_fem );
791  h1_dual_error = h1Norm( u_dual_error )/h1Norm( u_dual_fem );
792  }
793  }
794 
795  }//compute-fem-during-online
796 
797  if ( crb->errorType()==2 )
798  {
799  double ocrb = o.template get<0>();
800  std::vector<double> v = boost::assign::list_of( output_fem )( time_fem )( ocrb )( relative_estimated_error )( time_crb )( relative_error )( condition_number )( l2_error )( h1_error );
801  if( proc_number == Environment::worldComm().masterRank() )
802  {
803  std::cout << "output=" << ocrb << " with " << o.template get<1>() << " basis functions\n";
804  printEntry( file_summary_of_simulations, mu, v );
805  printEntry( ostr, mu, v );
806  //file_summary_of_simulations.close();
807 
808  if ( option(_name="crb.compute-stat").template as<bool>() && compute_fem )
809  {
810  relative_error_vector[curpar-1] = relative_error;
811  l2_error_vector[curpar-1] = l2_error;
812  h1_error_vector[curpar-1] = h1_error;
813  time_fem_vector[curpar-1] = time_fem;
814  time_crb_vector[curpar-1] = time_crb;
815  }
816 
817  std::ofstream res(option(_name="result-file").template as<std::string>() );
818  res << "output="<< ocrb << "\n";
819 
820  }
821 
822  }//end of crb->errorType==2
823  else
824  {
825  //if( ! boost::is_same< crbmodel_type , crbmodelbilinear_type >::value )
826  // throw std::logic_error( "ERROR TYPE must be 2 when using CRBTrilinear (no error estimation)" );
827  double ocrb = o.template get<0>();
828  std::vector<double> v = boost::assign::list_of( output_fem )( time_fem )( ocrb )( relative_estimated_error )( ti.elapsed() ) ( relative_error )( condition_number )( l2_error )( h1_error ) ;
829  if( proc_number == Environment::worldComm().masterRank() )
830  {
831  std::cout << "output=" << ocrb << " with " << o.template get<1>() << " basis functions (relative error estimation on this output : " << relative_estimated_error<<") \n";
832  //std::ofstream file_summary_of_simulations( ( boost::format( "summary_of_simulations_%d" ) % o.template get<2>() ).str().c_str() ,std::ios::out | std::ios::app );
833  printEntry( file_summary_of_simulations, mu, v );
834  printEntry( ostr, mu, v );
835  //file_summary_of_simulations.close();
836 
837  if ( option(_name="crb.compute-stat").template as<bool>() && compute_fem )
838  {
839  relative_error_vector[curpar-1] = relative_error;
840  l2_error_vector[curpar-1] = l2_error;
841  h1_error_vector[curpar-1] = h1_error;
842  time_fem_vector[curpar-1] = time_fem;
843  time_crb_vector[curpar-1] = time_crb;
844  relative_estimated_error_vector[curpar-1] = relative_estimated_error;
845  }
846  std::ofstream res(option(_name="result-file").template as<std::string>() );
847  res << "output="<< ocrb << "\n";
848  }//end of proc==master
849  }//end of else (errorType==2)
850 
851  if (option(_name="eim.cvg-study").template as<bool>() )
852  {
853  bool check_name = false;
854  std::string how_compute_unknown = option(_name=_o( this->about().appName(),"how-compute-unkown-for-eim" )).template as<std::string>();
855  if( how_compute_unknown == "CRB-with-ad")
856  {
857  LOG( INFO ) << "convergence eim with CRB-with-ad ";
858  check_name = true;
859  this->studyEimConvergence( mu , u_crb , curpar );
860  }
861  if( how_compute_unknown == "FEM-with-ad")
862  {
863  LOG( INFO ) << "convergence eim with FEM-with-ad ";
864  check_name = true;
865  //fem computed via solveFemUsingAffineDecomposition use the affine decomposition
866  auto fem_with_ad = model->solveFemUsingAffineDecompositionFixedPoint( mu );
867  this->studyEimConvergence( mu , fem_with_ad , curpar );
868  }
869  if( how_compute_unknown == "FEM-without-ad")
870  {
871  LOG( INFO ) << "convergence eim with FEM-without-ad ";
872  check_name = true;
873  auto fem_without_ad = model->solve( mu );
874  this->studyEimConvergence( mu , fem_without_ad , curpar );
875  }
876  if( ! check_name )
877  throw std::logic_error( "OpusApp error with option how-compute-unknown-for-eim, please use CRB-with-ad, FEM-with-ad or FEM-without-ad" );
878  }
879 
880  if (option(_name="crb.cvg-study").template as<bool>() && compute_fem )
881  {
882 
883  LOG(INFO) << "start convergence study...\n";
884  std::map<int, boost::tuple<double,double,double,double,double,double,double> > conver;
885 
886  int Nmax = crb->dimension();
887  for( int N = 1; N <= Nmax ; N++ )
888  {
889  auto o= crb->run( mu, option(_name="crb.online-tolerance").template as<double>() , N);
890  auto ocrb = o.template get<0>();
891  auto solutions=o.template get<2>();
892  auto u_crb = solutions.template get<0>();
893  auto u_crb_du = solutions.template get<1>();
894  int size = u_crb.size();
895  auto uN = crb->expansion( u_crb[size-1], N, WN );
896 
897  element_type uNdu;
898 
899  auto u_error = u_fem - uN;
900  auto u_dual_error = model->functionSpace()->element();
901  if( solve_dual_problem )
902  {
903  uNdu = crb->expansion( u_crb_du[0], N, WNdu );
904  u_dual_error = u_dual_fem - uNdu;
905  }
906 
907  auto all_upper_bounds = o.template get<6>();
908  output_estimated_error = all_upper_bounds.template get<0>();
909  solution_estimated_error = all_upper_bounds.template get<1>();
910  dual_solution_estimated_error = all_upper_bounds.template get<2>();
911 
912  //auto o = crb->run( mu, option(_name="crb.online-tolerance").template as<double>() , N);
913  double rel_err = std::abs( output_fem-ocrb ) /output_fem;
914 
915  double output_relative_estimated_error = output_estimated_error / output_fem;
916 
917  double primal_residual_norm = o.template get<4>();
918  double dual_residual_norm = o.template get<5>();
919 
920  double solution_error=0;
921  double dual_solution_error=0;
922  double ref_primal=0;
923  double ref_dual=0;
924  if( model->isSteady() )
925  {
926  //let ufem-ucrb = e
927  //||| e |||_mu = sqrt( a( e , e ; mu ) ) = solution_error
928  for(int q=0; q<model->Qa();q++)
929  {
930  for(int m=0; m<model->mMaxA(q); m++)
931  {
932  solution_error += ref_betaAqm[0][q][m]*model->Aqm(q,m,u_error,u_error) ;
933  ref_primal += ref_betaAqm[0][q][m]*model->Aqm(q,m,u_fem,u_fem);
934  }
935  }
936 
937  if( solve_dual_problem )
938  {
939  for(int q=0; q<model->Qa();q++)
940  {
941  for(int m=0; m<model->mMaxA(q); m++)
942  {
943  dual_solution_error += ref_betaAqm[0][q][m]*model->Aqm(q,m,u_dual_error,u_dual_error);
944  ref_dual += ref_betaAqm[0][q][m]*model->Aqm(q,m,u_dual_fem,u_dual_fem);
945  }
946  }
947  dual_solution_error = math::sqrt( dual_solution_error );
948  ref_dual = math::sqrt( ref_dual );
949  }
950 
951  solution_error = math::sqrt( solution_error );
952  ref_primal = math::sqrt( ref_primal );
953  //dual_solution_error = math::sqrt( model->scalarProduct( u_dual_error, u_dual_error ) );
954 
955  }
956  else
957  {
958  double dt = model->timeStep();
959  double ti = model->timeInitial();
960  double tf = model->timeFinal();
961  int K = ( tf - ti )/dt;
962 
963  for(int q=0; q<model->Qm();q++)
964  {
965  for(int m=0; m<model->mMaxM(q); m++)
966  {
967  solution_error += ref_betaMqm[q][m]*model->Mqm(q,m,u_error,u_error);
968  ref_primal += ref_betaMqm[q][m]*model->Mqm(q,m,u_fem,u_fem);
969  }
970  }
971  for(int time_index=0; time_index<K; time_index++)
972  {
973  double t=time_index*dt;
974  for(int q=0; q<model->Qa();q++)
975  {
976  for(int m=0; m<model->mMaxA(q); m++)
977  {
978  solution_error += ref_betaAqm[time_index][q][m]*model->Aqm(q,m,u_error,u_error) * dt;
979  ref_primal += ref_betaAqm[time_index][q][m]*model->Aqm(q,m,u_fem,u_fem) * dt;
980  }
981  }
982  }
983  solution_error = math::sqrt( solution_error );
984  ref_primal = math::sqrt( ref_primal );
985 
986  if( solve_dual_problem )
987  {
988  ti = model->timeFinal()+dt;
989  tf = model->timeInitial()+dt;
990  dt -= dt;
991 
992  for(int q=0; q<model->Qm();q++)
993  {
994  for(int m=0; m<model->mMaxM(q); m++)
995  {
996  dual_solution_error += ref_betaMqm[q][m]*model->Mqm(q,m,u_dual_error,u_dual_error);
997  ref_dual += ref_betaMqm[q][m]*model->Mqm(q,m,u_dual_fem,u_dual_fem);
998  }
999  }
1000 
1001  for(int time_index=0; time_index<K; time_index++)
1002  {
1003  double t=time_index*dt;
1004  for(int q=0; q<model->Qa();q++)
1005  {
1006  for(int m=0; m<model->mMaxA(q); m++)
1007  {
1008  dual_solution_error += ref_betaAqm[time_index][q][m]*model->Aqm(q,m,u_dual_error,u_dual_error) * dt;
1009  ref_dual += ref_betaAqm[time_index][q][m]*model->Aqm(q,m,u_dual_fem,u_dual_fem) * dt;
1010  }
1011  }
1012  }
1013 
1014  dual_solution_error = math::sqrt( dual_solution_error );
1015  ref_dual = math::sqrt( ref_dual );
1016 
1017  }//if solve-dual
1018 
1019  }//transient case
1020 
1021  double l2_error = l2Norm( u_error )/l2Norm( u_fem );
1022  double h1_error = h1Norm( u_error )/h1Norm( u_fem );
1023  double condition_number = o.template get<3>();
1024  double output_error_bound_efficiency = output_relative_estimated_error / rel_err;
1025 
1026  double relative_primal_solution_error = solution_error / ref_primal ;
1027  double relative_primal_solution_estimated_error = solution_estimated_error / ref_primal;
1028  double relative_primal_solution_error_bound_efficiency = relative_primal_solution_estimated_error / relative_primal_solution_error;
1029 
1030  if( relative_primal_solution_error_bound_efficiency < 1 )
1031  {
1032  vector_sampling_for_primal_efficiency_under_1[N-1]->push_back( mu , 1);
1033  }
1034 
1035  double relative_dual_solution_error = 1;
1036  double relative_dual_solution_estimated_error = 1;
1037  double relative_dual_solution_error_bound_efficiency = 1;
1038  if( solve_dual_problem )
1039  {
1040  relative_dual_solution_error = dual_solution_error / ref_dual ;
1041  relative_dual_solution_estimated_error = dual_solution_estimated_error / ref_dual;
1042  relative_dual_solution_error_bound_efficiency = relative_dual_solution_estimated_error / relative_dual_solution_error;
1043 
1044  if( relative_dual_solution_error_bound_efficiency < 1 )
1045  {
1046  vector_sampling_for_dual_efficiency_under_1[N-1]->push_back( mu , 0);
1047  }
1048 
1049  }
1050  conver[N]=boost::make_tuple( rel_err, l2_error, h1_error , relative_estimated_error, condition_number , output_error_bound_efficiency , relative_primal_solution_error_bound_efficiency );
1051 
1052  LOG(INFO) << "N=" << N << " " << rel_err << " " << l2_error << " " << h1_error << " " <<condition_number<<"\n";
1053  if ( proc_number == Environment::worldComm().masterRank() )
1054  std::cout << "N=" << N << " " << rel_err << " " << l2_error << " " << h1_error << " " <<relative_estimated_error<<" "<<condition_number<<std::endl;
1055  M_mapConvCRB["L2"][N-1](curpar - 1) = l2_error;
1056  M_mapConvCRB["H1"][N-1](curpar - 1) = h1_error;
1057  M_mapConvCRB["OutputError"][N-1](curpar - 1) = rel_err;
1058  M_mapConvCRB["OutputEstimatedError"][N-1](curpar - 1) = output_relative_estimated_error;
1059  M_mapConvCRB["OutputErrorBoundEfficiency"][N-1](curpar - 1) = output_error_bound_efficiency;
1060  M_mapConvCRB["SolutionErrorBoundEfficiency"][N-1](curpar - 1) = relative_primal_solution_error_bound_efficiency;
1061  M_mapConvCRB["SolutionError"][N-1](curpar - 1) = relative_primal_solution_error;
1062  M_mapConvCRB["SolutionErrorEstimated"][N-1](curpar - 1) = relative_primal_solution_estimated_error;
1063  M_mapConvCRB["SolutionDualErrorBoundEfficiency"][N-1](curpar - 1) = relative_dual_solution_error_bound_efficiency;
1064  M_mapConvCRB["SolutionDualError"][N-1](curpar - 1) = relative_dual_solution_error;
1065  M_mapConvCRB["SolutionDualErrorEstimated"][N-1](curpar - 1) = relative_dual_solution_estimated_error;
1066  M_mapConvCRB["PrimalResidualNorm"][N-1](curpar - 1) = primal_residual_norm;
1067  M_mapConvCRB["DualResidualNorm"][N-1](curpar - 1) = dual_residual_norm;
1068  LOG(INFO) << "N=" << N << " done.\n";
1069 
1070  if( relative_primal_solution_error_bound_efficiency < 1 )
1071  LOG( INFO ) << "efficiency of error estimation on primal solution is "<<relative_primal_solution_error_bound_efficiency<<" ( should be >= 1 )";
1072  if( relative_dual_solution_error_bound_efficiency < 1 )
1073  LOG( INFO ) << "efficiency of error estimation on dual solution is "<<relative_dual_solution_error_bound_efficiency<<" ( should be >= 1 )";
1074  if( output_error_bound_efficiency < 1 )
1075  LOG( INFO ) << "efficiency of error estimation on output is "<<output_error_bound_efficiency<<" ( should be >= 1 )";
1076 
1077  if ( option(_name="crb.check.residual").template as<bool>() && solve_dual_problem )
1078  {
1079  std::vector < std::vector<double> > primal_residual_coefficients = all_upper_bounds.template get<3>();
1080  std::vector < std::vector<double> > dual_residual_coefficients = all_upper_bounds.template get<4>();
1081  if( model->isSteady() )
1082  {
1083  crb->checkResidual( mu , primal_residual_coefficients, dual_residual_coefficients, uN, uNdu );
1084  }
1085  else
1086  {
1087  std::vector< element_type > uNelement;
1088  std::vector< element_type > uNelement_old;
1089  std::vector< element_type > uNelement_du;
1090  std::vector< element_type > uNelement_du_old;
1091 
1092  auto u_crb_old = solutions.template get<2>();
1093  auto u_crb_du_old = solutions.template get<3>();
1094 
1095  //size is the number of time step
1096  for(int t=0; t<size; t++)
1097  {
1098  uNelement.push_back( crb->expansion( u_crb[t], N, WN ) );
1099  uNelement_old.push_back( crb->expansion( u_crb_old[t], N, WN ) );
1100  uNelement_du.push_back( crb->expansion( u_crb_du[t], N, WNdu ) );
1101  uNelement_du_old.push_back( crb->expansion( u_crb_du_old[t], N, WNdu ) );
1102  }//loop over time step
1103 
1104  crb->compareResidualsForTransientProblems(N, mu ,
1105  uNelement, uNelement_old, uNelement_du, uNelement_du_old,
1106  primal_residual_coefficients, dual_residual_coefficients );
1107  }//transient case
1108  }//check residuals computations
1109 
1110  }//loop over basis functions ( N )
1111  if( proc_number == Environment::worldComm().masterRank() )
1112  {
1113  LOG(INFO) << "save in logfile\n";
1114  std::string mu_str;
1115  for ( int i=0; i<mu.size(); i++ )
1116  mu_str= mu_str + ( boost::format( "_%1%" ) %mu[i] ).str() ;
1117  std::string file_name = "convergence"+mu_str+".dat";
1118  std::ofstream conv( file_name );
1119  BOOST_FOREACH( auto en, conver )
1120  conv << en.first << "\t" << en.second.get<0>() << "\t" << en.second.get<1>() << "\t" << en.second.get<2>() <<
1121  "\t"<< en.second.get<3>() << "\t"<< en.second.get<4>()<< "\t" <<en.second.get<5>()<< "\n";
1122  }
1123  }//end of cvg-study
1124  }//case CRB
1125  break;
1126  case CRBModelMode::CRB_ONLINE:
1127  {
1128  std::cout << "CRB Online mode\n";
1129  boost::mpi::timer ti;
1130  ti.restart();
1131  auto o = crb->run( mu, option(_name="crb.online-tolerance").template as<double>() );
1132 
1133  if ( crb->errorType()==2 )
1134  {
1135  std::vector<double> v = boost::assign::list_of( o.template get<0>() )( ti.elapsed() );
1136  std::cout << "output=" << o.template get<0>() << " with " << o.template get<1>() << " basis functions\n";
1137  printEntry( ostr, mu, v );
1138  }
1139 
1140  else
1141  {
1142  auto all_upper_bounds = o.template get<6>();
1143  double output_estimated_error = all_upper_bounds.template get<0>();
1144  double relative_estimated_error = output_estimated_error / output_fem;
1145  std::vector<double> v = boost::assign::list_of( o.template get<0>() )( output_estimated_error )( ti.elapsed() );
1146  std::cout << "output=" << o.template get<0>() << " with " << o.template get<1>() <<
1147  " basis functions (relative error estimation on this output : " << relative_estimated_error<<") \n";
1148  printEntry( ostr, mu, v );
1149  }
1150 
1151  }
1152  break;
1153 
1154  case CRBModelMode::SCM:
1155  {
1156 
1157  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1158  std::cout << "SCM mode\n";
1159  int kmax = crb->scm()->KMax();
1160  auto o = crb->scm()->run( mu, kmax );
1161  printEntry( file_summary_of_simulations, mu, o );
1162  printEntry( ostr, mu, o );
1163 
1164  scm_relative_error[curpar-1] = o[6];
1165 
1166  if (option(_name="crb.scm.cvg-study").template as<bool>() )
1167  {
1168  LOG(INFO) << "start scm convergence study...\n";
1169  std::map<int, boost::tuple<double> > conver;
1170  for( int N = 1; N <= kmax; N++ )
1171  {
1172  auto o = crb->scm()->run( mu, N);
1173  double relative_error = o[6];
1174  conver[N]=boost::make_tuple( relative_error );
1175  if ( proc_number == Environment::worldComm().masterRank() )
1176  std::cout << "N=" << N << " " << relative_error <<std::endl;
1177  M_mapConvSCM["RelativeError"][N-1](curpar - 1) = relative_error;
1178  }
1179  if( proc_number == Environment::worldComm().masterRank() )
1180  {
1181  LOG(INFO) << "save in logfile\n";
1182  std::string mu_str;
1183  for ( int i=0; i<mu.size(); i++ )
1184  mu_str= mu_str + ( boost::format( "_%1%" ) %mu[i] ).str() ;
1185  std::string file_name = "convergence-scm-"+mu_str+".dat";
1186  std::ofstream conv( file_name );
1187  BOOST_FOREACH( auto en, conver )
1188  conv << en.first << "\t" << en.second.get<0>() ;
1189  }
1190  }//end of cvg-study
1191 
1192  }
1193  break;
1194 
1195  case CRBModelMode::SCM_ONLINE:
1196  {
1197  std::cout << "SCM online mode\n";
1198  auto o = crb->scm()->run( mu, crb->scm()->KMax() );
1199  printEntry( ostr, mu, o );
1200  }
1201  break;
1202 
1203  }
1204 
1205  LOG( INFO ) << "------------------------------------------------------------";
1206  }
1207  }
1208 
1209  //model->computationalTimeEimStatistics();
1210  if( export_solution )
1211  exporter->save();
1212 
1213  if( proc_number == Environment::worldComm().masterRank() ) std::cout << ostr.str() << "\n";
1214 
1215  if (option(_name="eim.cvg-study").template as<bool>() && M_mode==CRBModelMode::CRB)
1216  this->doTheEimConvergenceStat( Sampling->size() );
1217 
1218  if (option(_name="crb.cvg-study").template as<bool>() && compute_fem && M_mode==CRBModelMode::CRB )
1219  this->doTheCrbConvergenceStat( Sampling->size() );
1220 
1221  if (option(_name="crb.scm.cvg-study").template as<bool>() && M_mode==CRBModelMode::SCM )
1222  this->doTheScmConvergenceStat( Sampling->size() );
1223 
1224  if ( compute_stat && compute_fem && M_mode==CRBModelMode::CRB )
1225  {
1226  LOG( INFO ) << "compute statistics \n";
1227  Eigen::MatrixXf::Index index_max_l2;
1228  Eigen::MatrixXf::Index index_min_l2;
1229  Eigen::MatrixXf::Index index_max_h1;
1230  Eigen::MatrixXf::Index index_min_h1;
1231  Eigen::MatrixXf::Index index_max_time_fem;
1232  Eigen::MatrixXf::Index index_min_time_fem;
1233  Eigen::MatrixXf::Index index_max_time_crb;
1234  Eigen::MatrixXf::Index index_min_time_crb;
1235  Eigen::MatrixXf::Index index_max_estimated_error;
1236  Eigen::MatrixXf::Index index_min_estimated_error;
1237  Eigen::MatrixXf::Index index_max_output_error;
1238  Eigen::MatrixXf::Index index_min_output_error;
1239 
1240  double max_l2 = l2_error_vector.maxCoeff(&index_max_l2);
1241  double min_l2 = l2_error_vector.minCoeff(&index_min_l2);
1242  double mean_l2 = l2_error_vector.mean();
1243  double max_h1 = h1_error_vector.maxCoeff(&index_max_h1);
1244  double min_h1 = h1_error_vector.minCoeff(&index_min_h1);
1245  double mean_h1 = h1_error_vector.mean();
1246  double max_time_fem = time_fem_vector.maxCoeff(&index_max_time_fem);
1247  double min_time_fem = time_fem_vector.minCoeff(&index_min_time_fem);
1248  double mean_time_fem = time_fem_vector.mean();
1249  double max_time_crb = time_crb_vector.maxCoeff(&index_max_time_crb);
1250  double min_time_crb = time_crb_vector.minCoeff(&index_min_time_crb);
1251  double mean_time_crb = time_crb_vector.mean();
1252  double max_output_error = relative_error_vector.maxCoeff(&index_max_output_error);
1253  double min_output_error = relative_error_vector.minCoeff(&index_min_output_error);
1254  double mean_output_error = relative_error_vector.mean();
1255  double max_estimated_error = 0;
1256  double min_estimated_error = 0;
1257  double mean_estimated_error = 0;
1258 
1259  if( crb->errorType()!=2 )
1260  {
1261  max_estimated_error = relative_estimated_error_vector.maxCoeff(&index_max_estimated_error);
1262  min_estimated_error = relative_estimated_error_vector.minCoeff(&index_min_estimated_error);
1263  mean_estimated_error = relative_estimated_error_vector.mean();
1264  }
1265  if( proc_number == Environment::worldComm().masterRank() )
1266  {
1267  file_summary_of_simulations <<"\n\nStatistics\n";
1268  file_summary_of_simulations <<"max of output error : "<<max_output_error<<" at the "<<index_max_output_error+1<<"^th simulation\n";
1269  file_summary_of_simulations <<"min of output error : "<<min_output_error<<" at the "<<index_min_output_error+1<<"^th simulation\n";
1270  file_summary_of_simulations <<"mean of output error : "<<mean_output_error<<"\n\n";
1271  file_summary_of_simulations <<"max of estimated output error : "<<max_estimated_error<<" at the "<<index_max_estimated_error+1<<"^th simulation\n";
1272  file_summary_of_simulations <<"min of estimated output error : "<<min_estimated_error<<" at the "<<index_min_estimated_error+1<<"^th simulation\n";
1273  file_summary_of_simulations <<"mean of estimated output error : "<<mean_estimated_error<<"\n\n";
1274  file_summary_of_simulations <<"max of L2 error : "<<max_l2<<" at the "<<index_max_l2+1<<"^th simulation\n";
1275  file_summary_of_simulations <<"min of L2 error : "<<min_l2<<" at the "<<index_min_l2+1<<"^th simulation\n";
1276  file_summary_of_simulations <<"mean of L2 error : "<<mean_l2<<"\n\n";
1277  file_summary_of_simulations <<"max of H1 error : "<<max_h1<<" at the "<<index_max_h1+1<<"^th simulation\n";
1278  file_summary_of_simulations <<"min of H1 error : "<<min_h1<<" at the "<<index_min_h1+1<<"^th simulation\n";
1279  file_summary_of_simulations <<"mean of H1 error : "<<mean_h1<<"\n\n";
1280  file_summary_of_simulations <<"max of time FEM : "<<max_time_fem<<" at the "<<index_max_time_fem+1<<"^th simulation\n";
1281  file_summary_of_simulations <<"min of time FEM : "<<min_time_fem<<" at the "<<index_min_time_fem+1<<"^th simulation\n";
1282  file_summary_of_simulations <<"mean of time FEM : "<<mean_time_fem<<"\n\n";
1283  file_summary_of_simulations <<"max of time CRB : "<<max_time_crb<<" at the "<<index_max_time_crb+1<<"^th simulation\n";
1284  file_summary_of_simulations <<"min of time CRB : "<<min_time_crb<<" at the "<<index_min_time_crb+1<<"^th simulation\n";
1285  file_summary_of_simulations <<"mean of time CRB : "<<mean_time_crb<<"\n\n";
1286  }
1287  }//end of compute-stat CRB
1288  if( M_mode==CRBModelMode::SCM )
1289  {
1290  LOG( INFO ) << "compute statistics \n";
1291  Eigen::MatrixXf::Index index_max_error;
1292  Eigen::MatrixXf::Index index_min_error;
1293  double max_error = scm_relative_error.maxCoeff(&index_max_error);
1294  double min_error = scm_relative_error.minCoeff(&index_min_error);
1295  double mean_error = scm_relative_error.mean();
1296  if( proc_number == Environment::worldComm().masterRank() )
1297  {
1298  file_summary_of_simulations <<"\n\nStatistics\n";
1299  file_summary_of_simulations <<"max of relative error : "<<max_error<<" at the "<<index_max_error+1<<"^th simulation\n";
1300  file_summary_of_simulations <<"min of relative error : "<<min_error<<" at the "<<index_min_error+1<<"^th simulation\n";
1301  file_summary_of_simulations <<"mean of relative error : "<<mean_error<<"\n\n";
1302  }
1303 
1304  }//end of compute-stat SCM
1305 
1306 
1307  }
1308  void run( const double * X, unsigned long N,
1309  double * Y, unsigned long P )
1310  {
1311 
1312  switch ( M_mode )
1313  {
1314  case CRBModelMode::PFEM:
1315  {
1316  model->run( X, N, Y, P );
1317  }
1318  break;
1319 
1320  case CRBModelMode::SCM:
1321  case CRBModelMode::SCM_ONLINE:
1322  {
1323  crb->scm()->run( X, N, Y, P );
1324  }
1325  break;
1326 
1327  case CRBModelMode::CRB:
1328  case CRBModelMode::CRB_ONLINE:
1329  {
1330  crb->run( X, N, Y, P );
1331  }
1332  break;
1333  }
1334 
1335  fs::current_path( M_current_path );
1336  }
1337 
1338 private:
1339  int printParameterHdr( std::ostream& os, int N, std::vector<std::string> outputhdrs )
1340  {
1341  for ( int i = 0; i < N; ++i )
1342  {
1343  std::ostringstream s;
1344  s << "mu" << i;
1345  os << hdrmanip( prec+7 ) << s.str();
1346  }
1347 
1348  BOOST_FOREACH( auto output, outputhdrs )
1349  {
1350  os << hdrmanip( 15 ) << output;
1351  }
1352  os << "\n";
1353 
1354  return N*( prec+7 )+outputhdrs.size()*15;
1355  }
1356  void printEntry( std::ostream& os,
1357  typename ModelType::parameter_type const& mu,
1358  std::vector<double> const& outputs )
1359  {
1360  for ( int i = 0; i < mu.size(); ++i )
1361  os << std::right <<std::setw( prec+7 ) << dmanip << mu[i];
1362 
1363  BOOST_FOREACH( auto o, outputs )
1364  {
1365  os << tabmanip( 15 ) << o;
1366  }
1367  os << "\n";
1368  }
1369 
1370 
1371  //double l2Norm( typename ModelType::parameter_type const& mu, int N )
1372  double l2Norm( element_type const& u )
1373  {
1374  static const bool is_composite = functionspace_type::is_composite;
1375  return l2Norm( u, mpl::bool_< is_composite >() );
1376  }
1377  double l2Norm( element_type const& u, mpl::bool_<false> )
1378  {
1379  auto mesh = model->functionSpace()->mesh();
1380  return math::sqrt( integrate( elements(mesh), (vf::idv(u))*(vf::idv(u)) ).evaluate()(0,0) );
1381  }
1382  double l2Norm( element_type const& u, mpl::bool_<true>)
1383  {
1384  auto mesh = model->functionSpace()->mesh();
1385  auto uT = u.template element<1>();
1386  return math::sqrt( integrate( elements(mesh), (vf::idv(uT))*(vf::idv(uT)) ).evaluate()(0,0) );
1387  }
1388  //double h1Norm( typename ModelType::parameter_type const& mu, int N )
1389  double h1Norm( element_type const& u )
1390  {
1391  static const bool is_composite = functionspace_type::is_composite;
1392  return h1Norm( u, mpl::bool_< is_composite >() );
1393  }
1394  double h1Norm( element_type const& u, mpl::bool_<false> )
1395  {
1396  auto mesh = model->functionSpace()->mesh();
1397  double l22 = integrate( elements(mesh), (vf::idv(u))*(vf::idv(u)) ).evaluate()(0,0);
1398  double semih12 = integrate( elements(mesh), (vf::gradv(u))*trans(vf::gradv(u)) ).evaluate()(0,0);
1399  return math::sqrt( l22+semih12 );
1400  }
1401  double h1Norm( element_type const& u, mpl::bool_<true>)
1402  {
1403  auto mesh = model->functionSpace()->mesh();
1404  auto u_femT = u.template element<1>();
1405  double l22 = integrate( elements(mesh), (vf::idv(u_femT))*(vf::idv(u_femT)) ).evaluate()(0,0);
1406  double semih12 = integrate( elements(mesh), (vf::gradv(u_femT))*trans(vf::gradv(u_femT))).evaluate()(0,0);
1407  return math::sqrt( l22+semih12 );
1408  }
1409 
1410 
1411  void initializeConvergenceEimMap( int sampling_size )
1412  {
1413  auto eim_sc_vector = model->scalarContinuousEim();
1414  auto eim_sd_vector = model->scalarDiscontinuousEim();
1415 
1416  for(int i=0; i<eim_sc_vector.size(); i++)
1417  {
1418  auto eim = eim_sc_vector[i];
1419  M_mapConvEIM[eim->name()] = std::vector<vectorN_type>(eim->mMax());
1420  for(int j=0; j<eim->mMax(); j++)
1421  M_mapConvEIM[eim->name()][j].resize(sampling_size);
1422  }
1423 
1424  for(int i=0; i<eim_sd_vector.size(); i++)
1425  {
1426  auto eim = eim_sd_vector[i];
1427  M_mapConvEIM[eim->name()] = std::vector<vectorN_type>(eim->mMax());
1428  for(int j=0; j<eim->mMax(); j++)
1429  M_mapConvEIM[eim->name()][j].resize(sampling_size);
1430  }
1431  }
1432 
1433  void initializeConvergenceCrbMap( int sampling_size )
1434  {
1435  auto N = crb->dimension();
1436  M_mapConvCRB["L2"] = std::vector<vectorN_type>(N);
1437  M_mapConvCRB["H1"] = std::vector<vectorN_type>(N);
1438  M_mapConvCRB["OutputError"] = std::vector<vectorN_type>(N);//true error
1439  M_mapConvCRB["OutputEstimatedError"] = std::vector<vectorN_type>(N);//estimated error
1440  M_mapConvCRB["OutputErrorBoundEfficiency"] = std::vector<vectorN_type>(N);
1441  M_mapConvCRB["SolutionErrorBoundEfficiency"] = std::vector<vectorN_type>(N);
1442  M_mapConvCRB["SolutionError"] = std::vector<vectorN_type>(N);
1443  M_mapConvCRB["SolutionErrorEstimated"] = std::vector<vectorN_type>(N);
1444  M_mapConvCRB["SolutionDualErrorBoundEfficiency"] = std::vector<vectorN_type>(N);
1445  M_mapConvCRB["SolutionDualError"] = std::vector<vectorN_type>(N);
1446  M_mapConvCRB["SolutionDualErrorEstimated"] = std::vector<vectorN_type>(N);
1447  M_mapConvCRB["PrimalResidualNorm"] = std::vector<vectorN_type>(N);
1448  M_mapConvCRB["DualResidualNorm"] = std::vector<vectorN_type>(N);
1449 
1450  for(int j=0; j<N; j++)
1451  {
1452  M_mapConvCRB["L2"][j].resize(sampling_size);
1453  M_mapConvCRB["H1"][j].resize(sampling_size);
1454  M_mapConvCRB["OutputError"][j].resize(sampling_size);
1455  M_mapConvCRB["OutputEstimatedError"][j].resize(sampling_size);
1456  M_mapConvCRB["OutputErrorBoundEfficiency"][j].resize(sampling_size);
1457  M_mapConvCRB["SolutionErrorBoundEfficiency"][j].resize(sampling_size);
1458  M_mapConvCRB["SolutionError"][j].resize(sampling_size);
1459  M_mapConvCRB["SolutionErrorEstimated"][j].resize(sampling_size);
1460  M_mapConvCRB["SolutionDualErrorBoundEfficiency"][j].resize(sampling_size);
1461  M_mapConvCRB["SolutionDualError"][j].resize(sampling_size);
1462  M_mapConvCRB["SolutionDualErrorEstimated"][j].resize(sampling_size);
1463  M_mapConvCRB["PrimalResidualNorm"][j].resize( sampling_size );
1464  M_mapConvCRB["DualResidualNorm"][j].resize( sampling_size );
1465  }
1466  }
1467 
1468  void initializeConvergenceScmMap( int sampling_size )
1469  {
1470  auto N = crb->scm()->KMax();
1471  M_mapConvSCM["RelativeError"] = std::vector<vectorN_type>(N);
1472 
1473  for(int j=0; j<N; j++)
1474  {
1475  M_mapConvSCM["RelativeError"][j].resize(sampling_size);
1476  }
1477  }
1478 
1479  void studyEimConvergence( typename ModelType::parameter_type const& mu , element_type & model_solution , int mu_number)
1480  {
1481  auto eim_sc_vector = model->scalarContinuousEim();
1482  auto eim_sd_vector = model->scalarDiscontinuousEim();
1483 
1484  for(int i=0; i<eim_sc_vector.size(); i++)
1485  {
1486  std::vector<double> l2error;
1487  auto eim = eim_sc_vector[i];
1488  //take two parameters : the mu and the solution of the model
1489  l2error = eim->studyConvergence( mu , model_solution );
1490 
1491  for(int j=0; j<l2error.size(); j++)
1492  M_mapConvEIM[eim->name()][j][mu_number-1] = l2error[j];
1493  }
1494 
1495  for(int i=0; i<eim_sd_vector.size(); i++)
1496  {
1497  std::vector<double> l2error;
1498  auto eim = eim_sd_vector[i];
1499  l2error = eim->studyConvergence( mu , model_solution );
1500 
1501  for(int j=0; j<l2error.size(); j++)
1502  M_mapConvEIM[eim->name()][j][mu_number-1] = l2error[j];
1503  }
1504  }
1505 
1506  void doTheEimConvergenceStat( int sampling_size )
1507  {
1508  auto eim_sc_vector = model->scalarContinuousEim();
1509  auto eim_sd_vector = model->scalarDiscontinuousEim();
1510 
1511  for(int i=0; i<eim_sc_vector.size(); i++)
1512  {
1513  auto eim = eim_sc_vector[i];
1514 
1515  std::ofstream conv;
1516  std::string file_name = "cvg-eim-"+eim->name()+"-stats.dat";
1517 
1518  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1519  {
1520  conv.open(file_name, std::ios::app);
1521  conv << "NbBasis" << "\t" << "Min" << "\t" << "Max" << "\t" << "Mean" << "\t" << "Variance" << "\n";
1522  }
1523 
1524  for(int j=0; j<eim->mMax(); j++)
1525  {
1526  double mean = M_mapConvEIM[eim->name()][j].mean();
1527  double variance = 0.0;
1528  for( int k=0; k < sampling_size ; k++)
1529  {
1530  variance += (M_mapConvEIM[eim->name()][j](k) - mean)*(M_mapConvEIM[eim->name()][j](k) - mean)/sampling_size;
1531  }
1532 
1533  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1534  {
1535  conv << j+1 << "\t"
1536  << M_mapConvEIM[eim->name()][j].minCoeff() << "\t"
1537  << M_mapConvEIM[eim->name()][j].maxCoeff() << "\t"
1538  << mean << "\t" << variance << "\n";
1539  }
1540  }
1541  conv.close();
1542  }
1543 
1544  for(int i=0; i<eim_sd_vector.size(); i++)
1545  {
1546  auto eim = eim_sd_vector[i];
1547 
1548  std::ofstream conv;
1549  std::string file_name = "cvg-eim-"+eim->name()+"-stats.dat";
1550 
1551  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1552  {
1553  conv.open(file_name, std::ios::app);
1554  conv << "NbBasis" << "\t" << "Min" << "\t" << "Max" << "\t" << "Mean" << "\t" << "Variance" << "\n";
1555  }
1556 
1557  for(int j=0; j<eim->mMax(); j++)
1558  {
1559  double mean = M_mapConvEIM[eim->name()][j].mean();
1560  double variance = 0.0;
1561  for( int k=0; k < sampling_size; k++)
1562  variance += (M_mapConvEIM[eim->name()][j](k) - mean)*(M_mapConvEIM[eim->name()][j](k) - mean)/sampling_size;
1563 
1564  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1565  {
1566  conv << j+1 << "\t"
1567  << M_mapConvEIM[eim->name()][j].minCoeff() << "\t"
1568  << M_mapConvEIM[eim->name()][j].maxCoeff() << "\t"
1569  << mean << "\t" << variance << "\n";
1570  }
1571  }
1572  conv.close();
1573  }
1574  }
1575 
1576  void doTheCrbConvergenceStat( int sampling_size )
1577  {
1578  auto N = crb->dimension();
1579  //std::list<std::string> list_error_type;
1580  std::list<std::string> list_error_type = boost::assign::list_of("L2")("H1")("OutputError")("OutputEstimatedError")
1581  ("SolutionError")("SolutionDualError")("PrimalResidualNorm")("DualResidualNorm")("SolutionErrorEstimated")
1582  ("SolutionDualErrorEstimated")("SolutionErrorBoundEfficiency")("SolutionDualErrorBoundEfficiency")("OutputErrorBoundEfficiency");
1583 
1584 
1585  std::vector< Eigen::MatrixXf::Index > index_max_vector_solution_primal;
1586  std::vector< Eigen::MatrixXf::Index > index_max_vector_solution_dual;
1587  std::vector< Eigen::MatrixXf::Index > index_max_vector_output;
1588  std::vector< Eigen::MatrixXf::Index > index_min_vector_solution_primal;
1589  std::vector< Eigen::MatrixXf::Index > index_min_vector_solution_dual;
1590  std::vector< Eigen::MatrixXf::Index > index_min_vector_output;
1591  Eigen::MatrixXf::Index index_max;
1592  Eigen::MatrixXf::Index index_min;
1593 
1594  BOOST_FOREACH( auto error_name, list_error_type)
1595  {
1596  std::ofstream conv;
1597  std::string file_name = "cvg-crb-"+ error_name +"-stats.dat";
1598 
1599  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1600  {
1601  if( error_name=="SolutionErrorEstimated" || error_name=="SolutionDualErrorEstimated" || error_name=="OutputEstimatedError" ||
1602  error_name=="SolutionErrorBoundEfficiency" || error_name=="SolutionDualErrorBoundEfficiency" || error_name=="OutputErrorBoundEfficiency")
1603  {
1604  conv.open(file_name, std::ios::app);
1605  conv << "NbBasis" << "\t" << "Min2" << "\t" << "Max2" << "\t" << "Min" << "\t" << "Max\t" <<"Mean" << "\t"<< "Variance" << "\n";
1606  }
1607  else
1608  {
1609  conv.open(file_name, std::ios::app);
1610  conv << "NbBasis" << "\t" << "Min" << "\t" << "Max" << "\t" << "Mean" << "\t" << "Variance" << "\n";
1611  }
1612  }
1613 
1614  for(int j=0; j<N; j++)
1615  {
1616  double mean = M_mapConvCRB[error_name][j].mean();
1617  double variance = 0.0;
1618  for( int k=0; k < sampling_size; k++)
1619  variance += (M_mapConvCRB[error_name][j](k) - mean)*(M_mapConvCRB[error_name][j](k) - mean)/sampling_size;
1620 
1621  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1622  {
1623 
1624  if( error_name=="SolutionErrorEstimated" || error_name=="SolutionDualErrorEstimated" || error_name=="OutputEstimatedError" ||
1625  error_name=="SolutionErrorBoundEfficiency" || error_name=="SolutionDualErrorBoundEfficiency" || error_name=="OutputErrorBoundEfficiency")
1626  {
1627  double max2=0;
1628  double min2=0;
1629  if( error_name=="SolutionErrorEstimated" )
1630  {
1631  index_max = index_max_vector_solution_primal[j];
1632  index_min = index_min_vector_solution_primal[j];
1633  max2 = M_mapConvCRB[error_name][j]( index_max );
1634  min2 = M_mapConvCRB[error_name][j]( index_min );
1635  }
1636  if( error_name=="SolutionErrorBoundEfficiency" )
1637  {
1638  index_max = index_max_vector_solution_primal[j];
1639  index_min = index_min_vector_solution_primal[j];
1640  double max_estimated = M_mapConvCRB["SolutionErrorEstimated"][j]( index_max );
1641  double min_estimated = M_mapConvCRB["SolutionErrorEstimated"][j]( index_min );
1642  double max = M_mapConvCRB["SolutionError"][j]( index_max );
1643  double min = M_mapConvCRB["SolutionError"][j]( index_min );
1644  max2 = max_estimated / max;
1645  min2 = min_estimated / min;
1646  }
1647  if( error_name=="SolutionDualErrorEstimated" )
1648  {
1649  index_max = index_max_vector_solution_dual[j];
1650  index_min = index_min_vector_solution_dual[j];
1651  max2 = M_mapConvCRB[error_name][j]( index_max );
1652  min2 = M_mapConvCRB[error_name][j]( index_min );
1653  }
1654  if( error_name=="SolutionDualErrorBoundEfficiency" )
1655  {
1656  index_max = index_max_vector_solution_dual[j];
1657  index_min = index_min_vector_solution_dual[j];
1658  double max_estimated = M_mapConvCRB["SolutionDualErrorEstimated"][j]( index_max );
1659  double min_estimated = M_mapConvCRB["SolutionDualErrorEstimated"][j]( index_min );
1660  double max = M_mapConvCRB["SolutionDualError"][j]( index_max );
1661  double min = M_mapConvCRB["SolutionDualError"][j]( index_min );
1662  max2 = max_estimated / max;
1663  min2 = min_estimated / min;
1664  }
1665 
1666  if( error_name=="OutputEstimatedError" )
1667  {
1668  index_max = index_max_vector_output[j];
1669  index_min = index_min_vector_output[j];
1670  max2 = M_mapConvCRB[error_name][j]( index_max );
1671  min2 = M_mapConvCRB[error_name][j]( index_min );
1672  }
1673  if( error_name=="OutputErrorBoundEfficiency" )
1674  {
1675  index_max = index_max_vector_output[j];
1676  index_min = index_min_vector_output[j];
1677  double max_estimated = M_mapConvCRB["OutputEstimatedError"][j]( index_max );
1678  double min_estimated = M_mapConvCRB["OutputEstimatedError"][j]( index_min );
1679  double max = M_mapConvCRB["OutputError"][j]( index_max );
1680  double min = M_mapConvCRB["OutputError"][j]( index_min );
1681  max2 = max_estimated / max;
1682  min2 = min_estimated / min;
1683  }
1684  double min = M_mapConvCRB[error_name][j].minCoeff(&index_min) ;
1685  double max = M_mapConvCRB[error_name][j].maxCoeff(&index_max) ;
1686  conv << j+1 << "\t"
1687  << min2 << "\t" << max2 << "\t" << min << "\t" << max << "\t" << mean << "\t" << variance << "\n";
1688 
1689  }
1690  else
1691  {
1692  conv << j+1 << "\t"
1693  << M_mapConvCRB[error_name][j].minCoeff(&index_min) << "\t"
1694  << M_mapConvCRB[error_name][j].maxCoeff(&index_max) << "\t"
1695  << mean << "\t" << variance << "\n";
1696  if( error_name=="OutputError" )
1697  {
1698  index_max_vector_output.push_back( index_max );
1699  index_min_vector_output.push_back( index_min );
1700  }
1701  if( error_name=="SolutionError")
1702  {
1703  index_max_vector_solution_primal.push_back( index_max );
1704  index_min_vector_solution_primal.push_back( index_min );
1705  }
1706  if( error_name=="SolutionDualError")
1707  {
1708  index_max_vector_solution_dual.push_back( index_max );
1709  index_min_vector_solution_dual.push_back( index_min );
1710  }
1711  }
1712  }//master proc
1713  }//loop over number of RB elements
1714 
1715  conv.close();
1716  }
1717 
1718  int Nmax = vector_sampling_for_primal_efficiency_under_1.size();
1719  for(int N=0; N<Nmax; N++)
1720  {
1721  std::string file_name_primal = (boost::format("Sampling-Primal-Problem-Bad-Efficiency-N=%1%") %(N+1) ).str();
1722  if( vector_sampling_for_primal_efficiency_under_1[N]->size() > 1 )
1723  vector_sampling_for_primal_efficiency_under_1[N]->writeOnFile(file_name_primal);
1724  }
1725  Nmax = vector_sampling_for_dual_efficiency_under_1.size();
1726  for(int N=0; N<Nmax; N++)
1727  {
1728  std::string file_name_dual = (boost::format("Sampling-Dual-Problem-Bad-Efficiency-N=%1%") %(N+1) ).str();
1729  if( vector_sampling_for_dual_efficiency_under_1[N]->size() > 1 )
1730  vector_sampling_for_dual_efficiency_under_1[N]->writeOnFile(file_name_dual);
1731  }
1732 
1733 
1734  }
1735 
1736  void doTheScmConvergenceStat( int sampling_size )
1737  {
1738  auto N = crb->scm()->KMax();
1739  std::list<std::string> list_error_type = boost::assign::list_of("RelativeError");
1740  BOOST_FOREACH( auto error_name, list_error_type)
1741  {
1742  std::ofstream conv;
1743  std::string file_name = "cvg-scm-"+ error_name +"-stats.dat";
1744 
1745  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1746  {
1747  conv.open(file_name, std::ios::app);
1748  conv << "Nb_basis" << "\t" << "Min" << "\t" << "Max" << "\t" << "Mean" << "\t" << "Variance" << "\n";
1749  }
1750 
1751  for(int j=0; j<N; j++)
1752  {
1753  double mean = M_mapConvSCM[error_name][j].mean();
1754  double variance = 0.0;
1755  for( int k=0; k < sampling_size; k++)
1756  variance += (M_mapConvSCM[error_name][j](k) - mean)*(M_mapConvSCM[error_name][j](k) - mean)/sampling_size;
1757 
1758  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1759  {
1760  conv << j+1 << "\t"
1761  << M_mapConvSCM[error_name][j].minCoeff() << "\t"
1762  << M_mapConvSCM[error_name][j].maxCoeff() << "\t"
1763  << mean << "\t" << variance << "\n";
1764  }
1765  }
1766  conv.close();
1767  }
1768  }
1769 
1770 
1771  // Script write current mu in cfg => need to write it in SamplingForTest
1772  void buildSamplingFromCfg()
1773  {
1774  // Size of mu
1775  int mu_size = model->parameterSpace()->dimension();
1776 
1777  // Clear SamplingForTest is exists, and open a new one
1778  fs::path input_file ("SamplingForTest");
1779  if( fs::exists(input_file) )
1780  std::remove( "SamplingForTest" );
1781 
1782  std::ofstream input( "SamplingForTest" );
1783  input << "mu= [ ";
1784 
1785  // Check if cfg file is readable
1786  std::ifstream cfg_file( option(_name="config-file").template as<std::string>() );
1787  if(!cfg_file)
1788  std::cout << "[Script-mode] Config file cannot be read" << std::endl;
1789 
1790  // OT writes values of mu in config file => read it and copy in SamplingForTest with specific syntax
1791  for(int i=1; i<=mu_size; i++)
1792  {
1793  // convert i into string
1794  std::ostringstream oss;
1795  oss << i;
1796  std::string is = oss.str();
1797 
1798  // Read cfg file, collect line with current mu_i
1799  std::string cfg_line_mu, tmp_content;
1800  std::ifstream cfg_file( option(_name="config-file").template as<std::string>() );
1801  while(cfg_file)
1802  {
1803  std::getline(cfg_file, tmp_content);
1804  if(tmp_content.compare(0,2+is.size(),"mu"+is) == 0)
1805  cfg_line_mu += tmp_content;
1806  }
1807 
1808  //Regular expression : corresponds to one set in xml file (mu<i>=<value>)
1809  std::string expr_s = "mu"+is+"[[:space:]]*=[[:space:]]*([0-9]+(.?)[0-9]*(e(\\+|-)[0-9]+)?)[[:space:]]*";
1810  boost::regex expression( expr_s );
1811 
1812  //Match mu<i>=<value> in cfg file and copy to SamplingForTest
1813  boost::smatch what;
1814  auto is_match = boost::regex_match(cfg_line_mu, what, expression);
1815  //std::cout << "is match ?" << is_match << std::endl;
1816  if(is_match)
1817  {
1818  // what[0] is the complete string mu<i>=<value>
1819  // what[1] is the submatch <value>
1820  //std::cout << "what 1 = " << what[1] << std::endl;
1821  if( i!=mu_size )
1822  input << what[1] << " , ";
1823  else
1824  input << what[1] << " ]";
1825  }
1826  }
1827  }
1828 
1829 private:
1830  CRBModelMode M_mode;
1831  crbmodel_ptrtype model;
1832  crb_ptrtype crb;
1833 
1834  // For EIM convergence study
1835  std::map<std::string, std::vector<vectorN_type> > M_mapConvEIM;
1836  // For CRB convergence study
1837  std::map<std::string, std::vector<vectorN_type> > M_mapConvCRB;
1838  // For SCM convergence study
1839  std::map<std::string, std::vector<vectorN_type> > M_mapConvSCM;
1840 
1841  //vector of sampling to stock parameters for which the efficiency is under 1
1842  std::vector< sampling_ptrtype > vector_sampling_for_primal_efficiency_under_1;
1843  std::vector< sampling_ptrtype > vector_sampling_for_dual_efficiency_under_1;
1844 
1845  fs::path M_current_path;
1846 }; // OpusApp
1847 
1848 } // Feel
1849 
1850 #endif /* __OpusApp_H */
1851 
provides information about the Application
Definition: application.hpp:71
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
elements_type const & elements() const
Definition: elements.hpp:355
void setLogs()
Definition: application.cpp:687
void add(Simget *simget)
Definition: application.cpp:877
model_type::functionspace_type functionspace_type
function space type
Definition: opusapp.hpp:96
po::options_description opusapp_options(std::string const &prefix)
Definition: opusapp.cpp:33
Application & changeRepository(boost::format)
Definition: application.cpp:818
Certifed Reduced Basis class.
Definition: crb.hpp:95
Certified Reduced Basis application.
Definition: opusapp.hpp:84
Holds information needed by the &quot;About&quot; box and other classes.
Definition: about.hpp:173
ModelType model_type
model type
Definition: opusapp.hpp:92
Certified Reduced Basis Model class.
Definition: crbmodel.hpp:73
po::options_description const & optionsDescription() const
Definition: application.hpp:175
static boost::shared_ptr< Exporter< MeshType, N > > New(std::string const &exportername, std::string prefix=Environment::about().appName(), WorldComm const &worldComm=Environment::worldComm())
Definition: exporterimpl.hpp:127
FEELPP_DONT_INLINE void run()
Definition: opusapp.hpp:318
po::variables_map const & vm() const
Definition: application.hpp:186
AboutData const & about() const
Definition: application.hpp:198
void run(const double *X, unsigned long N, double *Y, unsigned long P)
Definition: opusapp.hpp:1308
std::string appName() const
Definition: about.cpp:225

Generated on Sun Dec 22 2013 13:11:10 for Feel++ by doxygen 1.8.5