Merge pull request #37 from memmett/feature/state

Add unit tests for AdvectionDiffusion example and change how state is passed to sweepers.

Author: torbjoernk

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -5,10 +5,11 @@
 #ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
 #define _ADVECTION_DIFFUSION_SWEEPER_HPP_
 
-#include <complex>
-#include <vector>
 #include <cassert>
+#include <complex>
+#include <map>
 #include <ostream>
+#include <vector>
 
 #include <pfasst/encap/imex_sweeper.hpp>
 
@@ -23,13 +24,16 @@ 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>
 {
     FFT fft;
 
     vector<complex<double>> ddx, lap;
+    error_map errors;
 
     double v  = 1.0;
     time   t0 = 1.0;
@@ -87,20 +91,35 @@ class AdvectionDiffusionSweeper
         double d = abs(qend[i] - qex[i]);
         if (d > max) { max = d; }
       }
-      cout << "err: " << scientific << max
+
+      auto n = this->get_controller()->get_step();
+      auto k = this->get_controller()->get_iteration();
+      cout << "err: " << n << " " << k << " " << scientific << max
            << " (" << qend.size() << ", " << predict << ")"
            << endl;
+
+      errors.insert(pair<pair<size_t, size_t>, double>
+		    (pair<size_t, size_t>(n, k), max));
+    }
+
+    error_map get_errors()
+    {
+      return errors;
     }
 
-    void predict(time t, time dt, bool initial)
+    void predict(bool initial)
     {
-      pfasst::encap::IMEXSweeper<time>::predict(t, dt, initial);
+      pfasst::encap::IMEXSweeper<time>::predict(initial);
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
       echo_error(t + dt, true);
     }
 
-    void sweep(time t, time dt)
+    void sweep()
     {
-      pfasst::encap::IMEXSweeper<time>::sweep(t, dt);
+      pfasst::encap::IMEXSweeper<time>::sweep();
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
       echo_error(t + dt);
     }
 

examples/advection_diffusion/serial_mlsdc.cpp

@@ -18,7 +18,7 @@
 using namespace pfasst;
 using namespace pfasst::encap;
 
-int main(int /*argc*/, char** /*argv*/)
+error_map run_serial_mlsdc()
 {
   MLSDC<> mlsdc;
 
@@ -70,8 +70,17 @@ int main(int /*argc*/, char** /*argv*/)
   /*
    * run mlsdc!
    */
-  mlsdc.set_duration(dt, nsteps, niters);
+  mlsdc.set_duration(0.0, nsteps*dt, dt, niters);
   mlsdc.run();
 
   fftw_cleanup();
+
+  return sweeper->get_errors();
+}
+
+#ifndef PFASST_UNIT_TESTING
+int main(int /*argc*/, char** /*argv*/)
+{
+  run_serial_mlsdc();
 }
+#endif

examples/advection_diffusion/serial_mlsdc_autobuild.cpp

@@ -65,7 +65,7 @@ int main(int /*argc*/, char** /*argv*/)
 
   auto_build(mlsdc, nodes, build_level);
   auto_setup(mlsdc, initial);
-  mlsdc.set_duration(dt, nsteps, niters);
+  mlsdc.set_duration(0.0, nsteps*dt, dt, niters);
   mlsdc.run();
 
   fftw_cleanup();

examples/advection_diffusion/vanilla_sdc.cpp

@@ -14,7 +14,7 @@
 
 #include "advection_diffusion_sweeper.hpp"
 
-int main(int /*argc*/, char** /*argv*/)
+error_map run_vanilla_sdc()
 {
   pfasst::SDC<> sdc;
 
@@ -32,7 +32,7 @@ int main(int /*argc*/, char** /*argv*/)
   sweeper->set_factory(factory);
 
   sdc.add_level(sweeper);
-  sdc.set_duration(dt, nsteps, niters);
+  sdc.set_duration(0.0, nsteps*dt, dt, niters);
   sdc.setup();
 
   auto q0 = sweeper->get_state(0);
@@ -41,4 +41,14 @@ int main(int /*argc*/, char** /*argv*/)
   sdc.run();
 
   fftw_cleanup();
+
+  return sweeper->get_errors();
+}
+
+
+#ifndef PFASST_UNIT_TESTING
+int main(int /*argc*/, char** /*argv*/)
+{
+  run_vanilla_sdc();
 }
+#endif

include/pfasst/controller.hpp

@@ -26,26 +26,30 @@ namespace pfasst
       deque<shared_ptr<ISweeper<time>>>  levels;
       deque<shared_ptr<ITransfer<time>>> transfer;
 
-      size_t nsteps, niters;
-      time   dt;
+      size_t step, iteration, max_iterations;
+      time t, dt, tend;
 
     public:
       //! @{
       void setup()
       {
         for (auto l = coarsest(); l <= finest(); ++l) {
+	  l.current()->set_controller(this);
           l.current()->setup();
         }
       }
 
-      void set_duration(time dt, size_t nsteps, size_t niters)
+      void set_duration(time t0, time tend, time dt, size_t niters)
       {
+        this->t = t0;
+        this->tend = tend;
         this->dt = dt;
-        this->nsteps = nsteps;
-        this->niters = niters;
+        this->step = 0;
+        this->iteration = 0;
+        this->max_iterations = niters;
       }
 
-      void add_level(shared_ptr<ISweeper<time>> swpr, 
+      void add_level(shared_ptr<ISweeper<time>> swpr,
                      shared_ptr<ITransfer<time>> trnsfr = shared_ptr<ITransfer<time>>(nullptr),
                      bool coarse = true)
       {
@@ -63,32 +67,26 @@ namespace pfasst
       template<typename R = ISweeper<time>>
       shared_ptr<R> get_level(size_t level)
       {
-        shared_ptr<R> r = dynamic_pointer_cast<R>(levels[level]);
-        assert(r);
+        shared_ptr<R> r = dynamic_pointer_cast<R>(levels[level]); assert(r);
         return r;
       }
 
       template<typename R = ISweeper<time>>
       shared_ptr<R> get_finest()
       {
-        shared_ptr<R> r = dynamic_pointer_cast<R>(levels.back());
-        assert(r);
-        return r;
+	return get_level<R>(nlevels()-1);
       }
 
       template<typename R = ISweeper<time>>
       shared_ptr<R> get_coarsest()
       {
-        shared_ptr<R> r = dynamic_pointer_cast<R>(levels.front());
-        assert(r);
-        return r;
+	return get_level<R>(0);
       }
 
       template<typename R = ITransfer<time>>
       shared_ptr<R> get_transfer(size_t level)
       {
-        shared_ptr<R> r = dynamic_pointer_cast<R>(transfer[level]);
-        assert(r);
+        shared_ptr<R> r = dynamic_pointer_cast<R>(transfer[level]); assert(r);
         return r;
       }
 
@@ -100,13 +98,13 @@ namespace pfasst
 
       /**
        * level (MLSDC/PFASST) iterator.
-       * 
+       *
        * This iterator is used to walk through the MLSDC/PFASST hierarchy of sweepers.
-       * It keeps track of the _current_ level, and has convenience routines to return the 
+       * It keeps track of the _current_ level, and has convenience routines to return the
        * LevelIter::current(), LevelIter::fine() (i.e. `current+1`), and LevelIter::coarse()
        * (`current-1`) sweepers.
-       * 
-       * Under the hood it satisfies the requirements of std::random_access_iterator_tag, thus 
+       *
+       * Under the hood it satisfies the requirements of std::random_access_iterator_tag, thus
        * implementing a `RandomAccessIterator`.
        */
       class LevelIter
@@ -179,6 +177,57 @@ namespace pfasst
       LevelIter finest()   { return LevelIter(nlevels() - 1, this); }
       LevelIter coarsest() { return LevelIter(0, this); }
       //! @}
+
+
+      /**
+       * Get current time step number.
+       */
+      size_t get_step()
+      {
+        return step;
+      }
+
+      time get_time_step()
+      {
+        return dt;
+      }
+
+      time get_time()
+      {
+        return t;
+      }
+
+      void advance_time(size_t nsteps=1)
+      {
+        step += nsteps;
+        t += nsteps*dt;
+      }
+
+      time get_end_time()
+      {
+        return tend;
+      }
+
+      size_t get_iteration()
+      {
+        return iteration;
+      }
+
+      void set_iteration(size_t iter)
+      {
+        this->iteration = iter;
+      }
+
+      void advance_iteration()
+      {
+        iteration++;
+      }
+
+      size_t get_max_iterations()
+      {
+        return max_iterations;
+      }
+
   };
 }
 

include/pfasst/encap/imex_sweeper.hpp

@@ -94,12 +94,15 @@ namespace pfasst
           }
         }
 
-        virtual void sweep(time t0, time dt)
+        virtual void sweep()
         {
           const auto   nodes  = this->get_nodes();
           const size_t nnodes = nodes.size();
           assert(nnodes >= 1);
 
+	  time dt = this->get_controller()->get_time_step();
+          time t  = this->get_controller()->get_time();
+
           // integrate
           S[0]->mat_apply(S, dt, SEmat, Fe, true);
           S[0]->mat_apply(S, dt, SImat, Fi, false);
@@ -112,7 +115,6 @@ namespace pfasst
           // sweep
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          time t = t0;
           for (size_t m = 0; m < nnodes - 1; m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
 
@@ -126,20 +128,22 @@ namespace pfasst
           }
         }
 
-        virtual void predict(time t0, time dt, bool initial)
+        virtual void predict(bool initial)
         {
           const auto   nodes  = this->get_nodes();
           const size_t nnodes = nodes.size();
           assert(nnodes >= 1);
 
+	  time dt = this->get_controller()->get_time_step();
+          time t  = this->get_controller()->get_time();
+
           if (initial) {
-            f1eval(Fe[0], Q[0], t0);
-            f2eval(Fi[0], Q[0], t0);
+            f1eval(Fe[0], Q[0], t);
+            f2eval(Fi[0], Q[0], t);
           }
 
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          time t = t0;
           for (size_t m = 0; m < nnodes - 1; m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
             rhs->copy(Q[m]);
@@ -165,19 +169,19 @@ namespace pfasst
           f2eval(Fi[m], Q[m], t);
         }
 
-        virtual void f1eval(shared_ptr<Encapsulation<time>> /*q*/, shared_ptr<Encapsulation<time>> /*f*/,
+        virtual void f1eval(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
 			    time /*t*/)
         {
           throw NotImplementedYet("imex (f1eval)");
         }
 
-        virtual void f2eval(shared_ptr<Encapsulation<time>> /*q*/, shared_ptr<Encapsulation<time>> /*f*/,
+        virtual void f2eval(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
 			    time /*t*/)
         {
           throw NotImplementedYet("imex (f2eval)");
         }
 
-        virtual void f2comp(shared_ptr<Encapsulation<time>> /*q*/, shared_ptr<Encapsulation<time>> /*f*/,
+        virtual void f2comp(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
                             time /*t*/, time /*dt*/,
                             shared_ptr<Encapsulation<time>> /*rhs*/)
         {