88#include < complex>
99#include < vector>
1010#include < cassert>
11+ #include < ostream>
1112
1213#include < pfasst/encap/imex_sweeper.hpp>
14+
1315#include " fft.hpp"
1416
1517#define PI 3.1415926535897932385
1820using namespace std ;
1921
2022template <typename time = pfasst::time_precision>
21- class AdvectionDiffusionSweeper : public pfasst ::encap::IMEXSweeper<time>
23+ class AdvectionDiffusionSweeper
24+ : public pfasst::encap::IMEXSweeper<time>
2225{
2326 typedef pfasst::encap::Encapsulation<time> Encapsulation;
2427 typedef pfasst::encap::VectorEncapsulation<double > DVectorT;
@@ -49,44 +52,46 @@ class AdvectionDiffusionSweeper : public pfasst::encap::IMEXSweeper<time>
4952 cout << " number of f1 evals: " << nf1evals << endl;
5053 }
5154
52- void exact (Encapsulation* q, time t)
55+ void exact (shared_ptr< Encapsulation> q, time t)
5356 {
54- DVectorT* q_cast = dynamic_cast <DVectorT* >(q);
55- assert (q_cast != nullptr );
56- this ->exact (* q_cast, t);
57+ shared_ptr< DVectorT> q_cast = dynamic_pointer_cast <DVectorT>(q);
58+ assert (q_cast);
59+ this ->exact (q_cast, t);
5760 }
5861
59- void exact (DVectorT& q, time t)
62+ void exact (shared_ptr< DVectorT> q, time t)
6063 {
61- size_t n = q. size ();
64+ size_t n = q-> size ();
6265 double a = 1.0 / sqrt (4 * PI * nu * (t + t0));
6366
6467 for (size_t i = 0 ; i < n; i++) {
65- q[i] = 0.0 ;
68+ q-> data () [i] = 0.0 ;
6669 }
6770
6871 for (int ii = -2 ; ii < 3 ; ii++) {
6972 for (size_t i = 0 ; i < n; i++) {
7073 double x = double (i) / n - 0.5 + ii - t * v;
71- q[i] += a * exp (-x * x / (4 * nu * (t + t0)));
74+ q-> data () [i] += a * exp (-x * x / (4 * nu * (t + t0)));
7275 }
7376 }
7477 }
7578
7679 void echo_error (time t, bool predict = false )
7780 {
78- DVectorT* qend = dynamic_cast <DVectorT* >(this ->get_state (this ->get_nodes ().size () - 1 ));
79- assert (qend != nullptr );
80- DVectorT qex = DVectorT (qend->size ());
81+ shared_ptr< DVectorT> qend = dynamic_pointer_cast <DVectorT>(this ->get_state (this ->get_nodes ().size () - 1 ));
82+ assert (qend);
83+ shared_ptr< DVectorT> qex = make_shared< DVectorT> (qend->size ());
8184
8285 exact (qex, t);
8386
8487 double max = 0.0 ;
8588 for (size_t i = 0 ; i < qend->size (); i++) {
86- double d = abs (qend->at (i) - qex[i]);
89+ double d = abs (qend->data ()[i] - qex-> data () [i]);
8790 if (d > max) { max = d; }
8891 }
89- cout << " err: " << scientific << max << " (" << qend->size () << " , " << predict << " )" << endl;
92+ cout << " err: " << scientific << max
93+ << " (" << qend->size () << " , " << predict << " )"
94+ << endl;
9095 }
9196
9297 void predict (time t, time dt, bool initial)
@@ -101,75 +106,76 @@ class AdvectionDiffusionSweeper : public pfasst::encap::IMEXSweeper<time>
101106 echo_error (t + dt);
102107 }
103108
104- void f1eval (Encapsulation* f, Encapsulation* q, time t)
109+ void f1eval (shared_ptr< Encapsulation> f, shared_ptr< Encapsulation> q, time t)
105110 {
106- DVectorT* f_cast = dynamic_cast <DVectorT* >(f);
107- assert (f_cast != nullptr );
108- DVectorT* q_cast = dynamic_cast <DVectorT* >(q);
109- assert (q_cast != nullptr );
111+ shared_ptr< DVectorT> f_cast = dynamic_pointer_cast <DVectorT>(f);
112+ assert (f_cast);
113+ shared_ptr< DVectorT> q_cast = dynamic_pointer_cast <DVectorT>(q);
114+ assert (q_cast);
110115
111116 this ->f1eval (f_cast, q_cast, t);
112117 }
113118
114- void f1eval (DVectorT* f, DVectorT* q, time t)
119+ void f1eval (shared_ptr< DVectorT> f, shared_ptr< DVectorT> q, time t)
115120 {
116121 double c = -v / double (q->size ());
117122
118- auto * z = fft.forward (* q);
123+ auto * z = fft.forward (q);
119124 for (size_t i = 0 ; i < q->size (); i++) {
120125 z[i] *= c * ddx[i];
121126 }
122- fft.backward (* f);
127+ fft.backward (f);
123128
124129 nf1evals++;
125130 }
126131
127- void f2eval (Encapsulation* f, Encapsulation* q, time t)
132+ void f2eval (shared_ptr< Encapsulation> f, shared_ptr< Encapsulation> q, time t)
128133 {
129- DVectorT* f_cast = dynamic_cast <DVectorT* >(f);
130- assert (f_cast != nullptr );
131- DVectorT* q_cast = dynamic_cast <DVectorT* >(q);
132- assert (q_cast != nullptr );
134+ shared_ptr< DVectorT> f_cast = dynamic_pointer_cast <DVectorT>(f);
135+ assert (f_cast);
136+ shared_ptr< DVectorT> q_cast = dynamic_pointer_cast <DVectorT>(q);
137+ assert (q_cast);
133138
134139 this ->f2eval (f_cast, q_cast, t);
135140 }
136141
137- void f2eval (DVectorT* f, DVectorT* q, time t)
142+ void f2eval (shared_ptr< DVectorT> f, shared_ptr< DVectorT> q, time t)
138143 {
139144 double c = nu / double (q->size ());
140145
141- auto * z = fft.forward (* q);
146+ auto * z = fft.forward (q);
142147 for (size_t i = 0 ; i < q->size (); i++) {
143148 z[i] *= c * lap[i];
144149 }
145- fft.backward (* f);
150+ fft.backward (f);
146151 }
147152
148- void f2comp (Encapsulation* f, Encapsulation* q, time t, time dt, Encapsulation* rhs)
153+ void f2comp (shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t, time dt,
154+ shared_ptr<Encapsulation> rhs)
149155 {
150- DVectorT* f_cast = dynamic_cast <DVectorT* >(f);
151- assert (f_cast != nullptr );
152- DVectorT* q_cast = dynamic_cast <DVectorT* >(q);
153- assert (q_cast != nullptr );
154- DVectorT* rhs_cast = dynamic_cast <DVectorT* >(rhs);
155- assert (rhs_cast != nullptr );
156+ shared_ptr< DVectorT> f_cast = dynamic_pointer_cast <DVectorT>(f);
157+ assert (f_cast);
158+ shared_ptr< DVectorT> q_cast = dynamic_pointer_cast <DVectorT>(q);
159+ assert (q_cast);
160+ shared_ptr< DVectorT> rhs_cast = dynamic_pointer_cast <DVectorT>(rhs);
161+ assert (rhs_cast);
156162
157163 this ->f2comp (f_cast, q_cast, t, dt, rhs_cast);
158164 }
159165
160- void f2comp (DVectorT* f, DVectorT* q, time t, time dt, DVectorT* rhs)
166+ void f2comp (shared_ptr<DVectorT> f, shared_ptr<DVectorT> q, time t, time dt,
167+ shared_ptr<DVectorT> rhs)
161168 {
162- auto * z = fft.forward (* rhs);
169+ auto * z = fft.forward (rhs);
163170 for (size_t i = 0 ; i < q->size (); i++) {
164171 z[i] /= (1.0 - nu * double (dt) * lap[i]) * double (q->size ());
165172 }
166- fft.backward (* q);
173+ fft.backward (q);
167174
168175 for (size_t i = 0 ; i < q->size (); i++) {
169- f->at (i) = (q->at (i) - rhs->at (i) ) / double (dt);
176+ f->data ()[i] = (q->data ()[i] - rhs->data ()[i] ) / double (dt);
170177 }
171178 }
172-
173179};
174180
175181#endif
0 commit comments