36 #include <boost/multi_array.hpp>
37 #include <boost/tuple/tuple.hpp>
38 #include "boost/tuple/tuple_io.hpp"
39 #include <boost/format.hpp>
40 #include <boost/foreach.hpp>
41 #include <boost/bimap.hpp>
42 #include <boost/bimap/support/lambda.hpp>
43 #include <boost/archive/text_oarchive.hpp>
44 #include <boost/archive/text_iarchive.hpp>
47 #include <feel/feelfilters/exporter.hpp>
48 #include <feel/feeldiscr/functionspace.hpp>
50 #include <boost/serialization/vector.hpp>
51 #include <boost/serialization/list.hpp>
52 #include <boost/serialization/string.hpp>
53 #include <boost/serialization/version.hpp>
54 #include <boost/serialization/split_member.hpp>
59 #include <Eigen/Dense>
64 #include <feel/feelcore/feel.hpp>
66 #include <feel/feelcore/parameter.hpp>
69 #include <feel/feelcrb/crbscm.hpp>
73 #include <feel/feeldiscr/bdf2.hpp>
90 template<
typename TruthModelType>
95 typedef TruthModelType truth_model_type;
96 typedef truth_model_type model_type;
97 typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
100 typedef typename model_type::value_type value_type;
104 typedef typename model_type::element_ptrtype element_ptrtype;
116 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
118 typedef typename model_type::space_type space_type;
122 typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
127 typedef boost::shared_ptr<export_type> export_ptrtype;
129 typedef std::vector<element_type> wn_type;
131 typedef std::vector<element_type> mode_set_type;
134 typedef Eigen::VectorXd vectorN_type;
135 typedef Eigen::MatrixXd matrixN_type;
137 typedef typename model_type::vector_ptrtype vector_ptrtype;
139 typedef typename model_type::backend_type backend_type;
140 typedef boost::shared_ptr<backend_type> backend_ptrtype;
145 M_store_pod_matrix ( false ),
146 M_store_pod_matrix_format_octave ( false ),
149 M_snapshots_matrix(),
156 POD( po::variables_map
const& vm,
const wn_type& WN,
const int Nm,
const int K,
const matrixN_type& SnapshotsMatrix )
158 M_store_pod_matrix( vm[
"pod.store-pod-matrix"].template as<bool>() ),
159 M_store_pod_matrix_format_octave( vm[
"pod.store-pod-matrix-format-octave"].template as<bool>() ),
162 M_snapshots_matrix( SnapshotsMatrix ),
173 M_store_pod_matrix( o.M_store_pod_matrix ),
174 M_store_pod_matrix_format_octave( o.M_store_pod_matrix_format_octave ),
176 M_pod_matrix( o.M_matrix ),
177 M_snapshots_matrix( o.M_snapshots_matrix ),
178 M_model( o.M_model ),
191 return M_store_pod_matrix;
197 return M_store_pod_matrix_format_octave;
209 return M_snapshots_matrix;
230 const double timeInitial()
232 return M_time_initial;
235 void setBdf( bdf_ptrtype& bdf )
240 void setSnapshotsMatrix ( matrixN_type& Matrix )
242 M_snapshots_matrix=Matrix;
245 void setWN ( wn_type& WN )
250 void setNm (
const int Nm )
255 void setModel ( truth_model_ptrtype Model )
260 void setTimeInitial(
double Ti )
272 int pod( mode_set_type& ModeSet,
bool is_primal );
274 void exportMode(
double time, element_ptrtype& mode );
276 void projectionOnPodSpace();
280 bool M_store_pod_matrix;
282 bool M_store_pod_matrix_format_octave;
286 matrixN_type M_pod_matrix;
288 matrixN_type M_snapshots_matrix;
290 truth_model_ptrtype M_model;
294 backend_ptrtype M_backend;
296 export_ptrtype exporter;
300 double M_time_initial;
304 template<
typename TruthModelType>
305 void POD<TruthModelType>::exportMode(
double time, element_ptrtype& mode )
308 LOG(INFO) <<
"exportResults starts\n";
310 functionspace_ptrtype function_space = M_model->functionSpace();
311 mesh_ptrtype mesh = function_space->mesh();
314 std::ofstream mode_file;
315 mode_file.open(
"mode.dat",std::ios::out );
316 std::map<double,double> data;
318 for (
auto it = mesh->beginElement( ),
319 en = mesh->endElement( );
322 for (
size_type i = 0; i < space_type::basis_type::nLocalDof; ++i )
324 value_type a = it->point( 0 ).node()[0];
325 value_type b = it->point( 1 ).node()[0];
335 x= a + ( i-1 )*( b-a )/( space_type::basis_type::nLocalDof-1 );
337 data[x] = mode->localToGlobal( it->id(), i, 0 );
342 BOOST_FOREACH(
auto d, data )
344 mode_file << d.first <<
" " << d.second <<
"\n";
355 template<
typename TruthModelType>
359 int K = M_bdf->timeValues().size()-1;
360 M_pod_matrix.resize( K,K );
362 auto bdfi = M_bdf->deepCopy();
363 auto bdfj = M_bdf->deepCopy();
365 bdfi->setRestart(
true );
366 bdfj->setRestart(
true );
367 bdfi->setTimeInitial( M_time_initial );
368 bdfj->setTimeInitial( M_time_initial );
369 bdfi->setRestartAtLastSave(
false);
370 bdfj->setRestartAtLastSave(
false);
371 for ( bdfi->restart(); !bdfi->isFinished(); bdfi->next() )
373 int i = bdfi->iteration()-1;
378 Environment::worldComm().barrier();
380 for ( bdfj->restart(); !bdfj->isFinished() && ( bdfj->iteration() < bdfi->iteration() ); bdfj->next() )
382 int j = bdfj->iteration()-1;
385 M_pod_matrix( i,j ) = M_model->scalarProductForPod( bdfj->unknown( 0 ), bdfi->unknown( 0 ) );
386 M_pod_matrix( j,i ) = M_pod_matrix( i,j );
389 M_pod_matrix( i,i ) = M_model->scalarProductForPod( bdfi->unknown( 0 ), bdfi->unknown( 0 ) );
397 template<
typename TruthModelType>
402 Eigen::SelfAdjointEigenSolver< matrixN_type > eigen_solver;
406 if ( M_store_pod_matrix )
408 std::ofstream matrix_file;
409 LOG(INFO)<<
"saving Pod matrix in a file \n";
410 matrix_file.open(
"PodMatrix",std::ios::out );
411 matrix_file<<M_pod_matrix.rows();
413 matrix_file<<M_pod_matrix.cols();
416 for (
int i=0; i<M_pod_matrix.rows(); i++ )
418 for (
int j=0; j<M_pod_matrix.cols(); j++ )
420 matrix_file<< std::setprecision( 16 ) <<M_pod_matrix( i,j )<<
" ";
426 LOG(INFO)<<
" matrix wrote in file named PodMatrix \n";
429 else if ( M_store_pod_matrix_format_octave )
431 std::ofstream matrix_file;
432 LOG(INFO)<<
"saving Pod matrix in a file \n";
433 matrix_file.open(
"PodMatrixOctave.mat",std::ios::out );
434 matrix_file<<
"# name: A\n";
435 matrix_file<<
"# type: matrix\n";
436 matrix_file<<
"# rows: "<<M_pod_matrix.rows()<<
"\n";
437 matrix_file<<
"# columns: "<<M_pod_matrix.cols()<<
"\n";
439 for (
int i=0; i<M_pod_matrix.rows(); i++ )
441 for (
int j=0; j<M_pod_matrix.cols(); j++ )
443 matrix_file<< std::setprecision( 16 ) <<M_pod_matrix( i,j )<<
" ";
449 LOG(INFO)<<
" matrix wrote in file named PodMatrix \n";
454 eigen_solver.compute( M_pod_matrix );
456 int number_of_eigenvalues = eigen_solver.eigenvalues().size();
457 LOG(INFO)<<
"Number of eigenvalues : "<<number_of_eigenvalues<<
"\n";
459 std::vector<double> eigen_values( number_of_eigenvalues );
461 int too_small_index = -1;
463 for (
int i=0; i<number_of_eigenvalues; i++ )
465 if ( imag( eigen_solver.eigenvalues()[i] )>1e-12 )
467 throw std::logic_error(
"[POD::pod] ERROR : complex eigenvalues were found" );
470 eigen_values[i]=real( eigen_solver.eigenvalues()[i] );
472 if ( eigen_values[i] < 1e-11 ) too_small_index=i+1;
475 int position_of_largest_eigenvalue=number_of_eigenvalues-1;
476 int number_of_good_eigenvectors = number_of_eigenvalues - too_small_index;
478 if ( M_Nm > number_of_good_eigenvectors && number_of_good_eigenvectors>0 && is_primal )
479 M_Nm=number_of_good_eigenvectors;
481 for (
int i=0; i<M_Nm; i++ )
483 element_ptrtype mode (
new element_type( M_model->functionSpace() ) );
486 M_bdf->setRestart(
true );
489 for ( M_bdf->restart(); !M_bdf->isFinished(); M_bdf->next() )
491 M_bdf->loadCurrent();
492 double psi_k = real( eigen_solver.eigenvectors().col( position_of_largest_eigenvalue )[index] );
493 M_bdf->unknown( 0 ).scale( psi_k );
494 mode->add( 1 , M_bdf->unknown( 0 ) );
498 --position_of_largest_eigenvalue;
499 ModeSet.push_back( *mode );
502 if ( M_store_pod_matrix_format_octave )
504 std::cout<<
"M_store_pod_matrix_format_octave = "<<M_store_pod_matrix_format_octave<<std::endl;
505 std::ofstream eigenvalue_file;
506 eigenvalue_file.open(
"eigen_values.mat",std::ios::out );
507 eigenvalue_file<<
"# name: E\n";
508 eigenvalue_file<<
"# type: matrix\n";
509 eigenvalue_file<<
"# rows: "<<number_of_eigenvalues<<
"\n";
510 eigenvalue_file<<
"# columns: 1\n";
512 for (
int i=0; i<number_of_eigenvalues; i++ )
514 eigenvalue_file<<std::setprecision( 16 )<<eigen_values[i]<<
"\n";
517 eigenvalue_file.close();
520 std::ofstream eigenvector_file( ( boost::format(
"eigen_vectors" ) ).str().c_str() );
522 for (
int j=0; j<M_pod_matrix.cols(); j++ )
524 for (
int i=0; i<number_of_eigenvalues; i++ )
526 eigenvector_file<<std::setprecision( 16 )<<eigen_solver.eigenvectors().col( i )[j] <<
" ";
529 eigenvector_file<<
"\n";
532 eigenvector_file.close();
535 for (
int i=0; i<M_Nm; i++ )
537 std::ofstream mode_file( ( boost::format(
"mode_%1%" ) %i ).str().c_str() );
540 for (
size_type j=0; j<e.size(); j++ ) mode_file<<std::setprecision( 16 )<<e( j )<<
"\n";