Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pod.hpp
Go to the documentation of this file.
1 /* -*- mode: c++ -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Stephane Veys <stephane.veys@imag.fr>
7  Date: 2009-09-30
8 
9  Copyright (C) 2011 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 2.1 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
32 #ifndef __POD_H
33 #define __POD_H 1
34 
35 
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>
45 #include <fstream>
46 
47 #include <feel/feelfilters/exporter.hpp>
48 #include <feel/feeldiscr/functionspace.hpp>
49 
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>
55 
56 
57 #include <Eigen/Core>
58 #include <Eigen/LU>
59 #include <Eigen/Dense>
60 
61 #include <vector>
62 
64 #include <feel/feelcore/feel.hpp>
66 #include <feel/feelcore/parameter.hpp>
68 #include <feel/feelcrb/crbdb.hpp>
69 #include <feel/feelcrb/crbscm.hpp>
71 
72 #include <feel/feelvf/vf.hpp>
73 #include <feel/feeldiscr/bdf2.hpp>
74 
75 namespace Feel
76 {
77 
78 
79 
90 template<typename TruthModelType>
91 class POD
92 {
93 public :
94 
95  typedef TruthModelType truth_model_type;
96  typedef truth_model_type model_type;
97  typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
98 
99 
100  typedef typename model_type::value_type value_type;
101 
103  typedef typename model_type::element_type element_type;
104  typedef typename model_type::element_ptrtype element_ptrtype;
105 
106 
107 
109  typedef typename model_type::mesh_type mesh_type;
110 
112  typedef typename model_type::mesh_ptrtype mesh_ptrtype;
113 
115  typedef typename model_type::functionspace_type functionspace_type;
116  typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
117 
118  typedef typename model_type::space_type space_type;
119 
122  typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
123 
124 
125  /* export */
127  typedef boost::shared_ptr<export_type> export_ptrtype;
128 
129  typedef std::vector<element_type> wn_type;
130 
131  typedef std::vector<element_type> mode_set_type;
132 
133 
134  typedef Eigen::VectorXd vectorN_type;
135  typedef Eigen::MatrixXd matrixN_type;
136 
137  typedef typename model_type::vector_ptrtype vector_ptrtype;
138 
139  typedef typename model_type::backend_type backend_type;
140  typedef boost::shared_ptr<backend_type> backend_ptrtype;
141 
142 
143  POD()
144  :
145  M_store_pod_matrix ( false ),
146  M_store_pod_matrix_format_octave ( false ),
147  M_Nm( 1 ),
148  M_pod_matrix(),
149  M_snapshots_matrix(),
150  M_model(),
151  M_WN()
152  {}
153 
154 
155 
156  POD( po::variables_map const& vm, const wn_type& WN, const int Nm, const int K, const matrixN_type& SnapshotsMatrix )
157  :
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>() ),
160  M_Nm( Nm ),
161  M_pod_matrix(),
162  M_snapshots_matrix( SnapshotsMatrix ),
163  M_model(),
164  M_WN( WN ),
165  M_backend( backend_type::build( vm ) )
166  {}
167 
171  POD( POD const & o )
172  :
173  M_store_pod_matrix( o.M_store_pod_matrix ),
174  M_store_pod_matrix_format_octave( o.M_store_pod_matrix_format_octave ),
175  M_Nm( o.M_Nm ),
176  M_pod_matrix( o.M_matrix ),
177  M_snapshots_matrix( o.M_snapshots_matrix ),
178  M_model( o.M_model ),
179  M_WN( o.M_WN )
180  {}
181 
184  {}
185 
186 
187 
189  bool storePodMatrix() const
190  {
191  return M_store_pod_matrix;
192  }
193 
196  {
197  return M_store_pod_matrix_format_octave;
198  }
199 
201  const matrixN_type& podMatrix()
202  {
203  return M_pod_matrix;
204  }
205 
207  const matrixN_type& snapshotsMatrix()
208  {
209  return M_snapshots_matrix;
210  }
211 
213  const int nm()
214  {
215  return M_Nm;
216  }
217 
219  const wn_type & wn()
220  {
221  return M_WN;
222  }
223 
225  const truth_model_ptrtype & model()
226  {
227  return M_model;
228  }
229 
230  const double timeInitial()
231  {
232  return M_time_initial;
233  }
234 
235  void setBdf( bdf_ptrtype& bdf )
236  {
237  M_bdf = bdf;
238  }
239 
240  void setSnapshotsMatrix ( matrixN_type& Matrix )
241  {
242  M_snapshots_matrix=Matrix;
243  }
244 
245  void setWN ( wn_type& WN )
246  {
247  M_WN=WN;
248  }
249 
250  void setNm ( const int Nm )
251  {
252  M_Nm=Nm;
253  }
254 
255  void setModel ( truth_model_ptrtype Model )
256  {
257  M_model=Model;
258  }
259 
260  void setTimeInitial( double Ti )
261  {
262  M_time_initial = Ti;
263  }
264 
266  void fillPodMatrix();
267 
272  int pod( mode_set_type& ModeSet, bool is_primal );
273 
274  void exportMode( double time, element_ptrtype& mode );
275 
276  void projectionOnPodSpace();
277 
278 private :
279 
280  bool M_store_pod_matrix;
281 
282  bool M_store_pod_matrix_format_octave;
283 
284  int M_Nm;
285 
286  matrixN_type M_pod_matrix;
287 
288  matrixN_type M_snapshots_matrix;
289 
290  truth_model_ptrtype M_model;
291 
292  wn_type M_WN;
293 
294  backend_ptrtype M_backend;
295 
296  export_ptrtype exporter;
297 
298  bdf_ptrtype M_bdf;
299 
300  double M_time_initial;
301 };//class POD
302 
303 
304 template<typename TruthModelType>
305 void POD<TruthModelType>::exportMode( double time, element_ptrtype& mode )
306 {
307 
308  LOG(INFO) << "exportResults starts\n";
309 
310  functionspace_ptrtype function_space = M_model->functionSpace();
311  mesh_ptrtype mesh = function_space->mesh();
312 
313  //1D
314  std::ofstream mode_file;
315  mode_file.open( "mode.dat",std::ios::out );
316  std::map<double,double> data;
317 
318  for ( auto it = mesh->beginElement( ),
319  en = mesh->endElement( );
320  it!=en; ++it )
321  {
322  for ( size_type i = 0; i < space_type::basis_type::nLocalDof; ++i )
323  {
324  value_type a = it->point( 0 ).node()[0];
325  value_type b = it->point( 1 ).node()[0];
326  value_type x = 0;
327 
328  if ( i == 0 )
329  x=a;
330 
331  else if ( i == 1 )
332  x=b;
333 
334  else
335  x= a + ( i-1 )*( b-a )/( space_type::basis_type::nLocalDof-1 );
336 
337  data[x] = mode->localToGlobal( it->id(), i, 0 );
338  }
339 
340  }
341 
342  BOOST_FOREACH( auto d, data )
343  {
344  mode_file << d.first << " " << d.second << "\n";
345  }
346  mode_file.close();
347 
348  //2D - 3D
349  //exporter->step(time)->setMesh( mesh );
350  //exporter->step(time)->add( "mode", *mode );
351  //exporter->save();
352 
353 }
354 
355 template<typename TruthModelType>
357 {
358  //M_bdf->setRestart( true );
359  int K = M_bdf->timeValues().size()-1;
360  M_pod_matrix.resize( K,K );
361 
362  auto bdfi = M_bdf->deepCopy();
363  auto bdfj = M_bdf->deepCopy();
364 
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() )
372  {
373  int i = bdfi->iteration()-1;
374  bdfi->loadCurrent();
375 
376  //here is a barrier. For now, if this barrier is removed, during the bdfj->restart() the program will wait and the memory increase.
377  //during the init step in bdf, the metadata are read, and it is during this step that the memory could increase dangerously.
378  Environment::worldComm().barrier();
379 
380  for ( bdfj->restart(); !bdfj->isFinished() && ( bdfj->iteration() < bdfi->iteration() ); bdfj->next() )
381  {
382  int j = bdfj->iteration()-1;
383  bdfj->loadCurrent();
384 
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 );
387  }
388 
389  M_pod_matrix( i,i ) = M_model->scalarProductForPod( bdfi->unknown( 0 ), bdfi->unknown( 0 ) );
390  }
391 }//fillPodMatrix
392 
393 
394 //the bool is_primal is used to determine the max number of modes ( M_Nm ) whithout taking eigenvectors
395 //associated with too small eigenvalues
396 //Since we always use M_Nm = 1 or 2 (so we never have to modify M_Nm) we only correct M_Nm in the primal case
397 template<typename TruthModelType>
398 int POD<TruthModelType>::pod( mode_set_type& ModeSet, bool is_primal )
399 {
400  M_backend = backend_type::build( BACKEND_PETSC );
401 
402  Eigen::SelfAdjointEigenSolver< matrixN_type > eigen_solver;
403 
404  fillPodMatrix();
405  //store the matrix
406  if ( M_store_pod_matrix )
407  {
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();
412  matrix_file<<"\n";
413  matrix_file<<M_pod_matrix.cols();
414  matrix_file<<"\n";
415 
416  for ( int i=0; i<M_pod_matrix.rows(); i++ )
417  {
418  for ( int j=0; j<M_pod_matrix.cols(); j++ )
419  {
420  matrix_file<< std::setprecision( 16 ) <<M_pod_matrix( i,j )<<" ";
421  }
422 
423  matrix_file<<"\n";
424  }
425 
426  LOG(INFO)<<" matrix wrote in file named PodMatrix \n";
427  }
428 
429  else if ( M_store_pod_matrix_format_octave )
430  {
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";
438 
439  for ( int i=0; i<M_pod_matrix.rows(); i++ )
440  {
441  for ( int j=0; j<M_pod_matrix.cols(); j++ )
442  {
443  matrix_file<< std::setprecision( 16 ) <<M_pod_matrix( i,j )<<" ";
444  }
445 
446  matrix_file<<"\n";
447  }
448 
449  LOG(INFO)<<" matrix wrote in file named PodMatrix \n";
450  matrix_file.close();
451  }
452 
453 
454  eigen_solver.compute( M_pod_matrix ); // solve M_pod_matrix psi = lambda psi
455 
456  int number_of_eigenvalues = eigen_solver.eigenvalues().size();
457  LOG(INFO)<<"Number of eigenvalues : "<<number_of_eigenvalues<<"\n";
458  //we copy eigenvalues in a std::vector beacause it's easier to manipulate it
459  std::vector<double> eigen_values( number_of_eigenvalues );
460 
461  int too_small_index = -1;
462 
463  for ( int i=0; i<number_of_eigenvalues; i++ )
464  {
465  if ( imag( eigen_solver.eigenvalues()[i] )>1e-12 )
466  {
467  throw std::logic_error( "[POD::pod] ERROR : complex eigenvalues were found" );
468  }
469 
470  eigen_values[i]=real( eigen_solver.eigenvalues()[i] );
471 
472  if ( eigen_values[i] < 1e-11 ) too_small_index=i+1;
473  }
474 
475  int position_of_largest_eigenvalue=number_of_eigenvalues-1;
476  int number_of_good_eigenvectors = number_of_eigenvalues - too_small_index;
477 
478  if ( M_Nm > number_of_good_eigenvectors && number_of_good_eigenvectors>0 && is_primal )
479  M_Nm=number_of_good_eigenvectors;
480 
481  for ( int i=0; i<M_Nm; i++ )
482  {
483  element_ptrtype mode ( new element_type( M_model->functionSpace() ) );
484  mode->zero();
485 
486  M_bdf->setRestart( true );
487  int index=0;
488 
489  for ( M_bdf->restart(); !M_bdf->isFinished(); M_bdf->next() )
490  {
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 ) );
495  index++;
496  }
497 
498  --position_of_largest_eigenvalue;
499  ModeSet.push_back( *mode );
500  }
501 
502  if ( M_store_pod_matrix_format_octave )
503  {
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";
511 
512  for ( int i=0; i<number_of_eigenvalues; i++ )
513  {
514  eigenvalue_file<<std::setprecision( 16 )<<eigen_values[i]<<"\n";
515  }
516 
517  eigenvalue_file.close();
518 
519 
520  std::ofstream eigenvector_file( ( boost::format( "eigen_vectors" ) ).str().c_str() );
521 
522  for ( int j=0; j<M_pod_matrix.cols(); j++ )
523  {
524  for ( int i=0; i<number_of_eigenvalues; i++ )
525  {
526  eigenvector_file<<std::setprecision( 16 )<<eigen_solver.eigenvectors().col( i )[j] <<" ";
527  }
528 
529  eigenvector_file<<"\n";
530  }
531 
532  eigenvector_file.close();
533 
534 
535  for ( int i=0; i<M_Nm; i++ )
536  {
537  std::ofstream mode_file( ( boost::format( "mode_%1%" ) %i ).str().c_str() );
538  element_type e = ModeSet[i];
539 
540  for ( size_type j=0; j<e.size(); j++ ) mode_file<<std::setprecision( 16 )<<e( j )<<"\n";
541 
542  mode_file.close();
543  }
544 
545  }
546 
547  return M_Nm;
548 
549 }
550 
551 
552 }//namespace Feel
553 
554 #endif /* __POD_H */
555 
556 

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4