external/mtl4/implicit_euler_mtl4.hpp

0001 /*
0002 [begin_description]
0003 Modification of the implicit Euler method, works with the MTL4 matrix library only.
0004 [end_description]
0005
0006 Copyright 2012-2013 Andreas Angelopoulos
0007 Copyright 2012-2013 Karsten Ahnert
0008 Copyright 2012-2013 Mario Mulansky
0009
0010 Distributed under the Boost Software License, Version 1.0.
0011 (See accompanying file LICENSE_1_0.txt or
0012 copy at http://www.boost.org/LICENSE_1_0.txt)
0013 */
0014
0015
0016 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
0017 #define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
0018
0019
0020 #include <utility>
0021
0022 #include <boost/numeric/odeint/util/bind.hpp>
0023 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0024 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0025
0026 #include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
0027
0028 #include <boost/numeric/mtl/mtl.hpp>
0029 #include <boost/numeric/itl/itl.hpp>
0030
0031
0032
0033
0034 namespace boost {
0035 namespace numeric {
0036 namespace odeint {
0037
0038
0039 template< class ValueType , class Resizer = initially_resizer >
0040 class implicit_euler_mtl4
0041 {
0042
0043 public:
0044
0045     typedef ValueType value_type;
0046     typedef value_type time_type;
0047     typedef mtl::dense_vector<value_type> state_type;
0048
0049     typedef state_wrapper< state_type > wrapped_state_type;
0050     typedef state_type deriv_type;
0051     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0052     typedef mtl::compressed2D< value_type > matrix_type;
0053     typedef state_wrapper< matrix_type > wrapped_matrix_type;
0054
0055     typedef Resizer resizer_type;
0056     typedef stepper_tag stepper_category;
0057
0058     typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type;
0059
0060
0061     implicit_euler_mtl4( const value_type epsilon = 1E-6 )
0062         : m_epsilon( epsilon ) , m_resizer() ,
0063           m_dxdt() , m_x() ,
0064           m_identity() , m_jacobi()
0065     { }
0066
0067
0068     template< class System >
0069     void do_step( System system , state_type &x , time_type t , time_type dt )
0070     {
0071         typedef typename odeint::unwrap_reference< System >::type system_type;
0072         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
0073         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
0074         system_type &sys = system;
0075         deriv_func_type &deriv_func = sys.first;
0076         jacobi_func_type &jacobi_func = sys.second;
0077
0078         m_resizer.adjust_size( x , detail::bind(
0079                                    &stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
0080
0081         m_identity.m_v = 1;
0082
0083         t += dt;
0084         m_x.m_v = x;
0085
0086         deriv_func( x , m_dxdt.m_v , t );
0087         jacobi_func( x , m_jacobi.m_v , t );
0088
0089
0090         m_dxdt.m_v *= -dt;
0091
0092         m_jacobi.m_v *= dt;
0093         m_jacobi.m_v -= m_identity.m_v ;
0094
0095
0096
0097         // using ilu_0 preconditioning -incomplete LU factorisation
0098         // itl::pc::diagonal<matrix_type,double> L(m_jacobi.m_v);
0099         itl::pc::ilu_0<matrix_type> L( m_jacobi.m_v );
0100
0101         solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L );
0102         x+= m_x.m_v;
0103
0104
0105     }
0106
0107
0108     template< class StateType >
0109     void adjust_size( const StateType &x )
0110     {
0111         resize_impl( x );
0112     }
0113
0114
0115 private:
0116
0117
0118     /*
0119       Applying approximate iterative linear solvers
0120       default solver is Biconjugate gradient stabilized method
0121       itl::bicgstab(A, x, b, L, iter);
0122     */
0123     template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner>
0124     void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
0125                const Preconditioner& L, int max_iteractions =500)
0126     {
0127         // Termination criterion: r < 1e-6 * b or N iterations
0128         itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 );
0129         itl::bicgstab( A , x , b , L , iter );
0130
0131     }
0132
0133
0134     template< class StateIn >
0135     bool resize_impl( const StateIn &x )
0136     {
0137         bool resized = false;
0138         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
0139         resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
0140         resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable<matrix_type>::type() );
0141         resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
0142         return resized;
0143     }
0144
0145
0146 private:
0147
0148     value_type m_epsilon;
0149     resizer_type m_resizer;
0150     wrapped_deriv_type m_dxdt;
0151     wrapped_state_type m_x;
0152     wrapped_matrix_type m_identity;
0153     wrapped_matrix_type m_jacobi;
0154 };
0155
0156
0157 } // odeint
0158 } // numeric
0159 } // boost
0160
0161
0162 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED