55#ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
66#define _ADVECTION_DIFFUSION_SWEEPER_HPP_
77
8+ #include < cstdlib>
89#include < cassert>
910#include < complex>
1011#include < map>
1112#include < ostream>
1213#include < vector>
14+ #include < utility>
1315
16+ #include < pfasst/globals.hpp>
1417#include < pfasst/encap/imex_sweeper.hpp>
1518
1619#include " fft.hpp"
@@ -24,42 +27,56 @@ using pfasst::encap::Encapsulation;
2427using pfasst::encap::as_vector;
2528
2629
30+ /* *
31+ * errors at different iterations and time nodes
32+ *
33+ * Mapping a pair of step/iteration indices onto the error of the solution.
34+ */
2735typedef map<pair<size_t , size_t >, double > error_map;
2836
2937template <typename time = pfasst::time_precision>
3038class AdvectionDiffusionSweeper
3139 : public pfasst::encap::IMEXSweeper<time>
3240{
41+ // ! @{
3342 FFT fft;
34-
3543 vector<complex <double >> ddx, lap;
44+ // ! @}
45+
46+ // ! @{
3647 error_map errors;
48+ // ! @}
3749
50+ // ! @{
3851 double v = 1.0 ;
3952 time t0 = 1.0 ;
4053 double nu = 0.02 ;
4154 size_t nf1evals = 0 ;
55+ // ! @}
4256
4357 public:
58+ // ! @{
4459 AdvectionDiffusionSweeper (size_t nvars)
4560 {
46- ddx.resize (nvars);
47- lap.resize (nvars);
61+ this -> ddx .resize (nvars);
62+ this -> lap .resize (nvars);
4863 for (size_t i = 0 ; i < nvars; i++) {
4964 double kx = 2 * PI * ((i <= nvars / 2 ) ? int (i) : int (i) - int (nvars));
50- ddx[i] = complex <double >(0.0 , 1.0 ) * kx;
51- lap[i] = (kx * kx < 1e-13 ) ? 0.0 : -kx * kx;
65+ this -> ddx [i] = complex <double >(0.0 , 1.0 ) * kx;
66+ this -> lap [i] = (kx * kx < 1e-13 ) ? 0.0 : -kx * kx;
5267 }
5368 }
5469
55- ~AdvectionDiffusionSweeper ()
70+ virtual ~AdvectionDiffusionSweeper ()
5671 {
57- cout << " number of f1 evals: " << nf1evals << endl;
72+ cout << " number of f1 evals: " << this -> nf1evals << endl;
5873 }
74+ // ! @}
5975
76+ // ! @{
6077 void exact (shared_ptr<Encapsulation<time>> q, time t)
6178 {
62- this ->exact (as_vector<double ,time>(q), t);
79+ this ->exact (as_vector<double , time>(q), t);
6380 }
6481
6582 void exact (DVectorT& q, time t)
@@ -81,10 +98,10 @@ class AdvectionDiffusionSweeper
8198
8299 void echo_error (time t, bool predict = false )
83100 {
84- auto & qend = as_vector<double ,time>(this ->get_end_state ());
101+ auto & qend = as_vector<double , time>(this ->get_end_state ());
85102 DVectorT qex (qend.size ());
86103
87- exact (qex, t);
104+ this -> exact (qex, t);
88105
89106 double max = 0.0 ;
90107 for (size_t i = 0 ; i < qend.size (); i++) {
@@ -98,79 +115,111 @@ class AdvectionDiffusionSweeper
98115 << " (" << qend.size () << " , " << predict << " )"
99116 << endl;
100117
101- errors.insert (pair<pair<size_t , size_t >, double >
102- (pair<size_t , size_t >(n, k), max));
118+ this ->errors .insert (pair<pair<size_t , size_t >, double >(pair<size_t , size_t >(n, k), max));
103119 }
104120
121+ /* *
122+ * retrieve errors at iterations and time nodes
123+ */
105124 error_map get_errors ()
106125 {
107- return errors;
126+ return this -> errors ;
108127 }
128+ // ! @}
109129
110- void predict (bool initial)
130+ // ! @{
131+ /* *
132+ * @copybrief pfasst::encap::IMEXSweeper::predict()
133+ */
134+ void predict (bool initial) override
111135 {
112136 pfasst::encap::IMEXSweeper<time>::predict (initial);
113137 time t = this ->get_controller ()->get_time ();
114138 time dt = this ->get_controller ()->get_time_step ();
115- echo_error (t + dt, true );
139+ this -> echo_error (t + dt, true );
116140 }
117141
118- void sweep ()
142+ /* *
143+ * @copybrief pfasst::encap::IMEXSweeper::sweep()
144+ */
145+ void sweep () override
119146 {
120147 pfasst::encap::IMEXSweeper<time>::sweep ();
121148 time t = this ->get_controller ()->get_time ();
122149 time dt = this ->get_controller ()->get_time_step ();
123- echo_error (t + dt);
150+ this -> echo_error (t + dt);
124151 }
125-
126- void f1eval (shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /* t*/ )
152+ // ! @}
153+
154+ // ! @{
155+ /* *
156+ * @copybrief pfasst::encap::IMEXSweeper::f_expl_eval()
157+ */
158+ void f_expl_eval (shared_ptr<Encapsulation<time>> f_expl_encap,
159+ shared_ptr<Encapsulation<time>> u_encap,
160+ time t) override
127161 {
128- auto & q = as_vector<double ,time>(_q);
129- auto & f = as_vector<double ,time>(_f);
162+ UNUSED (t);
163+ auto & u = as_vector<double , time>(u_encap);
164+ auto & f_expl = as_vector<double , time>(f_expl_encap);
130165
131- double c = -v / double (q .size ());
166+ double c = -v / double (u .size ());
132167
133- auto * z = fft.forward (q );
134- for (size_t i = 0 ; i < q .size (); i++) {
135- z[i] *= c * ddx[i];
168+ auto * z = this -> fft .forward (u );
169+ for (size_t i = 0 ; i < u .size (); i++) {
170+ z[i] *= c * this -> ddx [i];
136171 }
137- fft.backward (f );
172+ this -> fft .backward (f_expl );
138173
139- nf1evals++;
174+ this -> nf1evals ++;
140175 }
141176
142- void f2eval (shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /* t*/ )
177+ /* *
178+ * @copybrief pfasst::encap::IMEXSweeper::f_impl_eval()
179+ */
180+ void f_impl_eval (shared_ptr<Encapsulation<time>> f_impl_encap,
181+ shared_ptr<Encapsulation<time>> u_encap,
182+ time t) override
143183 {
144- auto & q = as_vector<double ,time>(_q);
145- auto & f = as_vector<double ,time>(_f);
184+ UNUSED (t);
185+ auto & u = as_vector<double , time>(u_encap);
186+ auto & f_impl = as_vector<double , time>(f_impl_encap);
146187
147- double c = nu / double (q .size ());
188+ double c = nu / double (u .size ());
148189
149- auto * z = fft.forward (q );
150- for (size_t i = 0 ; i < q .size (); i++) {
151- z[i] *= c * lap[i];
190+ auto * z = this -> fft .forward (u );
191+ for (size_t i = 0 ; i < u .size (); i++) {
192+ z[i] *= c * this -> lap [i];
152193 }
153- fft.backward (f );
194+ this -> fft .backward (f_impl );
154195 }
155196
156- void f2comp (shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /* t*/ , time dt,
157- shared_ptr<Encapsulation<time>> _rhs)
197+ /* *
198+ * @copybrief pfasst::encap::IMEXSweeper::impl_solve()
199+ */
200+ void impl_solve (shared_ptr<Encapsulation<time>> f_impl_encap,
201+ shared_ptr<Encapsulation<time>> u_encap,
202+ time t, time dt,
203+ shared_ptr<Encapsulation<time>> rhs_encap) override
158204 {
159- auto & q = as_vector<double ,time>(_q);
160- auto & f = as_vector<double ,time>(_f);
161- auto & rhs = as_vector<double ,time>(_rhs);
205+ UNUSED (t);
206+ auto & u = as_vector<double , time>(u_encap);
207+ auto & f_impl = as_vector<double , time>(f_impl_encap);
208+ auto & rhs = as_vector<double , time>(rhs_encap);
162209
163- auto * z = fft.forward (rhs);
164- for (size_t i = 0 ; i < q.size (); i++) {
165- z[i] /= (1.0 - nu * double (dt) * lap[i]) * double (q.size ());
166- }
167- fft.backward (q);
210+ double c = nu * double (dt);
168211
169- for (size_t i = 0 ; i < q.size (); i++) {
170- f[i] = (q[i] - rhs[i]) / double (dt);
212+ auto * z = this ->fft .forward (rhs);
213+ for (size_t i = 0 ; i < u.size (); i++) {
214+ z[i] /= (1.0 - c * this ->lap [i]) * double (u.size ());
171215 }
216+ this ->fft .backward (u);
172217
218+ for (size_t i = 0 ; i < u.size (); i++) {
219+ f_impl[i] = (u[i] - rhs[i]) / double (dt);
220+ }
173221 }
222+ // ! @}
174223
175224};
176225
0 commit comments