Merge branch 'matrix-pfasst'

Author: pancetta

.travis.yml

@@ -26,7 +26,7 @@ script:
   - "./docs/update_apidocs.sh"
   - which nosetests
   - nosetests --with-coverage --cover-erase --cover-inclusive --cover-package=pySDC/core,pySDC/implementations,pySDC/helpers,pySDC/tutorial --cover-html --cover-html-dir=target/doc/build/test_coverage pySDC/tests
-  - if [[ $PYTHON_VERSION == 3 ]]; then travis-sphinx -n build; fi
+  - if [[ $PYTHON_VERSION == 3 ]]; then travis-sphinx build; fi
 
 after_success:
   - if [[ $PYTHON_VERSION == 3 ]]; then travis-sphinx deploy; fi

docs/source/index.rst

@@ -84,7 +84,7 @@ Projects
    projects/node_failure.rst
    projects/fwsw.rst
    projects/RDC.rst
-..   projects/asymp_conv.rst
+   projects/asymp_conv.rst
 
 
 Tests

docs/source/projects/doc_matrixPFASST_matrix.rst

@@ -0,0 +1,7 @@
+Full code: `pySDC/projects/matrixPFASST/compare_to_matrixbased.py <https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/projects/matrixPFASST/compare_to_matrixbased.py>`_
+
+.. literalinclude:: ../../../pySDC/projects/matrixPFASST/compare_to_matrixbased.py
+
+Results:
+
+.. literalinclude:: ../../../comparison_matrix_vs_nomat_detail.txt

docs/source/projects/doc_matrixPFASST_propagator.rst

@@ -0,0 +1,7 @@
+Full code: `pySDC/projects/matrixPFASST/compare_to_propagator.py <https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/projects/matrixPFASST/compare_to_propagator.py>`_
+
+.. literalinclude:: ../../../pySDC/projects/matrixPFASST/compare_to_propagator.py
+
+Results:
+
+.. literalinclude:: ../../../comparison_matrix_vs_propagator_detail.txt

docs/source/projects/matrixPFASST.rst

@@ -0,0 +1 @@
+.. include:: /../../pySDC/projects/matrixPFASST/README.rst

pySDC/core/BaseTransfer.py

@@ -2,6 +2,7 @@
 import logging
 
 from pySDC.helpers.pysdc_helper import FrozenClass
+from scipy import interpolate
 import pySDC.helpers.transfer_helper as th
 from pySDC.core.Errors import UnlockError
 
@@ -49,33 +50,11 @@ def __init__(self, pars):
         self.fine = fine_level
         self.coarse = coarse_level
 
-        # for Q-based transfer, check if we need to add 0 to the list of nodes
-        if not self.fine.sweep.coll.left_is_node:
-            fine_grid = np.concatenate(([0], self.fine.sweep.coll.nodes))
-            coarse_grid = np.concatenate(([0], self.coarse.sweep.coll.nodes))
-        else:
-            fine_grid = self.fine.sweep.coll.nodes
-            coarse_grid = self.coarse.sweep.coll.nodes
-
-        if self.params.coll_iorder > len(coarse_grid):
-            self.logger.warning('requested order of Q-interpolation is not valid, resetting to %s' % len(coarse_grid))
-            self.params.coll_iorder = len(coarse_grid)
-        if self.params.coll_rorder != 1:
-            self.logger.warning('requested order of Q-restriction is != 1, can lead to weird behavior!')
-
-        # set up preliminary transfer matrices for Q-based coarsening
-        Pcoll = th.interpolation_matrix_1d(fine_grid, coarse_grid, k=self.params.coll_iorder, pad=0).toarray()
-        Rcoll = th.restriction_matrix_1d(fine_grid, coarse_grid, k=self.params.coll_rorder, pad=0).toarray()
-
-        # pad transfer matrices if necessary
-        if self.fine.sweep.coll.left_is_node:
-            self.Pcoll = np.zeros((self.fine.sweep.coll.num_nodes + 1, self.coarse.sweep.coll.num_nodes + 1))
-            self.Rcoll = np.zeros((self.coarse.sweep.coll.num_nodes + 1, self.fine.sweep.coll.num_nodes + 1))
-            self.Pcoll[1:, 1:] = Pcoll
-            self.Rcoll[1:, 1:] = Rcoll
-        else:
-            self.Pcoll = Pcoll
-            self.Rcoll = Rcoll
+        fine_grid = self.fine.sweep.coll.nodes
+        coarse_grid = self.coarse.sweep.coll.nodes
+
+        self.Pcoll = th.interpolation_matrix_1d(fine_grid, coarse_grid, k=self.params.coll_iorder, pad=0).toarray()
+        self.Rcoll = th.restriction_matrix_1d(fine_grid, coarse_grid, k=self.params.coll_rorder, pad=0).toarray()
 
         # set up spatial transfer
         self.space_transfer = space_transfer_class(fine_prob=self.fine.prob, coarse_prob=self.coarse.prob,
@@ -105,16 +84,16 @@ def restrict(self):
             raise UnlockError('fine level is still locked, cannot use data from there')
 
         # restrict fine values in space
-        tmp_u = [self.space_transfer.restrict(F.u[0])]
+        tmp_u = []
         for m in range(1, SF.coll.num_nodes + 1):
             tmp_u.append(self.space_transfer.restrict(F.u[m]))
 
         # restrict collocation values
-        G.u[0] = tmp_u[0]
+        G.u[0] = self.space_transfer.restrict(F.u[0])
         for n in range(1, SG.coll.num_nodes + 1):
-            G.u[n] = self.Rcoll[n, 0] * tmp_u[0]
-            for m in range(1, SF.coll.num_nodes + 1):
-                G.u[n] += self.Rcoll[n, m] * tmp_u[m]
+            G.u[n] = self.Rcoll[n - 1, 0] * tmp_u[0]
+            for m in range(1, SF.coll.num_nodes):
+                G.u[n] += self.Rcoll[n - 1, m] * tmp_u[m]
 
         # re-evaluate f on coarse level
         G.f[0] = PG.eval_f(G.u[0], G.time)
@@ -129,15 +108,15 @@ def restrict(self):
 
         # restrict fine level tau correction part in space
         tmp_tau = []
-        for m in range(0, SF.coll.num_nodes):
+        for m in range(SF.coll.num_nodes):
             tmp_tau.append(self.space_transfer.restrict(tauF[m]))
 
         # restrict fine level tau correction part in collocation
         tauFG = []
         for n in range(1, SG.coll.num_nodes + 1):
-            tauFG.append(self.Rcoll[n, 1] * tmp_tau[0])
+            tauFG.append(self.Rcoll[n - 1, 0] * tmp_tau[0])
             for m in range(1, SF.coll.num_nodes):
-                tauFG[-1] += self.Rcoll[n, m + 1] * tmp_tau[m]
+                tauFG[-1] += self.Rcoll[n - 1, m] * tmp_tau[m]
 
         # build tau correction
         for m in range(SG.coll.num_nodes):
@@ -146,14 +125,13 @@ def restrict(self):
         if F.tau is not None:
             # restrict possible tau correction from fine in space
             tmp_tau = []
-            for m in range(0, SF.coll.num_nodes):
+            for m in range(SF.coll.num_nodes):
                 tmp_tau.append(self.space_transfer.restrict(F.tau[m]))
 
             # restrict possible tau correction from fine in collocation
-            for n in range(1, SG.coll.num_nodes + 1):
-                G.tau[n - 1] += self.Rcoll[n, 1] * tmp_tau[0]
-                for m in range(1, SF.coll.num_nodes):
-                    G.tau[n - 1] += self.Rcoll[n, m + 1] * tmp_tau[m]
+            for n in range(SG.coll.num_nodes):
+                for m in range(SF.coll.num_nodes):
+                    G.tau[n] += self.Rcoll[n, m] * tmp_tau[m]
         else:
             pass
 
@@ -203,8 +181,8 @@ def prolong(self):
         # interpolate values in collocation
         F.u[0] += tmp_u[0]
         for n in range(1, SF.coll.num_nodes + 1):
-            for m in range(0, SG.coll.num_nodes + 1):
-                F.u[n] += self.Pcoll[n, m] * tmp_u[m]
+            for m in range(1, SG.coll.num_nodes + 1):
+                F.u[n] += self.Pcoll[n - 1, m - 1] * tmp_u[m]
 
         # re-evaluate f on fine level
         F.f[0] = PF.eval_f(F.u[0], F.time)
@@ -250,8 +228,8 @@ def prolong_f(self):
         F.u[0] += tmp_u[0]
         F.f[0] += tmp_f[0]
         for n in range(1, SF.coll.num_nodes + 1):
-            for m in range(0, SG.coll.num_nodes + 1):
-                F.u[n] += self.Pcoll[n, m] * tmp_u[m]
-                F.f[n] += self.Pcoll[n, m] * tmp_f[m]
+            for m in range(1, SG.coll.num_nodes + 1):
+                F.u[n] += self.Pcoll[n - 1, m - 1] * tmp_u[m]
+                F.f[n] += self.Pcoll[n - 1, m - 1] * tmp_f[m]
 
         return None

pySDC/implementations/controller_classes/allinclusive_multigrid_MPI.py

@@ -240,7 +240,7 @@ def pfasst(self, comm, num_procs):
             if len(self.S.levels) > 1 and self.params.predict:  # MLSDC or PFASST with predict
                 self.S.status.stage = 'PREDICT'
             else:
-                self.S.status.stage = 'IT_CHECK'
+                self.S.status.stage = 'IT_FINE'
 
         elif stage == 'PREDICT':
 
@@ -250,7 +250,7 @@ def pfasst(self, comm, num_procs):
 
             # update stage
             self.hooks.pre_iteration(step=self.S, level_number=0)
-            self.S.status.stage = 'IT_CHECK'
+            self.S.status.stage = 'IT_FINE'
 
         elif stage == 'IT_CHECK':
 
@@ -286,7 +286,10 @@ def pfasst(self, comm, num_procs):
                 # increment iteration count here (and only here)
                 self.S.status.iter += 1
                 self.hooks.pre_iteration(step=self.S, level_number=0)
-                self.S.status.stage = 'IT_FINE'
+                if len(self.S.levels) > 1:  # MLSDC or PFASST
+                    self.S.status.stage = 'IT_UP'
+                else:  # SDC
+                    self.S.status.stage = 'IT_FINE'
 
             else:
                 self.S.levels[0].sweep.compute_end_point()
@@ -295,20 +298,36 @@ def pfasst(self, comm, num_procs):
 
         elif stage == 'IT_FINE':
 
+            nsweeps = self.S.levels[0].params.nsweeps
+
             # do fine sweep
+            for k in range(nsweeps):
 
-            # standard sweep workflow: update nodes, compute residual, log progress
-            self.hooks.pre_sweep(step=self.S, level_number=0)
-            for k in range(self.S.levels[0].params.nsweeps):
+                self.hooks.pre_sweep(step=self.S, level_number=0)
                 self.S.levels[0].sweep.update_nodes()
-            self.S.levels[0].sweep.compute_residual()
-            self.hooks.post_sweep(step=self.S, level_number=0)
+
+                req_send = None
+                self.S.levels[0].sweep.compute_end_point()
+                if not self.S.status.last and self.params.fine_comm:
+                    self.logger.debug('send data: process %s, stage %s, time %s, target %s, tag %s, iter %s' %
+                                      (self.S.status.slot, self.S.status.stage, self.S.time, self.S.next,
+                                       0, self.S.status.iter))
+                    req_send = comm.isend(self.S.levels[0].uend, dest=self.S.next, tag=self.S.status.iter)
+
+                if not self.S.status.first and self.params.fine_comm:
+                    self.logger.debug('recv data: process %s, stage %s, time %s, source %s, tag %s, iter %s' %
+                                      (self.S.status.slot, self.S.status.stage, self.S.time, self.S.prev,
+                                       0, self.S.status.iter))
+                    self.recv(target=self.S.levels[0], source=self.S.prev, tag=self.S.status.iter, comm=comm)
+
+                if not self.S.status.last and self.params.fine_comm:
+                    req_send.wait()
+
+                self.S.levels[0].sweep.compute_residual()
+                self.hooks.post_sweep(step=self.S, level_number=0)
 
             # update stage
-            if len(self.S.levels) > 1:  # MLSDC or PFASST
-                self.S.status.stage = 'IT_UP'
-            else:  # SDC
-                self.S.status.stage = 'IT_CHECK'
+            self.S.status.stage = 'IT_CHECK'
 
         elif stage == 'IT_UP':
 
@@ -318,11 +337,33 @@ def pfasst(self, comm, num_procs):
 
             # sweep and send on middle levels (not on finest, not on coarsest, though)
             for l in range(1, len(self.S.levels) - 1):
-                self.hooks.pre_sweep(step=self.S, level_number=l)
-                for k in range(self.S.levels[l].params.nsweeps):
+
+                nsweeps = self.S.levels[l].params.nsweeps
+
+                for k in range(nsweeps):
+
+                    self.hooks.pre_sweep(step=self.S, level_number=l)
                     self.S.levels[l].sweep.update_nodes()
-                self.S.levels[l].sweep.compute_residual()
-                self.hooks.post_sweep(step=self.S, level_number=l)
+
+                    req_send = None
+                    self.S.levels[l].sweep.compute_end_point()
+                    if not self.S.status.last and self.params.fine_comm:
+                        self.logger.debug('send data: process %s, stage %s, time %s, target %s, tag %s, iter %s' %
+                                          (self.S.status.slot, self.S.status.stage, self.S.time, self.S.next,
+                                           l, self.S.status.iter))
+                        req_send = comm.isend(self.S.levels[l].uend, dest=self.S.next, tag=self.S.status.iter)
+
+                    if not self.S.status.first and self.params.fine_comm:
+                        self.logger.debug('recv data: process %s, stage %s, time %s, source %s, tag %s, iter %s' %
+                                          (self.S.status.slot, self.S.status.stage, self.S.time, self.S.prev,
+                                           l, self.S.status.iter))
+                        self.recv(target=self.S.levels[l], source=self.S.prev, tag=self.S.status.iter, comm=comm)
+
+                    if not self.S.status.last and self.params.fine_comm:
+                        req_send.wait()
+
+                    self.S.levels[l].sweep.compute_residual()
+                    self.hooks.post_sweep(step=self.S, level_number=l)
 
                 # transfer further up the hierarchy
                 self.S.transfer(source=self.S.levels[l], target=self.S.levels[l + 1])
@@ -343,8 +384,10 @@ def pfasst(self, comm, num_procs):
 
             # do the sweep
             self.hooks.pre_sweep(step=self.S, level_number=len(self.S.levels) - 1)
-            for k in range(self.S.levels[-1].params.nsweeps):
-                self.S.levels[-1].sweep.update_nodes()
+            assert self.S.levels[-1].params.nsweeps == 1, \
+                'ERROR: this controller can only work with one sweep on the coarse level, got %s' % \
+                self.S.levels[-1].params.nsweeps
+            self.S.levels[-1].sweep.update_nodes()
             self.S.levels[-1].sweep.compute_residual()
             self.hooks.post_sweep(step=self.S, level_number=len(self.S.levels) - 1)
             self.S.levels[-1].sweep.compute_end_point()
@@ -370,32 +413,52 @@ def pfasst(self, comm, num_procs):
                 self.S.transfer(source=self.S.levels[l], target=self.S.levels[l - 1])
                 self.S.levels[l - 1].sweep.compute_end_point()
 
-                # on middle levels: do sweep as usual
-                if l - 1 > 0:
-                    req_send = None
-                    if not self.S.status.last and self.params.fine_comm:
-                        self.logger.debug('send data: process %s, stage %s, time %s, target %s, tag %s, iter %s' %
-                                          (self.S.status.slot, self.S.status.stage, self.S.time, self.S.next,
-                                           l - 1, self.S.status.iter))
-                        req_send = comm.isend(self.S.levels[l - 1].uend, dest=self.S.next, tag=self.S.status.iter)
+                req_send = None
+                if not self.S.status.last and self.params.fine_comm:
+                    self.logger.debug('send data: process %s, stage %s, time %s, target %s, tag %s, iter %s' %
+                                      (self.S.status.slot, self.S.status.stage, self.S.time, self.S.next,
+                                       l - 1, self.S.status.iter))
+                    req_send = comm.isend(self.S.levels[l - 1].uend, dest=self.S.next, tag=self.S.status.iter)
 
-                    if not self.S.status.first and self.params.fine_comm:
-                        self.logger.debug('recv data: process %s, stage %s, time %s, source %s, tag %s, iter %s' %
-                                          (self.S.status.slot, self.S.status.stage, self.S.time, self.S.prev,
-                                           l - 1, self.S.status.iter))
-                        self.recv(target=self.S.levels[l - 1], source=self.S.prev, tag=self.S.status.iter, comm=comm)
+                if not self.S.status.first and self.params.fine_comm:
+                    self.logger.debug('recv data: process %s, stage %s, time %s, source %s, tag %s, iter %s' %
+                                      (self.S.status.slot, self.S.status.stage, self.S.time, self.S.prev,
+                                       l - 1, self.S.status.iter))
+                    self.recv(target=self.S.levels[l - 1], source=self.S.prev, tag=self.S.status.iter, comm=comm)
 
-                    if not self.S.status.last and self.params.fine_comm:
-                        req_send.wait()
+                if not self.S.status.last and self.params.fine_comm:
+                    req_send.wait()
+
+                # on middle levels: do sweep as usual
+                if l - 1 > 0:
 
-                    self.hooks.pre_sweep(step=self.S, level_number=l - 1)
                     for k in range(self.S.levels[l - 1].params.nsweeps):
+
+                        self.hooks.pre_sweep(step=self.S, level_number=l - 1)
                         self.S.levels[l - 1].sweep.update_nodes()
-                    self.S.levels[l - 1].sweep.compute_residual()
-                    self.hooks.post_sweep(step=self.S, level_number=l - 1)
+
+                        req_send = None
+                        if not self.S.status.last and self.params.fine_comm:
+                            self.logger.debug('send data: process %s, stage %s, time %s, target %s, tag %s, iter %s' %
+                                              (self.S.status.slot, self.S.status.stage, self.S.time, self.S.next,
+                                               l - 1, self.S.status.iter))
+                            req_send = comm.isend(self.S.levels[l - 1].uend, dest=self.S.next, tag=self.S.status.iter)
+
+                        if not self.S.status.first and self.params.fine_comm:
+                            self.logger.debug('recv data: process %s, stage %s, time %s, source %s, tag %s, iter %s' %
+                                              (self.S.status.slot, self.S.status.stage, self.S.time, self.S.prev,
+                                               l - 1, self.S.status.iter))
+                            self.recv(target=self.S.levels[l - 1], source=self.S.prev, tag=self.S.status.iter,
+                                      comm=comm)
+
+                        if not self.S.status.last and self.params.fine_comm:
+                            req_send.wait()
+
+                        self.S.levels[l - 1].sweep.compute_residual()
+                        self.hooks.post_sweep(step=self.S, level_number=l - 1)
 
             # update stage
-            self.S.status.stage = 'IT_CHECK'
+            self.S.status.stage = 'IT_FINE'
 
         else: