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
aitken.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): Goncalo Pena <goncalo.pena@epfl.ch>
6  Date: 15-07-2008
7 
8  Copyright (C) 2007-2008 Universite Joseph Fourier (Grenoble I)
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 */
24 
25 
26 #ifndef __AitkenExtrapolation
27 #define __AitkenExtrapolation 1
28 
29 #include <string>
30 #include <vector>
31 #include <map>
32 
33 #include <boost/assign/list_of.hpp>
34 #include <boost/assert.hpp>
35 #include <boost/timer.hpp>
36 #include <boost/tuple/tuple.hpp>
37 #include <feel/feeldiscr/functionspace.hpp>
38 #include <feel/feelalg/vector.hpp>
39 
40 
41 namespace Feel
42 {
82 template< typename fs_type >
83 class Aitken
84 {
85 
86 public:
87 
88  typedef Aitken<fs_type> self_type;
89 
90  typedef fs_type functionspace_type;
91  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
92 
93  typedef typename functionspace_type::element_type element_type;
94 
95  typedef typename functionspace_type::template Element<typename functionspace_type::value_type,
97 
104  typedef std::map<std::string, double> convergence_iteration_type;
105  typedef std::map<int, convergence_iteration_type> convergence_type;
106 
110  Aitken( functionspace_ptrtype _Xh,
111  AitkenType _aitkenType = AITKEN_STANDARD,
112  double _failsafeParameter = 1.0,
113  double _tol=1.0e-6,
114  double _minParam = 1e-4 )
115  :
116  M_Xh( _Xh ),
117  M_failsafeParameter( _failsafeParameter ),
118  M_previousParameter( _failsafeParameter ),
119  M_maxParameter( _failsafeParameter ),
120  M_minParameter( _minParam ),
121  M_previousResidual( M_Xh, "previous residual" ),
122  M_previousElement( M_Xh, "previous element" ),
123  M_currentResidual( M_Xh, "current residual" ),
124  M_currentElement( M_Xh, "current element" ),
125  M_cptIteration( 1 ),
126  M_aitkenType( _aitkenType ),
127  M_tolerance( _tol ),
128  M_residualConvergence( 1. ),
129  M_hasConverged( false )
130  {
131  }
132 
136  Aitken( Aitken const& tc )
137  :
138  M_Xh( tc.M_Xh ),
139  M_failsafeParameter( tc.M_failsafeParameter ),
140  M_previousParameter( tc.M_previousParameter ),
141  M_maxParameter( tc.M_maxParameter ),
142  M_minParameter( tc.M_minParameter ),
143  M_previousResidual( tc.M_previousResidual ),
144  M_previousElement( tc.M_previousElement ),
145  M_currentResidual( tc.M_currentResidual ),
146  M_currentElement( tc.M_currentElement ),
147  M_cptIteration( tc.M_cptIteration ),
148  M_aitkenType( tc.M_aitkenType ),
149  M_tolerance( tc.M_tolerance ),
150  M_residualConvergence( tc.M_residualConvergence ),
151  M_hasConverged( tc.M_hasConverged )
152  {
153  }
154 
158  ~Aitken() {}
159 
163  BOOST_PARAMETER_MEMBER_FUNCTION(
164  ( void ),
165  initialize,
166  tag,
167  ( required
168  ( residual, */*(element_type const&)*/ )
169  ( currentElt,*/*(element_type const& )*/ ) ) )
170  {
171  initializeimpl( residual,currentElt );
172  }
173 
177  BOOST_PARAMETER_MEMBER_FUNCTION(
178  ( element_type ),
179  apply,
180  tag,
181  ( required
182  ( residual, * )
183  ( currentElt,* )
184  ) //required
185  ( optional
186  ( forceRelaxation, ( bool ), false )
187  )//optional
188  )
189  {
190  element_type newElt( M_Xh );
191  applyimpl( newElt,residual,currentElt,forceRelaxation );
192  return newElt;
193  }
194 
198  template< typename eltType >
199  element_type operator()( element_type const& residual, eltType const& elem,bool _forceRelax=false )
200  {
201  element_type newElt( M_Xh );
202  applyimpl( newElt,residual,elem,_forceRelax );
203  return newElt;
204  }
205 
209  BOOST_PARAMETER_MEMBER_FUNCTION(
210  ( void ),
211  apply2,
212  tag,
213  ( required
214  ( newElt, * )
215  ( residual, * )
216  ( currentElt,* )
217  ) //required
218  ( optional
219  ( forceRelaxation, ( bool ), false )
220  ) //optional
221  )
222  {
223  applyimpl( newElt,residual,currentElt,forceRelaxation );
224  }
225 
227  void setType( AitkenType t )
228  {
229  M_aitkenType = t;
230  }
231 
233  AitkenType type( ) const
234  {
235  return M_aitkenType;
236  }
237 
241  void shiftRight();
242 
246  self_type & operator++();
247 
251  void restart();
252 
256  double theta()
257  {
258  return M_previousParameter;
259  }
260 
264  void setTheta(double v)
265  {
266  M_previousParameter=v;
267  }
268 
273  uint nIterations()
274  {
275  return M_cptIteration;
276  }
277 
278  bool isFinished()
279  {
280  return M_hasConverged;
281  }
282 
283  double residualNorm()
284  {
285  return M_residualConvergence;
286  }
287 
288  void printInfo();
289 
294  void saveConvergenceHistory( std::string const& fname ) const;
295 
296  void forceConvergence( bool b )
297  {
298  M_hasConverged=b;
299  }
300 
301 
305  void computeResidualNorm();
306 
310  void setElement( element_type const& residual, element_type const& elem );
311 
315  void setElement( element_type const& residual, element_range_type const& elem );
316 
320  void calculateParameter();
321 
325  //void relaxationStep( element_type& new_elem );
326  template< typename eltType >
327  void relaxationStep( eltType& new_elem );
328 
330  convergence_type const& convergenceHistory() const
331  {
332  return M_convergence;
333  }
334 
335 private:
339  void initializeimpl( element_type const& residual, element_type const& elem );
340 
344  void initializeimpl( element_type const& residual, element_range_type const& elem );
345 
349  void applyimpl( element_type & new_elem,element_type const& residual, element_type const& elem,bool _forceRelax=false );
350 
354  void applyimpl( element_range_type new_elem,element_type const& residual, element_range_type const& elem,bool _forceRelax=false );
355 
359  void calculateParameter( mpl::int_<AITKEN_STANDARD> );
360 
364  void calculateParameter( mpl::int_<AITKEN_METHOD_1> );
365 
366 
367 private:
368 
372  functionspace_ptrtype M_Xh;
373 
374  double M_failsafeParameter, M_previousParameter, M_maxParameter;
375  double M_minParameter;
376 
377  element_type M_previousResidual, M_previousElement, M_currentResidual, M_currentElement;
378 
379  uint M_cptIteration;
380  AitkenType M_aitkenType;
381  double M_tolerance;
382  double M_residualConvergence;
383  bool M_hasConverged;
384  convergence_type M_convergence;
385 
386 
387  mutable boost::timer M_timer;
388 };
389 
390 
391 //-----------------------------------------------------------------------------------------//
392 
393 template< typename fs_type >
394 void
395 Aitken<fs_type>::initializeimpl( element_type const& residual, element_type const& elem )
396 {
397  M_previousResidual = residual;
398  M_previousElement = elem;
399 }
400 
401 //-----------------------------------------------------------------------------------------//
402 
403 template< typename fs_type >
404 void
405 Aitken<fs_type>::initializeimpl( element_type const& residual, element_range_type const& elem )
406 {
407  M_previousResidual = residual;
408  M_previousElement.zero();
409  M_previousElement.add( 1.,elem );
410  /*M_previousElement = vf::project(M_previousElement.functionSpace(),
411  elements(M_previousElement.mesh()),
412  vf::idv(elem) );*/
413 }
414 
415 //-----------------------------------------------------------------------------------------//
416 
417 template< typename fs_type >
418 void
420 {
421  auto oldEltL2Norm = M_previousElement.l2Norm();
422 
423  if ( oldEltL2Norm > 1e-8 )
424  M_residualConvergence = M_currentResidual.l2Norm()/oldEltL2Norm;
425 
426  else
427  M_residualConvergence = M_currentResidual.l2Norm();
428 
429  if ( M_residualConvergence < M_tolerance ) M_hasConverged=true;
430 
431 }
432 
433 //-----------------------------------------------------------------------------------------//
434 
435 template< typename fs_type >
436 void
437 Aitken<fs_type>::setElement( element_type const& residual, element_type const& elem )
438 {
439  M_currentResidual = residual;
440  M_currentElement = elem;
441 }
442 
443 //-----------------------------------------------------------------------------------------//
444 
445 template< typename fs_type >
446 void
447 Aitken<fs_type>::setElement( element_type const& residual, element_range_type const& elem )
448 {
449  M_currentResidual = residual;
450  M_currentElement.zero();
451  M_currentElement.add( 1.,elem );
452  /*M_currentElement = vf::project(M_currentElement.functionSpace(),
453  elements(M_currentElement.mesh()),
454  vf::idv(elem) );*/
455 }
456 
457 //-----------------------------------------------------------------------------------------//
458 
459 template< typename fs_type >
460 void
461 Aitken<fs_type>::applyimpl( element_type & new_elem,element_type const& residual, element_type const& elem,bool _forceRelax )
462 {
463 
464  setElement( residual,elem );
465 
466  if ( M_cptIteration>=2 )
467  {
468  computeResidualNorm();
469 
470  if ( !M_hasConverged || _forceRelax )
471  {
472  calculateParameter();
473  }
474  }
475 
476  if ( !M_hasConverged || _forceRelax )
477  relaxationStep( new_elem );
478 
479 }
480 
481 //-----------------------------------------------------------------------------------------//
482 
483 template< typename fs_type >
484 void
485 Aitken<fs_type>::applyimpl( element_range_type new_elem,element_type const& residual, element_range_type const& elem,bool _forceRelax )
486 {
487 
488  setElement( residual,elem );
489 
490  if ( M_cptIteration>=2 )
491  {
492  computeResidualNorm();
493 
494  if ( !M_hasConverged || _forceRelax )
495  {
496  calculateParameter();
497  }
498  }
499 
500  if ( !M_hasConverged || _forceRelax )
501  relaxationStep( new_elem );
502 
503 }
504 
505 //-----------------------------------------------------------------------------------------//
506 template< typename fs_type >
507 template< typename eltType >
508 void
510 {
511 #if 0
512  new_elem = M_currentResidual;
513  //new_elem.scale( -M_previousParameter );
514  new_elem.scale( -( 1-M_previousParameter ) );
515 
516  //new_elem += M_currentElement;
517  new_elem.add( 1.,M_currentElement );
518 #else
519  new_elem = M_previousElement;
520  new_elem.add(M_previousParameter,M_currentResidual);
521 
522  M_currentElement = new_elem;
523 #endif
524 }
525 
526 //-----------------------------------------------------------------------------------------//
527 
528 template< typename fs_type >
529 void
531 {
532  // store convergence history
533  double rel_residual = 1;
534 
535  if ( M_cptIteration > 1 )
536  rel_residual = M_currentResidual.l2Norm()/M_convergence.find( 1 )->second.find( "absolute_residual" )->second;
537 
539  boost::assign::map_list_of
540  ( "absolute_residual", M_currentResidual.l2Norm() )
541  ( "relative_residual", rel_residual )
542  ( "time", M_timer.elapsed() )
543  ( "relaxation_parameter", M_previousParameter );
544  M_convergence.insert( std::make_pair( M_cptIteration, cit ) );
545  M_timer.restart();
546 
547  M_previousResidual = M_currentResidual;
548  M_previousElement = M_currentElement;
549 
550 
551  ++M_cptIteration;
552 }
553 
554 //-----------------------------------------------------------------------------------------//
555 
556 template< typename fs_type >
559 {
560  shiftRight();
561 
562  return *this;
563 }
564 
565 //-----------------------------------------------------------------------------------------//
566 
567 template< typename fs_type >
568 void
570 {
571  if ( !M_convergence.empty() )
572  {
573  std::string entry = "relaxation_parameter";
574  auto it = std::max_element( M_convergence.begin(), M_convergence.end(),
575  [entry] ( std::pair<int, convergence_iteration_type> const& x,
576  std::pair<int, convergence_iteration_type> const& y )
577  {
578  return x.second.find( entry )->second < y.second.find( entry )->second;
579  } );
580 
581  if ( it != M_convergence.end() )
582  M_maxParameter = std::max( M_maxParameter,it->second.find( entry )->second );
583 
584  M_previousParameter = M_maxParameter;
585  }
586 
587  M_cptIteration=1;
588  M_hasConverged=false;
589  M_convergence.clear();
590  M_timer.restart();
591 }
592 
593 //-----------------------------------------------------------------------------------------//
594 
595 template< typename fs_type >
596 void
598 {
599  std::cout << "[Aitken] iteration : "<< M_cptIteration
600  <<" theta=" << M_previousParameter
601  <<" residualNorm : " << M_residualConvergence
602  << "\n";
603 
604  LOG(INFO) << "[Aitken] iteration : "<< M_cptIteration
605  <<" theta=" << M_previousParameter
606  <<" residualNorm : " << M_residualConvergence
607  << "\n";
608 }
609 
610 template< typename fs_type >
611 void
612 Aitken<fs_type>::saveConvergenceHistory( std::string const& fname ) const
613 {
614  std::ofstream ofs( fname.c_str() );
615 
616  if ( ofs )
617  {
618  ofs.setf( std::ios::scientific );
619  BOOST_FOREACH( auto cit, M_convergence )
620  {
621  // iteration
622  ofs << std::setw( 6 ) << cit.first << " ";
623  BOOST_FOREACH( auto it, cit.second )
624  {
625  ofs << std::setw( 20 ) << std::setprecision( 16 ) << it.second << " ";
626  }
627  ofs << "\n";
628  }
629  }
630 
631  else
632  {
633  std::cerr << "[Aitken] convergence history filename " << fname << " could not be opened"
634  << std::endl;
635  }
636 }
637 //-----------------------------------------------------------------------------------------//
638 
639 template< typename fs_type >
640 void
642 {
643  if ( M_aitkenType==AITKEN_STANDARD )
644  calculateParameter( mpl::int_<AITKEN_STANDARD>() );
645 
646  if ( M_aitkenType==AITKEN_METHOD_1 )
647  calculateParameter( mpl::int_<AITKEN_METHOD_1>() );
648 
649  if ( M_aitkenType==FIXED_RELAXATION_METHOD )
650  M_previousParameter = M_failsafeParameter;
651 }
652 
653 //-----------------------------------------------------------------------------------------//
654 
655 template< typename fs_type >
656 void
657 Aitken<fs_type>::calculateParameter( mpl::int_<AITKEN_STANDARD> )
658 {
659 
660  element_type aux( M_Xh, "aux" );
661 
662  aux = M_currentResidual;
663  aux -= M_previousResidual;
664 
665  double scalar = inner_product( aux, aux );
666 
667  aux.scale( 1.0/scalar );
668 
669  element_type aux2( M_Xh, "aux2" );
670 
671  aux2 = M_currentElement;
672  aux2 -= M_previousElement;
673 
674  scalar = inner_product( aux2, aux );
675 
676  if ( scalar > 1 )
677  scalar = /*M_previousParameter*/M_failsafeParameter;
678 
679  if ( scalar < 0 )
680  scalar = /*M_previousParameter*/M_failsafeParameter;
681 
682  M_previousParameter = 1-scalar;
683 }
684 
685 //-----------------------------------------------------------------------------------------//
686 
687 template< typename fs_type >
688 void
689 Aitken<fs_type>::calculateParameter( mpl::int_<AITKEN_METHOD_1> )
690 {
691 
692  element_type aux( M_Xh, "aux" );
693 
694  aux = M_currentResidual;
695  aux -= M_previousResidual;
696 
697  double scalar = inner_product( aux, aux );
698 
699  aux.scale( 1.0/scalar );
700 
701  /*element_type aux2( M_Xh, "aux2");
702 
703  aux2 = M_currentElement;
704  aux2 -= M_previousElement;*/
705 
706  scalar = inner_product( M_previousResidual , aux );
707  scalar = -M_previousParameter*scalar;
708 
709 #if 1
710  if ( scalar > 1 )
711  scalar = /*M_previousParameter;*/M_failsafeParameter;
712 
713  if ( scalar < /*0*/M_minParameter )
714  scalar = /*M_previousParameter;*/M_failsafeParameter;
715 #endif
716  M_previousParameter = scalar;
717 
718 }
719 
720 
721 //-----------------------------------------------------------------------------------------//
722 //-----------------------------------------------------------------------------------------//
723 //-----------------------------------------------------------------------------------------//
724 
725 template <typename fs_type >
726 void
727 operator++( boost::shared_ptr<Aitken<fs_type> > & aitk )
728 {
729  aitk->shiftRight();
730 }
731 
732 //-----------------------------------------------------------------------------------------//
733 //-----------------------------------------------------------------------------------------//
734 //-----------------------------------------------------------------------------------------//
735 
736 template<typename SpaceType>
737 boost::shared_ptr<Aitken<SpaceType> >
738 aitkenNew( boost::shared_ptr<SpaceType> const& _space,
739  AitkenType _type,
740  double _init_theta,
741  double _tol,
742  double _minParam )
743 {
744 
745  boost::shared_ptr<Aitken<SpaceType> > Aitk( new Aitken<SpaceType>( _space,_type,_init_theta,_tol,_minParam ) );
746 
747  return Aitk;
748 }
749 
750 
751 
752 template<typename Args>
753 struct compute_aitken_return
754 {
755  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::space>::type>::type::element_type space_type;
756 
757  typedef Aitken<space_type> type;
758  typedef boost::shared_ptr<type> ptrtype;
759 };
760 
761 
762 BOOST_PARAMETER_FUNCTION(
763  ( typename compute_aitken_return<Args>::type ),
764  aitken,
765  tag,
766  ( required
767  ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
768  )//required
769  ( optional
770  ( type, ( AitkenType ), AITKEN_STANDARD )
771  ( initial_theta, *( boost::is_arithmetic<mpl::_> ), 1.0 )
772  ( min_theta, *( boost::is_arithmetic<mpl::_> ), 1e-4 )
773  ( tolerance, *( boost::is_arithmetic<mpl::_> ), 1.0e-6 )
774  )//optional
775 )
776 {
777  return *aitkenNew( space,type,initial_theta,tolerance,min_theta );
778 }
779 
780 BOOST_PARAMETER_FUNCTION(
781  ( typename compute_aitken_return<Args>::ptrtype ),
782  aitkenPtr,
783  tag,
784  ( required
785  ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
786  )//required
787  ( optional
788  ( type, ( AitkenType ), AITKEN_STANDARD )
789  ( initial_theta, *( boost::is_arithmetic<mpl::_> ), 1.0 )
790  ( min_theta, *( boost::is_arithmetic<mpl::_> ), 1e-4 )
791  ( tolerance, *( boost::is_arithmetic<mpl::_> ), 1.0e-6 )
792  )//optional
793 )
794 {
795  return aitkenNew( space,type,initial_theta,tolerance,min_theta );
796 }
797 
798 
799 
800 
801 } // End namespace Feel
802 
803 #endif

Generated on Fri Oct 25 2013 14:24:03 for Feel++ by doxygen 1.8.4