]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
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 > | |
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 > | |
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> | |
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 > | |
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 |