Merge pull request #59 from memmett/feature/fix-pf-interp

memmett · memmett · commit efba4be13831 · 2014-08-04T15:08:19.000-07:00
Fix PFASST predictor. Fixes #58.
diff --git a/examples/advection_diffusion/mpi_pfasst.cpp b/examples/advection_diffusion/mpi_pfasst.cpp
@@ -61,6 +61,7 @@ int main(int argc, char** argv)
 
   pf.set_comm(&comm);
   pf.set_duration(0.0, nsteps * dt, dt, niters);
+  pf.set_nsweeps({2, 1});
   pf.run();
 
   fftw_cleanup();
diff --git a/include/pfasst/encap/encap_sweeper.hpp b/include/pfasst/encap/encap_sweeper.hpp
@@ -25,6 +25,13 @@ namespace pfasst
 
       public:
 
+        virtual void spread()
+        {
+          for (size_t m = 1; m < nodes.size(); m++) {
+            this->get_state(m)->copy(this->get_state(0));
+          }
+        }
+
         virtual void post(ICommunicator* comm, int tag)
         {
           this->get_state(0)->post(comm, tag);
@@ -42,9 +49,9 @@ namespace pfasst
 
         virtual void broadcast(ICommunicator* comm)
         {
-	  if (comm->rank() == comm->size() - 1) {
-	    this->get_state(0)->copy(this->get_state(this->get_nodes().size() - 1));
-	  }
+          if (comm->rank() == comm->size() - 1) {
+            this->get_state(0)->copy(this->get_state(this->get_nodes().size() - 1));
+          }
           this->get_state(0)->broadcast(comm);
         }
 
@@ -112,6 +119,23 @@ namespace pfasst
         }
     };
 
+
+    template<typename time>
+    EncapSweeper<time>& as_encap_sweeper(shared_ptr<ISweeper<time>> x)
+    {
+      shared_ptr<EncapSweeper<time>> y = dynamic_pointer_cast<EncapSweeper<time>>(x);
+      assert(y);
+      return *y.get();
+    }
+
+    template<typename time>
+    const EncapSweeper<time>& as_encap_sweeper(shared_ptr<const ISweeper<time>> x)
+    {
+      shared_ptr<const EncapSweeper<time>> y = dynamic_pointer_cast<const EncapSweeper<time>>(x);
+      assert(y);
+      return *y.get();
+    }
+
   }
 
 }
diff --git a/include/pfasst/encap/poly_interp.hpp b/include/pfasst/encap/poly_interp.hpp
@@ -27,79 +27,74 @@ namespace pfasst
       public:
         virtual ~PolyInterpMixin() { }
 
-        virtual void interpolate(shared_ptr<ISweeper<time>> dst,
-                                 shared_ptr<const ISweeper<time>> src,
-                                 bool interp_delta_from_initial,
-                                 bool interp_initial,
-                                 bool interp_initial_only)
+        virtual void interpolate_initial(shared_ptr<ISweeper<time>> dst,
+                                         shared_ptr<const ISweeper<time>> src)
         {
-          shared_ptr<EncapSweeper<time>> fine = dynamic_pointer_cast<EncapSweeper<time>>(dst);
-          assert(fine);
-          shared_ptr<const EncapSweeper<time>> crse = dynamic_pointer_cast<const EncapSweeper<time>>(src);
-          assert(crse);
+          auto& fine = as_encap_sweeper(dst);
+          auto& crse = as_encap_sweeper(src);
+
+          auto crse_factory = crse.get_factory();
+          auto fine_factory = fine.get_factory();
+
+          auto crse_delta = crse_factory->create(solution);
+          this->restrict(crse_delta, fine.get_state(0));
+          crse_delta->saxpy(-1.0, crse.get_state(0));
+
+          auto fine_delta = fine_factory->create(solution);
+          this->interpolate(fine_delta, crse_delta);
+          fine.get_state(0)->saxpy(-1.0, fine_delta);
 
-          this->interpolate(fine, crse, interp_delta_from_initial, interp_initial, interp_initial_only);
+          fine.evaluate(0);
         }
 
-        virtual void interpolate(shared_ptr<EncapSweeper<time>> fine,
-                                 shared_ptr<const EncapSweeper<time>> crse,
-                                 bool interp_delta_from_initial,
-                                 bool interp_initial,
-                                 bool interp_initial_only)
+        virtual void interpolate(shared_ptr<ISweeper<time>> dst,
+                                 shared_ptr<const ISweeper<time>> src,
+                                 bool interp_initial)
         {
+          auto& fine = as_encap_sweeper(dst);
+          auto& crse = as_encap_sweeper(src);
+
           if (tmat.size1() == 0) {
-            tmat = pfasst::compute_interp<time>(fine->get_nodes(), crse->get_nodes());
+            tmat = pfasst::compute_interp<time>(fine.get_nodes(), crse.get_nodes());
           }
 
-          size_t nfine = fine->get_nodes().size();
-          size_t ncrse = crse->get_nodes().size();
-
-          auto crse_factory = crse->get_factory();
-          auto fine_factory = fine->get_factory();
-
-          if (interp_initial_only) {
-            auto crse_delta = crse_factory->create(solution);
-            this->restrict(crse_delta, fine->get_state(0));
-            crse_delta->saxpy(-1.0, crse->get_state(0));
+          size_t nfine = fine.get_nodes().size();
+          size_t ncrse = crse.get_nodes().size();
 
-            auto fine_delta = fine_factory->create(solution);
-            this->interpolate(fine_delta, crse_delta);
-            fine->get_state(0)->saxpy(-1.0, fine_delta);
-
-            fine->evaluate(0);
-            return;
-          }
+          auto crse_factory = crse.get_factory();
+          auto fine_factory = fine.get_factory();
 
           EncapVecT fine_state(nfine), fine_delta(ncrse);
 
-          for (size_t m = 0; m < nfine; m++) { fine_state[m] = fine->get_state(m); }
+          for (size_t m = 0; m < nfine; m++) { fine_state[m] = fine.get_state(m); }
           for (size_t m = 0; m < ncrse; m++) { fine_delta[m] = fine_factory->create(solution); }
 
-          if (interp_delta_from_initial) {
-            for (size_t m = 1; m < nfine; m++) {
-              fine_state[m]->copy(fine_state[0]);
-            }
-          }
+          // if (interp_delta_from_initial) {
+          //   for (size_t m = 1; m < nfine; m++) {
+          //     fine_state[m]->copy(fine_state[0]);
+          //   }
+          // }
 
           auto crse_delta = crse_factory->create(solution);
           size_t m0 = interp_initial ? 0 : 1;
           for (size_t m = m0; m < ncrse; m++) {
-            crse_delta->copy(crse->get_state(m));
-            if (interp_delta_from_initial) {
-              crse_delta->saxpy(-1.0, crse->get_state(0));
-            } else {
-              crse_delta->saxpy(-1.0, crse->get_saved_state(m));
-            }
+            crse_delta->copy(crse.get_state(m));
+            // if (interp_delta_from_initial) {
+            //   crse_delta->saxpy(-1.0, crse->get_saved_state(0));
+            //   // crse_delta->saxpy(-1.0, crse->get_state(0));
+            // } else {
+              crse_delta->saxpy(-1.0, crse.get_saved_state(m));
+            // }
             interpolate(fine_delta[m], crse_delta);
           }
 
           if (!interp_initial) {
             fine_delta[0]->zero();
           }
 
-          fine->get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
+          fine.get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
 
-          for (size_t m = m0; m < nfine; m++) { fine->evaluate(m); }
+          for (size_t m = m0; m < nfine; m++) { fine.evaluate(m); }
         }
 
         virtual void restrict(shared_ptr<ISweeper<time>> dst,
diff --git a/include/pfasst/interfaces.hpp b/include/pfasst/interfaces.hpp
@@ -112,29 +112,32 @@ namespace pfasst
       virtual void predict(bool initial) = 0;
 
       /**
-       * perform one SDC sweep/iteration.
-       * Compute a correction and update solution values.
-       * Note that this function can assume that valid function values exist from a previous
-       * pfasst::ISweeper::sweep() or pfasst::ISweeper::predict().
+       * Perform one SDC sweep/iteration.
+       *
+       * Compute a correction and update solution values.  Note that this function can assume that
+       * valid function values exist from a previous pfasst::ISweeper::sweep() or
+       * pfasst::ISweeper::predict().
        */
       virtual void sweep() = 0;
 
       /**
-       * advance from one time step to the next.
+       * Advance from one time step to the next.
        *
        * Essentially this means copying the solution and function values from the last node to the
        * first node.
        */
       virtual void advance() = 0;
 
       /**
-       * save solutions (and/or function values) at all nodes.
+       * Save states (and/or function values) at all nodes.
        *
-       * This is typically done in MLSDC/PFASST immediately after a call to restrict.
-       * The saved states are used to compute deltas during interpolation.
+       * This is typically done in MLSDC/PFASST immediately after a call to restrict.  The saved
+       * states are used to compute deltas during interpolation.
        */
       virtual void save(bool initial_only = false) { (void) initial_only; NotImplementedYet("mlsdc/pfasst"); }
 
+      virtual void spread() { NotImplementedYet("pfasst"); }
+
       virtual void post(ICommunicator* /*comm*/, int /*tag*/) { };
       virtual void send(ICommunicator* /*comm*/, int /*tag*/, bool /*blocking*/) { NotImplementedYet("pfasst"); }
       virtual void recv(ICommunicator* /*comm*/, int /*tag*/, bool /*blocking*/) { NotImplementedYet("pfasst"); }
@@ -151,24 +154,31 @@ namespace pfasst
   class ITransfer
   {
     public:
-      // XXX: pass level iterator to these routines as well
       virtual ~ITransfer() { }
 
+
+      /**
+       * Interpolate initial condition (in space) from the coarse sweeper to the fine sweeper.
+       */
+      virtual void interpolate_initial(shared_ptr<ISweeper<time>> dst,
+                                       shared_ptr<const ISweeper<time>> src)
+      {
+        NotImplementedYet("pfasst");
+      }
+
+
       /**
-       * interpolate, in time and space, from the coarse sweeper to the fine sweeper.
-       * @param[in] interp_delta_from_initial
-       *     `true` if the delta computed at each node should be relative to the initial condition.
+       * Interpolate, in time and space, from the coarse sweeper to the fine sweeper.
        * @param[in] interp_initial
        *     `true` if a delta for the initial condtion should also be computed (PFASST).
        */
       virtual void interpolate(shared_ptr<ISweeper<time>> dst,
                                shared_ptr<const ISweeper<time>> src,
-                               bool interp_delta_from_initial = false,
-                               bool interp_initial = false,
-                               bool interp_initial_only = false) = 0;
+                               bool interp_initial = false) = 0;
+
 
       /**
-       * restrict, in time and space, from the fine sweeper to the coarse sweeper.
+       * Restrict, in time and space, from the fine sweeper to the coarse sweeper.
        * @param[in] restrict_initial
        *     `true` if the initial condition should also be restricted.
        */
@@ -178,7 +188,7 @@ namespace pfasst
                             bool restrict_initial_only = false) = 0;
 
       /**
-       * compute FAS correction between the coarse and fine sweepers.
+       * Compute FAS correction between the coarse and fine sweepers.
        */
       virtual void fas(time dt, shared_ptr<ISweeper<time>> dst,
                        shared_ptr<const ISweeper<time>> src) = 0;
diff --git a/include/pfasst/mlsdc.hpp b/include/pfasst/mlsdc.hpp
@@ -46,11 +46,16 @@ namespace pfasst
         nsweeps.resize(this->nlevels());
         fill(nsweeps.begin(), nsweeps.end(), 1);
         for (auto leviter = this->coarsest(); leviter <= this->finest(); ++leviter) {
-	  leviter.current()->set_controller(this);
+          leviter.current()->set_controller(this);
           leviter.current()->setup(leviter != this->finest());
         }
       }
 
+      void set_nsweeps(vector<size_t> nsweeps)
+      {
+        this->nsweeps = nsweeps;
+      }
+
       /**
        * evolve ODE using MLSDC.
        *
@@ -64,7 +69,7 @@ namespace pfasst
           initial = this->get_step() == 0; // only evaluate node 0 functions on first step
 
           // iterate by performing v-cycles
-	  for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {
+          for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {
             cycle_v(this->finest());
           }
 
diff --git a/include/pfasst/pfasst.hpp b/include/pfasst/pfasst.hpp
@@ -19,7 +19,8 @@ namespace pfasst
 
       typedef typename pfasst::Controller<time>::LevelIter LevelIter;
 
-      bool predict, initial;
+      bool predict; //<! whether to use a 'predict' sweep
+      bool initial; //<! whether we're sweeping from a new initial condition
 
       void perform_sweeps(size_t level)
       {
@@ -46,13 +47,16 @@ namespace pfasst
        */
       void predictor()
       {
+        this->get_finest()->spread();
+
         // restrict fine initial condition
         for (auto l = this->finest() - 1; l >= this->coarsest(); --l) {
           auto crse = l.current();
           auto fine = l.fine();
           auto trns = l.transfer();
-          trns->restrict(crse, fine, false, true);
-          crse->save(true);
+          trns->restrict(crse, fine, true, true);
+          crse->spread();
+          crse->save();
         }
 
         // perform sweeps on the coarse level based on rank
@@ -63,7 +67,9 @@ namespace pfasst
           // XXX: set iteration?
 
           perform_sweeps(0);
-          crse->advance();
+          if (nstep < comm->rank()) {
+            crse->advance();
+          }
         }
 
         // return to finest level, sweeping as we go
@@ -72,7 +78,7 @@ namespace pfasst
           auto fine = l.current();
           auto trns = l.transfer();
 
-          trns->interpolate(fine, crse, true, true);
+          trns->interpolate(fine, crse, true);
           if (l < this->finest()) {
             perform_sweeps(l.level);
           }
@@ -128,9 +134,9 @@ namespace pfasst
 
             cycle_v(this->finest() - 1);
 
-            trns->interpolate(fine, crse, false, true, false);
+            trns->interpolate(fine, crse, true);
             fine->recv(comm, tag, false);
-            trns->interpolate(fine, crse, false, false, true);
+            trns->interpolate_initial(fine, crse);
             // XXX: call interpolate_q0(pf,F, G)
           }
 
@@ -176,12 +182,12 @@ namespace pfasst
         auto crse = l.coarse();
         auto trns = l.transfer();
 
-        trns->interpolate(fine, crse, false, true, false);
+        trns->interpolate(fine, crse, true);
 
         int tag = l.level * 10000 + this->get_iteration() + 10;
         fine->recv(comm, tag, false);
         // XXX          call interpolate_q0(pf,F, G)
-        trns->interpolate(fine, crse, false, false, true);
+        trns->interpolate_initial(fine, crse);
 
         if (l < this->finest()) {
           perform_sweeps(l.level);

Original file line number	Diff line number	Diff line change
`@@ -46,11 +46,16 @@ namespace pfasst`
`46`	`46`	`nsweeps.resize(this->nlevels());`
`47`	`47`	`fill(nsweeps.begin(), nsweeps.end(), 1);`
`48`	`48`	`for (auto leviter = this->coarsest(); leviter <= this->finest(); ++leviter) {`
`49`		`- leviter.current()->set_controller(this);`
	`49`	`+ leviter.current()->set_controller(this);`
`50`	`50`	`leviter.current()->setup(leviter != this->finest());`
`51`	`51`	`}`
`52`	`52`	`}`
`53`	`53`
	`54`	`+ void set_nsweeps(vector<size_t> nsweeps)`
	`55`	`+ {`
	`56`	`+ this->nsweeps = nsweeps;`
	`57`	`+ }`
	`58`	`+`
`54`	`59`	`/**`
`55`	`60`	`* evolve ODE using MLSDC.`
`56`	`61`	`*`
`@@ -64,7 +69,7 @@ namespace pfasst`
`64`	`69`	`initial = this->get_step() == 0; // only evaluate node 0 functions on first step`
`65`	`70`
`66`	`71`	`// iterate by performing v-cycles`
`67`		`- for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {`
	`72`	`+ for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {`
`68`	`73`	`cycle_v(this->finest());`
`69`	`74`	`}`
`70`	`75`