Merge pull request #35 from memmett/feature/by-reference

Go back to using references in AdvectionDiffusionSweeper.

Author: torbjoernk

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -14,18 +14,19 @@
 
 #include "fft.hpp"
 
-#define PI     3.1415926535897932385
-#define TWO_PI 6.2831853071795864769
+#ifndef PI
+#define PI 3.1415926535897932385
+#endif
 
 using namespace std;
+using pfasst::encap::Encapsulation;
+using pfasst::encap::as_vector;
+
 
 template<typename time = pfasst::time_precision>
-class AdvectionDiffusionSweeper 
+class AdvectionDiffusionSweeper
   : public pfasst::encap::IMEXSweeper<time>
 {
-    typedef pfasst::encap::Encapsulation<time> Encapsulation;
-    typedef pfasst::encap::VectorEncapsulation<double> DVectorT;
-
     FFT fft;
 
     vector<complex<double>> ddx, lap;
@@ -41,7 +42,7 @@ class AdvectionDiffusionSweeper
       ddx.resize(nvars);
       lap.resize(nvars);
       for (size_t i = 0; i < nvars; i++) {
-        double kx = TWO_PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
+        double kx = 2 * PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
         ddx[i] = complex<double>(0.0, 1.0) * kx;
         lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
       }
@@ -52,45 +53,42 @@ class AdvectionDiffusionSweeper
       cout << "number of f1 evals: " << nf1evals << endl;
     }
 
-    void exact(shared_ptr<Encapsulation> q, time t)
+    void exact(shared_ptr<Encapsulation<time>> q, time t)
     {
-      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
-      assert(q_cast);
-      this->exact(q_cast, t);
+      this->exact(as_vector<double,time>(q), t);
     }
 
-    void exact(shared_ptr<DVectorT> q, time t)
+    void exact(DVectorT& q, time t)
     {
-      size_t n = q->size();
+      size_t n = q.size();
       double a = 1.0 / sqrt(4 * PI * nu * (t + t0));
 
       for (size_t i = 0; i < n; i++) {
-        q->data()[i] = 0.0;
+        q[i] = 0.0;
       }
 
       for (int ii = -2; ii < 3; ii++) {
         for (size_t i = 0; i < n; i++) {
           double x = double(i) / n - 0.5 + ii - t * v;
-          q->data()[i] += a * exp(-x * x / (4 * nu * (t + t0)));
+          q[i] += a * exp(-x * x / (4 * nu * (t + t0)));
         }
       }
     }
 
     void echo_error(time t, bool predict = false)
     {
-      shared_ptr<DVectorT> qend = dynamic_pointer_cast<DVectorT>(this->get_state(this->get_nodes().size() - 1));
-      assert(qend);
-      shared_ptr<DVectorT> qex = make_shared<DVectorT>(qend->size());
+      auto& qend = as_vector<double,time>(this->get_end_state());
+      DVectorT qex(qend.size());
 
       exact(qex, t);
 
       double max = 0.0;
-      for (size_t i = 0; i < qend->size(); i++) {
-        double d = abs(qend->data()[i] - qex->data()[i]);
+      for (size_t i = 0; i < qend.size(); i++) {
+        double d = abs(qend[i] - qex[i]);
         if (d > max) { max = d; }
       }
-      cout << "err: " << scientific << max 
-           << " (" << qend->size() << ", " << predict << ")"
+      cout << "err: " << scientific << max
+           << " (" << qend.size() << ", " << predict << ")"
            << endl;
     }
 
@@ -106,76 +104,55 @@ class AdvectionDiffusionSweeper
       echo_error(t + dt);
     }
 
-    void f1eval(shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t)
+    void f1eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
     {
-      shared_ptr<DVectorT> f_cast = dynamic_pointer_cast<DVectorT>(f);
-      assert(f_cast);
-      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
-      assert(q_cast);
+      auto& q = as_vector<double,time>(_q);
+      auto& f = as_vector<double,time>(_f);
 
-      this->f1eval(f_cast, q_cast, t);
-    }
-
-    void f1eval(shared_ptr<DVectorT> f, shared_ptr<DVectorT> q, time t)
-    {
-      double c = -v / double(q->size());
+      double c = -v / double(q.size());
 
       auto* z = fft.forward(q);
-      for (size_t i = 0; i < q->size(); i++) {
+      for (size_t i = 0; i < q.size(); i++) {
         z[i] *= c * ddx[i];
       }
       fft.backward(f);
 
       nf1evals++;
     }
 
-    void f2eval(shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t)
+    void f2eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
     {
-      shared_ptr<DVectorT> f_cast = dynamic_pointer_cast<DVectorT>(f);
-      assert(f_cast);
-      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
-      assert(q_cast);
+      auto& q = as_vector<double,time>(_q);
+      auto& f = as_vector<double,time>(_f);
 
-      this->f2eval(f_cast, q_cast, t);
-    }
-
-    void f2eval(shared_ptr<DVectorT> f, shared_ptr<DVectorT> q, time t)
-    {
-      double c = nu / double(q->size());
+      double c = nu / double(q.size());
 
       auto* z = fft.forward(q);
-      for (size_t i = 0; i < q->size(); i++) {
+      for (size_t i = 0; i < q.size(); i++) {
         z[i] *= c * lap[i];
       }
       fft.backward(f);
     }
 
-    void f2comp(shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t, time dt, 
-                shared_ptr<Encapsulation> rhs)
+    void f2comp(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/, time dt,
+                shared_ptr<Encapsulation<time>> _rhs)
     {
-      shared_ptr<DVectorT> f_cast   = dynamic_pointer_cast<DVectorT>(f);
-      assert(f_cast);
-      shared_ptr<DVectorT> q_cast   = dynamic_pointer_cast<DVectorT>(q);
-      assert(q_cast);
-      shared_ptr<DVectorT> rhs_cast = dynamic_pointer_cast<DVectorT>(rhs);
-      assert(rhs_cast);
-
-      this->f2comp(f_cast, q_cast, t, dt, rhs_cast);
-    }
+      auto& q = as_vector<double,time>(_q);
+      auto& f = as_vector<double,time>(_f);
+      auto& rhs = as_vector<double,time>(_rhs);
 
-    void f2comp(shared_ptr<DVectorT> f, shared_ptr<DVectorT> q, time t, time dt, 
-                shared_ptr<DVectorT> rhs)
-    {
       auto* z = fft.forward(rhs);
-      for (size_t i = 0; i < q->size(); i++) {
-        z[i] /= (1.0 - nu * double(dt) * lap[i]) * double(q->size());
+      for (size_t i = 0; i < q.size(); i++) {
+        z[i] /= (1.0 - nu * double(dt) * lap[i]) * double(q.size());
       }
       fft.backward(q);
 
-      for (size_t i = 0; i < q->size(); i++) {
-        f->data()[i] = (q->data()[i] - rhs->data()[i]) / double(dt);
+      for (size_t i = 0; i < q.size(); i++) {
+        f[i] = (q[i] - rhs[i]) / double(dt);
       }
+
     }
+
 };
 
 #endif

examples/advection_diffusion/fft.hpp

@@ -15,10 +15,10 @@
 
 #include <fftw3.h>
 
+typedef pfasst::encap::VectorEncapsulation<double> DVectorT;
+
 class FFT
 {
-    typedef pfasst::encap::VectorEncapsulation<double> DVectorT;
-
     struct workspace {
       fftw_plan        ffft;
       fftw_plan        ifft;
@@ -36,7 +36,6 @@ class FFT
         fftw_free(wk->wk);
         fftw_destroy_plan(wk->ffft);
         fftw_destroy_plan(wk->ifft);
-//         delete wk;
       }
       workspaces.clear();
     }
@@ -55,22 +54,22 @@ class FFT
       return workspaces[ndofs];
     }
 
-    complex<double>* forward(shared_ptr<const DVectorT> x)
+    complex<double>* forward(const DVectorT& x)
     {
-      shared_ptr<workspace> wk = get_workspace(x->size());
-      for (size_t i = 0; i < x->size(); i++) {
-        wk->z[i] = x->data()[i];
+      shared_ptr<workspace> wk = get_workspace(x.size());
+      for (size_t i = 0; i < x.size(); i++) {
+        wk->z[i] = x[i];
       }
       fftw_execute_dft(wk->ffft, wk->wk, wk->wk);
       return wk->z;
     }
 
-    void backward(shared_ptr<DVectorT> x)
+    void backward(DVectorT& x)
     {
-      shared_ptr<workspace> wk = get_workspace(x->size());
+      shared_ptr<workspace> wk = get_workspace(x.size());
       fftw_execute_dft(wk->ifft, wk->wk, wk->wk);
-      for (size_t i = 0; i < x->size(); i++) {
-        x->data()[i] = real(wk->z[i]);
+      for (size_t i = 0; i < x.size(); i++) {
+        x[i] = real(wk->z[i]);
       }
     }
 

examples/advection_diffusion/serial_mlsdc.cpp

@@ -18,7 +18,7 @@
 using namespace pfasst;
 using namespace pfasst::encap;
 
-int main(int argc, char** argv)
+int main(int /*argc*/, char** /*argv*/)
 {
   MLSDC<> mlsdc;
 

examples/advection_diffusion/serial_mlsdc_autobuild.cpp

@@ -24,7 +24,7 @@ using namespace std;
 using namespace pfasst;
 using namespace pfasst::encap;
 
-int main(int argc, char** argv)
+int main(int /*argc*/, char** /*argv*/)
 {
   MLSDC<> mlsdc;
 

examples/advection_diffusion/spectral_transfer_1d.hpp

@@ -14,63 +14,49 @@
 #include "fft.hpp"
 
 template<typename time = pfasst::time_precision>
-class SpectralTransfer1D 
+class SpectralTransfer1D
   : public pfasst::encap::PolyInterpMixin<time>
 {
     typedef pfasst::encap::Encapsulation<double> Encapsulation;
-    typedef pfasst::encap::VectorEncapsulation<double> DVectorT;
 
     FFT fft;
 
   public:
+
     void interpolate(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
     {
-      shared_ptr<DVectorT> fine = dynamic_pointer_cast<DVectorT>(dst);
-      assert(fine);
-      shared_ptr<const DVectorT> crse = dynamic_pointer_cast<const DVectorT>(src);
-      assert(crse);
-
-      this->interpolate(fine, crse);
-    }
+      auto& fine = as_vector<double,time>(dst);
+      auto& crse = as_vector<double,time>(src);
 
-    void interpolate(shared_ptr<DVectorT> fine, shared_ptr<const DVectorT> crse)
-    {
       auto* crse_z = fft.forward(crse);
-      auto* fine_z = fft.get_workspace(fine->size())->z;
+      auto* fine_z = fft.get_workspace(fine.size())->z;
 
-      for (size_t i = 0; i < fine->size(); i++) {
+      for (size_t i = 0; i < fine.size(); i++) {
         fine_z[i] = 0.0;
       }
 
-      double c = 1.0 / crse->size();
+      double c = 1.0 / crse.size();
 
-      for (size_t i = 0; i < crse->size() / 2; i++) {
+      for (size_t i = 0; i < crse.size() / 2; i++) {
         fine_z[i] = c * crse_z[i];
       }
 
-      for (size_t i = 1; i < crse->size() / 2; i++) {
-        fine_z[fine->size() - crse->size() / 2 + i] = c * crse_z[crse->size() / 2 + i];
+      for (size_t i = 1; i < crse.size() / 2; i++) {
+        fine_z[fine.size() - crse.size() / 2 + i] = c * crse_z[crse.size() / 2 + i];
       }
 
       fft.backward(fine);
     }
 
     void restrict(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
     {
-      shared_ptr<DVectorT> crse = dynamic_pointer_cast<DVectorT>(dst);
-      assert(crse);
-      shared_ptr<const DVectorT> fine = dynamic_pointer_cast<const DVectorT>(src);
-      assert(fine);
+      auto& fine = as_vector<double,time>(src);
+      auto& crse = as_vector<double,time>(dst);
 
-      this->restrict(crse, fine);
-    }
-
-    void restrict(shared_ptr<DVectorT> crse, shared_ptr<const DVectorT> fine)
-    {
-      size_t xrat = fine->size() / crse->size();
+      size_t xrat = fine.size() / crse.size();
 
-      for (size_t i = 0; i < crse->size(); i++) {
-        crse->data()[i] = fine->at(xrat * i);
+      for (size_t i = 0; i < crse.size(); i++) {
+        crse[i] = fine[xrat*i];
       }
     }
 

examples/advection_diffusion/vanilla_sdc.cpp

@@ -14,7 +14,7 @@
 
 #include "advection_diffusion_sweeper.hpp"
 
-int main(int argc, char** argv)
+int main(int /*argc*/, char** /*argv*/)
 {
   pfasst::SDC<> sdc;