imex: Support virtual nodes (0 and 1).

Author: Matthew Emmett

include/pfasst/encap/imex_sweeper.hpp

@@ -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()
@@ -178,6 +201,12 @@ namespace pfasst
 
           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;
           for (size_t m = 0; m < nodes.size() - 1; m++) {
@@ -191,7 +220,12 @@ namespace pfasst
             if (coarse) {
               this->previous_u_state.push_back(this->get_factory()->create(pfasst::encap::solution));
             }
-            this->fs_expl.push_back(this->get_factory()->create(pfasst::encap::function));
+	    // XXX: can prolly save some f's here...
+	    // if ((m == 0) || this->get_is_proper()[m]) {
+	    this->fs_expl.push_back(this->get_factory()->create(pfasst::encap::function));
+	    // } else {
+	    //   this->fs_expl.push_back(nullptr);
+	    // }
             this->fs_impl.push_back(this->get_factory()->create(pfasst::encap::function));
           }
 
@@ -212,14 +246,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]);
@@ -228,6 +259,8 @@ namespace pfasst
 
             t += ds;
           }
+
+	  if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
         virtual void sweep()
@@ -251,7 +284,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]);
@@ -262,6 +295,8 @@ namespace pfasst
 
             t += ds;
           }
+
+	  if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
         virtual void advance() override
@@ -282,11 +317,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);
+	  }
         }
 
         /**