mlsdc: Don't interpolate delta of initial conditions.

Signed-off-by: Matthew Emmett <[email protected]>

Author: memmett

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -45,6 +45,11 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
     }
   }
 
+  ~AdvectionDiffusionSweeper()
+  {
+    cout << "number of f1 evals: " << nf1evals << endl;
+  }
+
   void exact(Encapsulation<scalar,time>* q, scalar t)
   {
     exact(*dynamic_cast<dvector*>(q), t);

include/pfasst/encap/encapsulation.hpp

@@ -134,7 +134,9 @@ namespace pfasst {
 
       virtual ~PolyInterpMixin() { }
 
-      virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial)
+      virtual void interpolate(ISweeper *dst, const ISweeper *src,
+			       bool interp_delta_from_initial,
+			       bool interp_initial)
       {
 	auto* fine = dynamic_cast<EncapSweeper<scalar,time>*>(dst);
 	auto* crse = dynamic_cast<const EncapSweeper<scalar,time>*>(src);
@@ -153,25 +155,29 @@ namespace pfasst {
 	for (int m=0; m<nfine; m++) fine_state[m] = fine->get_state(m);
 	for (int m=0; m<ncrse; m++) fine_delta[m] = fine_factory->create(solution);
 
-	if (initial)
+	if (interp_delta_from_initial)
 	  for (int m=1; m<nfine; m++)
 	    fine_state[m]->copy(fine_state[0]);
 
 	auto* crse_delta = crse_factory->create(solution);
-	for (int m=0; m<ncrse; m++) {
+	int m0 = interp_initial ? 0 : 1;
+	for (int m=m0; m<ncrse; m++) {
 	  crse_delta->copy(crse->get_state(m));
-	  if (initial)
+	  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));
 	  interpolate(fine_delta[m], crse_delta);
 	}
 	delete crse_delta;
 
+	if (! interp_initial)
+	  fine_delta[0]->setval(0.0);
+
 	fine->get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
 
 	for (int m=0; m<ncrse; m++) delete fine_delta[m];
-	for (int m=0; m<nfine; m++) fine->evaluate(m);
+	for (int m=m0; m<nfine; m++) fine->evaluate(m);
       }
 
       virtual void restrict(ISweeper *dst, const ISweeper *src, bool restrict_initial)

include/pfasst/interfaces.hpp

@@ -39,9 +39,13 @@ namespace pfasst {
   class ITransfer {
   public:
     // XXX: pass level iterator to these routines as well
-    virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial=false) = 0;
-    virtual void restrict(ISweeper *dst, const ISweeper *src, bool initial=false) = 0;
-    virtual void fas(double dt, ISweeper *dst, const ISweeper *src) = 0; // XXX: this needs to be templated
+    // XXX: these needs to be templated
+    virtual void interpolate(ISweeper *dst, const ISweeper *src,
+			     bool interp_delta_from_initial=false,
+			     bool interp_initial=false) = 0;
+    virtual void restrict(ISweeper *dst, const ISweeper *src,
+			  bool restrict_initial=false) = 0;
+    virtual void fas(double dt, ISweeper *dst, const ISweeper *src) = 0;
   };
 
   class ICommunicator {

include/pfasst/mlsdc.hpp

@@ -103,7 +103,7 @@ namespace pfasst {
       auto* crse = leviter.coarse();
       auto* trns = leviter.transfer();
 
-      trns->interpolate(fine, crse, false);
+      trns->interpolate(fine, crse);
 
       if (leviter < this->finest())
 	perform_sweeps(leviter, t, dt);