Merge branch 'development' into feature/cleanup-scalar

Signed-off-by: Torbjörn Klatt <[email protected]>

Conflicts:
	examples/scalar/scalar_sweeper.hpp

Author: torbjoernk

CHANGELOG.md

@@ -1,5 +1,40 @@
 # Changelog                                                                        {#page_changelog}
 
+## v0.2.0 -- MPI PFASST (2014/XX/XX)
+
+XXX: DOI: [10.5281/zenodo.11047](http://dx.doi.org/10.5281/zenodo.11047)
+
+### Notable Features
+
+* Addition of MPI based PFASST.
+
+### Details
+
+* Addition of MPI based PFASST algorithm using the standard predictor
+  stage in a block mode with fixed iterations.
+  ([#46][], [#57][], [#59][])
+
+* Addition of a simple scalar example.
+  ([#61][])
+
+* Various tidying
+  ()
+
+[#46]: https://github.com/Parallel-in-Time/PFASST/pull/46
+[#57]: https://github.com/Parallel-in-Time/PFASST/pull/56
+[#59]: https://github.com/Parallel-in-Time/PFASST/pull/59
+[#61]: https://github.com/Parallel-in-Time/PFASST/pull/61
+
+### Contributors
+
+* Matthew Emmett, Lawrence Berkeley National Laboratory ([memmett][])
+* Torbjörn Klatt, Jülich Supercomputing Centre ([torbjoernk][])
+* Daniel Ruprecht, Institute of Computational Science, University of Lugano ([danielru][])
+
+[memmett]: https://github.com/memmett
+[torbjoernk]: https://github.com/torbjoernk
+[danielru]: https://github.com/danielru
+
 ## v0.1.0 -- First Release (2014/07/25)
 
 DOI: [10.5281/zenodo.11047](http://dx.doi.org/10.5281/zenodo.11047)

README.md

@@ -1,18 +1,16 @@
 PFASST                                                                                   {#mainpage}
 ======
 
-The PFASST algorithm is a time-parallel algorithm for solving ODEs and PDEs.
-
-The PFASST project is a C++ implementation of the *parallel full approximation scheme in space and 
-time* (PFASST) algorithm.
-It also contains basic implementations of the *spectral deferred correction* (SDC) and 
-*multi-level spectral deferred correction* (MLSDC) algorithms.
+The PFASST project is a C++ implementation of the *parallel full approximation scheme in space and
+time* (PFASST) algorithm, which in turn is a time-parallel algorithm for solving ODEs and PDEs.  It
+also contains basic implementations of the *spectral deferred correction* (SDC) and *multi-level
+spectral deferred correction* (MLSDC) algorithms.
 
 
 News
 ----
 
-* July 25, 2014: PFASST v0.1.0 released.  Please see the [release notes](#releases) for more 
+* July 25, 2014: PFASST v0.1.0 released.  Please see the [release notes](#releases) for more
   information.
 
 
@@ -39,8 +37,8 @@ Releases
 
 * **v0.1.0** First Release (2014/07/25)
 
-  Initial release with basic implementations of SDC and MLSDC.  
-  DOI: [10.5281/zenodo.11047][DOI_v010]  
+  Initial release with basic implementations of SDC and MLSDC.
+  DOI: [10.5281/zenodo.11047][DOI_v010]
   See \subpage #page_changelog "the Changelog" for details.
 
 [DOI_v010]: http://dx.doi.org/10.5281/zenodo.11047

examples/scalar/scalar_sdc.cpp

@@ -1,3 +1,10 @@
+/*
+ * Scalar test equation example using an encapsulated IMEX sweeper.
+ * Solves Dahlquist's test equation
+ * y' = a y + i b y
+ * treating the real part implicitly and the imaginary part explicitly.
+ */
+ 
 #include<complex>
 
 #include<pfasst.hpp>
@@ -6,32 +13,45 @@
 
 #include "scalar_sweeper.hpp"
 
-int main(int, char**)
+/*
+ * A simple run routine with preset parameters
+ */
+
+double run_scalar_sdc(const size_t nsteps, const double dt, const size_t nnodes,
+  const size_t niters, const complex<double> lambda)
 {
+
   pfasst::SDC<> sdc;
 
-  const size_t nsteps = 2;
-  const double dt     = 1.0;
-  const size_t nnodes = 4;
-  const size_t niters = 6;
-  const complex<double> lambda = complex<double>(-1.0, 1.0);
   const complex<double> y0     = complex<double>(1.0, 0.0);
-
-  // so far, only lobatto nodes appear to be working
-  auto nodes   = pfasst::compute_nodes(nnodes, "gauss-lobatto");
+  auto nodes = pfasst::compute_nodes(nnodes, "gauss-lobatto");
   auto factory = make_shared<pfasst::encap::VectorFactory<complex<double>>>(1);
-  auto sweeper = make_shared<ScalarSweeper<>>(lambda, y0);
-
+  auto sweeper = make_shared<ScalarSweeper<>>(lambda,y0);
+  
   sweeper->set_nodes(nodes);
   sweeper->set_factory(factory);
-
+  
   sdc.add_level(sweeper);
-  sdc.set_duration(0.0, dt * nsteps, dt, niters);
+  sdc.set_duration(0.0, dt*nsteps, dt, niters);
   sdc.setup();
-
+  
   auto q0 = sweeper->get_state(0);
   sweeper->exact(q0, 0.0);
-
+  
   sdc.run();
+  
+  return sweeper->get_errors();
+}
 
+#ifndef PFASST_UNIT_TESTING
+int main(int /*argc*/, char** /*argv*/)
+{
+  const size_t nsteps = 2;
+  const double dt     = 1.0;
+  const size_t nnodes = 4;
+  const size_t niters = 6;
+  const complex<double> lambda = complex<double>(-1.0, 1.0);
+  
+  run_scalar_sdc(nsteps, dt, nnodes, niters, lambda);
 }
+#endif

examples/scalar/scalar_sweeper.hpp

@@ -24,19 +24,21 @@ class ScalarSweeper
     complex<double> _lambda, _u0;
     int _n_f_expl_eval, _n_f_impl_eval, _n_impl_solve;
     const complex<double> i_complex = complex<double>(0, 1);
+    double _error;
 
   public:
-
     ScalarSweeper(complex<double> lambda, complex<double> u0)
       :   _lambda(lambda)
         , _u0(u0)
         , _n_f_expl_eval(0)
         , _n_f_impl_eval(0)
         , _n_impl_solve(0)
+        , _error(0.0)
     {}
 
     virtual ~ScalarSweeper()
     {
+      cout << "Final error:               " << scientific << this->_error << endl;
       cout << "Number of explicit evaluations: " << this->_n_f_expl_eval << endl;
       cout << "Number of implicit evaluations: " << this->_n_f_impl_eval << endl;
       cout << "Number of implicit solves: " << this->_n_impl_solve << endl;
@@ -51,6 +53,12 @@ class ScalarSweeper
       this->exact(qex, t);
       double max_err = abs(qend[0] - qex[0]) / abs(qex[0]);
       cout << "err: " << scientific << max_err << endl;
+      this->_error = max_err;
+    }
+    
+    double get_errors()
+    {
+      return this->_error;    
     }
 
     void predict(bool initial) override

include/pfasst/encap/encap_sweeper.hpp

@@ -25,6 +25,7 @@ namespace pfasst
     {
         //! @{
         vector<time> nodes;
+        vector<bool> is_proper;
         shared_ptr<EncapFactory<time>> factory;
         //! @}
 
@@ -69,14 +70,21 @@ namespace pfasst
         //! @{
         void set_nodes(vector<time> nodes)
         {
-          this->nodes = nodes;
+          auto augmented = pfasst::augment_nodes(nodes);
+          this->nodes = get<0>(augmented);
+          this->is_proper = get<1>(augmented);
         }
 
         const vector<time> get_nodes() const
         {
           return nodes;
         }
 
+        const vector<bool> get_is_proper() const
+        {
+          return is_proper;
+        }
+
         void set_factory(shared_ptr<EncapFactory<time>> factory)
         {
           this->factory = shared_ptr<EncapFactory<time>>(factory);

include/pfasst/encap/imex_sweeper.hpp

@@ -30,9 +30,9 @@ namespace pfasst
      * explicitly and the stiff part implicitly.
      * Therefore, we define the splitting \\( F(t,u) = F_{expl}(t,u) + F_{impl}(t,u) \\).
      *
-     * This sweeper requires three interfaces to be implement for such ODEs: two routines to 
-     * evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a 
-     * given state, and one that solves (perhaps with an external solver) the backward-Euler 
+     * This sweeper requires three interfaces to be implement for such ODEs: two routines to
+     * evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a
+     * given state, and one that solves (perhaps with an external solver) the backward-Euler
      * equation \\( U^{n+1} - \\Delta t F_{\\rm impl}(U^{n+1}) = RHS \\) for \\( U^{n+1} \\).
      *
      * @tparam time precision type of the time dimension
@@ -86,6 +86,11 @@ namespace pfasst
          */
         matrix<time> s_mat;
 
+        /**
+         * quadrature matrix containing weights for 0-to-last node integration
+         */
+        matrix<time> b_mat;
+
         /**
          * quadrature matrix containing weights for node-to-node integration of explicit part
          *
@@ -101,6 +106,24 @@ namespace pfasst
         matrix<time> s_mat_impl;
         //! @}
 
+
+        /**
+        * Set end state to \\( U_0 + \\int F_{expl} + F_{expl} \\).
+        */
+        void integrate_end_state(time dt)
+        {
+          vector<shared_ptr<Encapsulation<time>>> dst = { this->u_state.back() };
+          dst[0]->copy(this->u_state.front());
+          dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_expl, false);
+          dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_impl, false);
+        }
+
+        inline bool last_node_is_virtual()
+        {
+          return !this->get_is_proper().back();
+        }
+
+
       public:
         //! @{
         virtual ~IMEXSweeper()
@@ -173,9 +196,16 @@ namespace pfasst
         virtual void setup(bool coarse) override
         {
           auto nodes = this->get_nodes();
+          auto is_proper = this->get_is_proper();
           assert(nodes.size() >= 1);
 
-          this->s_mat = compute_quadrature(nodes, nodes, 's');
+          this->s_mat = compute_quadrature(nodes, nodes, is_proper, 's');
+
+          auto q_mat = compute_quadrature(nodes, nodes, is_proper, 'q');
+          this->b_mat = matrix<time>(1, nodes.size());
+          for (size_t m = 0; m < nodes.size(); m++) {
+            this->b_mat(0, m) = q_mat(nodes.size() - 2, m);
+          }
 
           this->s_mat_expl = this->s_mat;
           this->s_mat_impl = this->s_mat;
@@ -190,6 +220,7 @@ namespace pfasst
             if (coarse) {
               this->previous_u_state.push_back(this->get_factory()->create(pfasst::encap::solution));
             }
+            // XXX: can prolly save some f's here for virtual nodes...
             this->fs_expl.push_back(this->get_factory()->create(pfasst::encap::function));
             this->fs_impl.push_back(this->get_factory()->create(pfasst::encap::function));
           }
@@ -211,14 +242,11 @@ namespace pfasst
           time dt = this->get_controller()->get_time_step();
           time t  = this->get_controller()->get_time();
 
-          if (initial) {
-            this->f_expl_eval(this->fs_expl[0], this->u_state[0], t);
-            this->f_impl_eval(this->fs_impl[0], this->u_state[0], t);
-          }
+          if (initial) { this->evaluate(0); }
 
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          for (size_t m = 0; m < nnodes - 1; m++) {
+          for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
             rhs->copy(this->u_state[m]);
             rhs->saxpy(ds, this->fs_expl[m]);
@@ -227,6 +255,8 @@ namespace pfasst
 
             t += ds;
           }
+
+          if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
         virtual void sweep()
@@ -250,7 +280,7 @@ namespace pfasst
           // sweep
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          for (size_t m = 0; m < nnodes - 1; m++) {
+          for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
 
             rhs->copy(this->u_state[m]);
@@ -261,6 +291,8 @@ namespace pfasst
 
             t += ds;
           }
+
+          if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
         virtual void advance() override
@@ -281,11 +313,23 @@ namespace pfasst
           }
         }
 
+        /**
+         * @copybrief EncapSweeper::evaluate()
+         *
+        * If the node `m` is virtual (not proper), we can save some
+        * implicit/explicit evaluations.
+        */
         virtual void evaluate(size_t m) override
         {
-          time t = this->get_nodes()[m]; // XXX
-          this->f_expl_eval(this->fs_expl[m], this->u_state[m], t);
-          this->f_impl_eval(this->fs_impl[m], this->u_state[m], t);
+          time t0 = this->get_controller()->get_time();
+          time dt = this->get_controller()->get_time_step();
+          time t =  t0 + dt * this->get_nodes()[m];
+          if ((m == 0) || this->get_is_proper()[m]) {
+            this->f_expl_eval(this->fs_expl[m], this->u_state[m], t);
+          }
+          if (this->get_is_proper()[m]) {
+            this->f_impl_eval(this->fs_impl[m], this->u_state[m], t);
+          }
         }
 
         /**