32 #include <feel/feel.hpp>
34 #include <boost/assign/std/vector.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>
45 #include <Eigen/Dense>
51 std::string _o( std::string
const& prefix, std::string
const& opt )
53 std::string o = prefix;
62 enum class SamplingMode
64 RANDOM = 0, EQUIDISTRIBUTED = 1, LOGEQUIDISTRIBUTED = 2
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
81 template<
typename ModelType,
82 template <
typename ReducedMethod >
class RM=CRB,
83 template <
typename ModelInterface >
class Model=CRBModel>
89 typedef double value_type;
93 typedef boost::shared_ptr<ModelType> model_ptrtype;
97 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
99 typedef typename model_type::element_type element_type;
101 typedef Eigen::VectorXd vectorN_type;
106 typedef boost::shared_ptr<crbmodel_type> crbmodel_ptrtype;
108 typedef boost::shared_ptr<crb_type> crb_ptrtype;
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;
118 typedef typename ModelType::parameter_type parameter_type;
119 typedef std::vector< parameter_type > vector_parameter_type;
121 typedef typename crb_type::sampling_ptrtype sampling_ptrtype;
126 M_mode( ( CRBModelMode )option(_name=_o( this->
about().appName(),
"run.mode" )).template as<int>() )
134 M_mode( ( CRBModelMode )option(_name=_o( this->
about().appName(),
"run.mode" )).template as<int>() )
139 OpusApp( AboutData
const& ad, po::options_description
const& od, CRBModelMode mode )
147 OpusApp(
int argc,
char** argv, AboutData
const& ad, po::options_description
const& od )
150 M_mode( ( CRBModelMode )option(_name=_o( this->
about().appName(),
"run.mode" )).template as<int>() )
154 OpusApp(
int argc,
char** argv, AboutData
const& ad, po::options_description
const& od, CRBModelMode mode )
165 M_current_path = fs::current_path();
167 std::srand( static_cast<unsigned>( std::time( 0 ) ) );
168 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
170 LOG(INFO) <<
"[OpusApp] constructor " << this->
about().
appName() <<
"\n";
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>();
178 results_repo_name =
"default_repo";
180 LOG(INFO) <<
"Name for results repo : " << results_repo_name <<
"\n";
182 % this->
about().appName()
186 LOG(INFO) <<
"[OpusApp] ch repo" <<
"\n";
188 LOG(INFO) <<
"[OpusApp] set Logs" <<
"\n";
189 LOG(INFO) <<
"[OpusApp] mode:" << ( int )M_mode <<
"\n";
191 model = crbmodel_ptrtype(
new crbmodel_type( this->
vm(),M_mode ) );
192 LOG(INFO) <<
"[OpusApp] get model done" <<
"\n";
194 crb = crb_ptrtype(
new crb_type( this->
about().appName(),
197 LOG(INFO) <<
"[OpusApp] get crb done" <<
"\n";
204 catch ( boost::bad_any_cast
const& e )
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";
212 void setMode( std::string
const& mode )
214 if ( mode ==
"pfem" ) M_mode = CRBModelMode::PFEM;
216 if ( mode ==
"crb" ) M_mode = CRBModelMode::CRB;
218 if ( mode ==
"scm" ) M_mode = CRBModelMode::SCM;
220 if ( mode ==
"scm_online" ) M_mode = CRBModelMode::SCM_ONLINE;
222 if ( mode ==
"crb_online" ) M_mode = CRBModelMode::CRB_ONLINE;
224 void setMode( CRBModelMode mode )
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 )
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);
247 if( NlogEquidistributedScm+NequidistributedScm > 0 )
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);
257 if ( M_mode == CRBModelMode::PFEM )
260 if ( !crb->scm()->isDBLoaded() || crb->scm()->rebuildDB() )
262 if ( M_mode == CRBModelMode::SCM )
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 );
269 crb->scm()->offline();
273 if ( !crb->isDBLoaded() || crb->rebuildDB() )
275 if ( M_mode == CRBModelMode::CRB )
278 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
279 std::cout <<
"No CRB DB available, do crb offline computations...\n";
280 crb->setOfflineStep(
true );
284 else if ( M_mode != CRBModelMode::SCM )
285 throw std::logic_error(
"CRB/SCM Database could not be loaded" );
288 if( crb->isDBLoaded() )
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);
297 if( sampling_size > current_dimension )
300 if( current_dimension < dimension_max && !crb_use_predefined )
310 crb->setOfflineStep(
true );
321 bool export_solution = option(_name=_o( this->
about().appName(),
"export-solution" )).template as<bool>();
323 int proc_number = Environment::worldComm().globalRank();
325 if ( this->
vm().count(
"help" ) )
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>();
337 typename crb_type::sampling_ptrtype Sampling(
new typename crb_type::sampling_type( model->parameterSpace() ) );
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>();
343 if( n_eval_computational_time > 0 )
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" );
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);
359 run_sampling_size = 0;
361 n_eval_computational_time = option(_name=
"crb.computational-time-neval").template as<int>();
362 if( n_eval_computational_time > 0 )
364 if( ! option(_name=
"crb.cvg-study").template as<bool>() )
367 run_sampling_size = 0;
370 crb->computationalTimeStatistics( appname );
374 switch ( run_sampling_type )
377 case SamplingMode::RANDOM:
378 Sampling->randomize( run_sampling_size );
381 case SamplingMode::EQUIDISTRIBUTED:
382 Sampling->equidistribute( run_sampling_size );
385 case SamplingMode::LOGEQUIDISTRIBUTED:
386 Sampling->logEquidistribute( run_sampling_size );
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;
407 if( crb->printErrorDuringOfflineStep() )
408 crb->printErrorsDuringRbConstruction();
409 if ( crb->showMuSelection() )
410 crb->printMuSelection();
414 if( export_solution )
415 exporter->step( 0 )->setMesh( model->functionSpace()->mesh() );
417 printParameterHdr( ostr, model->parameterSpace()->dimension(), hdrs[M_mode] );
419 int crb_error_type = option(_name=
"crb.error-type").template as<int>();
422 if( M_mode==CRBModelMode::CRB )
424 dim=crb->dimension();
426 Sampling = crb->wnmu();
428 if( option(_name=
"crb.run-on-scm-parameters").
template as<bool>() )
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." );
435 if( M_mode==CRBModelMode::SCM )
437 dim=crb->scm()->KMax();
438 if( option(_name=
"crb.scm.run-on-C").
template as<bool>() )
439 Sampling = crb->scm()->c();
442 std::ofstream file_summary_of_simulations( ( boost::format(
"summary_of_simulations_%d" ) %dim ).str().c_str() ,std::ios::out | std::ios::app );
461 if( option(_name=
"crb.script-mode").
template as<bool>() )
464 if( proc_number == Environment::worldComm().masterRank() )
465 buildSamplingFromCfg();
467 Environment::worldComm().barrier();
470 compute_stat =
false;
478 if( option(_name=
"crb.use-predefined-test-sampling").
template as<bool>() || option(_name=
"crb.script-mode").
template as<bool>() )
480 std::string file_name = ( boost::format(
"SamplingForTest") ).str();
481 std::ifstream file ( file_name );
485 Sampling->readFromFile( file_name ) ;
489 VLOG(2) <<
"proc number : " << proc_number <<
"can't find file \n";
490 throw std::logic_error(
"[OpusApp] file SamplingForTest was not found" );
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;
502 vectorN_type scm_relative_error;
504 bool solve_dual_problem = option(_name=
"crb.solve-dual-problem").template as<bool>();
506 if (option(_name=
"crb.cvg-study").template as<bool>() && compute_fem )
508 int Nmax = crb->dimension();
509 vector_sampling_for_primal_efficiency_under_1.resize(Nmax);
510 for(
int N=0; N<Nmax; N++)
512 sampling_ptrtype sampling_primal (
new typename crb_type::sampling_type( model->parameterSpace() ) );
513 vector_sampling_for_primal_efficiency_under_1[N]=sampling_primal;
515 if( solve_dual_problem )
517 vector_sampling_for_dual_efficiency_under_1.resize(Nmax);
518 for(
int N=0; N<Nmax; N++)
520 sampling_ptrtype sampling_dual (
new typename crb_type::sampling_type( model->parameterSpace() ) );
521 vector_sampling_for_dual_efficiency_under_1[N]=sampling_dual;
526 if( M_mode==CRBModelMode::CRB )
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() );
535 if( crb->errorType()!=2 )
536 relative_estimated_error_vector.resize( Sampling->size() );
538 crb->setOfflineStep(
false );
540 if (option(_name=
"eim.cvg-study").
template as<bool>() )
542 this->initializeConvergenceEimMap( Sampling->size() );
545 if (option(_name=
"crb.cvg-study").template as<bool>() )
546 this->initializeConvergenceCrbMap( Sampling->size() );
549 if( M_mode==CRBModelMode::SCM )
551 if (option(_name=
"crb.scm.cvg-study").template as<bool>() )
552 this->initializeConvergenceScmMap( Sampling->size() );
554 scm_relative_error.resize( Sampling->size() );
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>();
562 double output_fem = -1;
566 model->computeAffineDecomposition();
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++)
577 double time = time_index*dt;
578 ref_betaAqm.push_back( model->computeBetaQm( ref_mu , time ).template get<1>() );
580 auto ref_betaMqm = model->computeBetaQm( ref_mu , tf ).template get<0>() ;
583 BOOST_FOREACH(
auto mu, *Sampling )
585 int size = mu.size();
588 element_type u_crb_dual;
590 if( proc_number == Environment::worldComm().masterRank() )
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 ";
599 std::ostringstream mu_str;
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];
608 LOG(INFO) <<
"mu=" << mu <<
"\n";
611 if( option(_name=
"crb.script-mode").
template as<bool>() )
613 unsigned long N = mu.size() + 5;
615 std::vector<double> X( N );
616 std::vector<double> Y( P );
617 for(
int i=0; i<mu.size(); i++)
620 int N_dim = crb_dimension;
621 int N_dimMax = crb_dimension_max;
627 X[N-5] = output_index;
629 X[N-3] = crb_online_tolerance;
630 X[N-2] = crb_error_type;
633 bool compute_variance = crb_compute_variance;
634 if ( compute_variance )
637 this->
run( X.data(), X.size(), Y.data(), Y.size() );
640 std::ofstream res(option(_name=
"result-file").
template as<std::string>() );
641 res <<
"output="<< Y[0] <<
"\n";
647 case CRBModelMode::PFEM:
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;
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() );
660 LOG(INFO) <<
"compute output\n";
661 if( export_solution )
662 exporter->step(0)->add( u_fem.name(), u_fem );
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 );
668 std::ofstream res(option(_name=
"result-file").
template as<std::string>() );
669 res <<
"output="<< o[0] <<
"\n";
674 case CRBModelMode::CRB:
676 LOG(INFO) <<
"CRB mode\n";
677 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
678 std::cout <<
"CRB mode\n";
681 boost::mpi::timer ti;
684 LOG(INFO) <<
"solve crb\n";
685 google::FlushLogFiles(google::GLOG_INFO);
688 int N = option(_name=
"crb.dimension").template as<int>();
690 auto o = crb->run( mu, option(_name=
"crb.online-tolerance").
template as<double>() , N);
691 double time_crb = ti.elapsed();
694 auto WNdu = crb->wndu();
696 auto solutions=o.template get<2>();
697 auto uN = solutions.template get<0>();
698 auto uNdu = solutions.template get<1>();
700 int size = uN.size();
703 u_crb = crb->expansion( uN[size-1] , N , WN );
704 if( solve_dual_problem )
705 u_crb_dual = crb->expansion( uNdu[0] , N , WNdu );
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 );
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;
726 element_type u_dual_fem ;
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>();
735 bool use_newton = option(_name=
"crb.use-newton").template as<bool>();
738 LOG(INFO) <<
"solve u_fem\n";
743 if( boost::is_same< crbmodel_type , crbmodelbilinear_type >::value && ! use_newton )
745 u_fem = model->solveFemUsingAffineDecompositionFixedPoint( mu );
747 u_fem = model->solve( mu );
749 std::ostringstream u_fem_str;
750 u_fem_str <<
"u_fem(" << mu_str.str() <<
")";
751 u_fem.setName( u_fem_str.str() );
753 if( export_solution )
755 LOG(INFO) <<
"export u_fem \n";
756 exporter->step(0)->add( u_fem.name(), u_fem );
758 std::vector<double> ofem = boost::assign::list_of( model->output( output_index,mu, u_fem ) )( ti.elapsed() );
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];
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 );
779 output_fem = ofem[0];
782 if( boost::is_same< crbmodel_type , crbmodelbilinear_type >::value && ! use_newton )
784 if( solve_dual_problem )
786 u_dual_fem = model->solveFemDualUsingAffineDecompositionFixedPoint( mu );
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 );
797 if ( crb->errorType()==2 )
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() )
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 );
808 if ( option(_name=
"crb.compute-stat").
template as<bool>() && compute_fem )
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;
817 std::ofstream res(option(_name=
"result-file").
template as<std::string>() );
818 res <<
"output="<< ocrb <<
"\n";
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() )
831 std::cout <<
"output=" << ocrb <<
" with " << o.template get<1>() <<
" basis functions (relative error estimation on this output : " << relative_estimated_error<<
") \n";
833 printEntry( file_summary_of_simulations, mu, v );
834 printEntry( ostr, mu, v );
837 if ( option(_name=
"crb.compute-stat").
template as<bool>() && compute_fem )
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;
846 std::ofstream res(option(_name=
"result-file").
template as<std::string>() );
847 res <<
"output="<< ocrb <<
"\n";
851 if (option(_name=
"eim.cvg-study").
template as<bool>() )
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")
857 LOG( INFO ) <<
"convergence eim with CRB-with-ad ";
859 this->studyEimConvergence( mu , u_crb , curpar );
861 if( how_compute_unknown ==
"FEM-with-ad")
863 LOG( INFO ) <<
"convergence eim with FEM-with-ad ";
866 auto fem_with_ad = model->solveFemUsingAffineDecompositionFixedPoint( mu );
867 this->studyEimConvergence( mu , fem_with_ad , curpar );
869 if( how_compute_unknown ==
"FEM-without-ad")
871 LOG( INFO ) <<
"convergence eim with FEM-without-ad ";
873 auto fem_without_ad = model->solve( mu );
874 this->studyEimConvergence( mu , fem_without_ad , curpar );
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" );
880 if (option(_name=
"crb.cvg-study").
template as<bool>() && compute_fem )
883 LOG(INFO) <<
"start convergence study...\n";
884 std::map<int, boost::tuple<double,double,double,double,double,double,double> > conver;
886 int Nmax = crb->dimension();
887 for(
int N = 1; N <= Nmax ; N++ )
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 );
899 auto u_error = u_fem - uN;
900 auto u_dual_error = model->functionSpace()->element();
901 if( solve_dual_problem )
903 uNdu = crb->expansion( u_crb_du[0], N, WNdu );
904 u_dual_error = u_dual_fem - uNdu;
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>();
913 double rel_err = std::abs( output_fem-ocrb ) /output_fem;
915 double output_relative_estimated_error = output_estimated_error / output_fem;
917 double primal_residual_norm = o.template get<4>();
918 double dual_residual_norm = o.template get<5>();
920 double solution_error=0;
921 double dual_solution_error=0;
924 if( model->isSteady() )
928 for(
int q=0; q<model->Qa();q++)
930 for(
int m=0; m<model->mMaxA(q); m++)
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);
937 if( solve_dual_problem )
939 for(
int q=0; q<model->Qa();q++)
941 for(
int m=0; m<model->mMaxA(q); m++)
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);
947 dual_solution_error = math::sqrt( dual_solution_error );
948 ref_dual = math::sqrt( ref_dual );
951 solution_error = math::sqrt( solution_error );
952 ref_primal = math::sqrt( ref_primal );
958 double dt = model->timeStep();
959 double ti = model->timeInitial();
960 double tf = model->timeFinal();
961 int K = ( tf - ti )/dt;
963 for(
int q=0; q<model->Qm();q++)
965 for(
int m=0; m<model->mMaxM(q); m++)
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);
971 for(
int time_index=0; time_index<K; time_index++)
973 double t=time_index*dt;
974 for(
int q=0; q<model->Qa();q++)
976 for(
int m=0; m<model->mMaxA(q); m++)
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;
983 solution_error = math::sqrt( solution_error );
984 ref_primal = math::sqrt( ref_primal );
986 if( solve_dual_problem )
988 ti = model->timeFinal()+dt;
989 tf = model->timeInitial()+dt;
992 for(
int q=0; q<model->Qm();q++)
994 for(
int m=0; m<model->mMaxM(q); m++)
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);
1001 for(
int time_index=0; time_index<K; time_index++)
1003 double t=time_index*dt;
1004 for(
int q=0; q<model->Qa();q++)
1006 for(
int m=0; m<model->mMaxA(q); m++)
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;
1014 dual_solution_error = math::sqrt( dual_solution_error );
1015 ref_dual = math::sqrt( ref_dual );
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;
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;
1030 if( relative_primal_solution_error_bound_efficiency < 1 )
1032 vector_sampling_for_primal_efficiency_under_1[N-1]->push_back( mu , 1);
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 )
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;
1044 if( relative_dual_solution_error_bound_efficiency < 1 )
1046 vector_sampling_for_dual_efficiency_under_1[N-1]->push_back( mu , 0);
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 );
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";
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 )";
1077 if ( option(_name=
"crb.check.residual").
template as<bool>() && solve_dual_problem )
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() )
1083 crb->checkResidual( mu , primal_residual_coefficients, dual_residual_coefficients, uN, uNdu );
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;
1092 auto u_crb_old = solutions.template get<2>();
1093 auto u_crb_du_old = solutions.template get<3>();
1096 for(
int t=0; t<size; t++)
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 ) );
1104 crb->compareResidualsForTransientProblems(N, mu ,
1105 uNelement, uNelement_old, uNelement_du, uNelement_du_old,
1106 primal_residual_coefficients, dual_residual_coefficients );
1111 if( proc_number == Environment::worldComm().masterRank() )
1113 LOG(INFO) <<
"save in logfile\n";
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";
1126 case CRBModelMode::CRB_ONLINE:
1128 std::cout <<
"CRB Online mode\n";
1129 boost::mpi::timer ti;
1131 auto o = crb->run( mu, option(_name=
"crb.online-tolerance").
template as<double>() );
1133 if ( crb->errorType()==2 )
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 );
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 );
1154 case CRBModelMode::SCM:
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 );
1164 scm_relative_error[curpar-1] = o[6];
1166 if (option(_name=
"crb.scm.cvg-study").
template as<bool>() )
1168 LOG(INFO) <<
"start scm convergence study...\n";
1169 std::map<int, boost::tuple<double> > conver;
1170 for(
int N = 1; N <= kmax; N++ )
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;
1179 if( proc_number == Environment::worldComm().masterRank() )
1181 LOG(INFO) <<
"save in logfile\n";
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>() ;
1195 case CRBModelMode::SCM_ONLINE:
1197 std::cout <<
"SCM online mode\n";
1198 auto o = crb->scm()->run( mu, crb->scm()->KMax() );
1199 printEntry( ostr, mu, o );
1205 LOG( INFO ) <<
"------------------------------------------------------------";
1210 if( export_solution )
1213 if( proc_number == Environment::worldComm().masterRank() ) std::cout << ostr.str() <<
"\n";
1215 if (option(_name=
"eim.cvg-study").
template as<bool>() && M_mode==CRBModelMode::CRB)
1216 this->doTheEimConvergenceStat( Sampling->size() );
1218 if (option(_name=
"crb.cvg-study").template as<bool>() && compute_fem && M_mode==CRBModelMode::CRB )
1219 this->doTheCrbConvergenceStat( Sampling->size() );
1221 if (option(_name=
"crb.scm.cvg-study").template as<bool>() && M_mode==CRBModelMode::SCM )
1222 this->doTheScmConvergenceStat( Sampling->size() );
1224 if ( compute_stat && compute_fem && M_mode==CRBModelMode::CRB )
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;
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;
1259 if( crb->errorType()!=2 )
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();
1265 if( proc_number == Environment::worldComm().masterRank() )
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";
1288 if( M_mode==CRBModelMode::SCM )
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() )
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";
1308 void run(
const double * X,
unsigned long N,
1309 double * Y,
unsigned long P )
1314 case CRBModelMode::PFEM:
1316 model->run( X, N, Y, P );
1320 case CRBModelMode::SCM:
1321 case CRBModelMode::SCM_ONLINE:
1323 crb->scm()->run( X, N, Y, P );
1327 case CRBModelMode::CRB:
1328 case CRBModelMode::CRB_ONLINE:
1330 crb->run( X, N, Y, P );
1335 fs::current_path( M_current_path );
1339 int printParameterHdr( std::ostream& os,
int N, std::vector<std::string> outputhdrs )
1341 for (
int i = 0; i < N; ++i )
1343 std::ostringstream s;
1345 os << hdrmanip( prec+7 ) << s.str();
1348 BOOST_FOREACH(
auto output, outputhdrs )
1350 os << hdrmanip( 15 ) << output;
1354 return N*( prec+7 )+outputhdrs.size()*15;
1356 void printEntry( std::ostream& os,
1357 typename ModelType::parameter_type
const& mu,
1358 std::vector<double>
const& outputs )
1360 for (
int i = 0; i < mu.size(); ++i )
1361 os << std::right <<std::setw( prec+7 ) << dmanip << mu[i];
1363 BOOST_FOREACH(
auto o, outputs )
1365 os << tabmanip( 15 ) << o;
1372 double l2Norm( element_type
const& u )
1374 static const bool is_composite = functionspace_type::is_composite;
1375 return l2Norm( u, mpl::bool_< is_composite >() );
1377 double l2Norm( element_type
const& u, mpl::bool_<false> )
1379 auto mesh = model->functionSpace()->mesh();
1380 return math::sqrt(
integrate(
elements(mesh), (vf::idv(u))*(vf::idv(u)) ).evaluate()(0,0) );
1382 double l2Norm( element_type
const& u, mpl::bool_<true>)
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) );
1389 double h1Norm( element_type
const& u )
1391 static const bool is_composite = functionspace_type::is_composite;
1392 return h1Norm( u, mpl::bool_< is_composite >() );
1394 double h1Norm( element_type
const& u, mpl::bool_<false> )
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 );
1401 double h1Norm( element_type
const& u, mpl::bool_<true>)
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 );
1411 void initializeConvergenceEimMap(
int sampling_size )
1413 auto eim_sc_vector = model->scalarContinuousEim();
1414 auto eim_sd_vector = model->scalarDiscontinuousEim();
1416 for(
int i=0; i<eim_sc_vector.size(); i++)
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);
1424 for(
int i=0; i<eim_sd_vector.size(); i++)
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);
1433 void initializeConvergenceCrbMap(
int sampling_size )
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);
1439 M_mapConvCRB[
"OutputEstimatedError"] = std::vector<vectorN_type>(N);
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);
1450 for(
int j=0; j<N; j++)
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 );
1468 void initializeConvergenceScmMap(
int sampling_size )
1470 auto N = crb->scm()->KMax();
1471 M_mapConvSCM[
"RelativeError"] = std::vector<vectorN_type>(N);
1473 for(
int j=0; j<N; j++)
1475 M_mapConvSCM[
"RelativeError"][j].resize(sampling_size);
1479 void studyEimConvergence(
typename ModelType::parameter_type
const& mu , element_type & model_solution ,
int mu_number)
1481 auto eim_sc_vector = model->scalarContinuousEim();
1482 auto eim_sd_vector = model->scalarDiscontinuousEim();
1484 for(
int i=0; i<eim_sc_vector.size(); i++)
1486 std::vector<double> l2error;
1487 auto eim = eim_sc_vector[i];
1489 l2error = eim->studyConvergence( mu , model_solution );
1491 for(
int j=0; j<l2error.size(); j++)
1492 M_mapConvEIM[eim->name()][j][mu_number-1] = l2error[j];
1495 for(
int i=0; i<eim_sd_vector.size(); i++)
1497 std::vector<double> l2error;
1498 auto eim = eim_sd_vector[i];
1499 l2error = eim->studyConvergence( mu , model_solution );
1501 for(
int j=0; j<l2error.size(); j++)
1502 M_mapConvEIM[eim->name()][j][mu_number-1] = l2error[j];
1506 void doTheEimConvergenceStat(
int sampling_size )
1508 auto eim_sc_vector = model->scalarContinuousEim();
1509 auto eim_sd_vector = model->scalarDiscontinuousEim();
1511 for(
int i=0; i<eim_sc_vector.size(); i++)
1513 auto eim = eim_sc_vector[i];
1516 std::string file_name =
"cvg-eim-"+eim->name()+
"-stats.dat";
1518 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1520 conv.open(file_name, std::ios::app);
1521 conv <<
"NbBasis" <<
"\t" <<
"Min" <<
"\t" <<
"Max" <<
"\t" <<
"Mean" <<
"\t" <<
"Variance" <<
"\n";
1524 for(
int j=0; j<eim->mMax(); j++)
1526 double mean = M_mapConvEIM[eim->name()][j].mean();
1527 double variance = 0.0;
1528 for(
int k=0; k < sampling_size ; k++)
1530 variance += (M_mapConvEIM[eim->name()][j](k) - mean)*(M_mapConvEIM[eim->name()][j](k) - mean)/sampling_size;
1533 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1536 << M_mapConvEIM[eim->name()][j].minCoeff() <<
"\t"
1537 << M_mapConvEIM[eim->name()][j].maxCoeff() <<
"\t"
1538 << mean <<
"\t" << variance <<
"\n";
1544 for(
int i=0; i<eim_sd_vector.size(); i++)
1546 auto eim = eim_sd_vector[i];
1549 std::string file_name =
"cvg-eim-"+eim->name()+
"-stats.dat";
1551 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1553 conv.open(file_name, std::ios::app);
1554 conv <<
"NbBasis" <<
"\t" <<
"Min" <<
"\t" <<
"Max" <<
"\t" <<
"Mean" <<
"\t" <<
"Variance" <<
"\n";
1557 for(
int j=0; j<eim->mMax(); j++)
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;
1564 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1567 << M_mapConvEIM[eim->name()][j].minCoeff() <<
"\t"
1568 << M_mapConvEIM[eim->name()][j].maxCoeff() <<
"\t"
1569 << mean <<
"\t" << variance <<
"\n";
1576 void doTheCrbConvergenceStat(
int sampling_size )
1578 auto N = crb->dimension();
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");
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;
1594 BOOST_FOREACH(
auto error_name, list_error_type)
1597 std::string file_name =
"cvg-crb-"+ error_name +
"-stats.dat";
1599 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1601 if( error_name==
"SolutionErrorEstimated" || error_name==
"SolutionDualErrorEstimated" || error_name==
"OutputEstimatedError" ||
1602 error_name==
"SolutionErrorBoundEfficiency" || error_name==
"SolutionDualErrorBoundEfficiency" || error_name==
"OutputErrorBoundEfficiency")
1604 conv.open(file_name, std::ios::app);
1605 conv <<
"NbBasis" <<
"\t" <<
"Min2" <<
"\t" <<
"Max2" <<
"\t" <<
"Min" <<
"\t" <<
"Max\t" <<
"Mean" <<
"\t"<<
"Variance" <<
"\n";
1609 conv.open(file_name, std::ios::app);
1610 conv <<
"NbBasis" <<
"\t" <<
"Min" <<
"\t" <<
"Max" <<
"\t" <<
"Mean" <<
"\t" <<
"Variance" <<
"\n";
1614 for(
int j=0; j<N; j++)
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;
1621 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1624 if( error_name==
"SolutionErrorEstimated" || error_name==
"SolutionDualErrorEstimated" || error_name==
"OutputEstimatedError" ||
1625 error_name==
"SolutionErrorBoundEfficiency" || error_name==
"SolutionDualErrorBoundEfficiency" || error_name==
"OutputErrorBoundEfficiency")
1629 if( error_name==
"SolutionErrorEstimated" )
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 );
1636 if( error_name==
"SolutionErrorBoundEfficiency" )
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;
1647 if( error_name==
"SolutionDualErrorEstimated" )
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 );
1654 if( error_name==
"SolutionDualErrorBoundEfficiency" )
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;
1666 if( error_name==
"OutputEstimatedError" )
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 );
1673 if( error_name==
"OutputErrorBoundEfficiency" )
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;
1684 double min = M_mapConvCRB[error_name][j].minCoeff(&index_min) ;
1685 double max = M_mapConvCRB[error_name][j].maxCoeff(&index_max) ;
1687 << min2 <<
"\t" << max2 <<
"\t" << min <<
"\t" << max <<
"\t" << mean <<
"\t" << variance <<
"\n";
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" )
1698 index_max_vector_output.push_back( index_max );
1699 index_min_vector_output.push_back( index_min );
1701 if( error_name==
"SolutionError")
1703 index_max_vector_solution_primal.push_back( index_max );
1704 index_min_vector_solution_primal.push_back( index_min );
1706 if( error_name==
"SolutionDualError")
1708 index_max_vector_solution_dual.push_back( index_max );
1709 index_min_vector_solution_dual.push_back( index_min );
1718 int Nmax = vector_sampling_for_primal_efficiency_under_1.size();
1719 for(
int N=0; N<Nmax; N++)
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);
1725 Nmax = vector_sampling_for_dual_efficiency_under_1.size();
1726 for(
int N=0; N<Nmax; N++)
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);
1736 void doTheScmConvergenceStat(
int sampling_size )
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)
1743 std::string file_name =
"cvg-scm-"+ error_name +
"-stats.dat";
1745 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1747 conv.open(file_name, std::ios::app);
1748 conv <<
"Nb_basis" <<
"\t" <<
"Min" <<
"\t" <<
"Max" <<
"\t" <<
"Mean" <<
"\t" <<
"Variance" <<
"\n";
1751 for(
int j=0; j<N; j++)
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;
1758 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1761 << M_mapConvSCM[error_name][j].minCoeff() <<
"\t"
1762 << M_mapConvSCM[error_name][j].maxCoeff() <<
"\t"
1763 << mean <<
"\t" << variance <<
"\n";
1772 void buildSamplingFromCfg()
1775 int mu_size = model->parameterSpace()->dimension();
1778 fs::path input_file (
"SamplingForTest");
1779 if( fs::exists(input_file) )
1780 std::remove(
"SamplingForTest" );
1782 std::ofstream input(
"SamplingForTest" );
1786 std::ifstream cfg_file( option(_name=
"config-file").
template as<std::string>() );
1788 std::cout <<
"[Script-mode] Config file cannot be read" << std::endl;
1791 for(
int i=1; i<=mu_size; i++)
1794 std::ostringstream oss;
1796 std::string is = oss.str();
1799 std::string cfg_line_mu, tmp_content;
1800 std::ifstream cfg_file( option(_name=
"config-file").
template as<std::string>() );
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;
1809 std::string expr_s =
"mu"+is+
"[[:space:]]*=[[:space:]]*([0-9]+(.?)[0-9]*(e(\\+|-)[0-9]+)?)[[:space:]]*";
1810 boost::regex expression( expr_s );
1814 auto is_match = boost::regex_match(cfg_line_mu, what, expression);
1822 input << what[1] <<
" , ";
1824 input << what[1] <<
" ]";
1830 CRBModelMode M_mode;
1831 crbmodel_ptrtype model;
1835 std::map<std::string, std::vector<vectorN_type> > M_mapConvEIM;
1837 std::map<std::string, std::vector<vectorN_type> > M_mapConvCRB;
1839 std::map<std::string, std::vector<vectorN_type> > M_mapConvSCM;
1842 std::vector< sampling_ptrtype > vector_sampling_for_primal_efficiency_under_1;
1843 std::vector< sampling_ptrtype > vector_sampling_for_dual_efficiency_under_1;
1845 fs::path M_current_path;