[WIP] Quadrature Revamping

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -38,6 +38,7 @@ template<typename time = pfasst::time_precision>
 class AdvectionDiffusionSweeper
   : public pfasst::encap::IMEXSweeper<time>
 {
+  protected:
     //! @{
     FFT fft;
     vector<complex<double>> ddx, lap;
@@ -56,7 +57,7 @@ class AdvectionDiffusionSweeper
 
   public:
     //! @{
-    AdvectionDiffusionSweeper(size_t nvars)
+    explicit AdvectionDiffusionSweeper(size_t nvars)
     {
       this->ddx.resize(nvars);
       this->lap.resize(nvars);
@@ -67,12 +68,68 @@ class AdvectionDiffusionSweeper
       }
     }
 
+    AdvectionDiffusionSweeper() = default;
+
     virtual ~AdvectionDiffusionSweeper()
     {
       cout << "number of f1 evals: " << this->nf1evals << endl;
     }
     //! @}
 
+    //! @{
+    void exact(shared_ptr<Encapsulation<time>> q, time t)
+    {
+      this->exact(as_vector<double, time>(q), t);
+    }
+
+    void exact(DVectorT& q, time t)
+    {
+      size_t n = q.size();
+      double a = 1.0 / sqrt(4 * PI * nu * (t + t0));
+
+      for (size_t i = 0; i < n; i++) {
+        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[i] += a * exp(-x * x / (4 * nu * (t + t0)));
+        }
+      }
+    }
+
+    void echo_error(time t, bool predict = false)
+    {
+      auto& qend = as_vector<double, time>(this->get_u_end());
+      DVectorT qex(qend.size());
+
+      this->exact(qex, t);
+
+      double max = 0.0;
+      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 << ")"
+           << endl;
+
+      this->errors.insert(pair<pair<size_t, size_t>, double>(pair<size_t, size_t>(n, k), max));
+    }
+
+    /**
+     * retrieve errors at iterations and time nodes
+     */
+    error_map get_errors()
+    {
+      return this->errors;
+    }
+    //! @}
+
     //! @{
     /**
      * @copybrief pfasst::encap::IMEXSweeper::predict()
@@ -167,59 +224,6 @@ class AdvectionDiffusionSweeper
     }
     //! @}
 
-    //! @{
-    void exact(shared_ptr<Encapsulation<time>> q, time t)
-    {
-      this->exact(as_vector<double, time>(q), t);
-    }
-
-    void exact(DVectorT& q, time t)
-    {
-      size_t n = q.size();
-      double a = 1.0 / sqrt(4 * PI * nu * (t + t0));
-
-      for (size_t i = 0; i < n; i++) {
-        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[i] += a * exp(-x * x / (4 * nu * (t + t0)));
-        }
-      }
-    }
-
-    void echo_error(time t, bool predict = false)
-    {
-      auto& qend = as_vector<double, time>(this->get_end_state());
-      DVectorT qex(qend.size());
-
-      this->exact(qex, t);
-
-      double max = 0.0;
-      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 << ")"
-           << endl;
-
-      this->errors.insert(pair<pair<size_t, size_t>, double>(pair<size_t, size_t>(n, k), max));
-    }
-
-    /**
-     * retrieve errors at iterations and time nodes
-     */
-    error_map get_errors()
-    {
-      return this->errors;
-    }
-    //! @}
 };
 
 #endif

examples/advection_diffusion/vanilla_sdc.cpp

@@ -25,18 +25,18 @@ error_map run_vanilla_sdc()
   const size_t ndofs  = 64;
   const size_t niters = 4;
 
-  auto nodes   = pfasst::compute_nodes(nnodes, pfasst::QuadratureType::GaussLobatto);
+  auto quad    = pfasst::quadrature::quadrature_factory(nnodes, pfasst::quadrature::QuadratureType::GaussLobatto);
   auto factory = make_shared<pfasst::encap::VectorFactory<double>>(ndofs);
   auto sweeper = make_shared<AdvectionDiffusionSweeper<>>(ndofs);
 
-  sweeper->set_nodes(nodes);
+  sweeper->set_quadrature(quad);
   sweeper->set_factory(factory);
 
   sdc.add_level(sweeper);
   sdc.set_duration(0.0, nsteps*dt, dt, niters);
   sdc.setup();
 
-  auto q0 = sweeper->get_state(0);
+  auto q0 = sweeper->get_u_start();
   sweeper->exact(q0, 0.0);
 
   sdc.run();

include/pfasst/encap/automagic.hpp

@@ -23,15 +23,15 @@ namespace pfasst
                            >;
 
     template<typename ControllerT, typename BuildT, typename time = time_precision>
-    void auto_build(ControllerT& c, vector<pair<size_t, pfasst::QuadratureType>> nodes, BuildT build)
+    void auto_build(ControllerT& c, vector<pair<size_t, quadrature::QuadratureType>> nodes, BuildT build)
     {
       for (size_t l = 0; l < nodes.size(); l++) {
-        auto nds = pfasst::compute_nodes<time>(get<0>(nodes[l]), get<1>(nodes[l]));
+        auto quad = quadrature::quadrature_factory<time>(get<0>(nodes[l]), get<1>(nodes[l]));
         AutoBuildTuple<time> tpl = build(l);
         auto sweeper = get<0>(tpl);
         auto transfer = get<1>(tpl);
         auto factory = get<2>(tpl);
-        sweeper->set_nodes(nds);
+        sweeper->set_quadrature(quad);
         sweeper->set_factory(factory);
         c.add_level(sweeper, transfer, false);
       }

include/pfasst/encap/encap_sweeper.hpp

@@ -18,54 +18,77 @@ namespace pfasst
 {
   namespace encap
   {
-
     template<typename time = time_precision>
     class EncapSweeper
       : public ISweeper<time>
     {
-      public:
+      protected:
         //! @{
-        typedef Encapsulation<time> encap_type;
-        typedef EncapFactory<time> factory_type;
+        quadrature::IQuadrature<time>* quad;
+        shared_ptr<EncapFactory<time>> factory;
+        shared_ptr<Encapsulation<time>> u_start;
+        shared_ptr<Encapsulation<time>> u_end;
+        shared_ptr<Encapsulation<time>> u_end_old;
         //! @}
 
-      private:
+      public:
         //! @{
-        vector<time> nodes;
-        vector<bool> is_proper;
-        shared_ptr<factory_type> factory;
+        EncapSweeper()
+          : quad(nullptr)
+        {}
+
+        virtual ~EncapSweeper()
+        {
+          if (this->quad) delete this->quad;
+        }
         //! @}
 
-      public:
         //! @{
-        virtual ~EncapSweeper()
-        {}
+        virtual void spread() override
+        {
+          for (size_t m = 1; m < this->quad->get_num_nodes(); m++) {
+            this->get_state(m)->copy(this->u_start);
+          }
+        }
         //! @}
 
         //! @{
-        virtual void set_nodes(vector<time> nodes)
+        void set_quadrature(quadrature::IQuadrature<time>* quad)
+        {
+          this->quad = quad;
+        }
+
+        const quadrature::IQuadrature<time>* get_quadrature() const
         {
-          auto augmented = pfasst::augment_nodes(nodes);
-          this->nodes = get<0>(augmented);
-          this->is_proper = get<1>(augmented);
+          return this->quad;
         }
 
-        virtual const vector<time> get_nodes() const
+        shared_ptr<Encapsulation<time>> get_u_start() const
         {
-          return nodes;
+          return this->u_start;
         }
 
-        const vector<bool> get_is_proper() const
+        shared_ptr<Encapsulation<time>> get_u_end() const
         {
-          return is_proper;
+          return this->u_end;
         }
 
-        virtual void set_factory(shared_ptr<factory_type> factory)
+        shared_ptr<const Encapsulation<time>> get_u_end_old() const
+        {
+          return this->u_end_old;
+        }
+
+        const vector<time> get_nodes() const
+        {
+          return this->quad->get_nodes();
+        }
+
+        void set_factory(shared_ptr<EncapFactory<time>> factory)
         {
           this->factory = factory;
         }
 
-        virtual shared_ptr<factory_type> get_factory() const
+        virtual shared_ptr<EncapFactory<time>> get_factory() const
         {
           return factory;
         }
@@ -78,12 +101,22 @@ namespace pfasst
          *
          * @note This method must be implemented in derived sweepers.
          */
-        virtual void set_state(shared_ptr<const encap_type> u0, size_t m)
+        virtual void set_state(shared_ptr<const Encapsulation<time>> u0, size_t m)
         {
           UNUSED(u0); UNUSED(m);
           throw NotImplementedYet("sweeper");
         }
 
+        void set_u_start(shared_ptr<const Encapsulation<time>> u_start)
+        {
+          this->u_start = u_start;
+        }
+
+        void set_u_end(shared_ptr<const Encapsulation<time>> u_end)
+        {
+          this->u_end = u_end;
+        }
+
         /**
          * retrieve solution values of current iteration at time node index `m`
          *
@@ -127,8 +160,9 @@ namespace pfasst
         }
 
         virtual shared_ptr<Encapsulation<time>> get_end_state()
+        __attribute__((deprecated("use u_end")))
         {
-          return this->get_state(this->get_nodes().size() - 1);
+          return this->get_u_end();
         }
         //! @}
 
@@ -143,13 +177,6 @@ namespace pfasst
           throw NotImplementedYet("sweeper");
         }
 
-        virtual void spread() override
-        {
-          for (size_t m = 1; m < nodes.size(); m++) {
-            this->get_state(m)->copy(this->get_state(0));
-          }
-        }
-
         /**
          * evaluates the right hand side at given time node
          *
@@ -175,7 +202,7 @@ namespace pfasst
          *
          * @note This method must be implemented in derived sweepers.
          */
-        virtual void integrate(time dt, vector<shared_ptr<encap_type>> dst) const
+        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const
         {
           UNUSED(dt); UNUSED(dst);
           throw NotImplementedYet("sweeper");