1 /*
2 [begin_description]
3 Modification of the implicit Euler method, works with the MTL4 matrix library only.
4 [end_description]
5 
6 Copyright 2012-2013 Andreas Angelopoulos
7 Copyright 2012-2013 Karsten Ahnert
8 Copyright 2012-2013 Mario Mulansky
9 
10 Distributed under the Boost Software License, Version 1.0.
11 (See accompanying file LICENSE_1_0.txt or
12 copy at http://www.boost.org/LICENSE_1_0.txt)
13 */
14 
15 
16 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
17 #define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
18 
19 
20 #include <utility>
21 
22 #include <boost/numeric/odeint/util/bind.hpp>
23 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
24 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
25 
26 #include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
27 
28 #include <boost/numeric/mtl/mtl.hpp>
29 #include <boost/numeric/itl/itl.hpp>
30 
31 
32 
33 
34 namespace boost {
35 namespace numeric {
36 namespace odeint {
37 
38 
39 template< class ValueType , class Resizer = initially_resizer >
40 class implicit_euler_mtl4
41 {
42 
43 public:
44 
45     typedef ValueType value_type;
46     typedef value_type time_type;
47     typedef mtl::dense_vector<value_type> state_type;
48 
49     typedef state_wrapper< state_type > wrapped_state_type;
50     typedef state_type deriv_type;
51     typedef state_wrapper< deriv_type > wrapped_deriv_type;
52     typedef mtl::compressed2D< value_type > matrix_type;
53     typedef state_wrapper< matrix_type > wrapped_matrix_type;
54 
55     typedef Resizer resizer_type;
56     typedef stepper_tag stepper_category;
57 
58     typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type;
59 
60 
implicit_euler_mtl4(const value_type epsilon=1E-6)61     implicit_euler_mtl4( const value_type epsilon = 1E-6 )
62         : m_epsilon( epsilon ) , m_resizer() ,
63           m_dxdt() , m_x() ,
64           m_identity() , m_jacobi()
65     { }
66 
67 
68     template< class System >
do_step(System system,state_type & x,time_type t,time_type dt)69     void do_step( System system , state_type &x , time_type t , time_type dt )
70     {
71         typedef typename odeint::unwrap_reference< System >::type system_type;
72         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
73         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
74         system_type &sys = system;
75         deriv_func_type &deriv_func = sys.first;
76         jacobi_func_type &jacobi_func = sys.second;
77 
78         m_resizer.adjust_size( x , detail::bind(
79                                    &stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
80 
81         m_identity.m_v = 1;
82 
83         t += dt;
84         m_x.m_v = x;
85 
86         deriv_func( x , m_dxdt.m_v , t );
87         jacobi_func( x , m_jacobi.m_v , t );
88 
89 
90         m_dxdt.m_v *= -dt;
91 
92         m_jacobi.m_v *= dt;
93         m_jacobi.m_v -= m_identity.m_v ;
94 
95 
96 
97         // using ilu_0 preconditioning -incomplete LU factorisation
98         // itl::pc::diagonal<matrix_type,double> L(m_jacobi.m_v);
99         itl::pc::ilu_0<matrix_type> L( m_jacobi.m_v );
100 
101         solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L );
102         x+= m_x.m_v;
103 
104 
105     }
106 
107 
108     template< class StateType >
adjust_size(const StateType & x)109     void adjust_size( const StateType &x )
110     {
111         resize_impl( x );
112     }
113 
114 
115 private:
116 
117 
118     /*
119       Applying approximate iterative linear solvers
120       default solver is Biconjugate gradient stabilized method
121       itl::bicgstab(A, x, b, L, iter);
122     */
123     template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner>
solve(const LinearOperator & A,HilbertSpaceX & x,const HilbertSpaceB & b,const Preconditioner & L,int max_iteractions=500)124     void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
125                const Preconditioner& L, int max_iteractions =500)
126     {
127         // Termination criterion: r < 1e-6 * b or N iterations
128         itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 );
129         itl::bicgstab( A , x , b , L , iter );
130 
131     }
132 
133 
134     template< class StateIn >
resize_impl(const StateIn & x)135     bool resize_impl( const StateIn &x )
136     {
137         bool resized = false;
138         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
139         resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
140         resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable<matrix_type>::type() );
141         resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
142         return resized;
143     }
144 
145 
146 private:
147 
148     value_type m_epsilon;
149     resizer_type m_resizer;
150     wrapped_deriv_type m_dxdt;
151     wrapped_state_type m_x;
152     wrapped_matrix_type m_identity;
153     wrapped_matrix_type m_jacobi;
154 };
155 
156 
157 } // odeint
158 } // numeric
159 } // boost
160 
161 
162 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
163