22 * Sweeper for scalar test equation
33 */
44
5- #ifndef _SCALAR_SWEEPER_HPP_
6- #define _SCALAR_SWEEPER_HPP_
5+ #ifndef _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
6+ #define _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
77
88#include < complex>
99#include < vector>
10+
1011#include < pfasst/encap/imex_sweeper.hpp>
1112#include < pfasst/encap/vector.hpp>
1213
1314using namespace std ;
1415
1516template <typename time = pfasst::time_precision>
16- class ScalarSweeper : public pfasst ::encap::IMEXSweeper<time>
17+ class ScalarSweeper
18+ : public pfasst::encap::IMEXSweeper<time>
1719{
20+ private:
21+ typedef pfasst::encap::Encapsulation<time> encap_type;
22+ typedef pfasst::encap::VectorEncapsulation<complex <double >> complex_vector_type;
1823
19- typedef pfasst::encap::Encapsulation<time> Encapsulation;
20- typedef pfasst::encap::VectorEncapsulation<complex <double >> CVectorT;
24+ complex <double > lambda, u0;
25+ size_t n_f_expl_eval, n_f_impl_eval, n_impl_solve;
26+ const complex <double > i_complex = complex <double >(0 , 1 );
27+ double error;
2128
2229 public:
23-
24- ScalarSweeper (complex <double > lambda, complex <double > y0)
25- {
26- this ->_lambda = lambda;
27- this ->_y0 = y0;
28- this ->_nf1eval = 0 ;
29- this ->_nf2eval = 0 ;
30- this ->_nf2comp = 0 ;
31- }
30+ ScalarSweeper (complex <double > lambda, complex <double > u0)
31+ : lambda(lambda)
32+ , u0(u0)
33+ , n_f_expl_eval(0 )
34+ , n_f_impl_eval(0 )
35+ , n_impl_solve(0 )
36+ , error(0.0 )
37+ {}
3238
3339 virtual ~ScalarSweeper ()
3440 {
35- cout << " Final error: " << scientific << this ->_error << endl;
36- cout << " Number of calls to f1eval : " << this ->_nf1eval << endl;
37- cout << " Number of calls to f2eval : " << this ->_nf2eval << endl;
38- cout << " Number of calls to f2comp: " << this ->_nf2comp << endl;
41+ cout << " Final error: " << scientific << this ->error << endl;
42+ cout << " Number of explicit evaluations : " << this ->n_f_expl_eval << endl;
43+ cout << " Number of implicit evaluations : " << this ->n_f_impl_eval << endl;
44+ cout << " Number of implicit solves: " << this ->n_impl_solve << endl;
3945 }
4046
4147 void echo_error (time t)
4248 {
43- auto & qend = pfasst::encap::as_vector<complex <double >,time>(this ->get_end_state ());
49+ auto & qend = pfasst::encap::as_vector<complex <double >, time>(this ->get_end_state ());
4450
45- CVectorT qex (qend.size ());
51+ complex_vector_type qex (qend.size ());
4652
4753 this ->exact (qex, t);
4854 double max_err = abs (qend[0 ] - qex[0 ]) / abs (qex[0 ]);
4955 cout << " err: " << scientific << max_err << endl;
50- this ->_error = max_err;
56+ this ->error = max_err;
5157 }
52-
58+
5359 double get_errors ()
5460 {
55- return this ->_error ;
61+ return this ->error ;
5662 }
5763
5864 void predict (bool initial) override
@@ -71,64 +77,59 @@ class ScalarSweeper : public pfasst::encap::IMEXSweeper<time>
7177 this ->echo_error (t + dt);
7278 }
7379
74- void exact (CVectorT & q, time t)
80+ void exact (complex_vector_type & q, time t)
7581 {
76- q[0 ] = _y0 * exp (_lambda * t);
82+ q[0 ] = this -> u0 * exp (this -> lambda * t);
7783 }
7884
79- void exact (shared_ptr<Encapsulation > q_encap, time t)
85+ void exact (shared_ptr<encap_type > q_encap, time t)
8086 {
81- auto & q = pfasst::encap::as_vector<complex <double >,time>(q_encap);
82- exact (q, t);
87+ auto & q = pfasst::encap::as_vector<complex <double >, time>(q_encap);
88+ this -> exact (q, t);
8389 }
8490
85- void f_expl_eval (shared_ptr<Encapsulation > f_encap,
86- shared_ptr<Encapsulation > q_encap, time t) override
91+ void f_expl_eval (shared_ptr<encap_type > f_encap,
92+ shared_ptr<encap_type > q_encap, time t) override
8793 {
8894 UNUSED (t);
89- auto & f = pfasst::encap::as_vector<complex <double >,time>(f_encap);
90- auto & q = pfasst::encap::as_vector<complex <double >,time>(q_encap);
95+ auto & f = pfasst::encap::as_vector<complex <double >, time>(f_encap);
96+ auto & q = pfasst::encap::as_vector<complex <double >, time>(q_encap);
9197
92- // f1 = multiply with imaginary part of lambda
93- f[0 ] = i_complex * imag (this ->_lambda ) * q[0 ];
94- this ->_nf1eval ++;
98+ // f_expl = multiply with imaginary part of lambda
99+ f[0 ] = this ->i_complex * imag (this ->lambda ) * q[0 ];
100+
101+ this ->n_f_expl_eval ++;
95102 }
96103
97- void f_impl_eval (shared_ptr<Encapsulation > f_encap,
98- shared_ptr<Encapsulation > q_encap, time t) override
104+ void f_impl_eval (shared_ptr<encap_type > f_encap,
105+ shared_ptr<encap_type > q_encap, time t) override
99106 {
100107 UNUSED (t);
101- auto & f = pfasst::encap::as_vector<complex <double >,time>(f_encap);
102- auto & q = pfasst::encap::as_vector<complex <double >,time>(q_encap);
108+ auto & f = pfasst::encap::as_vector<complex <double >, time>(f_encap);
109+ auto & q = pfasst::encap::as_vector<complex <double >, time>(q_encap);
110+
111+ // f_impl = multiply with real part of lambda
112+ f[0 ] = real (this ->lambda ) * q[0 ];
103113
104- // f2 = multiply with real part of lambda
105- f[0 ] = real (this ->_lambda ) * q[0 ];
106- this ->_nf2eval ++;
114+ this ->n_f_impl_eval ++;
107115 }
108116
109- void impl_solve (shared_ptr<Encapsulation > f_encap,
110- shared_ptr<Encapsulation > q_encap, time t, time dt,
111- shared_ptr<Encapsulation > rhs_encap) override
117+ void impl_solve (shared_ptr<encap_type > f_encap,
118+ shared_ptr<encap_type > q_encap, time t, time dt,
119+ shared_ptr<encap_type > rhs_encap) override
112120 {
113121 UNUSED (t);
114- auto & f = pfasst::encap::as_vector<complex <double >,time>(f_encap);
115- auto & q = pfasst::encap::as_vector<complex <double >,time>(q_encap);
116- auto & rhs = pfasst::encap::as_vector<complex <double >,time>(rhs_encap);
122+ auto & f = pfasst::encap::as_vector<complex <double >, time>(f_encap);
123+ auto & q = pfasst::encap::as_vector<complex <double >, time>(q_encap);
124+ auto & rhs = pfasst::encap::as_vector<complex <double >, time>(rhs_encap);
117125
118- // invert f2= multiply with inverse of real part of lambda
119- double inv = 1 / (1 - double (dt) * real (this ->_lambda ));
126+ // invert f_impl = multiply with inverse of real part of lambda
127+ double inv = 1.0 / (1.0 - double (dt) * real (this ->lambda ));
120128 q[0 ] = inv * rhs[0 ];
121- f[0 ] = real (this ->_lambda ) * q[0 ];
122- this ->_nf2comp ++;
123- }
124-
125-
126- private:
127-
128- complex <double > _lambda, _y0;
129- int _nf1eval, _nf2eval, _nf2comp;
130- const complex <double > i_complex = complex <double >(0 , 1 );
131- double _error;
129+ f[0 ] = real (this ->lambda ) * q[0 ];
132130
131+ this ->n_impl_solve ++;
132+ }
133133};
134- #endif
134+
135+ #endif // _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
0 commit comments