33 #include <boost/multi_array.hpp>
34 #include <boost/tuple/tuple.hpp>
35 #include "boost/tuple/tuple_io.hpp"
36 #include <boost/format.hpp>
37 #include <boost/foreach.hpp>
38 #include <boost/bimap.hpp>
39 #include <boost/bimap/support/lambda.hpp>
40 #include <boost/archive/text_oarchive.hpp>
41 #include <boost/archive/text_iarchive.hpp>
42 #include <boost/math/special_functions/fpclassify.hpp>
46 #include <boost/serialization/vector.hpp>
47 #include <boost/serialization/list.hpp>
48 #include <boost/serialization/string.hpp>
49 #include <boost/serialization/version.hpp>
50 #include <boost/serialization/split_member.hpp>
56 #include <Eigen/Dense>
59 #include <feel/feelcore/feel.hpp>
61 #include <feel/feelcore/parameter.hpp>
64 #include <feel/feelcrb/crbscm.hpp>
67 #include <feel/feeldiscr/bdf2.hpp>
68 #include <feel/feelfilters/exporter.hpp>
83 enum CRBErrorType { CRB_RESIDUAL = 0, CRB_RESIDUAL_SCM=1, CRB_NO_RESIDUAL=2 , CRB_EMPIRICAL=3};
94 template<
typename TruthModelType>
112 typedef TruthModelType truth_model_type;
113 typedef truth_model_type model_type;
114 typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
116 typedef double value_type;
117 typedef boost::tuple<double,double> bounds_type;
120 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
122 typedef typename parameterspace_type::element_ptrtype parameter_ptrtype;
124 typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
126 typedef boost::tuple<double, parameter_type, size_type, double, double> relative_error_type;
127 typedef relative_error_type max_error_type;
129 typedef boost::tuple<double, std::vector< std::vector<double> > , std::vector< std::vector<double> >, double,
double > error_estimation_type;
130 typedef boost::tuple<double, std::vector<double> > residual_error_type;
132 typedef boost::bimap< int, boost::tuple<double,double,double> > convergence_type;
134 typedef typename convergence_type::value_type convergence;
140 typedef boost::shared_ptr<scm_type> scm_ptrtype;
144 typedef boost::shared_ptr<crb_elements_db_type> crb_elements_db_ptrtype;
148 typedef boost::shared_ptr<pod_type> pod_ptrtype;
152 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
157 typedef typename model_type::element_ptrtype element_ptrtype;
159 typedef typename model_type::backend_type backend_type;
160 typedef boost::shared_ptr<backend_type> backend_ptrtype;
161 typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
162 typedef typename model_type::vector_ptrtype vector_ptrtype;
163 typedef typename model_type::beta_vector_type beta_vector_type;
166 typedef Eigen::VectorXd y_type;
167 typedef std::vector<y_type> y_set_type;
168 typedef std::vector<boost::tuple<double,double> > y_bounds_type;
170 typedef std::vector<element_type> wn_type;
171 typedef boost::tuple< std::vector<wn_type> , std::vector<std::string> > export_vector_wn_type;
173 typedef std::vector<double> vector_double_type;
174 typedef boost::shared_ptr<vector_double_type> vector_double_ptrtype;
176 typedef Eigen::VectorXd vectorN_type;
177 typedef Eigen::MatrixXd matrixN_type;
179 typedef Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> > map_dense_matrix_type;
180 typedef Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, 1> > map_dense_vector_type;
182 typedef boost::tuple< std::vector<vectorN_type> , std::vector<vectorN_type> , std::vector<vectorN_type>, std::vector<vectorN_type> > solutions_tuple;
183 typedef boost::tuple< double,double,double , std::vector< std::vector< double > > , std::vector< std::vector< double > > > upper_bounds_tuple;
185 typedef std::vector<element_type> mode_set_type;
186 typedef boost::shared_ptr<mode_set_type> mode_set_ptrtype;
188 typedef boost::multi_array<value_type, 2> array_2_type;
189 typedef boost::multi_array<vectorN_type, 2> array_3_type;
190 typedef boost::multi_array<matrixN_type, 2> array_4_type;
194 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
198 typedef boost::shared_ptr<space_type> space_ptrtype;
202 typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
206 typedef boost::shared_ptr<export_type> export_ptrtype;
209 typedef boost::shared_ptr<preconditioner_type> preconditioner_ptrtype;
213 static const int nb_spaces = functionspace_type::nSpaces;
214 typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
215 typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
216 typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<4> >, fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3> >,
217 fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
218 >::type >::type >::type index_vector_type;
231 M_elements_database(),
232 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
238 M_error_type( CRB_NO_RESIDUAL ),
250 po::variables_map
const&
vm )
252 super( ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
254 ( boost::format(
"%1%-%2%-%3%" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
257 ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
259 ( boost::format(
"%1%-%2%-%3%-elements" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
261 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
263 M_backend( backend_type::build( vm ) ),
264 M_backend_primal( backend_type::build( vm ) ),
265 M_backend_dual( backend_type::build( vm ) ),
266 M_output_index( vm[
"crb.output-index"].template as<int>() ),
267 M_tolerance( vm[
"crb.error-max"].template as<double>() ),
268 M_iter_max( vm[
"crb.dimension-max"].template as<int>() ),
269 M_factor( vm[
"crb.factor"].template as<int>() ),
270 M_error_type(
CRBErrorType( vm[
"crb.error-type"].template as<int>() ) ),
275 M_scmA( new
scm_type( name+
"_a", vm ) ),
276 M_scmM( new
scm_type( name+
"_m", vm ) ),
288 po::variables_map
const&
vm,
289 truth_model_ptrtype
const & model )
291 super( ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
293 ( boost::format(
"%1%-%2%-%3%" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
296 ( boost::format(
"%1%" ) % vm[
"crb.error-type"].template as<int>() ).str(),
298 ( boost::format(
"%1%-%2%-%3%-elements" ) % name % vm[
"crb.output-index"].template as<int>() % vm[
"crb.error-type"].template as<int>() ).str(),
301 M_nlsolver(
SolverNonLinear<double>::build( SOLVERS_PETSC, Environment::worldComm() ) ),
303 M_backend( backend_type::build( vm ) ),
304 M_backend_primal( backend_type::build( vm ) ),
305 M_backend_dual( backend_type::build( vm ) ),
306 M_output_index( vm[
"crb.output-index"].template as<int>() ),
307 M_tolerance( vm[
"crb.error-max"].template as<double>() ),
308 M_iter_max( vm[
"crb.dimension-max"].template as<int>() ),
309 M_factor( vm[
"crb.factor"].template as<int>() ),
310 M_error_type(
CRBErrorType( vm[
"crb.error-type"].template as<int>() ) ),
315 M_scmA( new
scm_type( name+
"_a", vm , model , false ) ),
316 M_scmM( new
scm_type( name+
"_m", vm , model , true ) ),
321 LOG(INFO) <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
322 M_elements_database.setMN( M_N );
323 if( M_elements_database.loadDB() )
324 LOG(INFO) <<
"Element Database " << M_elements_database.lookForDB() <<
" available and loaded\n";
325 auto basis_functions = M_elements_database.wn();
326 M_WN = basis_functions.template get<0>();
327 M_WNdu = basis_functions.template get<1>();
329 M_preconditioner_primal = preconditioner(_pc=(
PreconditionerType) M_backend_primal->pcEnumType(),
330 _backend= M_backend_primal,
331 _pcfactormatsolverpackage=(MatSolverPackageType) M_backend_primal->matSolverPackageEnumType(),
332 _worldcomm=M_backend_primal->comm(),
333 _prefix=M_backend_primal->prefix() ,
335 M_preconditioner_dual = preconditioner(_pc=(
PreconditionerType) M_backend_dual->pcEnumType(),
336 _backend= M_backend_dual,
337 _pcfactormatsolverpackage=(MatSolverPackageType) M_backend_dual->matSolverPackageEnumType(),
338 _worldcomm=M_backend_dual->comm(),
339 _prefix=M_backend_dual->prefix() ,
348 M_elements_database( o.M_elements_database ),
349 M_nlsolver( o.M_nlsolver ),
350 M_output_index( o.M_output_index ),
351 M_tolerance( o.M_tolerance ),
352 M_iter_max( o.M_iter_max ),
353 M_factor( o.M_factor ),
354 M_error_type( o.M_error_type ),
355 M_maxerror( o.M_maxerror ),
359 M_WNmu_complement( o.M_WNmu_complement ),
360 M_C0_pr( o.M_C0_pr ),
361 M_C0_du( o.M_C0_du ),
362 M_Lambda_pr( o.M_Lambda_pr ),
363 M_Lambda_du( o.M_Lambda_du ),
364 M_Gamma_pr( o.M_Gamma_pr ),
365 M_Gamma_du( o.M_Gamma_du ),
366 M_Cmf_pr( o.M_Cmf_pr ),
367 M_Cmf_du( o.M_Cmf_du ),
368 M_Cmf_du_ini( o.M_Cmf_du_ini ),
369 M_Cma_pr( o.M_Cma_pr ),
370 M_Cma_du( o.M_Cma_du ),
371 M_Cmm_pr( o.M_Cmm_pr ),
372 M_Cmm_du( o.M_Cmm_du ),
373 M_coeff_pr_ini_online( o.M_coeff_pr_ini_online ),
374 M_coeff_du_ini_online( o.M_coeff_du_ini_online )
416 parameterspace_ptrtype
Dmu()
const
424 return M_output_index;
462 int M_prev_o = M_output_index;
463 M_output_index = oindex;
465 if ( (
size_type )M_output_index >= M_model->Nl() )
466 M_output_index = M_prev_o;
469 this->
setDBFilename( ( boost::format(
"%1%-%2%-%3%.crbdb" ) % this->
name() % M_output_index % M_error_type ).str() );
471 if ( M_output_index != M_prev_o )
481 M_error_type = error;
487 M_tolerance = tolerance;
494 M_Dmu = M_model->parameterSpace();
501 LOG(INFO) <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
529 ComputePhi( wn_type v)
535 template<
typename T >
537 operator()(
const T& t )
const
539 mesh_ptrtype mesh = t.functionSpace()->mesh();
540 double surface =
integrate( _range=
elements(mesh), _expr=vf::cst(1.) ).evaluate()(0,0);
541 double mean =
integrate( _range=
elements(mesh), _expr=vf::idv( t ) ).evaluate()(0,0);
545 int first_dof = t.start();
546 int nb_dof = t.functionSpace()->nLocalDof();
547 for(
int dof=0 ; dof<nb_dof; dof++)
548 M_vect[M_curs][first_dof + dof] = element_t[dof] - element_mean[dof];
558 mutable wn_type M_vect;
562 struct ComputeIntegrals
568 M_composite_e1 ( composite_e1 ),
569 M_composite_e2 ( composite_e2 )
572 template<
typename T >
574 operator()(
const T& t )
const
576 auto e1 = M_composite_e1.template element< T::value >();
577 auto e2 = M_composite_e2.template element< T::value >();
578 mesh_ptrtype mesh = e1.functionSpace()->mesh();
579 double integral =
integrate( _range=
elements(mesh) , _expr=trans( vf::idv( e1 ) ) * vf::idv( e2 ) ).evaluate()(0,0);
581 M_vect.push_back( integral );
584 std::vector< double > vectorIntegrals()
589 mutable std::vector< double > M_vect;
594 struct ComputeIntegralsSquare
597 ComputeIntegralsSquare(
element_type const composite_e1 ,
600 M_composite_e1 ( composite_e1 ),
601 M_composite_e2 ( composite_e2 )
604 template<
typename T >
606 operator()(
const T& t )
const
612 M_error.conservativeResize( i+1 );
614 auto e1 = M_composite_e1.template element< T::value >();
615 auto e2 = M_composite_e2.template element< T::value >();
617 auto Xh = M_composite_e1.functionSpace();
618 mesh_ptrtype mesh = Xh->mesh();
620 _expr=trans( vf::idv( e1 ) - vf::idv( e2 ) )
621 * ( vf::idv( e1 ) - vf::idv( e2 ) )
623 M_error(i) = math::sqrt( integral ) ;
626 vectorN_type vectorErrors()
631 mutable vectorN_type M_error;
648 void checkResidual( parameter_type
const& mu, std::vector< std::vector<double> >
const& primal_residual_coeffs,
651 void compareResidualsForTransientProblems(
int N, parameter_type
const& mu, std::vector<element_type>
const & Un,
652 std::vector<element_type>
const & Unold, std::vector<element_type>
const& Undu,
653 std::vector<element_type>
const & Unduold,
654 std::vector< std::vector<double> >
const& primal_residual_coeffs,
655 std::vector< std::vector<double> >
const& dual_residual_coeffs )
const ;
658 void buildFunctionFromRbCoefficients(
int N, std::vector< vectorN_type >
const & RBcoeff, wn_type
const & WN, std::vector<element_ptrtype> & FEMsolutions );
664 double checkOrthonormality(
int N,
const wn_type& wn )
const;
705 boost::tuple<double,double>
lb(
size_type N, parameter_type
const& mu, std::vector< vectorN_type >& uN, std::vector< vectorN_type >& uNdu , std::vector<vectorN_type> & uNold, std::vector<vectorN_type> & uNduold,
int K=0 )
const;
714 void updateJacobian(
const map_dense_vector_type& map_X, map_dense_matrix_type& map_J ,
const parameter_type & mu ,
int N)
const ;
723 void updateResidual(
const map_dense_vector_type& map_X, map_dense_vector_type& map_R ,
const parameter_type & mu,
int N )
const ;
731 void computeProjectionInitialGuess(
const parameter_type & mu,
int N , vectorN_type& initial_guess )
const ;
740 void newton(
size_type N, parameter_type
const& mu , vectorN_type & uN ,
double& condition_number,
double& output)
const ;
748 element_type offlineFixedPointPrimal( parameter_type
const& mu, sparse_matrix_ptrtype & A,
bool zero_iteration ) ;
757 element_type offlineFixedPointDual( parameter_type
const& mu , element_ptrtype & dual_initial_field,
const sparse_matrix_ptrtype & A,
const element_type & u,
bool zero_iteration ) ;
769 void fixedPointPrimal(
size_type N, parameter_type
const& mu, std::vector< vectorN_type > & uN, std::vector<vectorN_type> & uNold,
double& condition_number,
770 std::vector< double > & output_vector,
int K=0)
const;
781 void fixedPointDual(
size_type N, parameter_type
const& mu, std::vector< vectorN_type > & uNdu, std::vector<vectorN_type> & uNduold, std::vector< double > & output_vector,
int K=0)
const;
794 void fixedPoint(
size_type N, parameter_type
const& mu, std::vector< vectorN_type > & uN, std::vector< vectorN_type > & uNdu, std::vector<vectorN_type> & uNold, std::vector<vectorN_type> & uNduold,
795 double& condition_number, std::vector< double > & output_vector,
int K=0)
const;
814 boost::tuple<double,double>
lb( parameter_ptrtype
const& mu,
size_type N, std::vector< vectorN_type >& uN, std::vector< vectorN_type >& uNdu )
const
816 return lb( N, *mu, uN, uNdu );
833 auto o =
lb( N, mu, uN, uNdu );
834 auto e =
delta( N, mu, uN, uNdu );
835 return o.template get<0>() + e.template get<0>();
848 value_type
delta(
size_type N, parameter_ptrtype
const& mu, std::vector< vectorN_type >
const& uN, std::vector< vectorN_type >
const& uNdu, std::vector<vectorN_type>
const& uNold, std::vector<vectorN_type>
const& uNduold ,
int k=0 )
const
850 auto e =
delta( N, mu, uN, uNdu );
851 return e.template get<0>();
864 error_estimation_type
delta(
size_type N, parameter_type
const& mu, std::vector< vectorN_type >
const& uN, std::vector< vectorN_type >
const& uNdu, std::vector<vectorN_type>
const& uNold, std::vector<vectorN_type>
const& uNduold,
int k=0 )
const;
874 value_type
ub(
size_type K, parameter_ptrtype
const& mu, std::vector< vectorN_type >& uN, std::vector< vectorN_type >& uNdu )
const
876 return ub( K, *mu, uN, uNdu );
897 residual_error_type
transientPrimalResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Un, vectorN_type
const& Unold=vectorN_type(),
double time_step=1,
double time=1e30 )
const;
898 residual_error_type steadyPrimalResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Un,
double time=0 )
const;
901 residual_error_type transientDualResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Un, vectorN_type
const& Unold=vectorN_type(),
double time_step=1,
double time=1e30 )
const;
902 residual_error_type steadyDualResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Un,
double time=0 )
const;
905 value_type initialDualResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Uduini,
double time_step )
const ;
912 void offlineResidual(
int Ncur, mpl::bool_<true> ,
int number_of_added_elements=1 );
913 void offlineResidual(
int Ncur, mpl::bool_<false> ,
int number_of_added_elements=1 );
921 value_type empiricalErrorEstimation (
int Ncur, parameter_type
const& mu,
int k )
const ;
940 void checkInitialGuess(
const element_type expansion_uN , parameter_type
const& mu, vectorN_type & error )
const ;
941 void checkInitialGuess(
const element_type expansion_uN , parameter_type
const& mu, vectorN_type & error , mpl::bool_<true> )
const ;
942 void checkInitialGuess(
const element_type expansion_uN , parameter_type
const& mu, vectorN_type & error , mpl::bool_<false> )
const ;
950 boost::tuple<double,double, solutions_tuple, double, double, double, upper_bounds_tuple >
run( parameter_type
const& mu,
double eps = 1e-6,
int N = -1 );
955 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
956 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P, mpl::bool_<true> );
957 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P, mpl::bool_<false> );
964 M_Xi->randomize( N );
973 M_Xi->equidistribute( N );
977 sampling_ptrtype wnmu ( )
const
1013 bool use = this->
vm()[
"crb.run-on-WNmu"].template as<bool>();
1051 double correctionTerms(parameter_type
const& mu, std::vector< vectorN_type >
const & uN, std::vector< vectorN_type >
const & uNdu , std::vector<vectorN_type>
const & uNold,
int const K=0)
const;
1057 void buildVarianceMatrixPhi(
int const N);
1058 void buildVarianceMatrixPhi(
int const N , mpl::bool_<true> );
1059 void buildVarianceMatrixPhi(
int const N , mpl::bool_<false> );
1061 WorldComm
const& worldComm()
const {
return Environment::worldComm() ; }
1074 boost::shared_ptr<SolverNonLinear<double> > M_nlsolver;
1076 truth_model_ptrtype M_model;
1078 backend_ptrtype M_backend;
1079 backend_ptrtype M_backend_primal;
1080 backend_ptrtype M_backend_dual;
1095 parameterspace_ptrtype M_Dmu;
1098 sampling_ptrtype M_Xi;
1101 sampling_ptrtype M_WNmu;
1102 sampling_ptrtype M_WNmu_complement;
1109 export_ptrtype exporter;
1112 array_2_type M_C0_pr;
1113 array_2_type M_C0_du;
1114 array_3_type M_Lambda_pr;
1115 array_3_type M_Lambda_du;
1116 array_4_type M_Gamma_pr;
1117 array_4_type M_Gamma_du;
1118 array_3_type M_Cmf_pr;
1119 array_3_type M_Cmf_du;
1120 array_4_type M_Cma_pr;
1121 array_4_type M_Cma_du;
1122 array_4_type M_Cmm_pr;
1123 array_4_type M_Cmm_du;
1125 std::vector< std::vector< std::vector< std::vector< double > > > > M_C0_pr;
1126 std::vector< std::vector< std::vector< std::vector< double > > > > M_C0_du;
1127 std::vector< std::vector< std::vector< std::vector< vectorN_type > > > > M_Lambda_pr;
1128 std::vector< std::vector< std::vector< std::vector< vectorN_type > > > > M_Lambda_du;
1129 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Gamma_pr;
1130 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Gamma_du;
1131 std::vector< std::vector< std::vector< std::vector< vectorN_type > > > > M_Cmf_pr;
1132 std::vector< std::vector< std::vector< std::vector< vectorN_type > > > > M_Cmf_du;
1133 std::vector< std::vector< std::vector< std::vector< vectorN_type > > > > M_Cmf_du_ini;
1134 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Cma_pr;
1135 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Cma_du;
1136 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Cmm_pr;
1137 std::vector< std::vector< std::vector< std::vector< matrixN_type > > > > M_Cmm_du;
1140 std::vector<double> M_coeff_pr_ini_online;
1141 std::vector<double> M_coeff_du_ini_online;
1144 friend class boost::serialization::access;
1148 template<
class Archive>
1149 void save( Archive & ar,
const unsigned int version )
const;
1151 template<
class Archive>
1152 void load( Archive & ar,
const unsigned int version ) ;
1154 BOOST_SERIALIZATION_SPLIT_MEMBER()
1163 bool orthonormalize_primal;
1164 bool orthonormalize_dual;
1165 bool solve_dual_problem;
1167 convergence_type M_rbconv;
1170 bdf_ptrtype M_bdf_primal;
1171 bdf_ptrtype M_bdf_primal_save;
1172 bdf_ptrtype M_bdf_dual;
1173 bdf_ptrtype M_bdf_dual_save;
1176 std::vector < std::vector<matrixN_type> > M_Aqm_pr;
1177 std::vector < std::vector<matrixN_type> > M_Aqm_du;
1178 std::vector < std::vector<matrixN_type> > M_Aqm_pr_du;
1181 std::vector < std::vector<matrixN_type> > M_Jqm_pr;
1184 std::vector < std::vector<vectorN_type> > M_Rqm_pr;
1187 std::vector < std::vector<matrixN_type> > M_Mqm_pr;
1188 std::vector < std::vector<matrixN_type> > M_Mqm_du;
1189 std::vector < std::vector<matrixN_type> > M_Mqm_pr_du;
1192 std::vector < std::vector<vectorN_type> > M_Fqm_pr;
1193 std::vector < std::vector<vectorN_type> > M_Fqm_du;
1195 std::vector < std::vector<vectorN_type> > M_Lqm_pr;
1196 std::vector < std::vector<vectorN_type> > M_Lqm_du;
1199 std::vector < std::vector<vectorN_type> > M_InitialGuessV_pr;
1201 std::vector<
int> M_index;
1204 std::vector < matrixN_type > M_variance_matrix_phi;
1206 bool M_compute_variance;
1207 bool M_rbconv_contains_primal_and_dual_contributions;
1208 parameter_type M_current_mu;
1209 int M_no_residual_index;
1211 bool M_database_contains_variance_info;
1215 bool M_offline_step;
1217 preconditioner_ptrtype M_preconditioner;
1218 preconditioner_ptrtype M_preconditioner_primal;
1219 preconditioner_ptrtype M_preconditioner_dual;
1223 po::options_description crbOptions( std::
string const& prefix = "" );
1228 template<typename TruthModelType>
1230 CRB<TruthModelType>::offlineFixedPointPrimal(parameter_type const& mu, sparse_matrix_ptrtype & A,
bool zero_iteration )
1232 auto u = M_model->functionSpace()->element();
1234 sparse_matrix_ptrtype M = M_model->newMatrix();
1235 int nl = M_model->Nl();
1236 std::vector< vector_ptrtype > F( nl );
1237 for(
int l=0; l<nl; l++)
1238 F[l]=M_model->newVector();
1241 bool reuse_prec = this->
vm()[
"crb.reuse-prec"].template as<bool>() ;
1243 M_bdf_primal = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_primal" );
1244 M_bdf_primal_save = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_primal_save" );
1247 M_bdf_primal->setTimeInitial( M_model->timeInitial() );
1248 M_bdf_primal->setTimeStep( M_model->timeStep() );
1249 M_bdf_primal->setTimeFinal( M_model->timeFinal() );
1250 M_bdf_primal->setOrder( M_model->timeOrder() );
1252 M_bdf_primal_save->setTimeInitial( M_model->timeInitial() );
1253 M_bdf_primal_save->setTimeStep( M_model->timeStep() );
1254 M_bdf_primal_save->setTimeFinal( M_model->timeFinal() );
1255 M_bdf_primal_save->setOrder( M_model->timeOrder() );
1257 M_bdf_primal_save->setRankProcInNameOfFiles(
true );
1258 M_bdf_primal->setRankProcInNameOfFiles(
true );
1261 auto elementptr = M_model->functionSpace()->elementPtr();
1262 M_model->initializationField( elementptr, mu );
1266 auto Apr = M_model->newMatrix();
1269 int max_fixedpoint_iterations = option(_name=
"crb.max-fixedpoint-iterations").template as<int>();
1270 double increment_fixedpoint_tol = option(_name=
"crb.increment-fixedpoint-tol").template as<double>();
1271 double fixedpoint_critical_value = option(_name=
"crb.fixedpoint-critical-value").template as<double>();
1273 double increment_norm=1e3;
1275 if( zero_iteration )
1280 auto vec_bdf_poly = M_backend_primal->newVector( M_model->functionSpace() );
1283 if ( M_model->isSteady() )
1285 elementptr = M_model->assembleInitialGuess( mu ) ;
1289 auto uold = M_model->functionSpace()->element();
1291 element_ptrtype uproj(
new element_type( M_model->functionSpace() ) );
1293 vector_ptrtype Rhs( M_backend_primal->newVector( M_model->functionSpace() ) );
1295 for ( M_bdf_primal->start(u),M_bdf_primal_save->start(u);
1296 !M_bdf_primal->isFinished() , !M_bdf_primal_save->isFinished();
1297 M_bdf_primal->next() , M_bdf_primal_save->next() )
1300 bdf_coeff = M_bdf_primal->polyDerivCoefficient( 0 );
1302 auto bdf_poly = M_bdf_primal->polyDeriv();
1306 boost::tie( M, Apr, F) = M_model->update( mu , u, M_bdf_primal->time() );
1309 if( iteration == 0 )
1312 if ( ! M_model->isSteady() )
1314 Apr->addMatrix( bdf_coeff, M );
1316 *vec_bdf_poly = bdf_poly;
1317 Rhs->addVector( *vec_bdf_poly, *M );
1327 M_preconditioner_primal->setMatrix( Apr );
1330 auto ret = M_backend_primal->solve( _matrix=Apr, _solution=u, _rhs=Rhs, _prec=M_preconditioner_primal, _reuse_prec=( M_bdf_primal->iteration() >=2 ) );
1331 if ( !ret.template get<0>() )
1332 LOG(INFO)<<
"[CRB] WARNING : at time "<<M_bdf_primal->time()<<
" we have not converged ( nb_it : "<<ret.template get<1>()<<
" and residual : "<<ret.template get<2>() <<
" ) \n";
1337 auto ret = M_backend_primal->solve( _matrix=Apr, _solution=u, _rhs=Rhs , _prec=M_preconditioner_primal );
1338 if ( !ret.template get<0>() )
1339 LOG(INFO)<<
"[CRB] WARNING : at time "<<M_bdf_primal->time()<<
" we have not converged ( nb_it : "<<ret.template get<1>()<<
" and residual : "<<ret.template get<2>() <<
" ) \n";
1343 increment_norm = M_model->computeNormL2( u , uold );
1346 }
while( increment_norm > increment_fixedpoint_tol && iteration < max_fixedpoint_iterations );
1348 M_bdf_primal->shiftRight( u );
1350 if ( ! M_model->isSteady() )
1352 element_ptrtype projection (
new element_type ( M_model->functionSpace() ) );
1355 M_bdf_primal_save->shiftRight( *uproj );
1358 if( increment_norm > fixedpoint_critical_value )
1359 throw std::logic_error( (boost::format(
"[CRB::offlineFixedPointPrimal] at time %1% ERROR : increment > critical value " ) %M_bdf_primal->time() ).str() );
1361 for (
size_type l = 0; l < M_model->Nl(); ++l )
1364 element_ptrtype eltF(
new element_type( M_model->functionSpace() ) );
1366 LOG(INFO) <<
"u^T F[" << l <<
"]= " <<
inner_product( u, *eltF ) <<
" at time : "<<M_bdf_primal->time()<<
"\n";
1368 LOG(INFO) <<
"[CRB::offlineWithErrorEstimation] energy = " << A->energy( u, u ) <<
"\n";
1379 auto merge = M_model->update( mu );
1380 A = merge.template get<1>();
1381 F = merge.template get<2>();
1387 M_preconditioner_primal->setMatrix( A );
1388 M_backend_primal->solve( _matrix=A , _solution=un, _rhs=F[0] , _prec=M_preconditioner );
1391 increment_norm = M_model->computeNormL2( un , uold );
1394 }
while( increment_norm > increment_fixedpoint_tol && iteration < max_fixedpoint_iterations );
1396 element_ptrtype eltF(
new element_type( M_model->functionSpace() ) );
1398 LOG(INFO) <<
"[CRB::offlineFixedPoint] u^T F = " <<
inner_product( *u, *eltF ) ;
1399 LOG(INFO) <<
"[CRB::offlineFixedPoint] energy = " << A->energy( *u, *u ) ;
1404 template<
typename TruthModelType>
1405 typename CRB<TruthModelType>::element_type
1406 CRB<TruthModelType>::offlineFixedPointDual(parameter_type
const& mu, element_ptrtype & dual_initial_field,
const sparse_matrix_ptrtype & A,
const element_type & u,
bool zero_iteration )
1410 bool reuse_prec = this->
vm()[
"crb.reuse-prec"].template as<bool>() ;
1412 auto udu = M_model->functionSpace()->element();
1414 sparse_matrix_ptrtype M = M_model->newMatrix();
1415 int nl = M_model->Nl();
1416 std::vector< vector_ptrtype > F( nl );
1417 for(
int l=0; l<nl; l++)
1418 F[l]=M_model->newVector();
1421 M_bdf_dual = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_dual" );
1422 M_bdf_dual_save = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_dual_save" );
1424 M_bdf_dual->setTimeInitial( M_model->timeFinal()+M_model->timeStep() );
1426 M_bdf_dual->setTimeStep( -M_model->timeStep() );
1427 M_bdf_dual->setTimeFinal( M_model->timeInitial()+M_model->timeStep() );
1428 M_bdf_dual->setOrder( M_model->timeOrder() );
1430 M_bdf_dual_save->setTimeInitial( M_model->timeFinal()+M_model->timeStep() );
1431 M_bdf_dual_save->setTimeStep( -M_model->timeStep() );
1432 M_bdf_dual_save->setTimeFinal( M_model->timeInitial()+M_model->timeStep() );
1433 M_bdf_dual_save->setOrder( M_model->timeOrder() );
1435 M_bdf_dual_save->setRankProcInNameOfFiles(
true );
1436 M_bdf_dual->setRankProcInNameOfFiles(
true );
1438 auto Adu = M_model->newMatrix();
1439 auto Apr = M_model->newMatrix();
1441 double dt = M_model->timeStep();
1443 int max_fixedpoint_iterations = option(_name=
"crb.max-fixedpoint-iterations").template as<int>();
1444 double increment_fixedpoint_tol = option(_name=
"crb.increment-fixedpoint-tol").template as<double>();
1445 double fixedpoint_critical_value = option(_name=
"crb.fixedpoint-critical-value").template as<double>();
1447 double increment_norm=1e3;
1449 vector_ptrtype Rhs( M_backend_dual->newVector( M_model->functionSpace() ) );
1451 if( zero_iteration )
1456 auto vec_bdf_poly = M_backend_dual->newVector( M_model->functionSpace() );
1458 if ( M_model->isSteady() )
1462 boost::tie( M, Apr, F) = M_model->update( mu , M_bdf_dual->timeInitial() );
1465 Apr->addMatrix( 1./dt, M );
1466 Apr->transpose( Adu );
1467 *Rhs=*F[M_output_index];
1468 Rhs->scale( 1./dt );
1471 M_preconditioner_dual->setMatrix( Adu );
1472 M_backend_dual->solve( _matrix=Adu, _solution=dual_initial_field, _rhs=Rhs, _prec=M_preconditioner_dual );
1474 *Rhs=*F[M_output_index];
1477 M_preconditioner_dual->setMatrix( M );
1478 M_backend_dual->solve( _matrix=M, _solution=dual_initial_field, _rhs=Rhs, _prec=M_preconditioner_dual );
1480 udu=*dual_initial_field;
1483 auto uold = M_model->functionSpace()->element();
1485 element_ptrtype uproj(
new element_type( M_model->functionSpace() ) );
1488 for ( M_bdf_dual->start(udu),M_bdf_dual_save->start(udu);
1489 !M_bdf_dual->isFinished() , !M_bdf_dual_save->isFinished();
1490 M_bdf_dual->next() , M_bdf_dual_save->next() )
1493 bdf_coeff = M_bdf_dual->polyDerivCoefficient( 0 );
1495 auto bdf_poly = M_bdf_dual->polyDeriv();
1499 boost::tie( M, Apr, F) = M_model->update( mu , udu, M_bdf_dual->time() );
1501 if( ! M_model->isSteady() )
1503 Apr->addMatrix( bdf_coeff, M );
1505 *vec_bdf_poly = bdf_poly;
1506 Rhs->addVector( *vec_bdf_poly, *M );
1510 *Rhs = *F[M_output_index];
1515 if( option(
"crb.use-symmetric-matrix").
template as<bool>() )
1518 Apr->transpose( Adu );
1526 M_preconditioner_dual->setMatrix( Adu );
1529 auto ret = M_backend_dual->solve( _matrix=Adu, _solution=udu, _rhs=Rhs, _prec=M_preconditioner_dual, _reuse_prec=( M_bdf_primal->iteration() >=2 ) );
1530 if ( !ret.template get<0>() )
1531 LOG(INFO)<<
"[CRB] WARNING : at time "<<M_bdf_primal->time()<<
" we have not converged ( nb_it : "<<ret.template get<1>()<<
" and residual : "<<ret.template get<2>() <<
" ) \n";
1535 auto ret = M_backend_dual->solve( _matrix=Adu, _solution=udu, _rhs=Rhs , _prec=M_preconditioner_dual );
1536 if ( !ret.template get<0>() )
1537 LOG(INFO)<<
"[CRB] WARNING : at time "<<M_bdf_primal->time()<<
" we have not converged ( nb_it : "<<ret.template get<1>()<<
" and residual : "<<ret.template get<2>() <<
" ) \n";
1541 increment_norm = M_model->computeNormL2( udu , uold );
1544 }
while( increment_norm > increment_fixedpoint_tol && iteration < max_fixedpoint_iterations );
1546 M_bdf_dual->shiftRight( udu );
1549 double term1 = A->energy( udu, u );
1550 double term2 = Adu->energy( u, udu );
1551 double diff = math::abs( term1-term2 );
1552 LOG(INFO) <<
"< A u , udu > - < u , A* udu > = "<<diff<<
"\n";
1554 if ( ! M_model->isSteady() )
1556 element_ptrtype projection (
new element_type ( M_model->functionSpace() ) );
1559 M_bdf_dual_save->shiftRight( *uproj );
1562 if( increment_norm > fixedpoint_critical_value )
1563 throw std::logic_error( (boost::format(
"[CRB::offlineFixedPointDual] at time %1% ERROR : increment > critical value " ) %M_bdf_dual->time() ).str() );
1572 template<
typename TruthModelType>
1574 CRB<TruthModelType>::loadSCMDB()
1576 if ( M_error_type == CRB_RESIDUAL_SCM )
1578 M_scmA->setScmForMassMatrix(
false );
1579 if( ! M_scmA->loadDB() )
1580 std::vector<boost::tuple<double,double,double> > M_rbconv2 = M_scmA->offline();
1582 if ( ! M_model->isSteady() )
1584 M_scmM->setScmForMassMatrix(
true );
1585 if( ! M_scmM->loadDB() )
1586 std::vector<boost::tuple<double,double,double> > M_rbconv3 = M_scmM->offline();
1592 template<
typename TruthModelType>
1593 typename CRB<TruthModelType>::convergence_type
1596 M_database_contains_variance_info = this->
vm()[
"crb.save-information-for-variance"].template as<bool>();
1598 int proc_number = this->worldComm().globalRank();
1599 M_rbconv_contains_primal_and_dual_contributions =
true;
1601 bool rebuild_database = this->
vm()[
"crb.rebuild-database"].template as<bool>() ;
1602 orthonormalize_primal = this->
vm()[
"crb.orthonormalize-primal"].template as<bool>() ;
1603 orthonormalize_dual = this->
vm()[
"crb.orthonormalize-dual"].template as<bool>() ;
1604 solve_dual_problem = this->
vm()[
"crb.solve-dual-problem"].template as<bool>() ;
1606 M_use_newton = this->
vm()[
"crb.use-newton"].template as<bool>() ;
1611 if ( ! solve_dual_problem ) orthonormalize_dual=
false;
1613 M_Nm = this->
vm()[
"crb.Nm"].template as<int>() ;
1614 bool seek_mu_in_complement = this->
vm()[
"crb.seek-mu-in-complement"].template as<bool>() ;
1617 LOG(INFO)<<
"Offline CRB starts, this may take a while until Database is computed...\n";
1630 if ( rebuild_database || M_N == 0)
1638 LOG(INFO) <<
"[CRB::offline] compute random sampling\n";
1640 int sampling_size = this->
vm()[
"crb.sampling-size"].template as<int>();
1641 std::string file_name = ( boost::format(
"M_Xi_%1%") % sampling_size ).str();
1642 std::ifstream file ( file_name );
1646 M_Xi->randomize( sampling_size );
1648 M_Xi->writeOnFile(file_name);
1653 M_Xi->readFromFile(file_name);
1656 M_WNmu->setSuperSampling( M_Xi );
1659 if( proc_number == 0 ) std::cout<<
"[CRB offline] M_error_type = "<<M_error_type<<std::endl;
1661 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1662 std::cout <<
" -- sampling init done in " << ti.elapsed() <<
"s\n";
1665 if ( M_error_type == CRB_RESIDUAL || M_error_type == CRB_RESIDUAL_SCM )
1667 int __QLhs = M_model->Qa();
1668 int __QRhs = M_model->Ql( 0 );
1669 int __QOutput = M_model->Ql( M_output_index );
1670 int __Qm = M_model->Qm();
1675 M_C0_pr.resize( __QRhs );
1676 for(
int __q1=0; __q1< __QRhs; __q1++)
1678 int __mMaxQ1=M_model->mMaxF(0,__q1);
1679 M_C0_pr[__q1].resize( __mMaxQ1 );
1680 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1682 M_C0_pr[__q1][__m1].resize( __QRhs );
1683 for(
int __q2=0; __q2< __QRhs; __q2++)
1685 int __mMaxQ2=M_model->mMaxF(0,__q2);
1686 M_C0_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1691 M_C0_du.resize( __QOutput );
1692 for(
int __q1=0; __q1< __QOutput; __q1++)
1694 int __mMaxQ1=M_model->mMaxF(M_output_index,__q1);
1695 M_C0_du[__q1].resize( __mMaxQ1 );
1696 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1698 M_C0_du[__q1][__m1].resize( __QOutput );
1699 for(
int __q2=0; __q2< __QOutput; __q2++)
1701 int __mMaxQ2=M_model->mMaxF(M_output_index,__q2);
1702 M_C0_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1710 M_Lambda_pr.resize( __QLhs );
1711 for(
int __q1=0; __q1< __QLhs; __q1++)
1713 int __mMaxQ1=M_model->mMaxA(__q1);
1714 M_Lambda_pr[__q1].resize( __mMaxQ1 );
1715 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1717 M_Lambda_pr[__q1][__m1].resize( __QRhs );
1718 for(
int __q2=0; __q2< __QRhs; __q2++)
1720 int __mMaxQ2=M_model->mMaxF(0,__q2);
1721 M_Lambda_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1726 M_Lambda_du.resize( __QLhs );
1727 for(
int __q1=0; __q1< __QLhs; __q1++)
1729 int __mMaxQ1=M_model->mMaxA(__q1);
1730 M_Lambda_du[__q1].resize( __mMaxQ1 );
1731 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1733 M_Lambda_du[__q1][__m1].resize( __QOutput );
1734 for(
int __q2=0; __q2< __QOutput; __q2++)
1736 int __mMaxQ2=M_model->mMaxF(M_output_index,__q2);
1737 M_Lambda_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1745 M_Gamma_pr.resize( __QLhs );
1746 for(
int __q1=0; __q1< __QLhs; __q1++)
1748 int __mMaxQ1=M_model->mMaxA(__q1);
1749 M_Gamma_pr[__q1].resize( __mMaxQ1 );
1750 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1752 M_Gamma_pr[__q1][__m1].resize( __QLhs );
1753 for(
int __q2=0; __q2< __QLhs; __q2++)
1755 int __mMaxQ2=M_model->mMaxA(__q2);
1756 M_Gamma_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1761 M_Gamma_du.resize( __QLhs );
1762 for(
int __q1=0; __q1< __QLhs; __q1++)
1764 int __mMaxQ1=M_model->mMaxA(__q1);
1765 M_Gamma_du[__q1].resize( __mMaxQ1 );
1766 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1768 M_Gamma_du[__q1][__m1].resize( __QLhs );
1769 for(
int __q2=0; __q2< __QLhs; __q2++)
1771 int __mMaxQ2=M_model->mMaxA(__q2);
1772 M_Gamma_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1777 if ( model_type::is_time_dependent )
1781 M_Cmf_pr.resize( __Qm );
1782 for(
int __q1=0; __q1< __Qm; __q1++)
1784 int __mMaxQ1=M_model->mMaxM(__q1);
1785 M_Cmf_pr[__q1].resize( __mMaxQ1 );
1786 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1788 M_Cmf_pr[__q1][__m1].resize( __QRhs );
1789 for(
int __q2=0; __q2< __QRhs; __q2++)
1791 int __mMaxQ2=M_model->mMaxF(0,__q2);
1792 M_Cmf_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1797 M_Cmf_du.resize( __Qm );
1798 for(
int __q1=0; __q1< __Qm; __q1++)
1800 int __mMaxQ1=M_model->mMaxM(__q1);
1801 M_Cmf_du[__q1].resize( __mMaxQ1 );
1802 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1804 M_Cmf_du[__q1][__m1].resize( __QOutput );
1805 for(
int __q2=0; __q2< __QOutput; __q2++)
1807 int __mMaxQ2=M_model->mMaxF(M_output_index,__q2);
1808 M_Cmf_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1812 M_Cmf_du_ini.resize( __Qm );
1813 for(
int __q1=0; __q1< __Qm; __q1++)
1815 int __mMaxQ1=M_model->mMaxM(__q1);
1816 M_Cmf_du_ini[__q1].resize( __mMaxQ1 );
1817 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1819 M_Cmf_du_ini[__q1][__m1].resize( __QOutput );
1820 for(
int __q2=0; __q2< __QOutput; __q2++)
1822 int __mMaxQ2=M_model->mMaxF(M_output_index,__q2);
1823 M_Cmf_du_ini[__q1][__m1][__q2].resize( __mMaxQ2 );
1830 M_Cma_pr.resize( __Qm );
1831 for(
int __q1=0; __q1< __Qm; __q1++)
1833 int __mMaxQ1=M_model->mMaxM(__q1);
1834 M_Cma_pr[__q1].resize( __mMaxQ1 );
1835 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1837 M_Cma_pr[__q1][__m1].resize( __QLhs );
1838 for(
int __q2=0; __q2< __QLhs; __q2++)
1840 int __mMaxQ2=M_model->mMaxA(__q2);
1841 M_Cma_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1846 M_Cma_du.resize( __Qm );
1847 for(
int __q1=0; __q1< __Qm; __q1++)
1849 int __mMaxQ1=M_model->mMaxM(__q1);
1850 M_Cma_du[__q1].resize( __mMaxQ1 );
1851 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1853 M_Cma_du[__q1][__m1].resize( __QLhs );
1854 for(
int __q2=0; __q2< __QLhs; __q2++)
1856 int __mMaxQ2=M_model->mMaxA(__q2);
1857 M_Cma_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1864 M_Cmm_pr.resize( __Qm );
1865 for(
int __q1=0; __q1< __Qm; __q1++)
1867 int __mMaxQ1=M_model->mMaxM(__q1);
1868 M_Cmm_pr[__q1].resize( __mMaxQ1 );
1869 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1871 M_Cmm_pr[__q1][__m1].resize( __Qm );
1872 for(
int __q2=0; __q2< __Qm; __q2++)
1874 int __mMaxQ2=M_model->mMaxM(__q2);
1875 M_Cmm_pr[__q1][__m1][__q2].resize( __mMaxQ2 );
1880 M_Cmm_du.resize( __Qm );
1881 for(
int __q1=0; __q1< __Qm; __q1++)
1883 int __mMaxQ1=M_model->mMaxM(__q1);
1884 M_Cmm_du[__q1].resize( __mMaxQ1 );
1885 for(
int __m1=0; __m1< __mMaxQ1; __m1++)
1887 M_Cmm_du[__q1][__m1].resize( __Qm );
1888 for(
int __q2=0; __q2< __Qm; __q2++)
1890 int __mMaxQ2=M_model->mMaxM(__q2);
1891 M_Cmm_du[__q1][__m1][__q2].resize( __mMaxQ2 );
1898 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1899 std::cout <<
" -- residual data init done in " << ti.elapsed() << std::endl;
1905 if( M_error_type == CRB_NO_RESIDUAL )
1906 mu = M_Dmu->element();
1910 boost::tie( mu, index ) = M_Xi->min();
1913 int size = mu.size();
1914 if( proc_number == 0 )
1916 std::cout <<
" -- start with mu = [ ";
1917 for (
int i=0; i<size-1; i++ ) std::cout<<mu( i )<<
" ";
1918 std::cout<<mu( size-1 )<<
" ]"<<std::endl;
1933 LOG(INFO) <<
"[CRB::offline] allocate reduced basis data structures\n";
1934 M_Aqm_pr.resize( M_model->Qa() );
1935 M_Aqm_du.resize( M_model->Qa() );
1936 M_Aqm_pr_du.resize( M_model->Qa() );
1939 M_Jqm_pr.resize( M_model->Qa() );
1941 for(
int q=0; q<M_model->Qa(); q++)
1943 M_Aqm_pr[q].resize( M_model->mMaxA(q) );
1944 M_Aqm_du[q].resize( M_model->mMaxA(q) );
1945 M_Aqm_pr_du[q].resize( M_model->mMaxA(q) );
1948 M_Jqm_pr[q].resize( M_model->mMaxA(q) );
1951 M_Mqm_pr.resize( M_model->Qm() );
1952 M_Mqm_du.resize( M_model->Qm() );
1953 M_Mqm_pr_du.resize( M_model->Qm() );
1954 for(
int q=0; q<M_model->Qm(); q++)
1956 M_Mqm_pr[q].resize( M_model->mMaxM(q) );
1957 M_Mqm_du[q].resize( M_model->mMaxM(q) );
1958 M_Mqm_pr_du[q].resize( M_model->mMaxM(q) );
1961 int QInitialGuessV = M_model->QInitialGuess();
1962 M_InitialGuessV_pr.resize( QInitialGuessV );
1963 for(
int q=0; q<QInitialGuessV; q++)
1965 M_InitialGuessV_pr[q].resize( M_model->mMaxInitialGuess(q) );
1968 M_Fqm_pr.resize( M_model->Ql( 0 ) );
1969 M_Fqm_du.resize( M_model->Ql( 0 ) );
1970 M_Lqm_pr.resize( M_model->Ql( M_output_index ) );
1971 M_Lqm_du.resize( M_model->Ql( M_output_index ) );
1974 M_Rqm_pr.resize( M_model->Ql( 0 ) );
1976 for(
int q=0; q<M_model->Ql( 0 ); q++)
1978 M_Fqm_pr[q].resize( M_model->mMaxF( 0 , q) );
1979 M_Fqm_du[q].resize( M_model->mMaxF( 0 , q) );
1982 M_Rqm_pr[q].resize( M_model->mMaxF( 0 , q) );
1984 for(
int q=0; q<M_model->Ql( M_output_index ); q++)
1986 M_Lqm_pr[q].resize( M_model->mMaxF( M_output_index , q) );
1987 M_Lqm_du[q].resize( M_model->mMaxF( M_output_index , q) );
1996 if( proc_number == 0 )
1998 std::cout<<
"we are going to enrich the reduced basis"<<std::endl;
1999 std::cout<<
"there are "<<M_N<<
" elements in the database"<<std::endl;
2001 LOG(INFO) <<
"we are going to enrich the reduced basis"<<std::endl;
2002 LOG(INFO) <<
"there are "<<M_N<<
" elements in the database"<<std::endl;
2010 LOG(INFO) <<
"[CRB::offline] compute affine decomposition\n";
2011 std::vector< std::vector<sparse_matrix_ptrtype> > Aqm;
2012 std::vector< std::vector<sparse_matrix_ptrtype> > Mqm;
2017 std::vector< std::vector<element_ptrtype> > InitialGuessV;
2018 std::vector< std::vector<std::vector<vector_ptrtype> > > Fqm,Lqm;
2019 sparse_matrix_ptrtype Aq_transpose = M_model->newMatrix();
2021 sparse_matrix_ptrtype A = M_model->newMatrix();
2022 int nl = M_model->Nl();
2023 std::vector< vector_ptrtype > F( nl );
2024 for(
int l=0; l<nl; l++)
2025 F[l]=M_model->newVector();
2027 std::vector< std::vector<sparse_matrix_ptrtype> > Jqm;
2028 std::vector< std::vector<std::vector<vector_ptrtype> > > Rqm;
2030 InitialGuessV = M_model->computeInitialGuessVAffineDecomposition();
2034 boost::tie( Mqm , Jqm, Rqm ) = M_model->computeAffineDecomposition();
2037 if( option(
"crb.stock-matrices").
template as<bool>() )
2038 boost::tie( Mqm, Aqm, Fqm ) = M_model->computeAffineDecomposition();
2041 element_ptrtype dual_initial_field(
new element_type( M_model->functionSpace() ) );
2042 element_ptrtype uproj(
new element_type( M_model->functionSpace() ) );
2043 auto u = M_model->functionSpace()->element();
2044 auto udu = M_model->functionSpace()->element();
2048 LOG(INFO) <<
"[CRB::offline] starting offline adaptive loop\n";
2050 bool reuse_prec = this->
vm()[
"crb.reuse-prec"].template as<bool>() ;
2052 bool use_predefined_WNmu = this->
vm()[
"crb.use-predefined-WNmu"].template as<bool>() ;
2054 int N_log_equi = this->
vm()[
"crb.use-logEquidistributed-WNmu"].template as<int>() ;
2055 int N_equi = this->
vm()[
"crb.use-equidistributed-WNmu"].template as<int>() ;
2057 if( N_log_equi > 0 || N_equi > 0 )
2058 use_predefined_WNmu =
true;
2060 if ( use_predefined_WNmu )
2062 std::string file_name = ( boost::format(
"SamplingWNmu") ).str();
2063 std::ifstream file ( file_name );
2066 throw std::logic_error(
"[CRB::offline] ERROR the file SamplingWNmu doesn't exist so it's impossible to known which parameters you want to use to build the database" );
2071 int sampling_size = M_WNmu->readFromFile(file_name);
2072 M_iter_max = sampling_size;
2074 mu = M_WNmu->at( M_N );
2076 if( proc_number == this->worldComm().masterRank() )
2077 std::cout<<
"[CRB::offline] read WNmu ( sampling size : "<<M_iter_max<<
" )"<<std::endl;
2081 LOG(INFO) <<
"[CRB::offline] strategy "<< M_error_type <<
"\n";
2082 if( proc_number == this->worldComm().masterRank() ) std::cout <<
"[CRB::offline] strategy "<< M_error_type <<std::endl;
2084 if( M_error_type == CRB_NO_RESIDUAL || use_predefined_WNmu )
2090 while ( M_maxerror > M_tolerance && M_N < M_iter_max )
2093 boost::timer timer, timer2;
2094 LOG(INFO) <<
"========================================"<<
"\n";
2096 if ( M_error_type == CRB_NO_RESIDUAL )
2097 LOG(INFO) <<
"N=" << M_N <<
"/" << M_iter_max <<
" ( nb proc : "<<worldComm().globalSize()<<
")";
2099 LOG(INFO) <<
"N=" << M_N <<
"/" << M_iter_max <<
" maxerror=" << M_maxerror <<
" / " << M_tolerance <<
"( nb proc : "<<worldComm().globalSize()<<
")";
2102 u.setName( ( boost::format(
"fem-primal-N%1%-proc%2%" ) % (M_N) % proc_number ).str() );
2103 udu.setName( ( boost::format(
"fem-dual-N%1%-proc%2%" ) % (M_N) % proc_number ).str() );
2105 if ( M_model->isSteady() && ! M_use_newton )
2110 bool zero_iteration=
false;
2111 u = offlineFixedPointPrimal( mu , A , zero_iteration );
2112 if( solve_dual_problem )
2113 udu = offlineFixedPointDual( mu , dual_initial_field , A , u, zero_iteration );
2116 if ( M_model->isSteady() && M_use_newton )
2122 LOG(INFO) <<
"[CRB::offline] solving primal" <<
"\n";
2123 u = M_model->solve( mu );
2124 if( proc_number == this->worldComm().masterRank() )
2125 LOG(INFO) <<
" -- primal problem solved in " << timer2.elapsed() <<
"s";
2130 if( ! M_model->isSteady() )
2132 bool zero_iteration=
true;
2133 u = offlineFixedPointPrimal( mu, A , zero_iteration );
2134 if ( solve_dual_problem || M_error_type==CRB_RESIDUAL || M_error_type == CRB_RESIDUAL_SCM )
2135 udu = offlineFixedPointDual( mu , dual_initial_field , A , u, zero_iteration );
2139 if( ! use_predefined_WNmu )
2140 M_WNmu->push_back( mu, index );
2142 M_WNmu_complement = M_WNmu->complement();
2144 bool norm_zero =
false;
2146 if ( M_model->isSteady() )
2148 M_WN.push_back( u );
2149 M_WNdu.push_back( udu );
2154 if ( M_N == 0 && orthonormalize_primal==
false )
2164 element_ptrtype primal_initial_field (
new element_type ( M_model->functionSpace() ) );
2166 M_model->initializationField( primal_initial_field, mu );
2167 int norm_ini = M_model->scalarProduct( *primal_initial_field, *primal_initial_field );
2169 if ( norm_ini != 0 )
2171 M_WN.push_back( *primal_initial_field );
2172 M_WNdu.push_back( *dual_initial_field );
2182 LOG(INFO)<<
"[CRB::offline] start of POD \n";
2187 if ( M_mode_number == 1 )
2197 POD->setNm( M_mode_number*M_Nm );
2198 int size = mu.size();
2199 LOG(INFO)<<
"... CRB M_mode_number = "<<M_mode_number<<
"\n";
2200 LOG(INFO)<<
"for mu = [ ";
2202 for (
int i=0; i<size-1; i++ ) LOG(INFO)<<mu[i]<<
" , ";
2204 LOG(INFO)<<mu[ size-1 ];
2207 double Tf = M_model->timeFinal();
2208 double dt = M_model->timeStep();
2209 int nb_mode_max = Tf/dt;
2211 if ( M_mode_number>=nb_mode_max-1 )
2213 std::cout<<
"Error : we access to "<<M_mode_number<<
"^th mode"<<std::endl;
2214 std::cout<<
"parameter choosen : [ ";
2216 for (
int i=0; i<size-1; i++ ) std::cout<<mu[i]<<
" , ";
2218 std::cout<<mu[ size-1 ]<<
" ] "<<std::endl;
2219 throw std::logic_error(
"[CRB::offline] ERROR during the construction of the reduced basis, one parameter has been choosen too many times" );
2223 POD->setBdf( M_bdf_primal_save );
2224 POD->setModel( M_model );
2225 POD->setTimeInitial( M_model->timeInitial() );
2226 mode_set_type ModeSet;
2228 size_type number_max_of_mode = POD->pod( ModeSet,
true );
2230 if ( number_max_of_mode < M_Nm )
2232 std::cout<<
"With crb.Nm = "<<M_Nm<<
" there was too much too small eigenvalues so the value";
2233 std::cout<<
" of crb.Nm as been changed to "<<number_max_of_mode<<std::endl;
2234 M_Nm = number_max_of_mode;
2239 if ( !seek_mu_in_complement )
2242 M_WN.push_back( ModeSet[M_mode_number*M_Nm-1+i] ) ;
2249 M_WN.push_back( ModeSet[i] ) ;
2253 if ( solve_dual_problem || M_error_type==CRB_RESIDUAL || M_error_type == CRB_RESIDUAL_SCM )
2255 POD->setBdf( M_bdf_dual );
2256 POD->setTimeInitial( M_model->timeFinal()+M_model->timeStep() );
2257 mode_set_type ModeSetdu;
2258 POD->pod( ModeSetdu,
false );
2260 if ( !seek_mu_in_complement )
2263 M_WNdu.push_back( ModeSetdu[M_mode_number*M_Nm-1+i] ) ;
2269 M_WNdu.push_back( ModeSetdu[i] ) ;
2275 element_ptrtype element_zero (
new element_type ( M_model->functionSpace() ) );
2279 M_WNdu.push_back( *element_zero ) ;
2289 size_type number_of_added_elements = M_Nm + ( M_N==0 && orthonormalize_primal==
false && norm_zero==
false && !M_model->isSteady() );
2292 if( M_model->isSteady() )
2293 number_of_added_elements=1;
2295 M_N+=number_of_added_elements;
2297 double norm_max = option(_name=
"crb.orthonormality-tol").template as<double>();
2298 int max_iter = option(_name=
"crb.orthonormality-max-iter").template as<int>();
2299 if ( orthonormalize_primal )
2301 double norm = norm_max+1;
2304 while( norm >= norm_max && iter < max_iter)
2309 if( math::abs(old-norm) < norm_max )
2315 if ( orthonormalize_dual )
2317 double norm = norm_max+1;
2320 while( norm >= norm_max && iter < max_iter )
2324 if( math::abs(old-norm) < norm_max )
2330 if( ! M_use_newton )
2332 LOG(INFO) <<
"[CRB::offline] compute Aq_pr, Aq_du, Aq_pr_du" <<
"\n";
2334 for (
size_type q = 0; q < M_model->Qa(); ++q )
2336 for(
size_type m = 0; m < M_model->mMaxA(q); ++m )
2338 M_Aqm_pr[q][m].conservativeResize( M_N, M_N );
2339 M_Aqm_du[q][m].conservativeResize( M_N, M_N );
2340 M_Aqm_pr_du[q][m].conservativeResize( M_N, M_N );
2343 for (
size_type i = M_N-number_of_added_elements; i < M_N; i++ )
2347 M_Aqm_pr[q][m]( i, j ) = M_model->Aqm(q , m , M_WN[i], M_WN[j] );
2348 M_Aqm_du[q][m]( i, j ) = M_model->Aqm( q , m , M_WNdu[i], M_WNdu[j],
true );
2349 M_Aqm_pr_du[q][m]( i, j ) = M_model->Aqm(q , m , M_WNdu[i], M_WN[j] );
2353 for (
size_type j=M_N-number_of_added_elements; j < M_N; j++ )
2357 M_Aqm_pr[q][m]( i, j ) = M_model->Aqm(q , m , M_WN[i], M_WN[j] );
2358 M_Aqm_du[q][m]( i, j ) = M_model->Aqm(q , m , M_WNdu[i], M_WNdu[j],
true );
2359 M_Aqm_pr_du[q][m]( i, j ) = M_model->Aqm(q , m , M_WNdu[i], M_WN[j] );
2365 LOG(INFO) <<
"[CRB::offline] compute Mq_pr, Mq_du, Mq_pr_du" <<
"\n";
2368 LOG(INFO) <<
"[CRB::offline] compute Fq_pr, Fq_du" <<
"\n";
2370 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
2372 for(
size_type m = 0; m < M_model->mMaxF( 0, q ); ++m )
2374 M_Fqm_pr[q][m].conservativeResize( M_N );
2375 M_Fqm_du[q][m].conservativeResize( M_N );
2376 for (
size_type l = 1; l <= number_of_added_elements; ++l )
2379 M_Fqm_pr[q][m]( index ) = M_model->Fqm( 0, q, m, M_WN[index] );
2380 M_Fqm_du[q][m]( index ) = M_model->Fqm( 0, q, m, M_WNdu[index] );
2390 LOG(INFO) <<
"[CRB::offline] compute Jq_pr " <<
"\n";
2392 for (
size_type q = 0; q < M_model->Qa(); ++q )
2394 for(
size_type m = 0; m < M_model->mMaxA(q); ++m )
2396 M_Jqm_pr[q][m].conservativeResize( M_N, M_N );
2399 for (
size_type i = M_N-number_of_added_elements; i < M_N; i++ )
2403 M_Jqm_pr[q][m]( i, j ) = Jqm[q][m]->energy( M_WN[i], M_WN[j] );
2407 for (
size_type j=M_N-number_of_added_elements; j < M_N; j++ )
2411 M_Jqm_pr[q][m]( i, j ) = Jqm[q][m]->energy( M_WN[i], M_WN[j] );
2418 LOG(INFO) <<
"[CRB::offline] compute Rq_pr" <<
"\n";
2420 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
2422 for(
size_type m = 0; m < M_model->mMaxF( 0, q ); ++m )
2424 M_Rqm_pr[q][m].conservativeResize( M_N );
2425 for (
size_type l = 1; l <= number_of_added_elements; ++l )
2428 M_Rqm_pr[q][m]( index ) = M_model->Fqm( 0, q, m, M_WN[index] );
2436 LOG(INFO) <<
"[CRB::offline] compute MFqm" <<
"\n";
2437 int q_max = M_model->QInitialGuess();
2440 int m_max =M_model->mMaxInitialGuess(q);
2443 M_InitialGuessV_pr[q][m].conservativeResize( M_N );
2445 M_InitialGuessV_pr[q][m]( j ) =
inner_product( *InitialGuessV[q][m] , M_WN[j] );
2450 for (
size_type q = 0; q < M_model->Qm(); ++q )
2452 for(
size_type m = 0; m < M_model->mMaxM(q); ++m )
2454 M_Mqm_pr[q][m].conservativeResize( M_N, M_N );
2455 M_Mqm_du[q][m].conservativeResize( M_N, M_N );
2456 M_Mqm_pr_du[q][m].conservativeResize( M_N, M_N );
2459 for (
size_type i=M_N-number_of_added_elements ; i < M_N; i++ )
2463 M_Mqm_pr[q][m]( i, j ) = M_model->Mqm(q , m , M_WN[i], M_WN[j] );
2464 M_Mqm_du[q][m]( i, j ) = M_model->Mqm(q , m , M_WNdu[i], M_WNdu[j],
true );
2465 M_Mqm_pr_du[q][m]( i, j ) = M_model->Mqm( q , m , M_WNdu[i], M_WN[j] );
2468 for (
size_type j = M_N-number_of_added_elements; j < M_N ; j++ )
2472 M_Mqm_pr[q][m]( i, j ) = M_model->Mqm(q , m , M_WN[i], M_WN[j] );
2473 M_Mqm_du[q][m]( i, j ) = M_model->Mqm(q , m , M_WNdu[i], M_WNdu[j],
true );
2474 M_Mqm_pr_du[q][m]( i, j ) = M_model->Mqm(q , m , M_WNdu[i], M_WN[j] );
2480 LOG(INFO) <<
"[CRB::offline] compute Lq_pr, Lq_du" <<
"\n";
2482 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
2484 for(
size_type m = 0; m < M_model->mMaxF( M_output_index, q ); ++m )
2486 M_Lqm_pr[q][m].conservativeResize( M_N );
2487 M_Lqm_du[q][m].conservativeResize( M_N );
2489 for (
size_type l = 1; l <= number_of_added_elements; ++l )
2492 M_Lqm_pr[q][m]( index ) = M_model->Fqm( M_output_index, q, m, M_WN[index] );
2493 M_Lqm_du[q][m]( index ) = M_model->Fqm( M_output_index, q, m, M_WNdu[index] );
2499 LOG(INFO) <<
"compute coefficients needed for the initialization of unknown in the online step\n";
2501 element_ptrtype primal_initial_field (
new element_type ( M_model->functionSpace() ) );
2502 element_ptrtype projection (
new element_type ( M_model->functionSpace() ) );
2503 M_model->initializationField( primal_initial_field, mu );
2504 if ( model_type::is_time_dependent || !M_model->isSteady() )
2506 if ( orthonormalize_primal )
2508 for (
size_type elem=M_N-number_of_added_elements; elem<M_N; elem++ )
2511 double k = M_model->scalarProduct( *primal_initial_field, M_WN[elem] );
2512 M_coeff_pr_ini_online.push_back( k );
2516 else if ( !orthonormalize_primal )
2518 matrixN_type MN ( (
int )M_N, (
int )M_N ) ;
2519 vectorN_type FN ( (
int )M_N );
2526 MN( i,j ) = M_model->scalarProduct( M_WN[j] , M_WN[i] );
2527 MN( j,i ) = MN( i,j );
2530 MN( i,i ) = M_model->scalarProduct( M_WN[i] , M_WN[i] );
2531 FN( i ) = M_model->scalarProduct( *primal_initial_field,M_WN[i] );
2534 vectorN_type projectionN ( (
int ) M_N );
2535 projectionN = MN.lu().solve( FN );
2537 for (
size_type i=M_N-number_of_added_elements; i<M_N; i++ )
2539 M_coeff_pr_ini_online.push_back( projectionN( i ) );
2543 if ( solve_dual_problem )
2545 if ( orthonormalize_dual )
2547 for (
size_type elem=M_N-number_of_added_elements; elem<M_N; elem++ )
2549 double k = M_model->scalarProduct( *dual_initial_field, M_WNdu[elem] );
2550 M_coeff_du_ini_online.push_back( k );
2554 else if ( !orthonormalize_dual )
2556 matrixN_type MNdu ( (
int )M_N, (
int )M_N ) ;
2557 vectorN_type FNdu ( (
int )M_N );
2564 MNdu( i,j ) = M_model->scalarProduct( M_WNdu[j] , M_WNdu[i] );
2565 MNdu( j,i ) = MNdu( i,j );
2568 MNdu( i,i ) = M_model->scalarProduct( M_WNdu[i] , M_WNdu[i] );
2569 FNdu( i ) = M_model->scalarProduct( *dual_initial_field,M_WNdu[i] );
2572 vectorN_type projectionN ( (
int ) M_N );
2573 projectionN = MNdu.lu().solve( FNdu );
2575 for (
size_type i=M_N-number_of_added_elements; i<M_N; i++ )
2577 M_coeff_du_ini_online.push_back( projectionN( i ) );
2585 M_compute_variance = this->
vm()[
"crb.compute-variance"].template as<bool>();
2586 if ( M_database_contains_variance_info )
2587 throw std::logic_error(
"[CRB::offline] ERROR : build variance is not actived" );
2590 if ( M_error_type==CRB_RESIDUAL || M_error_type == CRB_RESIDUAL_SCM )
2592 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
2593 std::cout <<
" -- offlineResidual update starts\n";
2595 LOG(INFO)<<
"[CRB::offline] end of call offlineResidual and M_N = "<< M_N <<
"\n";
2596 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
2597 std::cout <<
" -- offlineResidual updated in " << timer2.elapsed() <<
"s\n";
2602 if ( M_error_type == CRB_NO_RESIDUAL && ! use_predefined_WNmu )
2610 already_exist=
false;
2612 mu = M_Dmu->element();
2614 BOOST_FOREACH(
auto _mu, *M_WNmu )
2620 while( already_exist );
2624 else if ( use_predefined_WNmu )
2627 if( M_N < M_iter_max )
2629 mu = M_WNmu->at( M_N );
2635 boost::tie( M_maxerror, mu, index , delta_pr , delta_du ) =
maxErrorBounds( M_N );
2637 M_index.push_back( index );
2640 int count = std::count( M_index.begin(),M_index.end(),index );
2641 M_mode_number = count;
2643 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
2644 std::cout <<
" -- max error bounds computed in " << timer2.elapsed() <<
"s\n";
2649 M_rbconv.insert( convergence( M_N, boost::make_tuple(M_maxerror,delta_pr,delta_du) ) );
2653 check( M_WNmu->size() );
2655 if ( this->
vm()[
"crb.check.rb"].template as<int>() == 1 )std::cout <<
" -- check reduced basis done in " << timer2.elapsed() <<
"s\n";
2658 LOG(INFO) <<
"time: " << timer.elapsed() <<
"\n";
2659 if( proc_number == 0 ) std::cout <<
"============================================================\n";
2660 LOG(INFO) <<
"========================================"<<
"\n";
2664 M_elements_database.setWn( boost::make_tuple( M_WN , M_WNdu ) );
2665 M_elements_database.saveDB();
2669 if( proc_number == 0 )
2670 std::cout<<
"number of elements in the reduced basis : "<<M_N<<
" ( nb proc : "<<worldComm().globalSize()<<
")"<<std::endl;
2671 LOG(INFO) <<
" index choosen : ";
2672 BOOST_FOREACH(
auto id, M_index )
2675 bool visualize_basis = this->
vm()[
"crb.visualize-basis"].template as<bool>() ;
2677 if ( visualize_basis )
2679 std::vector<wn_type> wn;
2680 std::vector<std::string> names;
2681 wn.push_back( M_WN );
2682 names.push_back(
"primal" );
2683 wn.push_back( M_WNdu );
2684 names.push_back(
"dual" );
2687 if ( orthonormalize_primal || orthonormalize_dual )
2689 std::cout<<
"[CRB::offline] Basis functions have been exported but warning elements have been orthonormalized"<<std::endl;
2702 template<
typename TruthModelType>
2706 static const bool is_composite = functionspace_type::is_composite;
2707 checkInitialGuess( expansion_uN, mu, error , mpl::bool_< is_composite >() );
2710 template<
typename TruthModelType>
2712 CRB<TruthModelType>::checkInitialGuess(
const element_type expansion_uN , parameter_type
const& mu, vectorN_type & error, mpl::bool_<false>)
const
2715 const element_ptrtype initial_guess = M_model->assembleInitialGuess( mu );
2716 auto Xh = expansion_uN.functionSpace();
2717 auto mesh = Xh->mesh();
2718 error(0) = math::sqrt(
2720 _expr=vf::idv( initial_guess ) * vf::idv( expansion_uN )
2721 * vf::idv( initial_guess ) * vf::idv( expansion_uN )
2728 template<
typename TruthModelType>
2730 CRB<TruthModelType>::checkInitialGuess(
const element_type expansion_uN , parameter_type
const& mu, vectorN_type & error, mpl::bool_<true>)
const
2733 index_vector_type index_vector;
2734 const element_ptrtype initial_guess = M_model->assembleInitialGuess( mu );
2735 ComputeIntegralsSquare compute_integrals_square( *initial_guess , expansion_uN );
2736 fusion::for_each( index_vector , compute_integrals_square );
2737 error.resize( functionspace_type::nSpaces );
2738 error = compute_integrals_square.vectorErrors();
2743 template<
typename TruthModelType>
2745 CRB<TruthModelType>::buildVarianceMatrixPhi(
int const N )
2747 static const bool is_composite = functionspace_type::is_composite;
2748 return buildVarianceMatrixPhi( N , mpl::bool_< is_composite >() );
2750 template<
typename TruthModelType>
2752 CRB<TruthModelType>::buildVarianceMatrixPhi(
int const N , mpl::bool_<true> )
2755 std::vector<std::string> s;
2756 static const int nb_spaces = functionspace_type::nSpaces;
2759 M_variance_matrix_phi.resize( nb_spaces );
2762 for(
int i=0;i<nb_spaces;i++)
2763 M_variance_matrix_phi[i].conservativeResize( M_N , M_N );
2766 int size_WN = M_WN.size();
2769 for(
int i = 0; i < size_WN; ++i)
2771 auto sub = subelements( M_WN[i] , s );
2772 element_type global_element = M_model->functionSpace()->element();
2773 wn_type v( nb_spaces , global_element);
2774 ComputePhi compute_phi(v);
2775 fusion::for_each( sub , compute_phi ) ;
2776 auto vect = compute_phi.vectorPhi();
2780 global_element.zero();
2781 BOOST_FOREACH(
auto element , vect )
2782 global_element += element;
2783 phi.push_back( global_element );
2787 index_vector_type index_vector;
2788 for(
int space=0; space<nb_spaces; space++)
2789 M_variance_matrix_phi[space].conservativeResize( N , N );
2792 for(
int i = 0; i < M_N; ++i)
2795 for(
int j = i+1; j < M_N; ++j)
2797 ComputeIntegrals compute_integrals( phi[i] , phi[j] );
2798 fusion::for_each( index_vector , compute_integrals );
2799 auto vect = compute_integrals.vectorIntegrals();
2800 for(
int space=0; space<nb_spaces; space++)
2801 M_variance_matrix_phi[space](i,j)=vect[space];
2803 ComputeIntegrals compute_integrals ( phi[i] , phi[i] );
2804 fusion::for_each( index_vector , compute_integrals );
2805 auto vect = compute_integrals.vectorIntegrals();
2806 for(
int space=0; space<nb_spaces; space++)
2807 M_variance_matrix_phi[space](i,i)=vect[space];
2813 template<
typename TruthModelType>
2815 CRB<TruthModelType>::buildVarianceMatrixPhi(
int const N , mpl::bool_<false> )
2817 std::vector<std::string> s;
2818 int nb_spaces = functionspace_type::nSpaces;
2821 M_variance_matrix_phi.resize( nb_spaces );
2823 M_variance_matrix_phi[0].conservativeResize( M_N , M_N );
2826 int size_WN = M_WN.size();
2830 mesh_ptrtype mesh = M_WN[0].functionSpace()->mesh();
2832 for(
int i = 0; i < size_WN; ++i)
2834 double surface =
integrate( _range=
elements(mesh), _expr=vf::cst(1.) ).evaluate()(0,0);
2835 double mean =
integrate( _range=
elements(mesh), _expr=vf::idv( M_WN[i] ) ).evaluate()(0,0);
2838 phi.push_back( M_WN[i] - element_mean );
2841 for(
int space=0; space<nb_spaces; space++)
2842 M_variance_matrix_phi[space].conservativeResize( N , N );
2844 for(
int i = 0; i < M_N; ++i)
2846 for(
int j = i+1; j < M_N; ++j)
2848 M_variance_matrix_phi[0]( i , j ) =
integrate( _range=
elements(mesh) , _expr=vf::idv( phi[i] ) * vf::idv( phi[j] ) ).evaluate()(0,0);
2849 M_variance_matrix_phi[0]( j , i ) = M_variance_matrix_phi[0]( i , j );
2852 M_variance_matrix_phi[0]( i , i ) =
integrate( _range=
elements(mesh) , _expr=vf::idv( phi[i] ) * vf::idv( phi[i] ) ).evaluate()(0,0);
2859 template<
typename TruthModelType>
2861 CRB<TruthModelType>::buildFunctionFromRbCoefficients(
int N, std::vector< vectorN_type >
const & RBcoeff, wn_type
const & WN, std::vector<element_ptrtype> & FEMsolutions )
2864 if( WN.size() == 0 )
2865 throw std::logic_error(
"[CRB::buildFunctionFromRbCoefficients] ERROR : reduced basis space is empty" );
2868 int nb_solutions = RBcoeff.size();
2870 for(
int i = 0; i < nb_solutions; i++ )
2872 element_ptrtype FEMelement (
new element_type( M_model->functionSpace() ) );
2873 FEMelement->setZero();
2874 for(
int j = 0; j < N; j++ )
2875 FEMelement->add( RBcoeff[i](j) , WN[j] );
2876 FEMsolutions.push_back( FEMelement );
2880 template<
typename TruthModelType>
2882 CRB<TruthModelType>::compareResidualsForTransientProblems(
int N, parameter_type
const& mu, std::vector<element_type>
const & Un, std::vector<element_type>
const & Unold,
2883 std::vector<element_type>
const& Undu, std::vector<element_type>
const & Unduold,
2884 std::vector< std::vector<double> >
const& primal_residual_coeffs, std::vector < std::vector<double> >
const& dual_residual_coeffs )
const
2887 LOG( INFO ) <<
"\n compareResidualsForTransientProblems \n";
2890 if ( M_model->isSteady() )
2892 throw std::logic_error(
"[CRB::compareResidualsForTransientProblems] ERROR : to check residual in a steady case, use checkResidual and not compareResidualsForTransientProblems" );
2895 sparse_matrix_ptrtype A,AM,M,Adu;
2897 std::vector<vector_ptrtype> F,L;
2899 vector_ptrtype Rhs( backend->newVector( M_model->functionSpace() ) );
2900 vector_ptrtype Aun( backend->newVector( M_model->functionSpace() ) );
2901 vector_ptrtype AduUn( backend->newVector( M_model->functionSpace() ) );
2903 vector_ptrtype Mun( backend->newVector( M_model->functionSpace() ) );
2904 vector_ptrtype Munold( backend->newVector( M_model->functionSpace() ) );
2905 vector_ptrtype Frhs( backend->newVector( M_model->functionSpace() ) );
2906 vector_ptrtype un( backend->newVector( M_model->functionSpace() ) );
2907 vector_ptrtype unold( backend->newVector( M_model->functionSpace() ) );
2908 vector_ptrtype undu( backend->newVector( M_model->functionSpace() ) );
2909 vector_ptrtype unduold( backend->newVector( M_model->functionSpace() ) );
2912 auto bdf_primal = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_primal_check_residual_transient" );
2913 bdf_primal->setTimeInitial( M_model->timeInitial() );
2914 bdf_primal->setTimeStep( M_model->timeStep() );
2915 bdf_primal->setTimeFinal( M_model->timeFinal() );
2916 bdf_primal->setOrder( M_model->timeOrder() );
2919 double time_step = bdf_primal->timeStep();
2921 for ( bdf_primal->start() ; !bdf_primal->isFinished(); bdf_primal->next() )
2924 auto bdf_poly = bdf_primal->polyDeriv();
2925 boost::tie( M, A, F) = M_model->update( mu , bdf_primal->time() );
2929 *un = Un[time_index];
2930 *unold = Unold[time_index];
2932 A->multVector( un, Aun );
2933 M->multVector( un, Mun );
2935 M->multVector( unold, Munold );
2940 vector_ptrtype __ef_pr( backend->newVector( M_model->functionSpace() ) );
2941 vector_ptrtype __ea_pr( backend->newVector( M_model->functionSpace() ) );
2942 vector_ptrtype __emu_pr( backend->newVector( M_model->functionSpace() ) );
2943 vector_ptrtype __emuold_pr( backend->newVector( M_model->functionSpace() ) );
2944 M_model->l2solve( __ef_pr, Frhs );
2945 M_model->l2solve( __ea_pr, Aun );
2946 M_model->l2solve( __emu_pr, Mun );
2947 M_model->l2solve( __emuold_pr, Munold );
2949 double check_Cff_pr = M_model->scalarProduct( __ef_pr,__ef_pr );
2950 double check_Caf_pr = 2*M_model->scalarProduct( __ef_pr,__ea_pr );
2951 double check_Caa_pr = M_model->scalarProduct( __ea_pr,__ea_pr );
2952 double check_Cmf_pr = 2./time_step*( M_model->scalarProduct( __emu_pr , __ef_pr )+M_model->scalarProduct( __emuold_pr , __ef_pr ) );
2953 double check_Cma_pr = 2./time_step*( M_model->scalarProduct( __emu_pr , __ea_pr )+M_model->scalarProduct( __emuold_pr , __ea_pr ) );
2954 double check_Cmm_pr = 1./(time_step*time_step)*( M_model->scalarProduct( __emu_pr , __emu_pr ) + 2*M_model->scalarProduct( __emu_pr , __emuold_pr ) + M_model->scalarProduct( __emuold_pr , __emuold_pr ) );
2956 double Cff_pr = primal_residual_coeffs[time_index][0];
2957 double Caf_pr = primal_residual_coeffs[time_index][1];
2958 double Caa_pr = primal_residual_coeffs[time_index][2];
2959 double Cmf_pr = primal_residual_coeffs[time_index][3];
2960 double Cma_pr = primal_residual_coeffs[time_index][4];
2961 double Cmm_pr = primal_residual_coeffs[time_index][5];
2963 LOG(INFO)<<
" --- time : "<<bdf_primal->time();
2964 LOG(INFO)<<
"Cff : "<< check_Cff_pr <<
" - "<<Cff_pr<<
" => "<<check_Cff_pr-Cff_pr;
2965 LOG(INFO)<<
"Caf : "<< check_Caf_pr <<
" - "<<Caf_pr<<
" => "<<check_Caf_pr-Caf_pr;
2966 LOG(INFO)<<
"Caa : "<< check_Caa_pr <<
" - "<<Caa_pr<<
" => "<<check_Caa_pr-Caa_pr;
2967 LOG(INFO)<<
"Cmf : "<< check_Cmf_pr <<
" - "<<Cmf_pr<<
" => "<<check_Cmf_pr-Cmf_pr;
2968 LOG(INFO)<<
"Cma : "<< check_Cma_pr <<
" - "<<Cma_pr<<
" => "<<check_Cma_pr-Cma_pr;
2969 LOG(INFO)<<
"Cmm : "<< check_Cmm_pr <<
" - "<<Cmm_pr<<
" => "<<check_Cmm_pr-Cmm_pr;
2975 bool solve_dual_problem = this->
vm()[
"crb.solve-dual-problem"].template as<bool>();
2980 if( solve_dual_problem )
2983 Adu = M_model->newMatrix();
2985 auto bdf_dual = bdf( _space=M_model->functionSpace(), _vm=this->
vm() , _name=
"bdf_dual_check_residual_transient" );
2987 bdf_dual->setTimeInitial( M_model->timeFinal()+M_model->timeStep() );
2988 bdf_dual->setTimeStep( -M_model->timeStep() );
2989 bdf_dual->setTimeFinal( M_model->timeInitial()+M_model->timeStep() );
2990 bdf_dual->setOrder( M_model->timeOrder() );
2994 LOG(INFO)<<
"**********dual problem************* "<<std::endl;
2996 boost::tie( M, A, F) = M_model->update( mu , bdf_dual->timeInitial() );
2998 vectorN_type dual_initial ( N );
2999 for(
int i=0; i<N; i++)
3000 dual_initial(i) = M_coeff_du_ini_online[i];
3001 auto dual_initial_field = this->
expansion( dual_initial , N , M_WNdu );
3003 *undu = dual_initial_field;
3004 M->multVector( undu, Mun );
3008 auto R = backend->newVector( M_model->functionSpace() );
3010 R->add( -1 , *Mun );
3013 vector_ptrtype __ef_du( backend->newVector( M_model->functionSpace() ) );
3014 vector_ptrtype __emu_du( backend->newVector( M_model->functionSpace() ) );
3015 M_model->l2solve( __ef_du, Frhs );
3016 M_model->l2solve( __emu_du, Mun );
3017 double check_Cff_du = M_model->scalarProduct( __ef_du,__ef_du );
3018 double check_Cmf_du = 2*M_model->scalarProduct( __ef_du,__emu_du );
3019 double check_Cmm_du = M_model->scalarProduct( __emu_du,__emu_du );
3020 double residual_final_condition = math::abs( check_Cff_du + check_Cmf_du + check_Cmm_du );
3023 time_step = bdf_dual->timeStep();
3026 for ( bdf_dual->start(); !bdf_dual->isFinished() ; bdf_dual->next() )
3028 auto bdf_poly = bdf_dual->polyDeriv();
3030 boost::tie( M, A, F ) = M_model->update( mu , bdf_dual->time() );
3031 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
3034 A->transpose( Adu );
3035 *undu = Undu[time_index];
3036 *unduold = Unduold[time_index];
3037 Adu->multVector( undu, AduUn );
3038 M->multVector( undu, Mun );
3039 M->multVector( unduold, Munold );
3041 Munold->scale( -1 );
3044 vector_ptrtype __ea_du( backend->newVector( M_model->functionSpace() ) );
3045 vector_ptrtype __emu_du( backend->newVector( M_model->functionSpace() ) );
3046 vector_ptrtype __emuold_du( backend->newVector( M_model->functionSpace() ) );
3047 M_model->l2solve( __ea_du, AduUn );
3048 M_model->l2solve( __emu_du, Mun );
3049 M_model->l2solve( __emuold_du, Munold );
3050 double check_Caa_du = M_model->scalarProduct( __ea_du,__ea_du );
3051 double check_Cma_du = 2./time_step*( M_model->scalarProduct( __emu_du , __ea_du )+M_model->scalarProduct( __emuold_du , __ea_du ) );
3052 double check_Cmm_du = 1./(time_step*time_step)*( M_model->scalarProduct( __emu_du , __emu_du ) + 2*M_model->scalarProduct( __emu_du , __emuold_du ) + M_model->scalarProduct( __emuold_du , __emuold_du ) );
3055 double Cff_du = dual_residual_coeffs[time_index][0];
3056 double Caf_du = dual_residual_coeffs[time_index][1];
3057 double Caa_du = dual_residual_coeffs[time_index][2];
3058 double Cmf_du = dual_residual_coeffs[time_index][3];
3059 double Cma_du = dual_residual_coeffs[time_index][4];
3060 double Cmm_du = dual_residual_coeffs[time_index][5];
3061 LOG(INFO)<<
" --- time : "<<bdf_dual->time()<<std::endl;
3062 LOG(INFO)<<
"Caa : "<< check_Caa_du <<
" - "<<Caa_du<<
" => "<<check_Caa_du-Caa_du<<std::endl;
3063 LOG(INFO)<<
"Cma : "<< check_Cma_du <<
" - "<<Cma_du<<
" => "<<check_Cma_du-Cma_du<<std::endl;
3064 LOG(INFO)<<
"Cmm : "<< check_Cmm_du <<
" - "<<Cmm_du<<
" => "<<check_Cmm_du-Cmm_du<<std::endl;
3070 sum += math::abs( check_Caa_du + check_Cma_du + check_Cmm_du );
3077 template<
typename TruthModelType>
3079 CRB<TruthModelType>::checkResidual( parameter_type
const& mu, std::vector< std::vector<double> >
const& primal_residual_coeffs,
3084 if ( orthonormalize_primal || orthonormalize_dual )
3086 throw std::logic_error(
"[CRB::checkResidual] ERROR : to check residual don't use orthonormalization" );
3089 if ( !M_model->isSteady() )
3091 throw std::logic_error(
"[CRB::checkResidual] ERROR : to check residual select steady state" );
3094 if ( M_error_type==CRB_NO_RESIDUAL || M_error_type==CRB_EMPIRICAL )
3096 throw std::logic_error(
"[CRB::checkResidual] ERROR : to check residual set option crb.error-type to 0 or 1" );
3099 int size = mu.size();
3103 std::cout<<
"[CRB::checkResidual] use mu = [";
3104 for (
int i=0; i<size-1; i++ ) std::cout<< mu[i] <<
" , ";
3105 std::cout<< mu[size-1]<<
" ]"<<std::endl;
3108 LOG( INFO ) <<
"[CRB::checkResidual] use mu = \n"<<mu;
3109 sparse_matrix_ptrtype A,At,M;
3110 std::vector<vector_ptrtype> F,L;
3117 vector_ptrtype U( M_backend->newVector( M_model->functionSpace() ) );
3118 vector_ptrtype Rhs( M_backend->newVector( M_model->functionSpace() ) );
3120 boost::timer timer, timer2;
3122 boost::tie( boost::tuples::ignore, A, F ) = M_model->update( mu );
3124 LOG(INFO) <<
" -- updated model for parameter in " << timer2.elapsed() <<
"s\n";
3127 LOG(INFO) <<
"[CRB::checkResidual] transpose primal matrix" <<
"\n";
3128 At = M_model->newMatrix();
3129 if( option(
"crb.use-symmetric-matrix").
template as<bool>() )
3140 *Rhs = *F[M_output_index];
3146 vector_ptrtype Aun( M_backend->newVector( M_model->functionSpace() ) );
3147 vector_ptrtype Atun( M_backend->newVector( M_model->functionSpace() ) );
3148 vector_ptrtype Un( M_backend->newVector( M_model->functionSpace() ) );
3149 vector_ptrtype Undu( M_backend->newVector( M_model->functionSpace() ) );
3150 vector_ptrtype Frhs( M_backend->newVector( M_model->functionSpace() ) );
3151 vector_ptrtype Lrhs( M_backend->newVector( M_model->functionSpace() ) );
3154 A->multVector( Un, Aun );
3155 At->multVector( Undu, Atun );
3159 *Lrhs = *F[M_output_index];
3160 LOG(INFO) <<
"[CRB::checkResidual] residual (f,f) " << M_N-1 <<
":=" << M_model->scalarProduct( Frhs, Frhs ) <<
"\n";
3161 LOG(INFO) <<
"[CRB::checkResidual] residual (f,A) " << M_N-1 <<
":=" << 2*M_model->scalarProduct( Frhs, Aun ) <<
"\n";
3162 LOG(INFO) <<
"[CRB::checkResidual] residual (A,A) " << M_N-1 <<
":=" << M_model->scalarProduct( Aun, Aun ) <<
"\n";
3164 LOG(INFO) <<
"[CRB::checkResidual] residual (l,l) " << M_N-1 <<
":=" << M_model->scalarProduct( Lrhs, Lrhs ) <<
"\n";
3165 LOG(INFO) <<
"[CRB::checkResidual] residual (l,At) " << M_N-1 <<
":=" << 2*M_model->scalarProduct( Lrhs, Atun ) <<
"\n";
3166 LOG(INFO) <<
"[CRB::checkResidual] residual (At,At) " << M_N-1 <<
":=" << M_model->scalarProduct( Atun, Atun ) <<
"\n";
3172 vector_ptrtype __ef_pr( M_backend->newVector( M_model->functionSpace() ) );
3173 vector_ptrtype __ea_pr( M_backend->newVector( M_model->functionSpace() ) );
3174 M_model->l2solve( __ef_pr, Frhs );
3175 M_model->l2solve( __ea_pr, Aun );
3176 double check_C0_pr = M_model->scalarProduct( __ef_pr,__ef_pr );
3177 double check_Lambda_pr = 2*M_model->scalarProduct( __ef_pr,__ea_pr );
3178 double check_Gamma_pr = M_model->scalarProduct( __ea_pr,__ea_pr );
3180 vector_ptrtype __ef_du( M_backend->newVector( M_model->functionSpace() ) );
3181 vector_ptrtype __ea_du( M_backend->newVector( M_model->functionSpace() ) );
3182 M_model->l2solve( __ef_du, Lrhs );
3183 M_model->l2solve( __ea_du, Atun );
3184 double check_C0_du = M_model->scalarProduct( __ef_du,__ef_du );
3185 double check_Lambda_du = 2*M_model->scalarProduct( __ef_du,__ea_du );
3186 double check_Gamma_du = M_model->scalarProduct( __ea_du,__ea_du );
3188 double primal_sum = check_C0_pr + check_Lambda_pr + check_Gamma_pr ;
3189 double dual_sum = check_C0_du + check_Lambda_du + check_Gamma_du ;
3196 double C0_pr = primal_residual_coeffs[time_index][0];
3197 double Lambda_pr = primal_residual_coeffs[time_index][1];
3198 double Gamma_pr = primal_residual_coeffs[time_index][2];
3200 double C0_du,Lambda_du,Gamma_du;
3201 C0_du = dual_residual_coeffs[time_index][0];
3202 Lambda_du = dual_residual_coeffs[time_index][1];
3203 Gamma_du = dual_residual_coeffs[time_index][2];
3205 double err_C0_pr = math::abs( C0_pr - check_C0_pr ) ;
3206 double err_Lambda_pr = math::abs( Lambda_pr - check_Lambda_pr ) ;
3207 double err_Gamma_pr = math::abs( Gamma_pr - check_Gamma_pr ) ;
3208 double err_C0_du = math::abs( C0_du - check_C0_du ) ;
3209 double err_Lambda_du = math::abs( Lambda_du - check_Lambda_du ) ;
3210 double err_Gamma_du = math::abs( Gamma_du - check_Gamma_du ) ;
3212 int start_dual_index = 6;
3214 LOG(INFO)<<
"[CRB::checkResidual]";
3215 LOG(INFO)<<
"====primal coefficients==== ";
3216 LOG(INFO)<<
" c0_pr \t\t lambda_pr \t\t gamma_pr";
3217 LOG(INFO)<<
"computed : ";
3218 LOG(INFO)<<std::setprecision( 16 )<<primal_residual_coeffs[time_index][0]<<
"\t"<<primal_residual_coeffs[time_index][1]<<
"\t"<<primal_residual_coeffs[time_index][2];
3221 LOG(INFO)<<
"true : ";
3222 LOG(INFO)<<std::setprecision( 16 )<<check_C0_pr<<
"\t"<<check_Lambda_pr<<
"\t"<<check_Gamma_pr;
3223 LOG(INFO)<<
"====dual coefficients==== ";
3224 LOG(INFO)<<
" c0_du \t\t lambda_du \t\t gamma_du";
3225 LOG(INFO)<<
"computed : ";
3227 LOG(INFO)<<std::setprecision( 16 )<<dual_residual_coeffs[time_index][0]<<
"\t"<<dual_residual_coeffs[time_index][1]<<
"\t"<<dual_residual_coeffs[time_index][2];
3229 LOG(INFO)<<
"\ntrue : ";
3230 LOG(INFO)<<std::setprecision( 16 )<<check_C0_du<<
"\t"<<check_Lambda_du<<
"\t"<<check_Gamma_du;
3231 LOG(INFO)<<std::setprecision( 16 )<<
"primal_true_sum = "<<check_C0_pr+check_Lambda_pr+check_Gamma_pr;
3232 LOG(INFO)<<std::setprecision( 16 )<<
"dual_true_sum = "<<check_C0_du+check_Lambda_du+check_Gamma_du;
3234 LOG(INFO)<<std::setprecision( 16 )<<
"primal_computed_sum = "<<C0_pr+Lambda_pr+Gamma_pr;
3235 LOG(INFO)<<std::setprecision( 16 )<<
"dual_computed_sum = "<<C0_du+Lambda_du+Gamma_du;
3238 LOG(INFO)<<
"errors committed on coefficients ";
3239 LOG(INFO)<<std::setprecision( 16 )<<
"C0_pr : "<<err_C0_pr<<
"\tLambda_pr : "<<err_Lambda_pr<<
"\tGamma_pr : "<<err_Gamma_pr;
3240 LOG(INFO)<<std::setprecision( 16 )<<
"C0_du : "<<err_C0_du<<
"\tLambda_du : "<<err_Lambda_du<<
"\tGamma_du : "<<err_Gamma_du;
3241 LOG(INFO)<<
"and now relative error : ";
3242 double errC0pr = err_C0_pr/check_C0_pr;
3243 double errLambdapr = err_Lambda_pr/check_Lambda_pr;
3244 double errGammapr = err_Gamma_pr/check_Gamma_pr;
3245 double errC0du = err_C0_du/check_C0_du;
3246 double errLambdadu = err_Lambda_pr/check_Lambda_pr;
3247 double errGammadu = err_Gamma_pr/check_Gamma_pr;
3248 LOG(INFO)<<std::setprecision( 16 )<<errC0pr<<
"\t"<<errLambdapr<<
"\t"<<errGammapr;
3249 LOG(INFO)<<std::setprecision( 16 )<<errC0du<<
"\t"<<errLambdadu<<
"\t"<<errGammadu;
3256 vector_ptrtype __e_pr( M_backend->newVector( M_model->functionSpace() ) );
3257 M_model->l2solve( __e_pr, Aun );
3258 double dual_norm_pr = math::sqrt ( M_model->scalarProduct( __e_pr,__e_pr ) );
3259 LOG(INFO)<<
"[CRB::checkResidual] dual norm of primal residual without isolate terms (c0,lambda,gamma) = "<<dual_norm_pr<<
"\n";
3261 vector_ptrtype __e_du( M_backend->newVector( M_model->functionSpace() ) );
3262 M_model->l2solve( __e_du, Atun );
3263 double dual_norm_du = math::sqrt ( M_model->scalarProduct( __e_du,__e_du ) );
3264 LOG(INFO) <<
"[CRB::checkResidual] dual norm of dual residual without isolate terms = "<<dual_norm_du<<
"\n";
3266 double err_primal = math::sqrt ( M_model->scalarProduct( Aun, Aun ) );
3267 double err_dual = math::sqrt ( M_model->scalarProduct( Atun, Atun ) );
3268 LOG(INFO) <<
"[CRB::checkResidual] true primal residual for reduced basis function " << M_N-1 <<
":=" << err_primal <<
"\n";
3269 LOG(INFO) <<
"[CRB::checkResidual] true dual residual for reduced basis function " << M_N-1 <<
":=" << err_dual <<
"\n";
3274 template<
typename TruthModelType>
3279 if ( this->
vm()[
"crb.check.rb"].
template as<int>() == 0)
3282 std::cout <<
" -- check reduced basis\n";
3285 LOG(INFO) <<
"----------------------------------------------------------------------\n";
3291 LOG(INFO) <<
"**********************************************************************\n";
3293 std::vector< vectorN_type > uN;
3294 std::vector< vectorN_type > uNdu;
3295 std::vector< vectorN_type > uNold;
3296 std::vector< vectorN_type > uNduold;
3297 auto tuple =
lb( N, mu, uN, uNdu, uNold, uNduold );
3298 double s = tuple.template get<0>();
3299 auto error_estimation =
delta( N, mu, uN, uNdu , uNold, uNduold );
3300 double err = error_estimation.template get<0>() ;
3306 std::cout <<
"[check] error bounds are not < 1e-10\n";
3307 std::cout <<
"[check] k = " << k <<
"\n";
3308 std::cout <<
"[check] mu = " << mu <<
"\n";
3309 std::cout <<
"[check] delta = " << err <<
"\n";
3310 std::cout <<
"[check] uN( " << k <<
" ) = " << uN( k ) <<
"\n";
3314 u_fem = M_model->solveFemUsingOfflineEim ( mu );
3315 double sfem = M_model->output( M_output_index, mu , u_fem , need_to_solve );
3316 int size = mu.size();
3317 std::cout<<
" o mu = [ ";
3319 for (
int i=0; i<size-1; i++ ) std::cout<< mu[i] <<
" , ";
3321 std::cout<< mu[size-1]<<
" ]"<<std::endl;
3323 LOG(INFO) <<
"[check] s= " << s <<
" +- " << err <<
" | sfem= " << sfem <<
" | abs(sfem-srb) =" << math::abs( sfem - s ) <<
"\n";
3324 std::cout <<
"[check] s = " << s <<
" +- " << err <<
" | sfem= " << sfem <<
" | abs(sfem-srb) =" << math::abs( sfem - s )<<
"\n";
3328 LOG(INFO) <<
"----------------------------------------------------------------------\n";
3332 template<
typename TruthModelType>
3336 std::vector< vectorN_type > uN;
3337 std::vector< vectorN_type > uNdu;
3340 sampling_ptrtype Sampling;
3343 Sampling->randomize( N );
3346 int RBsize = M_WNmu->size();
3348 y_type ei( Sampling->size() );
3350 for (
size_type k = 0; k < Sampling->size(); ++k )
3353 double s =
lb( RBsize, mu, uN, uNdu );
3355 double sfem = M_model->output( M_output_index, mu , u_fem , need_to_solve);
3356 double error_estimation =
delta( RBsize,mu,uN,uNdu );
3357 ei( k ) = error_estimation/math::abs( sfem-s );
3358 std::cout<<
" efficiency indicator = "<<ei( k )<<
" for parameters {";
3360 for (
int i=0; i<mu.size(); i++ ) std::cout<<mu[i]<<
" ";
3362 std::cout<<
"} -- |sfem - s| = "<<math::abs( sfem-s )<<
" and error estimation = "<<error_estimation<<std::endl;
3365 Eigen::MatrixXf::Index index_max_ei;
3366 max_ei = ei.array().abs().maxCoeff( &index_max_ei );
3367 Eigen::MatrixXf::Index index_min_ei;
3368 min_ei = ei.array().abs().minCoeff( &index_min_ei );
3370 std::cout<<
"[EI] max_ei = "<<max_ei<<
" and min_ei = "<<min_ei<<min_ei<<
" sampling size = "<<N<<std::endl;
3378 template<
typename TruthModelType>
3380 CRB<TruthModelType>::correctionTerms(parameter_type
const& mu, std::vector< vectorN_type >
const & uN, std::vector< vectorN_type >
const & uNdu, std::vector<vectorN_type>
const & uNold,
int const k )
const
3382 int N = uN[0].size();
3384 matrixN_type Aprdu ( (
int)N, (
int)N ) ;
3385 matrixN_type Mprdu ( (
int)N, (
int)N ) ;
3386 vectorN_type Fdu ( (
int)N );
3387 vectorN_type du ( (
int)N );
3388 vectorN_type pr ( (
int)N );
3389 vectorN_type oldpr ( (
int)N );
3394 beta_vector_type betaAqm;
3395 beta_vector_type betaMqm;
3396 std::vector<beta_vector_type> betaFqm;
3398 double correction=0;
3400 if( M_model->isSteady() )
3402 Aprdu.setZero( N , N );
3405 boost::tie( betaMqm, betaAqm, betaFqm ) = M_model->computeBetaQm( mu ,time);
3407 for(
size_type q = 0;q < M_model->Ql(0); ++q)
3409 for(
int m=0; m < M_model->mMaxF(0,q); m++)
3410 Fdu += betaFqm[0][q][m]*M_Fqm_du[q][m].head(N);
3412 for(
size_type q = 0;q < M_model->Qa(); ++q)
3414 for(
int m=0; m < M_model->mMaxA(q); m++)
3415 Aprdu += betaAqm[q][m]*M_Aqm_pr_du[q][m].block(0,0,N,N);
3420 correction = -( Fdu.dot( du ) - du.dot( Aprdu*pr ) );
3425 double dt = M_model->timeStep();
3426 double Tf = M_model->timeFinal();
3430 for(
int kp=1; kp<=k; kp++)
3433 Aprdu.setZero( N , N );
3434 Mprdu.setZero( N , N );
3437 time_index = K-k+kp;
3438 time = time_index*dt;
3440 boost::tie( betaMqm, betaAqm, betaFqm) = M_model->computeBetaQm( mu ,time);
3445 for(
size_type q = 0;q < M_model->Ql(0); ++q)
3447 for(
int m=0; m < M_model->mMaxF(0,q); m++)
3448 Fdu += betaFqm[0][q][m]*M_Fqm_du[q][m].head(N);
3452 for(
size_type q = 0;q < M_model->Qa(); ++q)
3454 for(
int m=0; m < M_model->mMaxA(q); m++)
3455 Aprdu += betaAqm[q][m]*M_Aqm_pr_du[q][m].block(0,0,N,N);
3459 for(
size_type q = 0;q < M_model->Qm(); ++q)
3461 for(
int m=0; m<M_model->mMaxM(q); m++)
3462 Mprdu += betaMqm[q][m]*M_Mqm_pr_du[q][m].block(0,0,N,N);
3466 du = uNdu[K-1-time_index];
3467 pr = uN[time_index];
3468 oldpr = uNold[time_index];
3469 correction += dt*( Fdu.dot( du ) - du.dot( Aprdu*pr ) ) - du.dot(Mprdu*pr) + du.dot(Mprdu*oldpr) ;
3478 template<
typename TruthModelType>
3480 CRB<TruthModelType>::computeProjectionInitialGuess(
const parameter_type & mu,
int N , vectorN_type& initial_guess )
const
3482 VLOG(2) <<
"Compute projection of initial guess\n";
3483 beta_vector_type betaMqm;
3484 beta_vector_type beta_initial_guess;
3486 matrixN_type Mass ( (
int )N, (
int )N ) ;
3487 vectorN_type F ( (
int )N );
3490 beta_initial_guess = M_model->computeBetaInitialGuess( mu );
3498 Mass.setZero( N,N );
3499 int q_max = M_Mqm_pr.size();
3502 int m_max = M_Mqm_pr[q].size();
3503 for(
int m=0; m<m_max; m++)
3505 Mass += 1*M_Mqm_pr[q][m].block( 0,0,N,N );
3510 q_max = M_InitialGuessV_pr.size();
3513 int m_max = M_InitialGuessV_pr[q].size();
3514 for(
int m=0; m<m_max; m++)
3515 F += beta_initial_guess[q][m]*M_InitialGuessV_pr[q][m].head( N );
3518 initial_guess = Mass.lu().solve( F );
3521 template<
typename TruthModelType>
3523 CRB<TruthModelType>::updateJacobian(
const map_dense_vector_type& map_X, map_dense_matrix_type& map_J ,
const parameter_type & mu ,
int N)
const
3527 beta_vector_type betaJqm;
3528 boost::tie( boost::tuples::ignore, betaJqm, boost::tuples::ignore ) = M_model->computeBetaQm( this->
expansion( map_X , N , M_WN ), mu , 0 );
3529 for (
size_type q = 0; q < M_model->Qa(); ++q )
3531 for(
int m=0; m<M_model->mMaxA(q); m++)
3532 map_J += betaJqm[q][m]*M_Jqm_pr[q][m].block( 0,0,N,N );
3536 template<
typename TruthModelType>
3538 CRB<TruthModelType>::updateResidual(
const map_dense_vector_type& map_X, map_dense_vector_type& map_R ,
const parameter_type & mu,
int N )
const
3541 std::vector< beta_vector_type > betaRqm;
3542 boost::tie( boost::tuples::ignore, boost::tuples::ignore, betaRqm ) = M_model->computeBetaQm( this->
expansion( map_X , N , M_WN ), mu , 0 );
3543 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
3545 for(
int m=0; m<M_model->mMaxF(0,q); m++)
3546 map_R += betaRqm[0][q][m]*M_Rqm_pr[q][m].head( N );
3550 template<
typename TruthModelType>
3552 CRB<TruthModelType>::newton(
size_type N, parameter_type
const& mu , vectorN_type & uN ,
double& condition_number,
double& output)
const
3554 matrixN_type J ( (
int )N, (
int )N ) ;
3555 vectorN_type R ( (
int )N );
3557 double *r_data = R.data();
3558 double *j_data = J.data();
3559 double *uN_data = uN.data();
3561 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , 1> > map_R ( r_data, N );
3562 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , 1> > map_uN ( uN_data, N );
3563 Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic , Eigen::Dynamic> > map_J ( j_data, N , N );
3565 computeProjectionInitialGuess( mu , N , uN );
3567 M_nlsolver->map_dense_jacobian = boost::bind( &self_type::updateJacobian, boost::ref( *
this ), _1, _2 , mu , N );
3568 M_nlsolver->map_dense_residual = boost::bind( &self_type::updateResidual, boost::ref( *
this ), _1, _2 , mu , N );
3569 M_nlsolver->solve( map_J , map_uN , map_R, 1e-12, 100);
3573 if( option(_name=
"crb.compute-conditioning").
template as<bool>() )
3578 vectorN_type L ( (
int )N );
3579 std::vector<beta_vector_type> betaFqm;
3580 boost::tie( boost::tuples::ignore, boost::tuples::ignore, betaFqm ) = M_model->computeBetaQm( this->
expansion( uN , N , M_WN ), mu , 0 );
3582 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
3584 for(
int m=0; m < M_model->mMaxF(M_output_index,q); m++)
3585 L += betaFqm[M_output_index][q][m]*M_Lqm_pr[q][m].head( N );
3588 output = L.dot( uN );
3593 template<
typename TruthModelType>
3597 Eigen::SelfAdjointEigenSolver< matrixN_type > eigen_solver;
3598 eigen_solver.compute( A );
3599 int number_of_eigenvalues = eigen_solver.eigenvalues().size();
3601 std::vector<double> eigen_values( number_of_eigenvalues );
3604 for (
int i=0; i<number_of_eigenvalues; i++ )
3606 if ( imag( eigen_solver.eigenvalues()[i] )>1e-12 )
3608 throw std::logic_error(
"[CRB::lb] ERROR : complex eigenvalues were found" );
3611 eigen_values[i]=real( eigen_solver.eigenvalues()[i] );
3614 int position_of_largest_eigenvalue=number_of_eigenvalues-1;
3615 int position_of_smallest_eigenvalue=0;
3616 double eig_max = eigen_values[position_of_largest_eigenvalue];
3617 double eig_min = eigen_values[position_of_smallest_eigenvalue];
3618 return eig_max / eig_min;
3623 template<
typename TruthModelType>
3628 double time_for_output;
3631 int number_of_time_step=1;
3634 if ( M_model->isSteady() )
3637 time_for_output = 1e30;
3645 time_step = M_model->timeStep();
3646 time_final = M_model->timeFinal();
3649 time_for_output = K * time_step;
3653 number_of_time_step = time_final / time_step;
3654 time_for_output = number_of_time_step * time_step;
3658 beta_vector_type betaAqm;
3659 beta_vector_type betaMqm;
3660 beta_vector_type betaMFqm;
3661 std::vector<beta_vector_type> betaFqm, betaLqm;
3663 matrixN_type Adu ( (
int )N, (
int )N ) ;
3664 matrixN_type Mdu ( (
int )N, (
int )N ) ;
3665 vectorN_type Fdu ( (
int )N );
3666 vectorN_type Ldu ( (
int )N );
3668 matrixN_type Aprdu( (
int )N, (
int )N );
3669 matrixN_type Mprdu( (
int )N, (
int )N );
3672 if ( M_model->isSteady() )
3676 boost::tie( betaMqm, betaAqm, betaFqm ) = M_model->computeBetaQm( mu ,time );
3680 for (
size_type q = 0; q < M_model->Qa(); ++q )
3682 for(
int m=0; m < M_model->mMaxA(q); m++)
3683 Adu += betaAqm[q][m]*M_Aqm_du[q][m].block( 0,0,N,N );
3686 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
3688 for(
int m=0; m < M_model->mMaxF(M_output_index,q); m++)
3689 Ldu += betaFqm[M_output_index][q][m]*M_Lqm_du[q][m].head( N );
3692 uNdu[0] = Adu.lu().solve( -Ldu );
3697 int time_index = number_of_time_step-1;
3700 double initial_dual_time = time_for_output+time_step;
3702 boost::tie( betaMqm, betaAqm, betaFqm, betaMFqm ) = M_model->computeBetaQm( mu ,initial_dual_time );
3705 for (
size_type q = 0; q < M_model->Qm(); ++q )
3707 for (
size_type m = 0; m < M_model->mMaxM( q ); ++m )
3709 for(
int m=0; m < M_model->mMaxM(q); m++)
3710 Mdu += betaMqm[q][m]*M_Mqm_du[q][m].block( 0,0,N,N );
3716 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
3718 for (
size_type m = 0; m < M_model->mMaxF( M_output_index , q ); ++m )
3720 for(
int m=0; m < M_model->mMaxF(M_output_index,q); m++)
3721 Ldu += betaFqm[M_output_index][q][m]*M_Lqm_du[q][m].head( N );
3735 uNduold[time_index]( n ) = M_coeff_du_ini_online[n];
3740 for ( time=time_for_output; time>=time_step; time-=time_step )
3743 boost::tie( betaMqm, betaAqm, betaFqm ) = M_model->computeBetaQm( mu ,time );
3746 for (
size_type q = 0; q < M_model->Qa(); ++q )
3748 for(
int m=0; m < M_model->mMaxA(q); m++)
3749 Adu += betaAqm[q][m]*M_Aqm_du[q][m].block( 0,0,N,N );
3757 for(
int m=0; m < M_model->mMaxM(q); m++)
3759 Adu += betaMqm[q][m]*M_Mqm_du[q][m].block( 0,0,N,N )/time_step;
3760 Fdu += betaMqm[q][m]*M_Mqm_du[q][m].block( 0,0,N,N )*uNduold[time_index]/time_step;
3764 uNdu[time_index] = Adu.lu().solve( Fdu );
3768 uNduold[time_index-1] = uNdu[time_index];
3778 template<
typename TruthModelType>
3780 CRB<TruthModelType>::fixedPointPrimal(
size_type N, parameter_type
const& mu, std::vector< vectorN_type > & uN, std::vector<vectorN_type> & uNold,
double& condition_number, std::vector< double > & output_vector ,
int K)
const
3783 double time_for_output;
3787 int number_of_time_step=1;
3792 if ( M_model->isSteady() )
3795 time_for_output = 1e30;
3803 time_step = M_model->timeStep();
3804 time_final = M_model->timeFinal();
3807 time_for_output = K * time_step;
3811 number_of_time_step = time_final / time_step;
3812 time_for_output = number_of_time_step * time_step;
3815 beta_vector_type betaAqm;
3816 beta_vector_type betaMqm;
3817 beta_vector_type betaMFqm;
3818 beta_vector_type beta_initial_guess;
3820 std::vector<beta_vector_type> betaFqm, betaLqm;
3822 matrixN_type A ( (
int )N, (
int )N ) ;
3823 vectorN_type F ( (
int )N );
3824 vectorN_type L ( (
int )N );
3826 if ( !M_model->isSteady() )
3829 uNold[0]( n ) = M_coeff_pr_ini_online[n];
3833 int max_fixedpoint_iterations = option(
"crb.max-fixedpoint-iterations").template as<int>();
3834 double increment_fixedpoint_tol = option(
"crb.increment-fixedpoint-tol").template as<double>();
3835 double output_fixedpoint_tol = option(
"crb.output-fixedpoint-tol").template as<double>();
3836 bool fixedpoint_verbose = option(
"crb.fixedpoint-verbose").template as<bool>();
3837 double fixedpoint_critical_value = option(_name=
"crb.fixedpoint-critical-value").template as<double>();
3838 for (
double time=time_step; time<time_for_output+time_step; time+=time_step )
3841 computeProjectionInitialGuess( mu , N , uN[time_index] );
3848 VLOG(2) <<
"lb: start fix point\n";
3850 vectorN_type previous_uN( M_N );
3858 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
3860 for(
int m=0; m < M_model->mMaxF(M_output_index,q); m++)
3862 L += betaFqm[M_output_index][q][m]*M_Lqm_pr[q][m].head( N );
3865 old_output = L.dot( uN[time_index] );
3871 boost::tie( betaMqm, betaAqm, betaFqm ) = M_model->computeBetaQm( this->
expansion( uN[time_index] , N , M_WN ), mu ,time );
3874 for (
size_type q = 0; q < M_model->Qa(); ++q )
3876 for(
int m=0; m<M_model->mMaxA(q); m++)
3877 A += betaAqm[q][m]*M_Aqm_pr[q][m].block( 0,0,N,N );
3881 for (
size_type q = 0; q < M_model->Ql( 0 ); ++q )
3883 for(
int m=0; m<M_model->mMaxF(0,q); m++)
3884 F += betaFqm[0][q][m]*M_Fqm_pr[q][m].head( N );
3889 for(
int m=0; m<M_model->mMaxM(q); m++)
3891 A += betaMqm[q][m]*M_Mqm_pr[q][m].block( 0,0,N,N )/time_step;
3892 F += betaMqm[q][m]*M_Mqm_pr[q][m].block( 0,0,N,N )*uNold[time_index]/time_step;
3897 previous_uN = uN[time_index];
3900 uN[time_index] = A.lu().solve( F );
3902 if ( time_index<number_of_time_step-1 )
3903 uNold[time_index+1] = uN[time_index];
3906 for (
size_type q = 0; q < M_model->Ql( M_output_index ); ++q )
3908 for(
int m=0; m < M_model->mMaxF(M_output_index,q); m++)
3910 L += betaFqm[M_output_index][q][m]*M_Lqm_pr[q][m].head( N );
3913 old_output = output;
3914 output = L.dot( uN[time_index] );
3917 output_vector[time_index] = output;
3918 DVLOG(2) <<
"iteration " << fi <<
" increment error: " << (uN[time_index]-previous_uN).norm() <<
"\n";
3921 if( fixedpoint_verbose && this->worldComm().globalRank()==this->worldComm().masterRank() )
3922 VLOG(2)<<
"[CRB::lb] fixedpoint iteration " << fi <<
" increment error: " << (uN[time_index]-previous_uN).norm()<<std::endl;
3924 double residual_norm = (A * uN[time_index] - F).norm() ;
3925 VLOG(2) <<
" residual_norm : "<<residual_norm;
3927 while ( (uN[time_index]-previous_uN).norm() > increment_fixedpoint_tol && fi<max_fixedpoint_iterations );
3930 if( (uN[time_index]-previous_uN).norm() > increment_fixedpoint_tol )
3931 DVLOG(2)<<
"[CRB::lb] fixed point, proc "<<this->worldComm().globalRank()
3932 <<
" fixed point has no converged : norm(uN-uNold) = "<<(uN[time_index]-previous_uN).norm()
3933 <<
" and tolerance : "<<increment_fixedpoint_tol<<
" so "<<max_fixedpoint_iterations<<
" iterations were done"<<std::endl;
3935 if( (uN[time_index]-previous_uN).norm() > fixedpoint_critical_value )
3936 throw std::logic_error(
"[CRB::lb] fixed point ERROR : norm(uN-uNold) > critical value " );
3938 if ( time_index<number_of_time_step-1 )
3942 condition_number = 0;
3943 if( option(_name=
"crb.compute-conditioning").template as<bool>() )
3947 template<
typename TruthModelType>
3949 CRB<TruthModelType>::fixedPoint(
size_type N,
parameter_type const& mu, std::vector< vectorN_type > & uN, std::vector< vectorN_type > & uNdu, std::vector<vectorN_type> & uNold, std::vector<vectorN_type> & uNduold,
double& condition_number, std::vector< double > & output_vector,
int K)
const
3952 double time_for_output;
3955 int number_of_time_step;
3957 if ( M_model->isSteady() )
3960 time_for_output = 1e30;
3961 number_of_time_step=1;
3965 time_step = M_model->timeStep();
3966 time_final = M_model->timeFinal();
3968 time_for_output = K * time_step;
3971 number_of_time_step = time_final / time_step;
3972 time_for_output = number_of_time_step * time_step;
3975 fixedPointPrimal( N, mu , uN , uNold , condition_number, output_vector, K ) ;
3977 int size=output_vector.size();
3978 double o =output_vector[size-1];
3979 bool solve_dual_problem = this->
vm()[
"crb.solve-dual-problem"].template as<bool>();
3983 if( solve_dual_problem )
3985 fixedPointDual( N, mu , uNdu , uNduold , output_vector , K ) ;
3989 for (
double time=time_step; time<=time_for_output; time+=time_step )
3991 int k = time_index+1;
3992 output_vector[time_index]+=correctionTerms(mu, uN , uNdu, uNold, k );
3999 template<
typename TruthModelType>
4000 boost::tuple<double,double>
4004 bool save_output_behavior = this->
vm()[
"crb.save-output-behavior"].template as<bool>();
4009 double time_for_output;
4013 int number_of_time_step = 0;
4016 if ( M_model->isSteady() )
4019 time_for_output = 1e30;
4021 number_of_time_step=1;
4027 time_step = M_model->timeStep();
4028 time_final = M_model->timeFinal();
4031 time_for_output = K * time_step;
4035 number_of_time_step = time_final / time_step;
4036 time_for_output = number_of_time_step * time_step;
4039 if ( N > M_N ) N = M_N;
4041 uN.resize( number_of_time_step );
4042 uNdu.resize( number_of_time_step );
4043 uNold.resize( number_of_time_step );
4044 uNduold.resize( number_of_time_step );
4047 BOOST_FOREACH(
auto elem, uN )
4049 uN[index].resize( N );
4050 uNdu[index].resize( N );
4051 uNold[index].resize( N );
4052 uNduold[index].resize( N );
4055 double condition_number;
4059 std::vector<double>output_vector;
4060 output_vector.resize( number_of_time_step );
4068 newton( N , mu , uN[0] , condition_number , output_vector[0] );
4070 fixedPoint( N , mu , uN , uNdu , uNold , uNduold , condition_number , output_vector , K );
4073 if( M_compute_variance )
4076 if( ! M_database_contains_variance_info )
4077 throw std::logic_error(
"[CRB::offline] ERROR there are no information available in the DataBase for variance computing, please set option crb.save-information-for-variance=true and rebuild the database" );
4079 int nb_spaces = functionspace_type::nSpaces;
4084 for (
double time=time_step; time<=time_for_output; time+=time_step )
4086 vectorN_type uNsquare = uN[time_index].array().pow(2);
4087 double first = uNsquare.dot( M_variance_matrix_phi[space].block(0,0,N,N).diagonal() );
4090 for(
int k = 1; k <= N-1; ++k)
4092 for(
int j = 1; j <= N-k; ++j)
4094 second += 2 * uN[time_index](k-1) * uN[time_index](k+j-1) * M_variance_matrix_phi[space](k-1 , j+k-1) ;
4097 output_vector[time_index] = first + second;
4104 if ( save_output_behavior )
4107 std::ofstream file_output;
4110 for (
int i=0; i<mu.size(); i++ )
4112 mu_str= mu_str + ( boost::format(
"_%1%" ) %mu[i] ).str() ;
4115 std::string
name =
"output_evolution" + mu_str;
4116 file_output.open( name.c_str(),std::ios::out );
4118 for (
double time=time_step; time<=time_for_output; time+=time_step )
4120 file_output<<time<<
"\t"<<output_vector[time_index]<<
"\n";
4124 file_output.close();
4127 int size=output_vector.size();
4128 return boost::make_tuple( output_vector[size-1], condition_number);
4132 template<
typename TruthModelType>
4133 typename CRB<TruthModelType>::error_estimation_type
4136 std::vector<vectorN_type>
const& uN,
4137 std::vector<vectorN_type>
const& uNdu,
4138 std::vector<vectorN_type>
const& uNold,
4139 std::vector<vectorN_type>
const& uNduold,
4143 std::vector< std::vector<double> > primal_residual_coeffs;
4144 std::vector< std::vector<double> > dual_residual_coeffs;
4148 if ( M_error_type == CRB_NO_RESIDUAL )
4149 return boost::make_tuple( -1,primal_residual_coeffs,dual_residual_coeffs,delta_pr,delta_du );
4151 else if ( M_error_type == CRB_EMPIRICAL )
4152 return boost::make_tuple( empiricalErrorEstimation ( N, mu , k ) , primal_residual_coeffs, dual_residual_coeffs , delta_pr, delta_du);
4157 double Tf = M_model->timeFinal();
4159 double dt = M_model->timeStep();
4162 double primal_sum=0;
4167 primal_residual_coeffs.resize( K );
4168 dual_residual_coeffs.resize( K );
4169 for (
double time=dt; time<=Tf; time+=dt )
4172 primal_sum += pr.template get<0>();
4173 primal_residual_coeffs[time_index].resize( pr.template get<1>().size() );
4174 primal_residual_coeffs[time_index] = pr.template get<1>() ;
4180 double dual_residual=0;
4182 if ( !M_model->isSteady() ) dual_residual = initialDualResidual( N,mu,uNduold[time_index],dt );
4183 bool solve_dual_problem = this->
vm()[
"crb.solve-dual-problem"].template as<bool>() ;
4185 if( solve_dual_problem )
4187 for (
double time=Tf; time>=dt; time-=dt )
4189 auto du = transientDualResidual( N, mu, uNdu[time_index], uNduold[time_index], dt, time );
4190 dual_sum += du.template get<0>();
4191 dual_residual_coeffs[time_index].resize( du.template get<1>().size() );
4192 dual_residual_coeffs[time_index] = du.template get<1>();
4198 bool show_residual = this->
vm()[
"crb.show-residual"].template as<bool>() ;
4199 if( ! M_offline_step && show_residual )
4202 bool seek_mu_in_complement = this->
vm()[
"crb.seek-mu-in-complement"].template as<bool>() ;
4203 LOG( INFO ) <<
" =========== Residual with "<<N<<
" basis functions - seek mu in complement of WNmu : "<<seek_mu_in_complement<<
"============ \n";
4205 for (
double time=dt; time<=Tf; time+=dt )
4209 sum+=pr.template get<0>();
4212 LOG(INFO) <<
"sum of primal residuals "<<sum<<std::endl;
4217 if (solve_dual_problem )
4219 for (
double time=Tf; time>=dt; time-=dt )
4221 auto du = transientDualResidual( N, mu, uNdu[time_index], uNduold[time_index], dt, time );
4223 sum += du.template get<0>();
4227 LOG(INFO) <<
"sum of dual residuals "<<sum<<std::endl;
4228 LOG( INFO ) <<
" ================================= \n";
4232 double alphaA=1,alphaM=1;
4234 if ( M_error_type == CRB_RESIDUAL_SCM )
4236 double alphaA_up, lbti;
4237 M_scmA->setScmForMassMatrix(
false );
4238 boost::tie( alphaA, lbti ) = M_scmA->lb( mu );
4239 if( option(_name=
"crb.scm.use-scm").
template as<bool>() )
4240 boost::tie( alphaA_up, lbti ) = M_scmA->ub( mu );
4243 if ( ! M_model->isSteady() )
4245 M_scmM->setScmForMassMatrix(
true );
4246 double alphaM_up, lbti;
4247 boost::tie( alphaM, lbti ) = M_scmM->lb( mu );
4248 if( option(_name=
"crb.scm.use-scm").
template as<bool>() )
4249 boost::tie( alphaM_up, lbti ) = M_scmM->ub( mu );
4254 double output_upper_bound;
4255 double solution_upper_bound;
4256 double solution_dual_upper_bound;
4259 if ( M_model->isSteady() )
4261 delta_pr = math::sqrt( primal_sum ) / math::sqrt( alphaA );
4262 if( solve_dual_problem )
4263 delta_du = math::sqrt( dual_sum ) / math::sqrt( alphaA );
4266 output_upper_bound = delta_pr * delta_du;
4272 delta_pr = math::sqrt( dt/alphaA * primal_sum );
4273 delta_du = math::sqrt( dt/alphaA * dual_sum + dual_residual/alphaM );
4274 output_upper_bound = delta_pr * delta_du;
4280 return boost::make_tuple( output_upper_bound, primal_residual_coeffs, dual_residual_coeffs , delta_pr, delta_du );
4285 template<
typename TruthModelType>
4286 typename CRB<TruthModelType>::max_error_type
4289 bool seek_mu_in_complement = this->
vm()[
"crb.seek-mu-in-complement"].template as<bool>() ;
4291 std::vector< vectorN_type > uN;
4292 std::vector< vectorN_type > uNdu;
4293 std::vector< vectorN_type > uNold;
4294 std::vector< vectorN_type > uNduold;
4296 y_type err( M_Xi->size() );
4297 y_type vect_delta_pr( M_Xi->size() );
4298 y_type vect_delta_du( M_Xi->size() );
4299 std::vector<double> check_err( M_Xi->size() );
4303 if ( M_error_type == CRB_EMPIRICAL )
4305 if ( M_WNmu->size()==1 )
4309 boost::tie ( mu,
id ) = M_Xi->max();
4310 return boost::make_tuple( 1e5, M_Xi->at(
id ), id , delta_pr, delta_du);
4315 err.resize( M_WNmu_complement->size() );
4316 check_err.resize( M_WNmu_complement->size() );
4318 for (
size_type k = 0; k < M_WNmu_complement->size(); ++k )
4322 auto error_estimation =
delta( N, mu, uN, uNdu, uNold, uNduold, k );
4323 double _err = error_estimation.template get<0>();
4324 vect_delta_pr( k ) = error_estimation.template get<3>();
4325 vect_delta_du( k ) = error_estimation.template get<4>();
4327 check_err[k] = _err;
4335 if ( seek_mu_in_complement )
4337 err.resize( M_WNmu_complement->size() );
4338 check_err.resize( M_WNmu_complement->size() );
4341 for (
size_type k = 0; k < M_WNmu_complement->size(); ++k )
4344 lb( N, mu, uN, uNdu , uNold ,uNduold );
4345 auto error_estimation =
delta( N, mu, uN, uNdu, uNold, uNduold, k );
4346 double _err = error_estimation.template get<0>();
4347 vect_delta_pr( k ) = error_estimation.template get<3>();
4348 vect_delta_du( k ) = error_estimation.template get<4>();
4350 check_err[k] = _err;
4357 for (
size_type k = 0; k < M_Xi->size(); ++k )
4362 lb( N, mu, uN, uNdu , uNold ,uNduold );
4366 auto error_estimation =
delta( N, mu, uN, uNdu, uNold, uNduold, k );
4367 double _err = error_estimation.template get<0>();
4368 vect_delta_pr( k ) = error_estimation.template get<3>();
4369 vect_delta_du( k ) = error_estimation.template get<4>();
4372 check_err[k] = _err;
4377 Eigen::MatrixXf::Index index;
4378 double maxerr = err.array().abs().maxCoeff( &index );
4379 delta_pr = vect_delta_pr( index );
4380 delta_du = vect_delta_du( index );
4383 std::vector<double>::iterator it = std::max_element( check_err.begin(), check_err.end() );
4384 int check_index = it - check_err.begin() ;
4385 double check_maxerr = *it;
4387 if ( index != check_index || maxerr != check_maxerr )
4389 std::cout<<
"[CRB::maxErrorBounds] index = "<<index<<
" / check_index = "<<check_index<<
" and maxerr = "<<maxerr<<
" / "<<check_maxerr<<std::endl;
4390 throw std::logic_error(
"[CRB::maxErrorBounds] index and check_index have different values" );
4395 if ( seek_mu_in_complement )
4397 LOG(INFO) <<
"[maxErrorBounds] WNmu_complement N=" << N <<
" max Error = " << maxerr <<
" at index = " << index <<
"\n";
4398 mu = M_WNmu_complement->at( index );
4399 _index = M_WNmu_complement->indexInSuperSampling( index );
4404 LOG(INFO) <<
"[maxErrorBounds] N=" << N <<
" max Error = " << maxerr <<
" at index = " << index <<
"\n";
4405 mu = M_Xi->at( index );
4409 int proc_number = this->worldComm().globalRank();
4410 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
4411 std::cout<< std::setprecision(15)<<
"[CRB maxerror] proc "<< proc_number<<
" delta_pr : "<<delta_pr<<
" - delta_du : "<<delta_du<<
" at index : "<<_index<<std::endl;
4412 lb( N, mu, uN, uNdu , uNold ,uNduold );
4414 return boost::make_tuple( maxerr, mu , _index , delta_pr, delta_du);
4418 template<
typename TruthModelType>
4422 int proc_number = this->worldComm().globalRank();
4423 if( proc_number == 0 ) std::cout <<
" -- orthonormalization (Gram-Schmidt)\n";
4424 DVLOG(2) <<
"[CRB::orthonormalize] orthonormalize basis for N=" << N <<
"\n";
4425 DVLOG(2) <<
"[CRB::orthonormalize] orthonormalize basis for WN="
4426 << wn.size() <<
"\n";
4427 DVLOG(2) <<
"[CRB::orthonormalize] starting ...\n";
4431 for (
size_type j = std::max( i+1,N-Nm ); j < N; ++j )
4433 value_type __rij_pr = M_model->scalarProduct( wn[i], wn[ j ] );
4434 wn[j].add( -__rij_pr, wn[i] );
4440 value_type __rii_pr = math::sqrt( M_model->scalarProduct( wn[i], wn[i] ) );
4441 wn[i].scale( 1./__rii_pr );
4444 DVLOG(2) <<
"[CRB::orthonormalize] finished ...\n";
4445 DVLOG(2) <<
"[CRB::orthonormalize] copying back results in basis\n";
4449 return checkOrthonormality( N , wn );
4452 template <
typename TruthModelType>
4459 throw std::logic_error(
"[CRB::checkOrthonormality] ERROR : size of wn is zero" );
4462 if ( orthonormalize_primal*orthonormalize_dual==0 )
4464 LOG(INFO)<<
"Warning : calling checkOrthonormality is called but ";
4465 LOG(INFO)<<
" orthonormalize_dual = "<<orthonormalize_dual;
4466 LOG(INFO)<<
" and orthonormalize_primal = "<<orthonormalize_primal;
4471 I.setIdentity( N, N );
4473 for (
int i = 0; i < N; ++i )
4475 for (
int j = 0; j < N; ++j )
4477 A( i, j ) = M_model->scalarProduct( wn[i], wn[j] );
4482 DVLOG(2) <<
"orthonormalization: " << A.norm() <<
"\n";
4483 if ( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
4485 LOG( INFO ) <<
" o check : " << A.norm() <<
" (should be 0)";
4493 template <
typename TruthModelType>
4499 std::vector<wn_type> vect_wn=export_vector_wn.template get<0>();
4500 std::vector<std::string> vect_names=export_vector_wn.template get<1>();
4502 if ( vect_wn.size()==0 )
4504 throw std::logic_error(
"[CRB::exportBasisFunctions] ERROR : there are no wn_type to export" );
4508 auto first_wn = vect_wn[0];
4509 auto first_element = first_wn[0];
4511 exporter->step( 0 )->setMesh( first_element.functionSpace()->mesh() );
4513 BOOST_FOREACH(
auto wn , vect_wn )
4518 throw std::logic_error(
"[CRB::exportBasisFunctions] ERROR : there are no element to export" );
4521 int element_number=0;
4524 BOOST_FOREACH(
auto element, wn )
4527 std::string basis_name = vect_names[basis_number];
4528 std::string number = ( boost::format(
"%1%_with_parameters" ) %element_number ).str() ;
4529 mu = M_WNmu->at( element_number );
4532 for (
int i=0; i<mu.size(); i++ )
4534 mu_str= mu_str + ( boost::format(
"_%1%" ) %mu[i] ).str() ;
4537 std::string
name = basis_name + number + mu_str;
4538 exporter->step( 0 )->add( name, element );
4550 template <
typename TruthModelType>
4551 typename CRB<TruthModelType>::value_type
4554 std::vector<vectorN_type> Un;
4555 std::vector<vectorN_type> Undu;
4556 std::vector<vectorN_type> Unold;
4557 std::vector<vectorN_type> Unduold;
4558 auto o =
lb( Nwn, mu, Un, Undu, Unold, Unduold );
4559 double sn = o.template get<0>();
4561 int nb_element =Nwn/M_factor*( M_factor>0 ) + ( Nwn+M_factor )*( M_factor<0 && ( int )Nwn>( -M_factor ) ) + 1*( M_factor<0 && ( int )Nwn<=( -M_factor ) ) ;
4563 std::vector<vectorN_type> Un2;
4564 std::vector<vectorN_type> Undu2;
4567 auto tuple =
lb( nb_element, mu, Un2, Undu2, Unold, Unduold );
4568 double output_smaller_basis = tuple.template get<0>();
4570 double error_estimation = math::abs( sn-output_smaller_basis );
4572 return error_estimation;
4579 template<
typename TruthModelType>
4580 typename CRB<TruthModelType>::value_type
4581 CRB<TruthModelType>::initialDualResidual(
int Ncur, parameter_type
const& mu, vectorN_type
const& Unduini,
double time_step )
const
4583 beta_vector_type betaAqm;
4584 beta_vector_type betaMqm;
4585 beta_vector_type betaMFqm;
4586 std::vector<beta_vector_type> betaFqm;
4587 double time = M_model->timeFinal();
4588 boost::tie( betaMqm, betaAqm, betaFqm) = M_model->computeBetaQm( mu, time );
4590 int __QLhs = M_model->Qa();
4591 int __QOutput = M_model->Ql( M_output_index );
4592 int __Qm = M_model->Qm();
4596 value_type __c0_du = 0.0;
4598 for (
int __q1 = 0; __q1 < __QOutput; ++__q1 )
4600 for (
int __m1 = 0; __m1 < M_model->mMaxF(M_output_index,__q1); ++__m1 )
4602 value_type fq1 = betaFqm[M_output_index][__q1][__m1];
4604 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
4606 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
4608 value_type fq2 = betaFqm[M_output_index][__q2][__m2];
4611 __c0_du += M_C0_du[__q1][__m1][__q2][__m2]*fq1*fq2;
4618 value_type __Caf_du = 0.0;
4619 value_type __Caa_du = 0.0;
4621 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
4623 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1) ; ++__m1 )
4626 value_type a_q1 = betaAqm[__q1][__m1];
4628 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
4631 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
4633 value_type f_q2 = betaFqm[M_output_index][__q2][__m2];
4634 __Caf_du += 1./time_step * a_q1*f_q2*M_Lambda_du[__q1][__m1][__q2][__m2].dot( Unduini );
4639 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
4641 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
4643 value_type a_q2 = betaAqm[__q2][__m2];
4644 auto m = M_Gamma_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Unduini;
4645 __Caa_du += a_q1 * a_q2 * Unduini.dot( m );
4652 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
4654 for (
int __m2=0 ; __m2< M_model->mMaxA(__q2); ++__m2 )
4656 value_type a_q2 = betaAqm[__q2][__m2];
4657 auto m = M_Cma_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Unduini;
4659 __Cma_du += m_q1 * a_q2 * Unduini.dot( m );
4665 value_type __Cmf_du=0;
4666 value_type __Cmm_du=0;
4667 value_type __Cma_du=0;
4669 for (
int __q1=0 ; __q1<__Qm; ++__q1 )
4671 for (
int __m1=0 ; __m1< M_model->mMaxM(__q1); ++__m1 )
4674 value_type m_q1 = betaMqm[__q1][__m1];
4676 for (
int __q2=0 ; __q2<__QOutput; ++__q2 )
4678 for (
int __m2=0 ; __m2< M_model->mMaxF(M_output_index,__q2); ++__m2 )
4680 value_type f_q2 = betaFqm[M_output_index][__q2][__m2];
4682 __Cmf_du += m_q1 * f_q2 * M_Cmf_du_ini[__q1][__m1][__q2][__m2].head( __N ).dot( Unduini );
4686 for (
int __q2 = 0; __q2 < __Qm; ++__q2 )
4688 for (
int __m2=0 ; __m2< M_model->mMaxM(__q2); ++__m2 )
4690 value_type m_q2 = betaMqm[__q2][__m2];
4691 auto m = M_Cmm_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Unduini;
4693 __Cmm_du += m_q1 * m_q2 * Unduini.dot( m );
4701 return math::abs( __c0_du+__Cmf_du+__Cmm_du ) ;
4708 template<
typename TruthModelType>
4709 typename CRB<TruthModelType>::residual_error_type
4713 beta_vector_type betaAqm;
4714 beta_vector_type betaMqm;
4715 beta_vector_type betaMFqm;
4716 std::vector<beta_vector_type> betaFqm, betaLqm;
4718 boost::tie( betaMqm, betaAqm, betaFqm ) = M_model->computeBetaQm( mu, time );
4727 residual_error_type steady_residual_contribution = steadyPrimalResidual( Ncur, mu, Un, time );
4728 std::vector<double> steady_coeff_vector = steady_residual_contribution.template get<1>();
4729 value_type __c0_pr = steady_coeff_vector[0];
4730 value_type __lambda_pr = steady_coeff_vector[1];
4731 value_type __gamma_pr = steady_coeff_vector[2];
4733 value_type __Cmf_pr=0;
4734 value_type __Cma_pr=0;
4735 value_type __Cmm_pr=0;
4738 if ( ! M_model->isSteady() )
4741 for (
size_type __q1=0 ; __q1<__Qm; ++__q1 )
4743 for (
size_type __m1=0 ; __m1< M_model->mMaxM(__q1); ++__m1 )
4745 value_type m_q1 = betaMqm[__q1][__m1];
4747 for (
size_type __q2=0 ; __q2<__QRhs; ++__q2 )
4749 for (
size_type __m2=0 ; __m2< M_model->mMaxF(0,__q2); ++__m2 )
4751 value_type f_q2 = betaFqm[0][__q2][__m2];
4752 __Cmf_pr += 1./time_step * m_q1 * f_q2 * M_Cmf_pr[__q1][__m1][__q2][__m2].head( __N ).dot( Un );
4753 __Cmf_pr -= 1./time_step * m_q1 * f_q2 * M_Cmf_pr[__q1][__m1][__q2][__m2].head( __N ).dot( Unold );
4757 for (
size_type __q2 = 0; __q2 < __QLhs; ++__q2 )
4759 for (
size_type __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
4761 value_type a_q2 = betaAqm[__q2][__m2];
4762 auto m = M_Cma_pr[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Un;
4763 __Cma_pr += 1./time_step * m_q1 * a_q2 * Un.dot( m );
4764 __Cma_pr -= 1./time_step * m_q1 * a_q2 * Unold.dot( m );
4768 for (
size_type __q2 = 0; __q2 < __Qm; ++__q2 )
4770 for (
size_type __m2 = 0; __m2 < M_model->mMaxM(__q2); ++__m2 )
4772 value_type m_q2 = betaMqm[__q2][__m2];
4773 auto m1 = M_Cmm_pr[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Un;
4774 auto m2 = M_Cmm_pr[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Unold;
4775 __Cmm_pr += 1./( time_step*time_step ) * m_q1 * m_q2 * Un.dot( m1 );
4776 __Cmm_pr -= 1./( time_step*time_step ) * m_q1 * m_q2 * Un.dot( m2 );
4777 __Cmm_pr -= 1./( time_step*time_step ) * m_q1 * m_q2 * Unold.dot( m1 );
4778 __Cmm_pr += 1./( time_step*time_step ) * m_q1 * m_q2 * Unold.dot( m2 );
4786 std::cout<<
"[transientResidual] time "<<time<<std::endl;
4787 std::cout<<
"Un = \n"<<Un<<std::endl;
4788 std::cout<<
"Unold = \n"<<Unold<<std::endl;
4789 std::cout<<
"__C0_pr = "<<__c0_pr<<std::endl;
4790 std::cout<<
"__lambda_pr = "<<__lambda_pr<<std::endl;
4791 std::cout<<
"__gamma_pr = "<<__gamma_pr<<std::endl;
4792 std::cout<<
"__Cmf_pr = "<<__Cmf_pr<<std::endl;
4793 std::cout<<
"__Cma_pr = "<<__Cma_pr<<std::endl;
4794 std::cout<<
"__Cmm_pr = "<<__Cmm_pr<<std::endl;
4795 std::cout<<
"__N = "<<__N<<std::endl;
4800 value_type delta_pr = math::abs( __c0_pr+__lambda_pr+__gamma_pr+__Cmf_pr+__Cma_pr+__Cmm_pr ) ;
4801 std::vector<double> transient_coeffs_vector;
4802 transient_coeffs_vector.push_back( __c0_pr );
4803 transient_coeffs_vector.push_back( __lambda_pr );
4804 transient_coeffs_vector.push_back( __gamma_pr );
4805 transient_coeffs_vector.push_back( __Cmf_pr );
4806 transient_coeffs_vector.push_back( __Cma_pr );
4807 transient_coeffs_vector.push_back( __Cmm_pr );
4809 return boost::make_tuple( delta_pr , transient_coeffs_vector ) ;
4814 template<
typename TruthModelType>
4815 typename CRB<TruthModelType>::residual_error_type
4819 beta_vector_type betaAqm;
4820 std::vector<beta_vector_type> betaFqm, betaLqm;
4821 boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu , time );
4823 int __QLhs = M_model->Qa();
4824 int __QRhs = M_model->Ql( 0 );
4828 value_type __c0_pr = 0.0;
4830 for (
int __q1 = 0; __q1 < __QRhs; ++__q1 )
4832 for (
int __m1 = 0; __m1 < M_model->mMaxF(0,__q1); ++__m1 )
4834 value_type fq1 = betaFqm[0][__q1][__m1];
4835 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
4837 for (
int __m2 = 0; __m2 < M_model->mMaxF(0,__q2); ++__m2 )
4839 value_type fq2 = betaFqm[0][__q2][__m2];
4840 __c0_pr += M_C0_pr[__q1][__m1][__q2][__m2]*fq1*fq2;
4846 value_type __lambda_pr = 0.0;
4847 value_type __gamma_pr = 0.0;
4850 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
4852 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1); ++__m1 )
4854 value_type a_q1 = betaAqm[__q1][__m1];
4856 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
4858 for (
int __m2 = 0; __m2 < M_model->mMaxF(0,__q2); ++__m2 )
4860 value_type f_q2 = betaFqm[0][__q2][__m2];
4861 __lambda_pr += a_q1*f_q2*M_Lambda_pr[__q1][__m1][__q2][__m2].head( __N ).dot( Un );
4865 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
4867 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
4869 value_type a_q2 = betaAqm[__q2][__m2];
4870 auto m = M_Gamma_pr[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Un;
4871 __gamma_pr += a_q1 * a_q2 * Un.dot( m );
4878 value_type delta_pr = math::abs( __c0_pr+__lambda_pr+__gamma_pr ) ;
4881 if ( !boost::math::isfinite( delta_pr ) )
4883 std::cout <<
"delta_pr is not finite:\n"
4884 <<
" - c0_pr = " << __c0_pr <<
"\n"
4885 <<
" - lambda_pr = " << __lambda_pr <<
"\n"
4886 <<
" - gamma_pr = " << __gamma_pr <<
"\n";
4887 std::cout <<
" - betaAqm = " << betaAqm <<
"\n"
4888 <<
" - betaFqm = " << betaFqm <<
"\n";
4889 std::cout <<
" - Un : " << Un <<
"\n";
4895 std::vector<double> coeffs_vector;
4896 coeffs_vector.push_back( __c0_pr );
4897 coeffs_vector.push_back( __lambda_pr );
4898 coeffs_vector.push_back( __gamma_pr );
4900 return boost::make_tuple( delta_pr,coeffs_vector );
4904 template<
typename TruthModelType>
4905 typename CRB<TruthModelType>::residual_error_type
4906 CRB<TruthModelType>::transientDualResidual(
int Ncur,parameter_type
const& mu, vectorN_type
const& Undu ,vectorN_type
const& Unduold ,
double time_step,
double time )
const
4908 int __QLhs = M_model->Qa();
4909 int __QOutput = M_model->Ql( M_output_index );
4910 int __Qm = M_model->Qm();
4913 beta_vector_type betaAqm;
4914 beta_vector_type betaMqm;
4915 beta_vector_type betaMFqm;
4916 std::vector<beta_vector_type> betaFqm;
4917 boost::tie( betaMqm, betaAqm, betaFqm) = M_model->computeBetaQm( mu, time );
4919 residual_error_type steady_residual_contribution = steadyDualResidual( Ncur, mu, Undu, time );
4920 std::vector<double> steady_coeff_vector = steady_residual_contribution.template get<1>();
4921 value_type __c0_du = steady_coeff_vector[0];
4922 value_type __lambda_du = steady_coeff_vector[1];
4923 value_type __gamma_du = steady_coeff_vector[2];
4926 value_type __Cmf_du=0;
4927 value_type __Cma_du=0;
4928 value_type __Cmm_du=0;
4932 if ( ! M_model->isSteady() )
4935 for (
int __q1=0 ; __q1<__Qm; ++__q1 )
4937 for (
int __m1=0 ; __m1< M_model->mMaxM(__q1); ++__m1 )
4939 value_type m_q1 = betaMqm[__q1][__m1];
4941 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
4943 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
4945 value_type a_q2 = betaAqm[__q2][__m2];
4946 auto m = M_Cma_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Undu;
4947 __Cma_du += 1./time_step * m_q1 * a_q2 * Undu.dot( m );
4948 __Cma_du -= 1./time_step * m_q1 * a_q2 * Unduold.dot( m );
4952 for (
int __q2 = 0; __q2 < __Qm; ++__q2 )
4954 for (
int __m2 = 0; __m2 < M_model->mMaxM(__q2); ++__m2 )
4956 value_type m_q2 = betaMqm[__q2][__m2];
4957 auto m1 = M_Cmm_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Undu;
4958 auto m2 = M_Cmm_du[__q1][__m1][__q2][__m2].block( 0,0,__N,__N )*Unduold;
4959 __Cmm_du += 1./( time_step*time_step ) * m_q1 * m_q2 * Undu.dot( m1 );
4960 __Cmm_du -= 1./( time_step*time_step ) * m_q1 * m_q2 * Undu.dot( m2 );
4961 __Cmm_du -= 1./( time_step*time_step ) * m_q1 * m_q2 * Unduold.dot( m1 );
4962 __Cmm_du += 1./( time_step*time_step ) * m_q1 * m_q2 * Unduold.dot( m2 );
4973 value_type delta_du;
4975 if ( M_model->isSteady() )
4976 delta_du = math::abs( __c0_du+__lambda_du+__gamma_du ) ;
4978 delta_du = math::abs( __gamma_du+__Cma_du+__Cmm_du ) ;
4983 std::cout<<
"[dualN2Q2] time "<<time<<std::endl;
4984 std::cout<<
"Undu = \n"<<Undu<<std::endl;
4985 std::cout<<
"Unduold = \n"<<Unduold<<std::endl;
4986 std::cout<<
"__c0_du = "<<__c0_du<<std::endl;
4987 std::cout<<
"__lambda_du = "<<__lambda_du<<std::endl;
4988 std::cout<<
"__gamma_du = "<<__gamma_du<<std::endl;
4989 std::cout<<
"__Cma_du = "<<__Cma_du<<std::endl;
4990 std::cout<<
"__Cmm_du = "<<__Cmm_du<<std::endl;
4994 std::vector<double> coeffs_vector;
4995 coeffs_vector.push_back( __c0_du );
4996 coeffs_vector.push_back( __lambda_du );
4997 coeffs_vector.push_back( __gamma_du );
4998 coeffs_vector.push_back( __Cmf_du );
4999 coeffs_vector.push_back( __Cma_du );
5000 coeffs_vector.push_back( __Cmm_du );
5002 return boost::make_tuple( delta_du,coeffs_vector );
5008 template<
typename TruthModelType>
5009 typename CRB<TruthModelType>::residual_error_type
5010 CRB<TruthModelType>::steadyDualResidual(
int Ncur,parameter_type
const& mu, vectorN_type
const& Undu,
double time )
const
5013 int __QLhs = M_model->Qa();
5014 int __QOutput = M_model->Ql( M_output_index );
5017 beta_vector_type betaAqm;
5018 beta_vector_type betaMqm;
5019 std::vector<beta_vector_type> betaFqm, betaLqm;
5020 boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu , time );
5022 value_type __c0_du = 0.0;
5024 for (
int __q1 = 0; __q1 < __QOutput; ++__q1 )
5026 for (
int __m1 = 0; __m1 < M_model->mMaxF(M_output_index,__q1); ++__m1 )
5028 value_type fq1 = betaFqm[M_output_index][__q1][__m1];
5029 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
5031 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
5033 value_type fq2 = betaFqm[M_output_index][__q2][__m2];
5034 __c0_du += M_C0_du[__q1][__m1][__q2][__m2]*fq1*fq2;
5040 value_type __lambda_du = 0.0;
5041 value_type __gamma_du = 0.0;
5043 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
5045 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1); ++__m1 )
5047 value_type a_q1 = betaAqm[__q1][__m1];
5048 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
5050 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
5052 value_type a_q2 = betaFqm[M_output_index][__q2][__m2]*a_q1;
5053 __lambda_du += a_q2 * M_Lambda_du[__q1][__m1][__q2][__m2].head( __N ).dot( Undu );
5057 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
5059 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
5061 value_type a_q2 = betaAqm[__q2][__m2]*a_q1;
5062 auto m = M_Gamma_du[ __q1][ __m1][ __q2][ __m2].block( 0,0,__N,__N )*Undu;
5063 __gamma_du += a_q2*Undu.dot( m );
5069 value_type delta_du = math::abs( __c0_du+__lambda_du+__gamma_du );
5073 if ( !boost::math::isfinite( delta_du ) )
5075 std::cout <<
"delta_du is not finite:\n"
5076 <<
" - c0_du = " << __c0_du <<
"\n"
5077 <<
" - lambda_du = " << __lambda_du <<
"\n"
5078 <<
" - gamma_du = " << __gamma_du <<
"\n";
5079 std::cout <<
" - betaAqm = " << betaAqm <<
"\n"
5080 <<
" - betaFqm = " << betaFqm <<
"\n";
5081 std::cout <<
" - Undu : " << Undu <<
"\n";
5086 std::vector<double> coeffs_vector;
5087 coeffs_vector.push_back( __c0_du );
5088 coeffs_vector.push_back( __lambda_du );
5089 coeffs_vector.push_back( __gamma_du );
5090 return boost::make_tuple( delta_du,coeffs_vector );
5093 template<
typename TruthModelType>
5097 return offlineResidual( Ncur, mpl::bool_<model_type::is_time_dependent>(), number_of_added_elements );
5100 template<
typename TruthModelType>
5106 int __QLhs = M_model->Qa();
5107 int __QRhs = M_model->Ql( 0 );
5108 int __QOutput = M_model->Ql( M_output_index );
5109 int __Qm = M_model->Qm();
5112 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5113 std::cout <<
" o N=" << Ncur <<
" QLhs=" << __QLhs
5114 <<
" QRhs=" << __QRhs <<
" Qoutput=" << __QOutput
5115 <<
" Qm=" << __Qm <<
"\n";
5117 vector_ptrtype __X( M_backend->newVector( M_model->functionSpace() ) );
5118 vector_ptrtype __Fdu( M_backend->newVector( M_model->functionSpace() ) );
5119 vector_ptrtype __Z1( M_backend->newVector( M_model->functionSpace() ) );
5120 vector_ptrtype __Z2( M_backend->newVector( M_model->functionSpace() ) );
5121 vector_ptrtype __W( M_backend->newVector( M_model->functionSpace() ) );
5122 namespace ublas = boost::numeric::ublas;
5124 std::vector< std::vector<sparse_matrix_ptrtype> > Aqm;
5125 std::vector< std::vector<sparse_matrix_ptrtype> > Mqm;
5126 std::vector< std::vector<std::vector<vector_ptrtype> > > Fqm;
5127 std::vector<std::vector<vector_ptrtype> > MFqm;
5128 boost::tie( Mqm, Aqm, Fqm ) = M_model->computeAffineDecomposition();
5133 ublas::vector<value_type> mu( P );
5135 for (
int q = 0; q < P; ++q )
5142 offlineResidual( Ncur, mpl::bool_<false>(), number_of_added_elements );
5145 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5147 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5150 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
5152 for (
int __m2 = 0; __m2 < M_model->mMaxF(0,__q2); ++__m2 )
5155 M_Cmf_pr[__q1][__m1][__q2][__m2].conservativeResize( __N );
5157 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5160 Mqm[__q1][__m1]->multVector( __X, __W );
5162 M_model->l2solve( __Z1, __W );
5164 M_model->l2solve( __Z2, Fqm[0][__q2][__m2] );
5166 M_Cmf_pr[ __q1][ __m1][ __q2][ __m2]( elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5173 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5174 std::cout <<
" o M_Cmf_pr updated in " << ti.elapsed() <<
"s\n";
5178 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5180 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5183 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
5185 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
5188 M_Cma_pr[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5190 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5193 Mqm[__q1][__m1]->multVector( __X, __W );
5195 M_model->l2solve( __Z1, __W );
5197 for (
int __l = 0; __l < ( int )__N; ++__l )
5200 Aqm[__q2][__m2]->multVector( __X, __W );
5202 M_model->l2solve( __Z2, __W );
5204 M_Cma_pr[ __q1][ __m1][ __q2][ __m2]( elem,__l ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5208 for (
int __j = 0; __j < ( int )__N; ++__j )
5211 Mqm[__q1][ __m1]->multVector( __X, __W );
5213 M_model->l2solve( __Z1, __W );
5215 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5218 Aqm[__q2][ __m2]->multVector( __X, __W );
5220 M_model->l2solve( __Z2, __W );
5222 M_Cma_pr[ __q1][ __m1][ __q2][ __m2]( __j,elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5231 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5232 std::cout <<
" o M_Cma_pr updated in " << ti.elapsed() <<
"s\n";
5235 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5237 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5240 for (
int __q2 = 0; __q2 < __Qm; ++__q2 )
5243 for (
int __m2 = 0; __m2 < M_model->mMaxM(__q2); ++__m2 )
5246 M_Cmm_pr[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5248 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5251 Mqm[__q1][__m1]->multVector( __X, __W );
5253 M_model->l2solve( __Z1, __W );
5255 for (
int __l = 0; __l < ( int )__N; ++__l )
5258 Mqm[__q2][__m2]->multVector( __X, __W );
5260 M_model->l2solve( __Z2, __W );
5262 M_Cmm_pr[ __q1][ __m1][ __q2][ __m2]( elem,__l ) = M_model->scalarProduct( __Z1, __Z2 );
5266 for (
int __j = 0; __j < ( int )__N; ++__j )
5269 Mqm[__q1][ __m1]->multVector( __X, __W );
5271 M_model->l2solve( __Z1, __W );
5273 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5276 Mqm[__q2][ __m2]->multVector( __X, __W );
5278 M_model->l2solve( __Z2, __W );
5280 M_Cmm_pr[ __q1][ __m1][ __q2][ __m2]( __j,elem ) = M_model->scalarProduct( __Z1, __Z2 );
5290 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5291 std::cout <<
" o M_Cmm_pr updated in " << ti.elapsed() <<
"s\n";
5295 sparse_matrix_ptrtype Atq1 = M_model->newMatrix();
5296 sparse_matrix_ptrtype Atq2 = M_model->newMatrix();
5297 sparse_matrix_ptrtype Mtq1 = M_model->newMatrix();
5298 sparse_matrix_ptrtype Mtq2 = M_model->newMatrix();
5304 LOG(INFO) <<
"[offlineResidual] Cmf_du Cma_du Cmm_du\n";
5306 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5309 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5312 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5313 Mtq1 = Mqm[__q1][__m1];
5315 Mqm[__q1][__m1]->transpose( Mtq1 );
5317 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
5319 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
5322 M_Cmf_du[__q1][__m1][__q2][__m2].conservativeResize( __N );
5323 M_Cmf_du_ini[__q1][__m1][__q2][__m2].conservativeResize( __N );
5325 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5328 Mtq1->multVector( __X, __W );
5331 M_model->l2solve( __Z1, __W );
5333 *__Fdu = *Fqm[M_output_index][__q2][__m2];
5335 __Fdu->scale( -1.0 );
5336 M_model->l2solve( __Z2, __Fdu );
5338 M_Cmf_du[ __q1][__m1][ __q2][__m2]( elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5340 *__Fdu = *Fqm[M_output_index][__q2][__m2];
5342 M_model->l2solve( __Z2, __Fdu );
5343 M_Cmf_du_ini[ __q1][__m1][ __q2][__m2]( elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5351 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5353 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5356 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5357 Mtq1 = Mqm[__q1][__m1];
5359 Mqm[__q1][__m1]->transpose( Mtq1 );
5361 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
5363 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
5366 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5367 Atq2 = Aqm[__q2][__m2];
5369 Aqm[__q2][__m2]->transpose( Atq2 );
5371 M_Cma_du[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5373 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5377 Mtq1->multVector( __X, __W );
5379 M_model->l2solve( __Z1, __W );
5381 for (
int __l = 0; __l < ( int )__N; ++__l )
5384 Atq2->multVector( __X, __W );
5386 M_model->l2solve( __Z2, __W );
5388 M_Cma_du[ __q1][__m1][__q2][ __m2]( elem,__l ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5392 for (
int __j = 0; __j < ( int )__N; ++__j )
5395 Mtq1->multVector( __X, __W );
5397 M_model->l2solve( __Z1, __W );
5399 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5402 Atq2->multVector( __X, __W );
5404 M_model->l2solve( __Z2, __W );
5406 M_Cma_du[ __q1][__m1][__q2][ __m2]( __j,elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5415 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5416 std::cout <<
" o Cma_du updated in " << ti.elapsed() <<
"s\n";
5422 for (
int __q1 = 0; __q1 < __Qm; ++__q1 )
5424 for (
int __m1 = 0; __m1 < M_model->mMaxM(__q1); ++__m1 )
5426 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5427 Mtq1 = Mqm[__q1][__m1];
5429 Mqm[__q1][__m1]->transpose( Mtq1 );
5431 for (
int __q2 = 0; __q2 < __Qm; ++__q2 )
5433 for (
int __m2 = 0; __m2 < M_model->mMaxM(__q2); ++__m2 )
5436 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5437 Mtq2 = Mqm[__q2][__m2];
5439 Mqm[__q2][__m2]->transpose( Mtq2 );
5441 M_Cmm_du[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5443 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5446 Mtq1->multVector( __X, __W );
5448 M_model->l2solve( __Z1, __W );
5450 for (
int __l = 0; __l < ( int )__N; ++__l )
5453 Mtq2->multVector( __X, __W );
5455 M_model->l2solve( __Z2, __W );
5457 M_Cmm_du[ __q1][ __m1][ __q2][ __m2]( elem,__l ) = M_model->scalarProduct( __Z1, __Z2 );
5461 for (
int __j = 0; __j < ( int )__N; ++__j )
5464 Mtq1->multVector( __X, __W );
5466 M_model->l2solve( __Z1, __W );
5468 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5471 Mtq2->multVector( __X, __W );
5473 M_model->l2solve( __Z2, __W );
5475 M_Cmm_du[ __q1][ __m1][ __q2][ __m2]( __j,elem ) = M_model->scalarProduct( __Z1, __Z2 );
5484 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5485 std::cout <<
" o Cmm_du updated in " << ti.elapsed() <<
"s\n";
5488 LOG(INFO) <<
"[offlineResidual] Done.\n";
5494 template<
typename TruthModelType>
5499 int __QLhs = M_model->Qa();
5500 int __QRhs = M_model->Ql( 0 );
5501 int __QOutput = M_model->Ql( M_output_index );
5504 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5505 std::cout <<
" o N=" << Ncur <<
" QLhs=" << __QLhs
5506 <<
" QRhs=" << __QRhs <<
" Qoutput=" << __QOutput <<
"\n";
5508 vector_ptrtype __X( M_backend->newVector( M_model->functionSpace() ) );
5509 vector_ptrtype __Fdu( M_backend->newVector( M_model->functionSpace() ) );
5510 vector_ptrtype __Z1( M_backend->newVector( M_model->functionSpace() ) );
5511 vector_ptrtype __Z2( M_backend->newVector( M_model->functionSpace() ) );
5512 vector_ptrtype __W( M_backend->newVector( M_model->functionSpace() ) );
5513 namespace ublas = boost::numeric::ublas;
5515 std::vector< std::vector<sparse_matrix_ptrtype> > Aqm,Mqm;
5516 std::vector< std::vector<vector_ptrtype> > MFqm;
5517 std::vector< std::vector<std::vector<vector_ptrtype> > > Fqm,Lqm;
5518 boost::tie( Mqm, Aqm, Fqm ) = M_model->computeAffineDecomposition();
5523 ublas::vector<value_type> mu( P );
5525 for (
int q = 0; q < P; ++q )
5536 LOG(INFO) <<
"[offlineResidual] Compute Primal residual data\n";
5537 LOG(INFO) <<
"[offlineResidual] C0_pr\n";
5540 for (
int __q1 = 0; __q1 < __QRhs; ++__q1 )
5542 for (
int __m1 = 0; __m1 < M_model->mMaxF(0,__q1); ++__m1 )
5545 M_model->l2solve( __Z1, Fqm[0][__q1][__m1] );
5547 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
5549 for (
int __m2 = 0; __m2 < M_model->mMaxF(0,__q2); ++__m2 )
5552 M_model->l2solve( __Z2, Fqm[0][__q2][__m2] );
5554 M_C0_pr[__q1][__m1][__q2][__m2] = M_model->scalarProduct( __Z1, __Z2 );
5565 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5566 std::cout <<
" o initialize offlineResidual in " << ti.elapsed() <<
"s\n";
5570 parameter_type
const& mu = M_WNmu->at( 0 );
5572 beta_vector_type betaAqm;
5573 std::vector<beta_vector_type> betaFqm;
5574 boost::tie( boost::tuples::ignore, betaAqm, betaFqm ) = M_model->computeBetaQm( mu );
5575 value_type __c0_pr = 0.0;
5577 for (
int __q1 = 0; __q1 < __QRhs; ++__q1 )
5579 for (
int __m1 = 0; __m1 < M_model->mMaxF(0,__q1); ++__m1 )
5581 value_type fq1 = betaFqm[0][__q1][__m1];
5583 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
5585 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
5587 value_type fq2 = betaFqm[0][__q2][__m2];
5588 __c0_pr += M_C0_pr[__q1][__m1][__q2][__m2]*fq1*fq2;
5597 std::vector< sparse_matrix_ptrtype > A;
5598 std::vector< std::vector<vector_ptrtype> > F;
5599 boost::tie( A, F ) = M_model->update( mu );
5600 M_model->l2solve( __X, F[0] );
5607 LOG(INFO) <<
"[offlineResidual] Lambda_pr, Gamma_pr\n";
5609 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5613 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
5615 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1) ; ++__m1 )
5617 Aqm[__q1][__m1]->multVector( __X, __W );
5620 M_model->l2solve( __Z1, __W );
5621 for (
int __q2 = 0; __q2 < __QRhs; ++__q2 )
5623 for (
int __m2 = 0; __m2 < M_model->mMaxF(0, __q2); ++__m2 )
5625 M_Lambda_pr[__q1][__m1][__q2][__m2].conservativeResize( __N );
5628 M_model->l2solve( __Z2, Fqm[0][__q2][__m2] );
5629 M_Lambda_pr[ __q1][ __m1][ __q2][ __m2]( elem ) = 2.0*M_model->scalarProduct( __Z1, __Z2 );
5638 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5639 std::cout <<
" o Lambda_pr updated in " << ti.elapsed() <<
"s\n";
5643 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
5645 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1); ++__m1 )
5648 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
5651 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
5654 M_Gamma_pr[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5656 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5659 Aqm[__q1][__m1]->multVector( __X, __W );
5662 M_model->l2solve( __Z1, __W );
5665 for (
int __l = 0; __l < ( int )__N; ++__l )
5668 Aqm[__q2][__m2]->multVector( __X, __W );
5671 M_model->l2solve( __Z2, __W );
5672 M_Gamma_pr[ __q1][ __m1][ __q2][ __m2]( elem,__l ) = M_model->scalarProduct( __Z1, __Z2 );
5679 for (
int __j = 0; __j < ( int )__N; ++__j )
5682 Aqm[__q1][__m1]->multVector( __X, __W );
5685 M_model->l2solve( __Z1, __W );
5690 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5693 Aqm[__q2][__m2]->multVector( __X, __W );
5696 M_model->l2solve( __Z2, __W );
5697 M_Gamma_pr[ __q1][ __m1][ __q2][ __m2]( __j,elem ) = M_model->scalarProduct( __Z1, __Z2 );
5708 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5709 std::cout <<
" o Gamma_pr updated in " << ti.elapsed() <<
"s\n";
5711 sparse_matrix_ptrtype Atq1 = M_model->newMatrix();
5712 sparse_matrix_ptrtype Atq2 = M_model->newMatrix();
5718 if( solve_dual_problem )
5722 LOG(INFO) <<
"[offlineResidual] Compute Dual residual data\n";
5723 LOG(INFO) <<
"[offlineResidual] C0_du\n";
5725 for (
int __q1 = 0; __q1 < __QOutput; ++__q1 )
5727 for (
int __m1 = 0; __m1 < M_model->mMaxF(M_output_index,__q1); ++__m1 )
5729 *__Fdu = *Fqm[M_output_index][__q1][__m1];
5731 __Fdu->scale( -1.0 );
5732 M_model->l2solve( __Z1, __Fdu );
5734 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
5736 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2); ++__m2 )
5738 *__Fdu = *Fqm[M_output_index][__q2][__m2];
5740 __Fdu->scale( -1.0 );
5741 M_model->l2solve( __Z2, __Fdu );
5742 M_C0_du[__q1][__m1][__q2][__m2] = M_model->scalarProduct( __Z1, __Z2 );
5750 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5751 std::cout <<
" o C0_du updated in " << ti.elapsed() <<
"s\n";
5755 LOG(INFO) <<
"[offlineResidual] Lambda_du, Gamma_du\n";
5757 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5761 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
5763 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1); ++__m1 )
5766 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5767 Atq1 = Aqm[__q1][__m1];
5769 Aqm[__q1][__m1]->transpose( Atq1 );
5771 Atq1->multVector( __X, __W );
5775 M_model->l2solve( __Z1, __W );
5777 for (
int __q2 = 0; __q2 < __QOutput; ++__q2 )
5779 for (
int __m2 = 0; __m2 < M_model->mMaxF(M_output_index,__q2) ; ++__m2 )
5781 M_Lambda_du[__q1][__m1][__q2][__m2].conservativeResize( __N );
5783 *__Fdu = *Fqm[M_output_index][__q2][__m2];
5785 __Fdu->scale( -1.0 );
5786 M_model->l2solve( __Z2, __Fdu );
5787 M_Lambda_du[ __q1][ __m1][ __q2][ __m2]( elem ) = 2.0*M_model->scalarProduct( __Z2, __Z1 );
5796 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5797 std::cout <<
" o Lambda_du updated in " << ti.elapsed() <<
"s\n";
5800 for (
int __q1 = 0; __q1 < __QLhs; ++__q1 )
5802 for (
int __m1 = 0; __m1 < M_model->mMaxA(__q1); ++__m1 )
5805 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5806 Atq1=Aqm[__q1][__m1];
5808 Aqm[__q1][__m1]->transpose( Atq1 );
5810 for (
int __q2 = 0; __q2 < __QLhs; ++__q2 )
5812 for (
int __m2 = 0; __m2 < M_model->mMaxA(__q2); ++__m2 )
5815 if( option(
"crb.use-symmetric-matrix").template as<bool>() )
5816 Atq2 = Aqm[__q2][__m2];
5818 Aqm[__q2][__m2]->transpose( Atq2 );
5820 M_Gamma_du[__q1][__m1][__q2][__m2].conservativeResize( __N, __N );
5822 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5826 Atq1->multVector( __X, __W );
5829 M_model->l2solve( __Z1, __W );
5833 for (
int __l = 0; __l < ( int )__N; ++__l )
5836 Atq2->multVector( __X, __W );
5839 M_model->l2solve( __Z2, __W );
5840 M_Gamma_du[ __q1][ __m1][ __q2][ __m2]( elem,__l ) = M_model->scalarProduct( __Z1, __Z2 );
5848 for (
int __j = 0; __j < ( int )__N; ++__j )
5851 Atq1->multVector( __X, __W );
5854 M_model->l2solve( __Z1, __W );
5858 for (
int elem=__N-number_of_added_elements; elem<__N; elem++ )
5861 Atq2->multVector( __X, __W );
5864 M_model->l2solve( __Z2, __W );
5865 M_Gamma_du[ __q1][ __m1][ __q2][ __m2]( __j,elem ) = M_model->scalarProduct( __Z1, __Z2 );
5876 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
5877 std::cout <<
" o Gamma_du updated in " << ti.elapsed() <<
"s\n";
5879 LOG(INFO) <<
"[offlineResidual] Done.\n";
5884 template<
typename TruthModelType>
5889 if( M_rbconv_contains_primal_and_dual_contributions )
5891 typedef convergence_type::left_map::const_iterator iterator;
5892 LOG(INFO)<<
"\nMax error during offline stage\n";
5893 for(iterator it = M_rbconv.left.begin(); it != M_rbconv.left.end(); ++it)
5894 LOG(INFO)<<
"N : "<<it->first<<
" - maxerror : "<<it->second.template get<0>()<<
"\n";
5896 LOG(INFO)<<
"\nPrimal contribution\n";
5897 for(iterator it = M_rbconv.left.begin(); it != M_rbconv.left.end(); ++it)
5898 LOG(INFO)<<
"N : "<<it->first<<
" - delta_pr : "<<it->second.template get<1>()<<
"\n";
5900 LOG(INFO)<<
"\nDual contribution\n";
5901 for(iterator it = M_rbconv.left.begin(); it != M_rbconv.left.end(); ++it)
5902 LOG(INFO)<<
"N : "<<it->first<<
" - delta_du : "<<it->second.template get<2>()<<
"\n";
5906 throw std::logic_error(
"[CRB::printErrorsDuringRbConstruction] ERROR, the database is too old to print the error during offline step, use the option rebuild-database = true" );
5912 template<
typename TruthModelType>
5916 LOG(INFO)<<
" List of parameter selectionned during the offline algorithm \n";
5917 for(
int k=0;k<M_WNmu->size();k++)
5919 std::cout<<
" mu "<<k<<
" = [ ";
5920 LOG(INFO)<<
" mu "<<k<<
" = [ ";
5922 for(
int i=0; i<_mu.size()-1; i++ )
5924 LOG(INFO)<<_mu(i)<<
" , ";
5925 std::cout<<_mu(i)<<
" , ";
5927 LOG(INFO)<<_mu( _mu.size()-1 )<<
" ] \n";
5928 std::cout<<_mu( _mu.size()-1 )<<
" ] "<<std::endl;
5933 template<
typename TruthModelType>
5944 std::vector<vectorN_type> uN;
5945 std::vector<vectorN_type> uNdu;
5946 std::vector<vectorN_type> uNold;
5947 std::vector<vectorN_type> uNduold;
5949 auto o =
lb( Nwn, mu, uN, uNdu , uNold, uNduold );
5950 int size = uN.size();
5951 FEELPP_ASSERT( N <= M_WN.size() )( N )( M_WN.size() ).error(
"invalid expansion size ( N and M_WN ) ");
5953 if( time_index == -1 )
5954 ucrb = Feel::expansion( M_WN, uN[size-1] , Nwn);
5957 CHECK( time_index < size )<<
" call crb::expansion with a wrong value of time index : "<<time_index<<
" or size of uN vector is only "<<size;
5958 ucrb = Feel::expansion( M_WN, uN[time_index] , Nwn);
5964 template<
typename TruthModelType>
5976 FEELPP_ASSERT( Nwn <= WN.size() )( Nwn )( WN.size() ).error(
"invalid expansion size ( Nwn and WN ) ");
5977 FEELPP_ASSERT( Nwn <= u.size() )( Nwn )( WN.size() ).error(
"invalid expansion size ( Nwn and u ) ");
5979 return Feel::expansion( WN, u, N );
5983 template<
typename TruthModelType>
5984 boost::tuple<double,double, typename CRB<TruthModelType>::solutions_tuple, double, double, double,
typename CRB<TruthModelType>::upper_bounds_tuple >
5988 M_compute_variance = this->
vm()[
"crb.compute-variance"].template as<bool>();
5991 int Nwn_max =
vm()[
"crb.dimension-max"].template as<int>();
5993 if ( M_error_type!=CRB_EMPIRICAL )
5995 auto lo = M_rbconv.right.range( boost::bimaps::unbounded, boost::bimaps::_key <= eps );
5997 for (
auto it = lo.first; it != lo.second; ++it )
5999 std::cout <<
"rbconv[" << it->first <<
"]=" << it->second <<
"\n";
6002 auto it = M_rbconv.project_left( lo.first );
6004 std::cout <<
"Nwn = "<< Nwn <<
" error = "<< it->second.template get<0>() <<
" eps=" << eps <<
"\n";
6007 if ( this->
vm()[
"crb.check.residual"].template as<bool>() )
6009 std::vector< std::vector<double> > primal_residual_coefficients = error_estimation.template get<1>();
6010 std::vector< std::vector<double> > dual_residual_coefficients = error_estimation.template get<2>();
6011 std::vector<element_ptrtype> Un,Unold,Undu,Unduold;
6012 buildFunctionFromRbCoefficients(Nwn, uN, M_WN, Un );
6013 buildFunctionFromRbCoefficients(Nwn, uNold, M_WN, Unold );
6014 buildFunctionFromRbCoefficients(Nwn, uNdu, M_WNdu, Undu );
6015 buildFunctionFromRbCoefficients(Nwn, uNduold, M_WNdu, Unduold );
6016 compareResidualsForTransientProblems(Nwn, mu , Un, Unold, Undu, Unduold, primal_residual_coefficients, dual_residual_coefficients );
6020 std::vector<vectorN_type> uN;
6021 std::vector<vectorN_type> uNdu;
6022 std::vector<vectorN_type> uNold;
6023 std::vector<vectorN_type> uNduold;
6037 auto o =
lb( Nwn, mu, uN, uNdu , uNold, uNduold );
6038 double output = o.template get<0>();
6040 auto error_estimation =
delta( Nwn, mu, uN, uNdu , uNold, uNduold );
6041 double output_upper_bound = error_estimation.template get<0>();
6042 double condition_number = o.template get<1>();
6043 auto primal_coefficients = error_estimation.template get<1>();
6044 auto dual_coefficients = error_estimation.template get<2>();
6047 int nb_dt = primal_coefficients.size();
6048 int final_time_index = nb_dt-1;
6049 double primal_residual_norm = 0;
6050 double dual_residual_norm = 0;
6052 if ( M_error_type == CRB_RESIDUAL || M_error_type == CRB_RESIDUAL_SCM )
6054 int nb_coeff = primal_coefficients[final_time_index].size();
6055 for(
int i=0 ; i<nb_coeff ; i++)
6056 primal_residual_norm += primal_coefficients[final_time_index][i] ;
6058 bool solve_dual_problem = this->
vm()[
"crb.solve-dual-problem"].template as<bool>() ;
6059 if( solve_dual_problem )
6061 if ( M_model->isSteady() )
6062 dual_residual_norm = math::abs( dual_coefficients[0][0]+dual_coefficients[0][1]+dual_coefficients[0][2] ) ;
6064 dual_residual_norm = math::abs( dual_coefficients[0][2]+dual_coefficients[0][4]+dual_coefficients[0][5] ) ;
6066 primal_residual_norm = math::sqrt( math::abs(primal_residual_norm) );
6067 dual_residual_norm = math::sqrt( math::abs(dual_residual_norm) );
6069 double delta_pr = error_estimation.template get<3>();
6070 double delta_du = error_estimation.template get<4>();
6071 int size = uN.size();
6073 auto upper_bounds = boost::make_tuple(output_upper_bound , delta_pr, delta_du , primal_coefficients , dual_coefficients );
6074 auto solutions = boost::make_tuple( uN , uNdu, uNold, uNduold);
6075 return boost::make_tuple( output , Nwn , solutions, condition_number , primal_residual_norm , dual_residual_norm, upper_bounds );
6079 template<
typename TruthModelType>
6083 return run( X, N, Y, P, mpl::bool_<model_type::is_time_dependent>() );
6088 template<
typename TruthModelType>
6094 for (
unsigned long p= 0; p < N-5; ++p ) std::cout<<
"mu["<<p<<
"] = "<<X[p]<<std::endl;
6097 parameter_type mu( M_Dmu );
6100 for (
unsigned long p= 0; p < N-5; ++p )
6119 int maxerror = X[N-3];
6124 M_compute_variance = X[N-1];
6130 if ( M_error_type!=CRB_EMPIRICAL )
6132 auto lo = M_rbconv.right.range( boost::bimaps::unbounded,boost::bimaps::_key <= maxerror );
6134 for (
auto it = lo.first; it != lo.second; ++it )
6136 std::cout <<
"rbconv[" << it->first <<
"]=" << it->second <<
"\n";
6139 auto it = M_rbconv.project_left( lo.first );
6145 std::vector<vectorN_type> uN;
6146 std::vector<vectorN_type> uNdu;
6147 std::vector<vectorN_type> uNold ;
6148 std::vector<vectorN_type> uNduold;
6150 FEELPP_ASSERT( P == 2 )( P ).warn(
"invalid number of outputs" );
6151 auto o =
lb( Nwn, mu, uN, uNdu , uNold, uNduold );
6152 auto e =
delta( Nwn, mu, uN, uNdu , uNold, uNduold );
6153 Y[0] = o.template get<0>();
6154 Y[1] = e.template get<0>();
6158 template<
typename TruthModelType>
6163 parameter_type mu( M_Dmu );
6165 for (
unsigned long p= 0; p < N-5; ++p )
6168 std::cout <<
"mu( " << p <<
" ) = " << mu( p ) << std::endl;
6180 int maxerror = X[N-3];
6185 M_compute_variance = X[N-1];
6189 auto lo = M_rbconv.right.range( boost::bimaps::unbounded,boost::bimaps::_key <= maxerror );
6192 for (
auto it = lo.first; it != lo.second; ++it )
6194 std::cout <<
"rbconv[" << it->first <<
"]=" << it->second <<
"\n";
6197 auto it = M_rbconv.project_left( lo.first );
6202 std::vector<vectorN_type> uN;
6203 std::vector<vectorN_type> uNdu;
6204 std::vector<vectorN_type> uNold ;
6205 std::vector<vectorN_type> uNduold;
6207 FEELPP_ASSERT( P == 2 )( P ).warn(
"invalid number of outputs" );
6208 auto o =
lb( Nwn, mu, uN, uNdu , uNold, uNduold );
6209 auto e =
delta( Nwn, mu, uN, uNdu , uNold, uNduold );
6210 Y[0] = o.template get<0>();
6211 Y[1] = e.template get<0>();
6219 template<
typename TruthModelType>
6226 if ( name_of_space==
"dual" )
6228 if ( orthonormalize_dual )
6231 BOOST_FOREACH(
auto du, M_WNdu )
6235 double k = M_model->scalarProduct( u, e );
6237 projection->add( 1 , e );
6243 matrixN_type MN ( (
int )M_N, (
int )M_N ) ;
6244 vectorN_type FN ( (
int )M_N );
6250 MN( i,j ) = M_model->scalarProduct( M_WNdu[j] , M_WNdu[i] );
6251 MN( j,i ) = MN( i,j );
6254 MN( i,i ) = M_model->scalarProduct( M_WNdu[i] , M_WNdu[i] );
6255 FN( i ) = M_model->scalarProduct( u,M_WNdu[i] );
6258 vectorN_type projectionN ( (
int ) M_N );
6259 projectionN = MN.lu().solve( FN );
6261 BOOST_FOREACH(
auto du, M_WNdu )
6265 double k = projectionN( index );
6267 projection->add( 1 , e );
6275 if ( orthonormalize_primal )
6277 BOOST_FOREACH(
auto pr, M_WN )
6279 auto e = pr.functionSpace()->element();
6281 double k = M_model->scalarProduct( u, e );
6283 projection->add( 1 , e );
6289 matrixN_type MN ( (
int )M_N, (
int )M_N ) ;
6290 vectorN_type FN ( (
int )M_N );
6296 MN( i,j ) = M_model->scalarProduct( M_WN[j] , M_WN[i] );
6297 MN( j,i ) = MN( i,j );
6300 MN( i,i ) = M_model->scalarProduct( M_WN[i] , M_WN[i] );
6301 FN( i ) = M_model->scalarProduct( u,M_WN[i] );
6304 vectorN_type projectionN ( (
int ) M_N );
6305 projectionN = MN.lu().solve( FN );
6307 BOOST_FOREACH(
auto pr, M_WN )
6311 double k = projectionN( index );
6313 projection->add( 1 , e );
6323 template<
typename TruthModelType>
6328 double min=0,max=0,mean=0,standard_deviation=0;
6330 int n_eval = option(_name=
"crb.computational-time-neval").template as<int>();
6332 Eigen::Matrix<double, Eigen::Dynamic, 1> time_crb;
6333 time_crb.resize( n_eval );
6336 Sampling->logEquidistribute( n_eval );
6338 bool cvg = option(_name=
"crb.cvg-study").template as<bool>();
6340 double tol = option(_name=
"crb.online-tolerance").template as<double>();
6347 int proc_number = Environment::worldComm().globalRank();
6348 int master = Environment::worldComm().masterRank();
6351 std::string file_name =
"cvg-timing-crb.dat";
6354 if( proc_number == master )
6356 conv.open(file_name, std::ios::app);
6357 conv <<
"NbBasis" <<
"\t" <<
"min" <<
"\t"<<
"max" <<
"\t"<<
"mean"<<
"\t"<<
"standard_deviation" <<
"\n";
6365 BOOST_FOREACH(
auto mu, *Sampling )
6367 boost::mpi::timer tcrb;
6368 auto o = this->
run( mu, tol , N);
6369 time_crb( mu_number ) = tcrb.elapsed() ;
6373 auto stat = M_model->computeStatistics( time_crb , appname );
6378 standard_deviation=stat(3);
6380 if( proc_number == master )
6381 conv << N <<
"\t" << min <<
"\t" << max<<
"\t"<< mean<<
"\t"<< standard_deviation<<
"\n";
6387 template<
typename TruthModelType>
6388 template<
class Archive>
6392 int proc_number = this->worldComm().globalRank();
6394 LOG(INFO) <<
"[CRB::save] version : "<<version<<std::endl;
6396 auto mesh = mesh_type::New();
6397 auto is_mesh_loaded = mesh->load( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
6399 if ( ! is_mesh_loaded )
6401 auto first_element = M_WN[0];
6402 mesh = first_element.functionSpace()->mesh() ;
6403 mesh->save( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
6407 ar & boost::serialization::base_object<super>( *this );
6408 ar & BOOST_SERIALIZATION_NVP( M_output_index );
6409 ar & BOOST_SERIALIZATION_NVP( M_N );
6410 ar & BOOST_SERIALIZATION_NVP( M_rbconv );
6411 ar & BOOST_SERIALIZATION_NVP( M_error_type );
6412 ar & BOOST_SERIALIZATION_NVP( M_Xi );
6413 ar & BOOST_SERIALIZATION_NVP( M_WNmu );
6414 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr );
6415 ar & BOOST_SERIALIZATION_NVP( M_Aqm_du );
6416 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr_du );
6417 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
6418 ar & BOOST_SERIALIZATION_NVP( M_Fqm_du );
6419 ar & BOOST_SERIALIZATION_NVP( M_Lqm_pr );
6420 ar & BOOST_SERIALIZATION_NVP( M_Lqm_du );
6422 ar & BOOST_SERIALIZATION_NVP( M_C0_pr );
6423 ar & BOOST_SERIALIZATION_NVP( M_C0_du );
6424 ar & BOOST_SERIALIZATION_NVP( M_Lambda_pr );
6425 ar & BOOST_SERIALIZATION_NVP( M_Lambda_du );
6426 ar & BOOST_SERIALIZATION_NVP( M_Gamma_pr );
6427 ar & BOOST_SERIALIZATION_NVP( M_Gamma_du );
6430 ar & BOOST_SERIALIZATION_NVP( M_Mqm_pr );
6431 ar & BOOST_SERIALIZATION_NVP( M_Mqm_du );
6432 ar & BOOST_SERIALIZATION_NVP( M_Mqm_pr_du );
6434 if ( model_type::is_time_dependent )
6437 ar & BOOST_SERIALIZATION_NVP( M_coeff_pr_ini_online );
6438 ar & BOOST_SERIALIZATION_NVP( M_coeff_du_ini_online );
6439 ar & BOOST_SERIALIZATION_NVP( M_Cmf_pr );
6440 ar & BOOST_SERIALIZATION_NVP( M_Cma_pr );
6441 ar & BOOST_SERIALIZATION_NVP( M_Cmm_pr );
6442 ar & BOOST_SERIALIZATION_NVP( M_Cmf_du );
6443 ar & BOOST_SERIALIZATION_NVP( M_Cmf_du_ini );
6444 ar & BOOST_SERIALIZATION_NVP( M_Cma_du );
6445 ar & BOOST_SERIALIZATION_NVP( M_Cmm_du );
6447 ar & BOOST_SERIALIZATION_NVP ( M_database_contains_variance_info );
6448 if( M_database_contains_variance_info )
6449 ar & BOOST_SERIALIZATION_NVP( M_variance_matrix_phi );
6451 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
6452 ar & BOOST_SERIALIZATION_NVP( M_InitialGuessV_pr );
6454 ar & BOOST_SERIALIZATION_NVP( M_current_mu );
6455 ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
6458 for(
int i=0; i<M_N; i++)
6459 ar & BOOST_SERIALIZATION_NVP( M_WN[i] );
6460 for(
int i=0; i<M_N; i++)
6461 ar & BOOST_SERIALIZATION_NVP( M_WNdu[i] );
6464 ar & BOOST_SERIALIZATION_NVP( M_maxerror );
6465 ar & BOOST_SERIALIZATION_NVP( M_use_newton );
6466 ar & BOOST_SERIALIZATION_NVP( M_Jqm_pr );
6467 ar & BOOST_SERIALIZATION_NVP( M_Rqm_pr );
6471 template<
typename TruthModelType>
6472 template<
class Archive>
6474 CRB<TruthModelType>::load( Archive & ar,
const unsigned int version )
6479 int proc_number = this->worldComm().globalRank();
6481 LOG(INFO) <<
"[CRB::load] version"<< version <<
"\n";
6489 LOG(INFO) <<
"[load] model not initialized, loading fdb files...\n";
6490 mesh = mesh_type::New();
6491 bool is_mesh_loaded = mesh->load( _name=
"mymesh",_path=this->
dbLocalPath(),_type=
"binary" );
6492 Xh = space_type::New( mesh );
6493 LOG(INFO) <<
"[load] loading fdb files done.\n";
6497 LOG(INFO) <<
"[load] get mesh/Xh from model...\n";
6498 mesh = M_model->functionSpace()->mesh();
6499 Xh = M_model->functionSpace();
6500 LOG(INFO) <<
"[load] get mesh/Xh from model done.\n";
6504 typedef boost::bimap< int, double > old_convergence_type;
6505 ar & boost::serialization::base_object<super>( *this );
6506 ar & BOOST_SERIALIZATION_NVP( M_output_index );
6507 ar & BOOST_SERIALIZATION_NVP( M_N );
6509 ar & BOOST_SERIALIZATION_NVP( M_rbconv );
6511 ar & BOOST_SERIALIZATION_NVP( M_error_type );
6512 ar & BOOST_SERIALIZATION_NVP( M_Xi );
6513 ar & BOOST_SERIALIZATION_NVP( M_WNmu );
6514 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr );
6515 ar & BOOST_SERIALIZATION_NVP( M_Aqm_du );
6516 ar & BOOST_SERIALIZATION_NVP( M_Aqm_pr_du );
6517 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
6518 ar & BOOST_SERIALIZATION_NVP( M_Fqm_du );
6519 ar & BOOST_SERIALIZATION_NVP( M_Lqm_pr );
6520 ar & BOOST_SERIALIZATION_NVP( M_Lqm_du );
6521 ar & BOOST_SERIALIZATION_NVP( M_C0_pr );
6522 ar & BOOST_SERIALIZATION_NVP( M_C0_du );
6523 ar & BOOST_SERIALIZATION_NVP( M_Lambda_pr );
6524 ar & BOOST_SERIALIZATION_NVP( M_Lambda_du );
6525 ar & BOOST_SERIALIZATION_NVP( M_Gamma_pr );
6526 ar & BOOST_SERIALIZATION_NVP( M_Gamma_du );
6528 ar & BOOST_SERIALIZATION_NVP( M_Mqm_pr );
6529 ar & BOOST_SERIALIZATION_NVP( M_Mqm_du );
6530 ar & BOOST_SERIALIZATION_NVP( M_Mqm_pr_du );
6532 if ( model_type::is_time_dependent )
6534 ar & BOOST_SERIALIZATION_NVP( M_coeff_pr_ini_online );
6535 ar & BOOST_SERIALIZATION_NVP( M_coeff_du_ini_online );
6536 ar & BOOST_SERIALIZATION_NVP( M_Cmf_pr );
6537 ar & BOOST_SERIALIZATION_NVP( M_Cma_pr );
6538 ar & BOOST_SERIALIZATION_NVP( M_Cmm_pr );
6539 ar & BOOST_SERIALIZATION_NVP( M_Cmf_du );
6540 ar & BOOST_SERIALIZATION_NVP( M_Cmf_du_ini );
6541 ar & BOOST_SERIALIZATION_NVP( M_Cma_du );
6542 ar & BOOST_SERIALIZATION_NVP( M_Cmm_du );
6545 ar & BOOST_SERIALIZATION_NVP ( M_database_contains_variance_info );
6546 if( M_database_contains_variance_info )
6547 ar & BOOST_SERIALIZATION_NVP( M_variance_matrix_phi );
6548 ar & BOOST_SERIALIZATION_NVP( M_Fqm_pr );
6549 ar & BOOST_SERIALIZATION_NVP( M_InitialGuessV_pr );
6551 ar & BOOST_SERIALIZATION_NVP( M_current_mu );
6552 ar & BOOST_SERIALIZATION_NVP( M_no_residual_index );
6554 ar & BOOST_SERIALIZATION_NVP( M_maxerror );
6555 ar & BOOST_SERIALIZATION_NVP( M_use_newton );
6556 ar & BOOST_SERIALIZATION_NVP( M_Jqm_pr );
6557 ar & BOOST_SERIALIZATION_NVP( M_Rqm_pr );
6559 if( this->
vm()[
"crb.use-newton"].
template as<bool>() != M_use_newton )
6562 throw std::logic_error(
"[CRB::loadDB] ERROR in the database used the option use-newton=true and it's not the case in your option" );
6564 throw std::logic_error(
"[CRB::loadDB] ERROR in the database used the option use-newton=false and it's not the case in your option" );
6568 std::cout <<
"[loadDB] output index : " << M_output_index <<
"\n"
6569 <<
"[loadDB] N : " << M_N <<
"\n"
6570 <<
"[loadDB] error type : " << M_error_type <<
"\n";
6572 for (
auto it = M_rbconv.begin(), en = M_rbconv.end();it != en; ++it )
6573 std::cout <<
"[loadDB] convergence: (" << it->left <<
"," << it->right <<
")\n";
6578 M_WNdu.resize( M_N );
6580 for(
int i = 0 ; i < M_N ; i++ )
6582 temp.setName( (boost::format(
"fem-primal-%1%" ) % ( i ) ).str() );
6583 ar & BOOST_SERIALIZATION_NVP( temp );
6587 for(
int i = 0 ; i < M_N ; i++ )
6589 temp.setName( (boost::format(
"fem-dual-%1%" ) % ( i ) ).str() );
6590 ar & BOOST_SERIALIZATION_NVP( temp );
6595 LOG(INFO) <<
"[CRB::load] end of load function" << std::endl;
6599 template<
typename TruthModelType>
6603 bool print = this->
vm()[
"crb.print-error-during-rb-construction"].template as<bool>();
6607 template<
typename TruthModelType>
6611 bool rebuild = this->
vm()[
"crb.rebuild-database"].template as<bool>();
6615 template<
typename TruthModelType>
6619 bool show = this->
vm()[
"crb.show-mu-selection"].template as<bool>();
6624 template<
typename TruthModelType>
6632 boost::archive::text_oarchive oa( ofs );
6639 template<
typename TruthModelType>
6654 if ( !fs::exists( db ) )
6658 fs::ifstream ifs( db );
6662 boost::archive::text_iarchive ia( ifs );
6666 this->setIsLoaded(
true );
6678 namespace serialization
6680 template<
typename T>
6681 struct version< Feel::
CRB<T> >
6686 typedef mpl::int_<0> type;
6687 typedef mpl::integral_c_tag tag;
6688 static const unsigned int value = version::type::value;
6690 template<
typename T>
const unsigned int version<Feel::CRB<T> >::value;
void computeErrorEstimationEfficiencyIndicator(parameterspace_ptrtype const &Dmu, double &max_ei, double &min_ei, int N=4)
Definition: crb.hpp:3334
void setDBFilename(std::string const &filename)
set the DB filename
Definition: crbdb.hpp:157
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
sampling_ptrtype trainSampling() const
Definition: crb.hpp:434
CRBErrorType
Definition: crb.hpp:83
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
brief description
Definition: crbdb.hpp:50
int outputIndex() const
Definition: crb.hpp:422
model_type::mesh_type mesh_type
mesh type
Definition: crb.hpp:193
PreconditionerType
Definition: feelalg/enums.hpp:114
std::string const & name() const
Definition: crbdb.hpp:112
boost::tuple< double, double, solutions_tuple, double, double, double, upper_bounds_tuple > run(parameter_type const &mu, double eps=1e-6, int N=-1)
Definition: crb.hpp:5985
~CRB()
destructor
Definition: crb.hpp:378
export Feel generated data to some file formatsUse the visitor and factory pattern.
Definition: exporter.hpp:82
void offlineResidual(int Ncur, int number_of_added_elements=1)
Definition: crb.hpp:5095
sampling_type randomSampling(int N)
Definition: crb.hpp:962
Non linear solver base interface.
Definition: solvernonlinear.hpp:62
base class for preconditioner
Definition: preconditioner.hpp:54
void setOfflineStep(bool b)
set boolean indicates if we are in offline_step or not
Definition: crb.hpp:521
model_type::space_type space_type
space type
Definition: crb.hpp:197
bool printErrorDuringOfflineStep()
Definition: crb.hpp:6601
bool isDBLoaded() const
Definition: crbdb.hpp:145
value_type delta(size_type N, parameter_ptrtype const &mu, std::vector< vectorN_type > const &uN, std::vector< vectorN_type > const &uNdu, std::vector< vectorN_type > const &uNold, std::vector< vectorN_type > const &uNduold, int k=0) const
Definition: crb.hpp:848
residual_error_type transientPrimalResidual(int Ncur, parameter_type const &mu, vectorN_type const &Un, vectorN_type const &Unold=vectorN_type(), double time_step=1, double time=1e30) const
Definition: crb.hpp:4710
Certifed Reduced Basis class.
Definition: crb.hpp:95
fs::path lookForDB() const
Definition: crbdb.cpp:125
CRB(std::string name, po::variables_map const &vm, truth_model_ptrtype const &model)
constructor from command line options
Definition: crb.hpp:287
Backward differencing formula time discretization.
Definition: bdf.hpp:91
void setTolerance(double tolerance)
set offline tolerance
Definition: crb.hpp:485
scm_ptrtype scm() const
Definition: crb.hpp:446
double computeConditioning(matrixN_type &A) const
Definition: crb.hpp:3595
void setOutputIndex(uint16_type oindex)
set the output index
Definition: crb.hpp:460
bool showMuSelection()
Definition: crb.hpp:6617
void setTruthModel(truth_model_ptrtype const &model)
set the truth offline model
Definition: crb.hpp:491
boost::tuple< double, double > lb(parameter_ptrtype const &mu, size_type N, std::vector< vectorN_type > &uN, std::vector< vectorN_type > &uNdu) const
Definition: crb.hpp:814
CRB(std::string name, po::variables_map const &vm)
constructor from command line options
Definition: crb.hpp:249
void computationalTimeStatistics(std::string appname)
Definition: crb.hpp:6325
void check(size_type N) const
Definition: crb.hpp:3276
std::string const & dbFilename() const
Definition: crbdb.hpp:118
CRBSCM< truth_model_type > scm_type
scm
Definition: crb.hpp:139
static backend_ptrtype build(BackendType=BACKEND_GMM, WorldComm const &worldComm=Environment::worldComm())
Definition: backend.cpp:167
int dimension() const
Definition: crb.hpp:428
type_traits< T >::real_type inner_product(Vector< T > const &v1, Vector< T > const &v2)
Definition: vector.hpp:579
model_type::functionspace_type functionspace_type
function space type
Definition: crb.hpp:151
void printErrorsDuringRbConstruction(void)
Definition: crb.hpp:5886
Parameter space sampling class.
Definition: parameterspace.hpp:219
value_type ub(size_type K, parameter_ptrtype const &mu, std::vector< vectorN_type > &uN, std::vector< vectorN_type > &uNdu) const
Definition: crb.hpp:874
element_type expansion(parameter_type const &mu, int N=-1, int time_index=-1)
Definition: crb.hpp:5935
parameterspace_ptrtype Dmu() const
Definition: crb.hpp:416
POD< truth_model_type > pod_type
POD.
Definition: crb.hpp:147
void setFactor(int Factor)
set factor
Definition: crb.hpp:515
void printMuSelection(void)
Definition: crb.hpp:5914
sampling_type equidistributedSampling(int N)
Definition: crb.hpp:971
boost::tuple< double, double > lb(size_type N, parameter_type const &mu, std::vector< vectorN_type > &uN, std::vector< vectorN_type > &uNdu, std::vector< vectorN_type > &uNold, std::vector< vectorN_type > &uNduold, int K=0) const
Definition: crb.hpp:4001
CRB(CRB const &o)
copy constructor
Definition: crb.hpp:345
max_error_type maxErrorBounds(size_type N) const
Retuns maximum value of the relative error.
Definition: crb.hpp:4287
void saveDB()
Definition: crb.hpp:6626
POD class.
Definition: pod.hpp:91
void exportBasisFunctions(const export_vector_wn_type &wn) const
Definition: crb.hpp:4495
CRBErrorType errorType() const
Definition: crb.hpp:440
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Parameter space class.
Definition: parameterspace.hpp:58
double orthonormalize(size_type N, wn_type &wn, int Nm=1)
Definition: crb.hpp:4420
SCM algorithm.
Definition: crbscm.hpp:69
CRB & operator=(CRB const &o)
copy operator
Definition: crb.hpp:389
model_type::element_type element_type
element of the functionspace type
Definition: crb.hpp:156
int maxIter() const
Definition: crb.hpp:410
convergence_type offline()
Definition: crb.hpp:1594
CRB()
default constructor
Definition: crb.hpp:228
po::variables_map vm()
Definition: crbdb.hpp:133
value_type ub(size_type N, parameter_type const &mu, std::vector< vectorN_type > &uN, std::vector< vectorN_type > &uNdu) const
Definition: crb.hpp:831
bool loadDB()
Definition: crb.hpp:6641
int factor() const
return factor
Definition: crb.hpp:405
void fixedPoint(size_type N, parameter_type const &mu, std::vector< vectorN_type > &uN, std::vector< vectorN_type > &uNdu, std::vector< vectorN_type > &uNold, std::vector< vectorN_type > &uNduold, double &condition_number, std::vector< double > &output_vector, int K=0) const
Definition: crb.hpp:3949
CRBElementsDB< truth_model_type > crb_elements_db_type
elements database
Definition: crb.hpp:143
Bdf< space_type > bdf_type
time discretization
Definition: crb.hpp:201
element of a parameter space
Definition: parameterspace.hpp:84
fs::path dbLocalPath() const
Definition: crbdb.cpp:95
bool rebuildDB()
Definition: crb.hpp:6609
void projectionOnPodSpace(const element_type &u, element_ptrtype &projection, const std::string &name_of_space="primal")
Definition: crb.hpp:6221
void setMaxIter(int K)
set max iteration number
Definition: crb.hpp:509
void setCRBErrorType(CRBErrorType error)
set the crb error type
Definition: crb.hpp:479