1313#include < vector>
1414#include < utility>
1515
16+ #include < pfasst/globals.hpp>
1617#include < pfasst/encap/imex_sweeper.hpp>
1718
1819#include " fft.hpp"
@@ -26,23 +27,35 @@ using pfasst::encap::Encapsulation;
2627using pfasst::encap::as_vector;
2728
2829
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+ */
2935typedef map<pair<size_t , size_t >, double > error_map;
3036
3137template <typename time = pfasst::time_precision>
3238class AdvectionDiffusionSweeper
3339 : public pfasst::encap::IMEXSweeper<time>
3440{
41+ // ! @{
3542 FFT fft;
36-
3743 vector<complex <double >> ddx, lap;
44+ // ! @}
45+
46+ // ! @{
3847 error_map errors;
48+ // ! @}
3949
50+ // ! @{
4051 double v = 1.0 ;
4152 time t0 = 1.0 ;
4253 double nu = 0.02 ;
4354 size_t nf1evals = 0 ;
55+ // ! @}
4456
4557 public:
58+ // ! @{
4659 AdvectionDiffusionSweeper (size_t nvars)
4760 {
4861 this ->ddx .resize (nvars);
@@ -54,11 +67,13 @@ class AdvectionDiffusionSweeper
5467 }
5568 }
5669
57- ~AdvectionDiffusionSweeper ()
70+ virtual ~AdvectionDiffusionSweeper ()
5871 {
5972 cout << " number of f1 evals: " << this ->nf1evals << endl;
6073 }
74+ // ! @}
6175
76+ // ! @{
6277 void exact (shared_ptr<Encapsulation<time>> q, time t)
6378 {
6479 this ->exact (as_vector<double , time>(q), t);
@@ -103,11 +118,19 @@ class AdvectionDiffusionSweeper
103118 this ->errors .insert (pair<pair<size_t , size_t >, double >(pair<size_t , size_t >(n, k), max));
104119 }
105120
121+ /* *
122+ * retrieve errors at iterations and time nodes
123+ */
106124 error_map get_errors ()
107125 {
108126 return this ->errors ;
109127 }
128+ // ! @}
110129
130+ // ! @{
131+ /* *
132+ * @copybrief pfasst::encap::IMEXSweeper::predict()
133+ */
111134 void predict (bool initial) override
112135 {
113136 pfasst::encap::IMEXSweeper<time>::predict (initial);
@@ -116,66 +139,87 @@ class AdvectionDiffusionSweeper
116139 this ->echo_error (t + dt, true );
117140 }
118141
142+ /* *
143+ * @copybrief pfasst::encap::IMEXSweeper::sweep()
144+ */
119145 void sweep () override
120146 {
121147 pfasst::encap::IMEXSweeper<time>::sweep ();
122148 time t = this ->get_controller ()->get_time ();
123149 time dt = this ->get_controller ()->get_time_step ();
124150 this ->echo_error (t + dt);
125151 }
126-
127- void f_eval_expl (shared_ptr<Encapsulation<time>> _f,
128- shared_ptr<Encapsulation<time>> _q,
129- time /* t*/ ) override
152+ // ! @}
153+
154+ // ! @{
155+ /* *
156+ * @copybrief pfasst::encap::IMEXSweeper::f_expl_eval()
157+ */
158+ void f_expl_eval (shared_ptr<Encapsulation<time>> f_expl,
159+ shared_ptr<Encapsulation<time>> u,
160+ time t) override
130161 {
131- auto & q = as_vector<double , time>(_q);
132- auto & f = as_vector<double , time>(_f);
162+ UNUSED (t);
163+ auto & u_cast = as_vector<double , time>(u);
164+ auto & f_expl_cast = as_vector<double , time>(f_expl);
133165
134- double c = -v / double (q .size ());
166+ double c = -v / double (u_cast .size ());
135167
136- auto * z = this ->fft .forward (q );
137- for (size_t i = 0 ; i < q .size (); i++) {
168+ auto * z = this ->fft .forward (u_cast );
169+ for (size_t i = 0 ; i < u_cast .size (); i++) {
138170 z[i] *= c * this ->ddx [i];
139171 }
140- this ->fft .backward (f );
172+ this ->fft .backward (f_expl_cast );
141173
142174 this ->nf1evals ++;
143175 }
144176
145- void f_eval_impl (shared_ptr<Encapsulation<time>> _f,
146- shared_ptr<Encapsulation<time>> _q,
147- time /* t*/ ) override
177+ /* *
178+ * @copybrief pfasst::encap::IMEXSweeper::f_impl_eval()
179+ */
180+ void f_impl_eval (shared_ptr<Encapsulation<time>> f_impl,
181+ shared_ptr<Encapsulation<time>> u,
182+ time t) override
148183 {
149- auto & q = as_vector<double , time>(_q);
150- auto & f = as_vector<double , time>(_f);
184+ UNUSED (t);
185+ auto & u_cast = as_vector<double , time>(u);
186+ auto & f_impl_cast = as_vector<double , time>(f_impl);
151187
152- double c = nu / double (q .size ());
188+ double c = nu / double (u_cast .size ());
153189
154- auto * z = this ->fft .forward (q );
155- for (size_t i = 0 ; i < q .size (); i++) {
190+ auto * z = this ->fft .forward (u_cast );
191+ for (size_t i = 0 ; i < u_cast .size (); i++) {
156192 z[i] *= c * this ->lap [i];
157193 }
158- this ->fft .backward (f );
194+ this ->fft .backward (f_impl_cast );
159195 }
160196
161- void impl_solve (shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q,
162- time /* t*/ , time dt, shared_ptr<Encapsulation<time>> _rhs) override
197+ /* *
198+ * @copybrief pfasst::encap::IMEXSweeper::impl_solve()
199+ */
200+ void impl_solve (shared_ptr<Encapsulation<time>> f_impl,
201+ shared_ptr<Encapsulation<time>> u,
202+ time t, time dt,
203+ shared_ptr<Encapsulation<time>> rhs) override
163204 {
164- auto & q = as_vector<double , time>(_q);
165- auto & f = as_vector<double , time>(_f);
166- auto & rhs = as_vector<double , time>(_rhs);
205+ UNUSED (t);
206+ auto & u_cast = as_vector<double , time>(u);
207+ auto & f_impl_cast = as_vector<double , time>(f_impl);
208+ auto & rhs_cast = as_vector<double , time>(rhs);
167209
168- auto * z = this ->fft .forward (rhs);
169- for (size_t i = 0 ; i < q.size (); i++) {
170- z[i] /= (1.0 - nu * double (dt) * this ->lap [i]) * double (q.size ());
171- }
172- this ->fft .backward (q);
210+ double c = nu * double (dt);
173211
174- for (size_t i = 0 ; i < q.size (); i++) {
175- f[i] = (q[i] - rhs[i]) / double (dt);
212+ auto * z = this ->fft .forward (rhs_cast);
213+ for (size_t i = 0 ; i < u_cast.size (); i++) {
214+ z[i] /= (1.0 - c * this ->lap [i]) * double (u_cast.size ());
176215 }
216+ this ->fft .backward (u_cast);
177217
218+ for (size_t i = 0 ; i < u_cast.size (); i++) {
219+ f_impl_cast[i] = (u_cast[i] - rhs_cast[i]) / double (dt);
220+ }
178221 }
222+ // ! @}
179223
180224};
181225
0 commit comments