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
bdf2.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2006-12-30
7 
8  Copyright (C) 2006-2008 Université Joseph Fourier (Grenoble)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef _BDF_H
30 #define _BDF_H
31 
32 #include <string>
33 #include <iostream>
34 #include <sstream>
35 #include <algorithm>
36 
37 #include <boost/timer.hpp>
38 #include <boost/shared_array.hpp>
39 #include <boost/lambda/lambda.hpp>
40 #include <boost/lambda/bind.hpp>
41 #include <boost/utility.hpp>
42 
43 #include <boost/foreach.hpp>
44 #include <boost/filesystem.hpp>
45 #include <boost/filesystem/fstream.hpp>
46 #include <boost/serialization/vector.hpp>
47 #include <boost/serialization/array.hpp>
48 #include <boost/serialization/base_object.hpp>
49 #include <boost/archive/text_oarchive.hpp>
50 #include <boost/archive/text_iarchive.hpp>
51 #include <boost/archive/binary_oarchive.hpp>
52 #include <boost/archive/binary_iarchive.hpp>
53 
54 #include <boost/numeric/ublas/vector.hpp>
55 #include <boost/numeric/ublas/vector_proxy.hpp>
56 
57 #include <boost/parameter.hpp>
58 #include <feel/feelcore/parameter.hpp>
59 
60 #include <feel/feelcore/feel.hpp>
61 #include <feel/feelalg/glas.hpp>
62 #include <feel/feeldiscr/functionspace.hpp>
63 
64 namespace Feel
65 {
66 namespace ublas = boost::numeric::ublas;
67 namespace fs = boost::filesystem;
68 
69 enum BDFState { BDF_UNITIALIZED = 0, BDF_RUNNING, BDF_STOPPED };
70 
71 enum BDFTimeScheme { BDF_ORDER_ONE=1, BDF_ORDER_TWO, BDF_ORDER_THREE, BDF_ORDER_FOUR, BDF_MAX_ORDER = 4 };
72 enum BDFStragegy { BDF_STRATEGY_DT_CONSTANT,
73  BDF_STRATEGY_DT_ADAPTATIVE
74  };
75 
76 class BdfBase
77 {
78  friend class boost::serialization::access;
79 public:
80  typedef std::vector<double>::iterator time_iterator;
81  typedef std::vector<double>::const_iterator time_const_iterator;
82  typedef std::vector<double> time_values_map_type;
83 
84  BdfBase()
85  :
86  M_order( 1 ),
87  M_order_cur( 1 ),
88  M_name( "bdf" ),
89  M_time( 0.0 ),
90  M_iteration( 0 ),
91  M_Ti( 0.0 ),
92  M_Tf( 1.0 ),
93  M_dt( 1.0 ),
94  M_strategy( BDF_STRATEGY_DT_CONSTANT ),
95  M_state( BDF_UNITIALIZED ),
96  M_n_restart( 0 ),
97  M_restart( false ),
98  M_restartPath( "" ),
99  M_restartAtLastSave( false ),
100  M_alpha( BDF_MAX_ORDER ),
101  M_beta( BDF_MAX_ORDER ),
102  M_saveInFile( true ),
103  M_saveFreq( 1 ),
104  M_rankProcInNameOfFiles( false ),
105  M_worldComm( Environment::worldComm() )
106  {}
107 
108 #if 0
109  template <class ArgumentPack>
110  BdfBase( ArgumentPack const& args )
111  :
112  M_order( args[_order | 1] ),
113  M_name( args[_name | "bdf"] ),
114  M_time( args[_initial_time | 0] ),
115  M_iteration( 0 ),
116  M_Ti( args[_initial_time | 0] ),
117  M_Tf( args[_final_time | 1] ),
118  M_dt( args[_time_step | 0.1] ),
119  M_iterations_between_order_change( 1 ),
120  M_strategy( ( BDFStragegy )args[_strategy | BDF_STRATEGY_DT_CONSTANT] ),
121  M_state( BDF_UNITIALIZED ),
122  M_n_restart( 0 ),
123  M_restart( args[_restart | false] ),
124  M_restartPath( args[_restart_path | ""] ),
125  M_restartAtLastSave( args[_restart_at_last_save | false] )
126  M_alpha( BDF_MAX_ORDER ),
127  M_beta( BDF_MAX_ORDER ),
128  M_saveInFile( true ),
129  M_saveFreq( 1 ),
130  M_rankProcInNameOfFiles( false ),
131  M_worldComm( Environment::worldComm() )
132  {
133 
134  }
135 #endif
136  BdfBase( po::variables_map const& vm, std::string name, std::string const& prefix, WorldComm const& worldComm )
137  :
138  M_order( vm[prefixvm( prefix, "bdf.order" )].as<int>() ),
139  M_order_cur( M_order ),
140  M_name( name ),
141  M_time( vm[prefixvm( prefix, "bdf.time-initial" )].as<double>() ),
142  M_iteration( 0 ),
143  M_iterations_between_order_change( vm[prefixvm( prefix, "bdf.iterations-between-order-change" )].as<int>() ),
144  M_Ti( vm[prefixvm( prefix, "bdf.time-initial" )].as<double>() ),
145  M_Tf( vm[prefixvm( prefix, "bdf.time-final" )].as<double>() ),
146  M_dt( vm[prefixvm( prefix, "bdf.time-step" )].as<double>() ),
147  M_strategy( ( BDFStragegy )vm[prefixvm( prefix, "bdf.strategy" )].as<int>() ),
148  M_state( BDF_UNITIALIZED ),
149  M_n_restart( 0 ),
150  M_restart( vm[prefixvm( prefix, "bdf.restart" )].as<bool>() ),
151  M_restartPath( vm[prefixvm( prefix, "bdf.restart.path" )].as<std::string>() ),
152  M_restartAtLastSave( vm[prefixvm( prefix, "bdf.restart.at-last-save" )].as<bool>() ),
153  M_alpha( BDF_MAX_ORDER ),
154  M_beta( BDF_MAX_ORDER ),
155  M_saveInFile( vm[prefixvm( prefix, "bdf.save" )].as<bool>() ),
156  M_saveFreq( vm[prefixvm( prefix, "bdf.save.freq" )].as<int>() ),
157  M_rankProcInNameOfFiles( vm[prefixvm( prefix, "bdf.rank-proc-in-files-name" )].as<bool>() ),
158  M_worldComm( worldComm )
159  {
160  }
161  BdfBase( std::string name, WorldComm const& worldComm )
162  :
163  M_order( 1 ),
164  M_order_cur( 1 ),
165  M_name( name ),
166  M_time( 0. ),
167  M_iteration( 0 ),
168  M_iterations_between_order_change( 1 ),
169  M_Ti( 0. ),
170  M_Tf( 1.0 ),
171  M_dt( 1.0 ),
172  M_strategy( ( BDFStragegy )0 ),
173  M_state( BDF_UNITIALIZED ),
174  M_n_restart( 0 ),
175  M_restart( false ),
176  M_restartPath( "" ),
177  M_restartAtLastSave( false ),
178  M_alpha( BDF_MAX_ORDER ),
179  M_beta( BDF_MAX_ORDER ),
180  M_saveInFile( true ),
181  M_saveFreq( 1 ),
182  M_worldComm( worldComm )
183  {
184  }
185 
186  BdfBase( BdfBase const& b )
187  :
188  M_order( b.M_order ),
189  M_order_cur( b.M_order_cur ),
190  M_name( b.M_name ),
191  M_time( b.M_time ),
192  M_iteration( b.M_iteration ),
193  M_iterations_between_order_change( b.M_iterations_between_order_change ),
194  M_Ti( b.M_Ti ),
195  M_Tf( b.M_Tf ),
196  M_dt( b.M_dt ),
197  M_strategy( b.M_strategy ),
198  M_state( b.M_state ),
199  M_n_restart( b.M_n_restart ),
200  M_restart( b.M_restart ),
201  M_restartPath( b.M_restartPath ),
202  M_restartAtLastSave( b.M_restartAtLastSave ),
203  M_time_values_map( b.M_time_values_map ),
204  M_alpha( b.M_alpha ),
205  M_beta( b.M_beta ),
206  M_saveInFile( b.M_saveInFile ),
207  M_saveFreq( b.M_saveFreq ),
208  M_rankProcInNameOfFiles( b.M_rankProcInNameOfFiles ),
209  M_worldComm( b.M_worldComm )
210  {}
211 
212  virtual ~BdfBase() {}
213 
214  double polyCoefficient( int i ) const
215  {
216  FEELPP_ASSERT( i >=0 && i < BDF_MAX_ORDER-1 ).error( "[BDF] invalid index" );
217  return M_beta[this->timeOrder()-1][i];
218  }
219  double polyDerivCoefficient( int i ) const
220  {
221  FEELPP_ASSERT( i >=0 && i <= BDF_MAX_ORDER ).error( "[BDF] invalid index" );
222  return M_alpha[this->timeOrder()-1][i]/math::abs( this->timeStep() );
223  //return M_alpha[this->timeOrder()-1][i]/this->timeStep();
224  }
225 
226  BdfBase& operator=( BdfBase const& b )
227  {
228  if ( this != &b )
229  {
230  M_order = b.M_order;
231  M_order_cur = b.M_order_cur;
232  M_name = b.M_name;
233  M_time = b.M_time;
234  M_iteration = b.M_iteration;
235  M_Ti = b.M_Ti;
236  M_Tf = b.M_Tf;
237  M_dt = b.M_dt;
238  M_iterations_between_order_change = b.M_iterations_between_order_change;
239  M_n_restart = b.M_n_restart;
240  M_restart = b.M_restart;
241  M_restartPath = b.M_restartPath;
242  M_restartAtLastSave = b.M_restartAtLastSave;
243  M_strategy = b.M_strategy;
244  M_state = b.M_state;
245 
246  M_alpha = b.M_alpha;
247  M_beta = b.M_beta;
248  M_saveInFile = b.M_saveInFile;
249  M_rankProcInNameOfFiles = b.M_rankProcInNameOfFiles;
250 
251  M_time_values_map = b.M_time_values_map;
252  M_worldComm = b.M_worldComm;
253  }
254 
255  return *this;
256  }
257 
259  template<class Archive>
260  void serialize( Archive & ar, const unsigned int version )
261  {
262  DVLOG(2) << "[BDF::serialize] serialize BDFBase\n";
263 #if 0
264  ar & M_order;
265  ar & M_name;
266  ar & M_time;
267 
268  ar & M_n_restart;
269 
270  ar & M_Tf;
271 #endif
272  //ar & M_time_orders;
273  ar & boost::serialization::make_nvp( "time_values", M_time_values_map );
274 
275  DVLOG(2) << "[BDF::serialize] time orders size: " << M_time_orders.size() << "\n";
276  DVLOG(2) << "[BDF::serialize] time values size: " << M_time_values_map.size() << "\n";
277 
278  for ( auto it = M_time_values_map.begin(), en = M_time_values_map.end(); it!=en; ++it )
279  {
280  //LOG(INFO) << "[Bdf] order " << i << "=" << M_time_orders[i] << "\n";
281  DVLOG(2) << "[Bdf::serialize] value " << *it << "\n";
282 
283  }
284 
285  DVLOG(2) << "[BDF::serialize] serialize BDFBase done\n";
286  }
287 
289  int timeOrder() const
290  {
291  return M_order_cur;
292  }
293 
295  double timeInitial() const
296  {
297  return M_Ti;
298  }
299 
301  double timeFinal() const
302  {
303  return M_Tf;
304  }
305 
307  double timeStep() const
308  {
309  return M_dt;
310  }
311 
313  BDFStragegy strategy() const
314  {
315  return M_strategy;
316  }
317 
319  int numberOfIterationsBetweenOrderChange() const
320  {
321  return M_iterations_between_order_change;
322  }
323 
325  int numberOfIterationsSinceOrderChange() const
326  {
327  return M_iteration-M_last_iteration_since_order_change;
328  }
329 
331  int nRestart() const
332  {
333  return M_n_restart;
334  }
335 
337  bool isRestart() const
338  {
339  return M_restart;
340  }
341 
343  bool doRestartAtLastSave() const
344  {
345  return M_restartAtLastSave;
346  }
347 
349  double time() const
350  {
351  return M_time;
352  }
353 
355  int iteration() const
356  {
357  return M_iteration;
358  }
359 
361  double realTimePerIteration() const
362  {
363  FEELPP_ASSERT( state() == BDF_RUNNING ).error( "invalid BDF state" );
364  M_real_time_per_iteration = M_timer.elapsed();
365  return M_real_time_per_iteration;
366  }
367 
369  time_values_map_type timeValues() const
370  {
371  return M_time_values_map;
372  }
373 
375  double start()
376  {
377  M_state = BDF_RUNNING;
378  M_timer.restart();
379  M_iteration = 1;
380  M_time = M_Ti+this->timeStep();
381  // warning: this needs to be fixed wrt restart
382  M_last_iteration_since_order_change = 1;
383  M_order_cur = 1;
384  return M_Ti;
385  }
386 
388  double restart()
389  {
390  M_state = BDF_RUNNING;
391  M_timer.restart();
392  M_time = M_Ti+this->timeStep();
393  M_last_iteration_since_order_change = 1;
394  M_order_cur = 1;
395  ++M_iteration;
396 
397  for ( int i = 2; i<=M_iteration; ++i )
398  {
399  if ( ( ( i - M_last_iteration_since_order_change ) == M_iterations_between_order_change )&&
400  M_order_cur < M_order )
401  {
402  M_last_iteration_since_order_change = i;
403  ++M_order_cur;
404  }
405  }
406 
407  return M_Ti;
408  }
409 
410 
412  bool isFinished() const
413  {
414  bool finished=false;
415 
416  if ( M_Tf>M_Ti )
417  {
418  if ( M_time > M_Tf )
419  {
420  M_state = BDF_STOPPED;
421  finished=true;
422  }
423  }
424 
425  else
426  {
427  if ( M_time < M_Tf )
428  {
429  M_state = BDF_STOPPED;
430  finished=true;
431  }
432  }
433 
434  return finished;
435  }
436 
443  double next() const
444  {
445  FEELPP_ASSERT( state() == BDF_RUNNING ).error( "invalid BDF state" );
446  M_real_time_per_iteration = M_timer.elapsed();
447  M_timer.restart();
448  M_time += M_dt;
449  ++M_iteration;
450 
451  if ( ( ( M_iteration - M_last_iteration_since_order_change ) == M_iterations_between_order_change )&&
452  M_order_cur < M_order )
453  {
454  M_last_iteration_since_order_change = M_iteration;
455  ++M_order_cur;
456  }
457 
458  return M_time;
459  }
460 
461  virtual void shiftRight()
462  {
463  // create and open a character archive for output
464  std::ostringstream ostr;
465 
466  if( M_rankProcInNameOfFiles )
467  ostr << M_name << "-" << M_iteration<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
468  else
469  ostr << M_name << "-" << M_iteration;
470 
471  DVLOG(2) << "[BdfBase::shiftRight] solution name " << ostr.str() << "\n";
472 
473  //M_time_values_map.insert( std::make_pair( M_iteration, this->time() ) );
474  M_time_values_map.push_back( this->time() );
475  //update();
476  }
477  virtual void update()
478  {
479  double tn = M_time_values_map[M_iteration];
480  double tn1 = M_time_values_map[M_iteration-1];
481  double tn2= 0;
482 
483  if ( M_iteration >= 2 )
484  tn2 = M_time_values_map[M_iteration-2];
485 
486  // order 3 and 4 not yet done
487 #if 0
488  double tn3= 0;
489 
490  if ( M_iteration >= 3 )
491  tn3 = M_time_values_map[M_iteration-3];
492 
493 #endif
494 
495  for ( int i = 0; i < BDF_MAX_ORDER; ++i )
496  {
497  if ( i == 0 ) // BDF_ORDER_ONE:
498  {
499  M_alpha[i][ 0 ] = 1./( tn-tn1 ); // Backward Euler
500  M_alpha[i][ 1 ] = 1./( tn-tn1 );
501  M_beta[i][ 0 ] = 1.; // u^{n+1} \approx u^n
502  }
503 
504  if ( i == 1 ) // BDF_ORDER_TWO:
505  {
506  double denom = ( tn-tn1 )*( tn1-tn2 );
507  M_alpha[i][ 0 ] = 0.5*( tn1-2*tn2+tn )/denom;
508  M_alpha[i][ 1 ] = ( tn2-tn )/denom;
509  M_alpha[i][ 2 ] = 0.5/( tn1-tn2 );
510  M_beta[i][ 0 ] = 1+denom;
511  M_beta[i][ 1 ] = denom;
512  }
513  }
514  }
516  BDFState state() const
517  {
518  return M_state;
519  }
520 
522  fs::path path()
523  {
524  return M_path_save;
525  }
526  fs::path restartPath()
527  {
528  return M_restartPath;
529  }
530 
531  bool saveInFile() const
532  {
533  return M_saveInFile;
534  }
535 
536  int saveFreq() const
537  {
538  return M_saveFreq;
539  }
540 
541  bool rankProcInNameOfFiles() const
542  {
543  return M_rankProcInNameOfFiles ;
544  }
545 
546  WorldComm const& worldComm() const
547  {
548  return M_worldComm;
549  }
550 
551  void setOrder( int order )
552  {
553  M_order = order;
554  }
555  void setTimeInitial( double ti )
556  {
557  M_Ti = ti;
558  }
559  void setTimeStep( double dt )
560  {
561  M_dt = dt;
562  }
563  void setTimeFinal( double T )
564  {
565  M_Tf = T;
566  }
567  void setStrategy( int strategy )
568  {
569  M_strategy = ( BDFStragegy )strategy;
570  }
571  void setSteady( bool steady = true )
572  {
573  if ( steady )
574  {
575  M_dt=1e30;
576  M_Tf=1e30;
577  }
578  }
579  void setRestart( bool doRestart )
580  {
581  M_restart=doRestart;
582  }
583  void setRestartPath( std::string s )
584  {
585  M_restartPath=s;
586  }
587  void setRestartAtLastSave( bool b )
588  {
589  M_restartAtLastSave=b;
590  }
591  void setSaveInFile( bool b )
592  {
593  M_saveInFile = b;
594  }
595  void setSaveFreq(int v)
596  {
597  M_saveFreq=v;
598  }
599  void setRankProcInNameOfFiles( bool b )
600  {
601  M_rankProcInNameOfFiles = b;
602  }
603 
604  void print() const
605  {
606  LOG(INFO) << "============================================================\n";
607  LOG(INFO) << "BDF Information\n";
608  LOG(INFO) << " time step : " << this->timeStep() << "\n";
609  LOG(INFO) << "time initial : " << this->timeInitial() << "\n";
610  LOG(INFO) << " time final : " << this->timeFinal() << "\n";
611  LOG(INFO) << " time order : " << this->timeOrder() << "\n";
612  }
613 protected:
615  int M_order;
616  mutable int M_order_cur;
617 
619  std::string M_name;
620 
622  mutable double M_time;
623  mutable int M_iteration;
624  mutable int M_last_iteration_since_order_change;
625  int M_iterations_between_order_change;
626 
628  double M_Ti;
629 
631  double M_Tf;
632 
634  double M_dt;
635 
637  BDFStragegy M_strategy;
638 
640  mutable BDFState M_state;
641 
642  int M_n_restart;
643 
645  bool M_restart;
646 
648  fs::path M_restartPath;
649 
651  bool M_restartAtLastSave;
652 
654  mutable boost::timer M_timer;
655 
657  mutable double M_real_time_per_iteration;
658 
660  std::vector<int> M_time_orders;
661 
663  time_values_map_type M_time_values_map;
664 
666  fs::path M_path_save;
667 
669  std::vector<ublas::vector<double> > M_alpha;
670 
672  std::vector<ublas::vector<double> > M_beta;
673 
675  bool M_saveInFile;
676 
678  int M_saveFreq;
679 
681  bool M_rankProcInNameOfFiles;
682 
684  WorldComm M_worldComm;
685 
686 protected:
687  void init()
688  {
689 
690  for ( int i = 0; i < BDF_MAX_ORDER; ++i )
691  {
692  M_alpha[ i ].resize( i+2 );
693  M_beta[ i ].resize( i+1 );
694  }
695 
696  for ( int i = 0; i < BDF_MAX_ORDER; ++i )
697  {
698  if ( i == 0 ) // BDF_ORDER_ONE:
699  {
700  M_alpha[i][ 0 ] = 1.; // Backward Euler
701  M_alpha[i][ 1 ] = 1.;
702  M_beta[i][ 0 ] = 1.; // u^{n+1} \approx u^n
703  }
704 
705  else if ( i == 1 ) // BDF_ORDER_TWO:
706  {
707  M_alpha[i][ 0 ] = 3. / 2.;
708  M_alpha[i][ 1 ] = 2.;
709  M_alpha[i][ 2 ] = -1. / 2.;
710  M_beta[i][ 0 ] = 2.;
711  M_beta[i][ 1 ] = -1.;
712  }
713 
714  else if ( i == 2 ) // BDF_ORDER_THREE:
715  {
716  M_alpha[i][ 0 ] = 11. / 6.;
717  M_alpha[i][ 1 ] = 3.;
718  M_alpha[i][ 2 ] = -3. / 2.;
719  M_alpha[i][ 3 ] = 1. / 3.;
720  M_beta[i][ 0 ] = 3.;
721  M_beta[i][ 1 ] = -3.;
722  M_beta[i][ 2 ] = 1.;
723  }
724 
725  else if ( i == 3 )
726  {
727  M_alpha[i][ 0 ] = 25. / 12.;
728  M_alpha[i][ 1 ] = 4.;
729  M_alpha[i][ 2 ] = -3.;
730  M_alpha[i][ 3 ] = 4. / 3.;
731  M_alpha[i][ 4 ] = -1. / 4.;
732  M_beta[i][ 0 ] = 4.;
733  M_beta[i][ 1 ] = -6.;
734  M_beta[i][ 2 ] = 4.;
735  M_beta[i][ 3 ] = -1.;
736  }
737  }
738 
739  std::ostringstream ostr;
740  ostr << "bdf_o_" << M_order << "_dt_" << M_dt;
741  //ostr << "bdf/" << M_name << "/o_" << M_order << "/dt_" << M_dt;
742  M_path_save = ostr.str();
743 
744  // if directory does not exist, create it
745  if ( !fs::exists( M_path_save ) && this->saveInFile() )
746  fs::create_directories( M_path_save );
747 
748  if ( M_restart )
749  {
750  //fs::ifstream ifs;
751  fs::path thepath;
752 
753  if ( this->restartPath().empty() ) thepath=this->path()/"metadata";// ifs.open(this->path()/"metadata");
754 
755  else thepath=this->restartPath()/this->path()/"metadata"; //ifs.open(this->restartPath()/this->path()/"metadata");
756 
757  // read the saved bdf data
758  if ( fs::exists( thepath/* this->restartPath() / this->path() / "metadata" )*/ ) )
759  {
760  DVLOG(2) << "[Bdf] loading metadata from " << M_path_save.string() << "\n";
761 
762  //fs::ifstream ifs( this->restartPath() / this->path() / "metadata")
763  fs::ifstream ifs( thepath );
764 
765 
766  boost::archive::text_iarchive ia( ifs );
767  ia >> BOOST_SERIALIZATION_NVP( *this );
768  DVLOG(2) << "[Bdf::init()] metadata loaded\n";
769  //BdfBaseMetadata bdfloader( *this );
770  //bdfloader.load();
771 
772  // modify Ti with last saved time
773  if (this->doRestartAtLastSave()) M_Ti = M_time_values_map.back();
774 
775  M_iteration = 0;
776  // look for M_ti in the time values
777  bool found = false;
778  BOOST_FOREACH( auto time, M_time_values_map )
779  {
780  if ( math::abs( time-M_Ti ) < 1e-10 )
781  {
782  //M_iteration = time.first;
783  //std::cout << "time found " << time << std::endl;
784  found = true;
785  break;
786  }
787 
788  ++M_iteration;
789  }
790  //std::cout << "M_iteration " << M_iteration << std::endl;
791 
792  if ( !found )
793  {
794  DVLOG(2) << "[Bdf] intial time " << M_Ti << " not found\n";
795  M_Ti = 0.0;
796  M_iteration = 0;
797  M_time_values_map.clear();
798  return;
799  }
800 
801  else
802  {
803  if (this->saveFreq()==1)
804  {
805  M_time_values_map.resize( M_iteration+1 );
806  }
807  else
808  {
809  int nItBack = M_iteration % this->saveFreq();
810  M_iteration-=nItBack;
811  M_time_values_map.resize( M_iteration+1 );
812  M_Ti = M_time_values_map.back();
813  }
814  }
815 
816  DVLOG(2) << "[Bdf] initial time is Ti=" << M_Ti << "\n";
817 
818  DVLOG(2) << "[Bdf::init()] file index: " << M_iteration << "\n";
819  }
820 
821  else
822  {
823  M_Ti = 0.0;
824  M_time_values_map.clear();
825  }
826  }
827 
828  } // init
829 };
830 class BdfBaseMetadata
831 {
832 public:
833  BdfBaseMetadata( BdfBase& bdf )
834  :
835  M_bdf( bdf )
836  {
837  }
838 
839  void load()
840  {
841  fs::ifstream ifs;
842 
843  if ( M_bdf.restartPath().empty() ) ifs.open( M_bdf.path()/"metadata" );
844 
845  else ifs.open( M_bdf.restartPath()/M_bdf.path()/"metadata" );
846 
847  //fs::ifstream ifs( M_bdf.path() / "metadata");
848 
849  boost::archive::text_iarchive ia( ifs );
850  ia >> BOOST_SERIALIZATION_NVP( M_bdf );
851  DVLOG(2) << "[Bdf::init()] metadata loaded\n";
852  }
853 
854  void save()
855  {
856  if ( !M_bdf.saveInFile() || M_bdf.worldComm().globalRank()!=M_bdf.worldComm().masterRank() ) return;
857 
858  fs::ofstream ofs( M_bdf.path() / "metadata" );
859 
860 
861  boost::archive::text_oarchive oa( ofs );
862  oa << BOOST_SERIALIZATION_NVP( ( BdfBase const& )M_bdf );
863  DVLOG(2) << "[Bdf::init()] metadata saved\n";
864  }
865 
866 private:
867  BdfBase& M_bdf;
868 };
905 template<typename SpaceType>
906 class Bdf : public BdfBase
907 {
908  friend class boost::serialization::access;
909  typedef BdfBase super;
910 public:
911  typedef Bdf<SpaceType> bdf_type;
912  typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
913  typedef SpaceType space_type;
914  typedef boost::shared_ptr<space_type> space_ptrtype;
915  typedef typename space_type::element_type element_type;
916  typedef typename space_type::return_type return_type;
917  typedef typename element_type::value_type value_type;
918  //typedef boost::numeric::ublas::vector< element_type > unknowns_type;
919  typedef boost::shared_ptr<element_type> unknown_type;
920  typedef std::vector< unknown_type > unknowns_type;
921  typedef typename node<value_type>::type node_type;
922 
923  typedef typename super::time_iterator time_iterator;
924 #if 0
925  BOOST_PARAMETER_CONSTRUCTOR(
926  Bdf, ( BdfBase ), tag,
927  ( required ( final_time,* ) )
928  ( optional
929  ( name,* )
930  ( order,* )
931  ( initial_time,* )
932  ( time_step,* )
933  ( strategy,* ) ) )
934 #endif
935 
941  Bdf( po::variables_map const& vm, space_ptrtype const& space, std::string const& name, std::string const& prefix="" );
942 
948  Bdf( space_ptrtype const& space, std::string const& name );
949 
951  Bdf( Bdf const& b )
952  :
953  super( b ),
954  M_space( b.M_space ),
955  M_unknowns( b.M_unknowns )
956  {}
957 
958  ~Bdf();
959 
961  bdf_ptrtype deepCopy() const
962  {
963  auto b = bdf_ptrtype( new bdf_type( *this ) );
964 
965  for ( auto it = b->M_unknowns.begin(), en = b->M_unknowns.end(); it != en; ++ it )
966  {
967  *it = unknown_type( new element_type( M_space ) );
968  }
969 
970  return b;
971  }
972 
977  void initialize( element_type const& u0 );
978 
983  void initialize( unknowns_type const& uv0 );
984 
988  double start();
989  double start( element_type const& u0 );
990  double start( unknowns_type const& uv0 );
991 
995  double restart();
996 
1002  template<typename container_type>
1003  void shiftRight( typename space_type::template Element<value_type, container_type> const& u_curr );
1004 
1007  element_type polyDeriv() const;
1008 
1011  element_type poly() const;
1012 
1014  unknowns_type const& unknowns() const;
1015 
1017  element_type& unknown( int i );
1018 
1019  template<typename container_type>
1020  void setUnknown( int i, typename space_type::template Element<value_type, container_type> const& e )
1021  {
1022  M_unknowns[i]->assign( e );
1023  }
1024 
1025  void showMe( std::ostream& __out = std::cout ) const;
1026 
1027  void loadCurrent();
1028 private:
1029  void init();
1030 
1031 
1032  void saveCurrent();
1033 
1035  template<class Archive>
1036  void serialize( Archive & ar, const unsigned int version )
1037  {
1038  DVLOG(2) << "[BDF::serialize] saving/loading archive\n";
1039  ar & boost::serialization::base_object<BdfBase>( *this );
1040  }
1041 
1042 private:
1043 
1045  space_ptrtype M_space;
1046 
1048  unknowns_type M_unknowns;
1049 };
1050 
1051 template <typename SpaceType>
1052 Bdf<SpaceType>::Bdf( po::variables_map const& vm,
1053  space_ptrtype const& __space,
1054  std::string const& name,
1055  std::string const& prefix )
1056  :
1057  super( vm, name, prefix, __space->worldComm() ),
1058  M_space( __space )
1059 {
1060  M_unknowns.resize( BDF_MAX_ORDER );
1061 
1062  for ( uint8_type __i = 0; __i < ( uint8_type )BDF_MAX_ORDER; ++__i )
1063  {
1064  M_unknowns[__i] = unknown_type( new element_type( M_space ) );
1065  M_unknowns[__i]->zero();
1066  }
1067 
1068 }
1069 
1070 template <typename SpaceType>
1071 Bdf<SpaceType>::Bdf( space_ptrtype const& __space,
1072  std::string const& name )
1073  :
1074  super( name, __space->worldComm() ),
1075  M_space( __space )
1076 {
1077  M_unknowns.resize( BDF_MAX_ORDER );
1078 
1079  for ( uint8_type __i = 0; __i < ( uint8_type )BDF_MAX_ORDER; ++__i )
1080  {
1081  M_unknowns[__i] = unknown_type( new element_type( M_space ) );
1082  M_unknowns[__i]->zero();
1083  }
1084 
1085 }
1086 
1087 template <typename SpaceType>
1088 void
1090 {
1091  super::init();
1092 
1093  if ( this->isRestart() )
1094  {
1095  for ( int p = 0; p < std::min( M_order, M_iteration+1 ); ++p )
1096  {
1097  // create and open a character archive for output
1098  std::ostringstream ostr;
1099 
1100  if( M_rankProcInNameOfFiles )
1101  ostr << M_name << "-" << M_iteration-p<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
1102  else
1103  ostr << M_name << "-" << M_iteration-p;
1104 
1105  DVLOG(2) << "[Bdf::init()] load file: " << ostr.str() << "\n";
1106 
1107  fs::ifstream ifs;
1108 
1109  if ( this->restartPath().empty() ) ifs.open( this->path()/ostr.str() );
1110 
1111  else ifs.open( this->restartPath()/this->path()/ostr.str() );
1112 
1113  //fs::ifstream ifs (this->restartPath() / this->path() / ostr.str(), std::ios::binary);
1114 
1115  // load data from archive
1116  boost::archive::binary_iarchive ia( ifs );
1117  ia >> *M_unknowns[p];
1118  }
1119  }
1120 }
1121 
1122 
1123 template <typename SpaceType>
1124 Bdf<SpaceType>::~Bdf()
1125 {}
1126 
1127 
1128 template <typename SpaceType>
1129 void
1130 Bdf<SpaceType>::initialize( element_type const& u0 )
1131 {
1132  M_time_values_map.clear();
1133  std::ostringstream ostr;
1134 
1135  if( M_rankProcInNameOfFiles )
1136  ostr << M_name << "-" << 0<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
1137  else
1138  ostr << M_name << "-" << 0;
1139  //M_time_values_map.insert( std::make_pair( 0, boost::make_tuple( 0, ostr.str() ) ) );
1140  //M_time_values_map.push_back( 0 );
1141  M_time_values_map.push_back( M_Ti );
1142  std::for_each( M_unknowns.begin(), M_unknowns.end(), *boost::lambda::_1 = u0 );
1143  this->saveCurrent();
1144 }
1145 
1146 template <typename SpaceType>
1147 void
1148 Bdf<SpaceType>::initialize( unknowns_type const& uv0 )
1149 {
1150  M_time_values_map.clear();
1151  std::ostringstream ostr;
1152 
1153  if( M_rankProcInNameOfFiles )
1154  ostr << M_name << "-" << 0<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
1155  else
1156  ostr << M_name << "-" << 0;
1157  //M_time_values_map.insert( std::make_pair( 0, boost::make_tuple( 0, ostr.str() ) ) );
1158  //M_time_values_map.push_back( 0);
1159  M_time_values_map.push_back( M_Ti );
1160 
1161  std::copy( uv0.begin(), uv0.end(), M_unknowns.begin() );
1162  this->saveCurrent();
1163 }
1164 
1165 template <typename SpaceType>
1166 double
1168 {
1169  this->init();
1170 
1171  return super::start();
1172 }
1173 
1174 template <typename SpaceType>
1175 double
1176 Bdf<SpaceType>::start( element_type const& u0 )
1177 {
1178  this->init();
1179  this->initialize( u0 );
1180  auto res = super::start();
1181  return res;
1182 }
1183 
1184 template <typename SpaceType>
1185 double
1186 Bdf<SpaceType>::start( unknowns_type const& uv0 )
1187 {
1188  this->init();
1189  this->initialize( uv0 );
1190  auto res = super::start();
1191  return res;
1192 }
1193 
1194 template <typename SpaceType>
1195 double
1197 {
1198  this->init();
1199 
1200  return super::restart();
1201 }
1202 
1203 template <typename SpaceType>
1204 const
1205 typename Bdf<SpaceType>::unknowns_type&
1207 {
1208  return M_unknowns;
1209 }
1210 
1211 template <typename SpaceType>
1212 typename Bdf<SpaceType>::element_type&
1213 Bdf<SpaceType>::unknown( int i )
1214 {
1215  DVLOG(2) << "[Bdf::unknown] id: " << i << " l2norm = " << M_unknowns[i]->l2Norm() << "\n";
1216  return *M_unknowns[i];
1217 }
1218 
1219 
1220 template <typename SpaceType>
1221 void
1222 Bdf<SpaceType>::saveCurrent()
1223 {
1224  if (!this->saveInFile()) return;
1225 
1226  bool doSave=false;
1227  for ( uint8_type i = 0; i < this->timeOrder() && !doSave; ++i )
1228  {
1229  int iterTranslate = M_iteration + this->timeOrder()-(i+1);
1230  if (iterTranslate % this->saveFreq()==0) doSave=true;
1231  }
1232 
1233  if (!doSave) return;
1234 
1235  BdfBaseMetadata bdfsaver( *this );
1236  bdfsaver.save();
1237 
1238  {
1239  std::ostringstream ostr;
1240 
1241  if( M_rankProcInNameOfFiles )
1242  ostr << M_name << "-" << M_iteration<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
1243  else
1244  ostr << M_name << "-" << M_iteration;
1245 
1246  fs::ofstream ofs( M_path_save / ostr.str() );
1247 
1248 
1249  // load data from archive
1250  boost::archive::binary_oarchive oa( ofs );
1251  oa << *M_unknowns[0];
1252  }
1253 }
1254 
1255 template <typename SpaceType>
1256 void
1257 Bdf<SpaceType>::loadCurrent()
1258 {
1259  BdfBaseMetadata bdfsaver( *this );
1260  bdfsaver.save();
1261 
1262  {
1263  std::ostringstream ostr;
1264 
1265  if( M_rankProcInNameOfFiles )
1266  ostr << M_name << "-" << M_iteration<<"-proc"<<Environment::worldComm().globalRank()<<"on"<<Environment::worldComm().globalSize();
1267  else
1268  ostr << M_name << "-" << M_iteration;
1269 
1270  fs::ifstream ifs( M_path_save / ostr.str() );
1271 
1272  // load data from archive
1273  boost::archive::binary_iarchive ia( ifs );
1274  ia >> *M_unknowns[0];
1275  }
1276 }
1277 
1278 template <typename SpaceType>
1279 template<typename container_type>
1280 void
1281 Bdf<SpaceType>::shiftRight( typename space_type::template Element<value_type, container_type> const& __new_unk )
1282 {
1283  DVLOG(2) << "shiftRight: inserting time " << this->time() << "s\n";
1284  super::shiftRight();
1285 
1286  // shift all previously stored bdf data
1287  using namespace boost::lambda;
1288  typename unknowns_type::reverse_iterator __it = boost::next( M_unknowns.rbegin() );
1289  std::for_each( M_unknowns.rbegin(), boost::prior( M_unknowns.rend() ),
1290  ( *lambda::_1 = *( *lambda::var( __it ) ), ++lambda::var( __it ) ) );
1291  // u(t^{n}) coefficient is in M_unknowns[0]
1292  *M_unknowns[0] = __new_unk;
1293  int i = 0;
1294  BOOST_FOREACH( boost::shared_ptr<element_type>& t, M_unknowns )
1295  {
1296  DVLOG(2) << "[Bdf::shiftright] id: " << i << " l2norm = " << t->l2Norm() << "\n";
1297  ++i;
1298  }
1299 
1300  // save newly stored bdf data
1301  this->saveCurrent();
1302 }
1303 
1304 
1305 template <typename SpaceType>
1306 typename Bdf<SpaceType>::element_type
1308 {
1309  element_type __t( M_space );
1310  __t.zero();
1311 
1312  FEELPP_ASSERT( __t.size() == M_space->nDof() )( __t.size() )( M_space->nDof() ).error( "invalid space element size" );
1313  FEELPP_ASSERT( __t.size() == M_unknowns[0]->size() )( __t.size() )( M_unknowns[0]->size() ).error( "invalid space element size" );
1314 
1315  for ( uint8_type i = 0; i < this->timeOrder(); ++i )
1316  __t.add( this->polyDerivCoefficient( i+1 ), *M_unknowns[i] );
1317 
1318  return __t;
1319 }
1320 
1321 template <typename SpaceType>
1322 typename Bdf<SpaceType>::element_type
1324 {
1325  element_type __t( M_space );
1326  __t.zero();
1327 
1328  FEELPP_ASSERT( __t.size() == M_space->nDof() )( __t.size() )( M_space->nDof() ).error( "invalid space element size" );
1329  FEELPP_ASSERT( __t.size() == M_unknowns[0]->size() )( __t.size() )( M_unknowns[0]->size() ).error( "invalid space element size" );
1330 
1331  for ( uint8_type i = 0; i < this->timeOrder(); ++i )
1332  __t.add( this->polyCoefficient( i ), *M_unknowns[ i ] );
1333 
1334  return __t;
1335 }
1336 
1337 
1338 
1339 BOOST_PARAMETER_FUNCTION(
1340  ( boost::shared_ptr<Bdf<typename meta::remove_all<typename parameter::binding<Args, tag::space>::type>::type::element_type> > ),
1341  bdf, tag,
1342  ( required
1343  ( space,*( boost::is_convertible<mpl::_,boost::shared_ptr<Feel::FunctionSpaceBase> > ) ) )
1344  ( optional
1345  ( vm,*, Environment::vm() )
1346  ( prefix,*,"" )
1347  ( name,*,"bdf" )
1348  ( order,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.order" )].template as<int>() )
1349  ( initial_time,*( boost::is_floating_point<mpl::_> ),vm[prefixvm( prefix,"bdf.time-initial" )].template as<double>() )
1350  ( final_time,*( boost::is_floating_point<mpl::_> ),vm[prefixvm( prefix,"bdf.time-final" )].template as<double>() )
1351  ( time_step,*( boost::is_floating_point<mpl::_> ),vm[prefixvm( prefix,"bdf.time-step" )].template as<double>() )
1352  ( strategy,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.strategy" )].template as<int>() )
1353  ( steady,*( bool ),vm[prefixvm( prefix,"bdf.steady" )].template as<bool>() )
1354  ( restart,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.restart" )].template as<bool>() )
1355  ( restart_path,*,vm[prefixvm( prefix,"bdf.restart.path" )].template as<std::string>() )
1356  ( restart_at_last_save,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.restart.at-last-save" )].template as<bool>() )
1357  ( save,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.save" )].template as<bool>() )
1358  ( freq,*(boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.save.freq" )].template as<int>() )
1359  ( rank_proc_in_files_name,*( boost::is_integral<mpl::_> ),vm[prefixvm( prefix,"bdf.rank-proc-in-files-name" )].template as<bool>() )
1360  ) )
1361 {
1362  typedef typename meta::remove_all<space_type>::type::element_type _space_type;
1363  auto thebdf = boost::shared_ptr<Bdf<_space_type> >( new Bdf<_space_type>( vm,space,name,prefix ) );
1364  thebdf->setTimeInitial( initial_time );
1365  thebdf->setTimeFinal( final_time );
1366  thebdf->setTimeStep( time_step );
1367  thebdf->setOrder( order );
1368  thebdf->setSteady( steady );
1369  thebdf->setStrategy( strategy );
1370  thebdf->setRestart( restart );
1371  thebdf->setRestartPath( restart_path );
1372  thebdf->setRestartAtLastSave( restart_at_last_save );
1373  thebdf->setSaveInFile( save );
1374  thebdf->setSaveFreq( freq );
1375  thebdf->setRankProcInNameOfFiles( rank_proc_in_files_name );
1376  return thebdf;
1377 }
1378 
1379 
1380 }
1381 #endif

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4