Merge branch 'release/v0.1.0' into feature/state

Signed-off-by: Matthew Emmett <[email protected]>

Conflicts:
	examples/advection_diffusion/advection_diffusion_sweeper.hpp
	examples/advection_diffusion/serial_mlsdc.cpp
	examples/advection_diffusion/vanilla_sdc.cpp
	include/pfasst/encap/imex_sweeper.hpp

Author: memmett

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -15,20 +15,21 @@
 
 #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;
+
 
 typedef map<pair<size_t, size_t>, double> error_map;
 
 template<typename time = pfasst::time_precision>
 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;
@@ -45,7 +46,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;
       }
@@ -56,48 +57,45 @@ 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; }
       }
 
       auto n = this->get_controller()->get_step();
       auto k = this->get_controller()->get_iteration();
       cout << "err: " << n << " " << k << " " << scientific << max
-           << " (" << qend->size() << ", " << predict << ")"
+           << " (" << qend.size() << ", " << predict << ")"
            << endl;
 
       errors.insert(pair<pair<size_t, size_t>, double>
@@ -125,76 +123,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

@@ -79,7 +79,7 @@ error_map run_serial_mlsdc()
 }
 
 #ifndef PFASST_UNIT_TESTING
-int main(int argc, char** argv)
+int main(int /*argc*/, char** /*argv*/)
 {
   run_serial_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];
       }
     }
 

include/pfasst/encap/encap_sweeper.hpp

@@ -17,7 +17,7 @@ namespace pfasst
   {
 
     template<typename time = time_precision>
-    class EncapSweeper 
+    class EncapSweeper
       : public ISweeper<time>
     {
         vector<time> nodes;
@@ -45,24 +45,24 @@ namespace pfasst
           return factory;
         }
 
-        virtual void set_state(shared_ptr<const Encapsulation<time>> q0, size_t m)
+        virtual void set_state(shared_ptr<const Encapsulation<time>> /*q0*/, size_t /*m*/)
         {
           throw NotImplementedYet("sweeper");
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_state(size_t m) const
+        virtual shared_ptr<Encapsulation<time>> get_state(size_t /*m*/) const
         {
           throw NotImplementedYet("sweeper");
           return NULL;
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_tau(size_t m) const
+        virtual shared_ptr<Encapsulation<time>> get_tau(size_t /*m*/) const
         {
           throw NotImplementedYet("sweeper");
           return NULL;
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_saved_state(size_t m) const
+        virtual shared_ptr<Encapsulation<time>> get_saved_state(size_t /*m*/) const
         {
           throw NotImplementedYet("sweeper");
           return NULL;
@@ -73,7 +73,7 @@ namespace pfasst
           return this->get_state(this->get_nodes().size() - 1);
         }
 
-        virtual void evaluate(size_t m)
+      virtual void evaluate(size_t /*m*/)
         {
           throw NotImplementedYet("sweeper");
         }
@@ -83,7 +83,7 @@ namespace pfasst
           throw NotImplementedYet("sweeper");
         }
 
-        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const
+        virtual void integrate(time /*dt*/, vector<shared_ptr<Encapsulation<time>>> /*dst*/) const
         {
           throw NotImplementedYet("sweeper");
         }