many fixes for transfer

Author: pancetta

pySDC/implementations/problem_classes/TestEquation_0D.py

@@ -35,11 +35,11 @@ def __init__(self, problem_params, dtype_u, dtype_f):
 
         assert not any(isinstance(i, list) for i in problem_params['lambdas']), \
             'ERROR: expect flat list here, got %s' % problem_params['lambdas']
-        nvars = len(problem_params['lambdas'])
-        assert nvars > 0, 'ERROR: expect at least one lambda parameter here'
+        problem_params['nvars'] = len(problem_params['lambdas'])
+        assert problem_params['nvars'] > 0, 'ERROR: expect at least one lambda parameter here'
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
-        super(testequation0d, self).__init__(init=nvars, dtype_u=dtype_u, dtype_f=dtype_f,
+        super(testequation0d, self).__init__(init=problem_params['nvars'], dtype_u=dtype_u, dtype_f=dtype_f,
                                              params=problem_params)
 
         self.A = self.__get_A(self.params.lambdas)

pySDC/implementations/transfer_classes/TransferMesh_NoCoarse.py

@@ -1,5 +1,7 @@
 from __future__ import division
 
+import scipy.sparse as sp
+
 from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh
 from pySDC.core.SpaceTransfer import space_transfer
 from pySDC.core.Errors import TransferError
@@ -29,6 +31,9 @@ def __init__(self, fine_prob, coarse_prob, params):
         # invoke super initialization
         super(mesh_to_mesh, self).__init__(fine_prob, coarse_prob, params)
 
+        self.Rspace = sp.eye(self.coarse_prob.params.nvars)
+        self.Pspace = sp.eye(self.fine_prob.params.nvars)
+
     def restrict(self, F):
         """
         Restriction implementation

pySDC/projects/matrixPFASST/allinclusive_matrix_nonMPI.py

@@ -1,6 +1,4 @@
 import numpy as np
-import scipy.sparse as sp
-import scipy.sparse.linalg as spla
 
 from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
 
@@ -55,16 +53,16 @@ def __init__(self, num_procs, controller_params, description):
         assert self.nlevels <= 2, 'ERROR: cannot use matrix-PFASST with more than 2 levels'  # TODO: fixme
         assert [level.dt for step in self.MS for level in step.levels].count(self.dt) == self.nlevels * self.nsteps, \
             'ERROR: dt must be equal for all steps and all levels'
-        assert [level.sweep.coll.num_nodes for step in self.MS for level in step.levels].count(self.nnodes) == \
-            self.nlevels * self.nsteps, 'ERROR: nnodes must be equal for all steps and all levels'
+        # assert [level.sweep.coll.num_nodes for step in self.MS for level in step.levels].count(self.nnodes) == \
+        #     self.nlevels * self.nsteps, 'ERROR: nnodes must be equal for all steps and all levels'
         assert [type(level.prob) for step in self.MS for level in step.levels].count(type(prob)) == \
             self.nlevels * self.nsteps, 'ERROR: all probem classes have to be the same'
 
         assert self.params.predict is False, 'ERROR: no predictor for matrix controller yet'  # TODO: fixme
 
         assert hasattr(prob, 'A'), 'ERROR: need system matrix A for this (and linear problems!)'
 
-        A = prob.A
+        A = prob.A.todense()
         Q = self.MS[0].levels[0].sweep.coll.Qmat[1:, 1:]
         Qd = self.MS[0].levels[0].sweep.QI[1:, 1:]
 
@@ -74,26 +72,33 @@ def __init__(self, num_procs, controller_params, description):
         N = np.zeros((self.nnodes, self.nnodes))
         N[:, -1] = 1
 
-        self.C = sp.eye(self.nsteps * self.nnodes * self.nspace) - \
-            self.dt * sp.kron(sp.eye(self.nsteps), sp.kron(Q, A)) - sp.kron(E, sp.kron(N, sp.eye(self.nspace)))
-        self.P = sp.eye(self.nsteps * self.nnodes * self.nspace) - \
-            self.dt * sp.kron(sp.eye(self.nsteps), sp.kron(Qd, A))
+        self.C = np.array(np.eye(self.nsteps * self.nnodes * self.nspace) - \
+            self.dt * np.kron(np.eye(self.nsteps), np.kron(Q, A)) - np.kron(E, np.kron(N, np.eye(self.nspace))))
+        self.P = np.array(np.eye(self.nsteps * self.nnodes * self.nspace) - \
+            self.dt * np.kron(np.eye(self.nsteps), np.kron(Qd, A)))
 
         if self.nlevels > 1:
             prob_c = self.MS[0].levels[1].prob
             self.nspace_c = prob_c.init
 
-            Ac = prob_c.A
+            Ac = prob_c.A.todense()
             Qdc = self.MS[0].levels[1].sweep.QI[1:, 1:]
 
-            TcfA = self.MS[0].base_transfer.space_transfer.Pspace
-            TfcA = self.MS[0].base_transfer.space_transfer.Rspace
+            nnodesc = self.MS[0].levels[1].sweep.coll.num_nodes
+            Nc = np.zeros((nnodesc, nnodesc))
+            Nc[:, -1] = 1
 
-            self.Tcf = sp.kron(sp.eye(self.nsteps * self.nnodes), TcfA)
-            self.Tfc = sp.kron(sp.eye(self.nsteps * self.nnodes), TfcA)
+            TcfA = self.MS[0].base_transfer.space_transfer.Pspace.todense()
+            TfcA = self.MS[0].base_transfer.space_transfer.Rspace.todense()
+            TcfQ = self.MS[0].base_transfer.Pcoll
+            TfcQ = self.MS[0].base_transfer.Rcoll
 
-            self.Pc = sp.eye(self.nsteps * self.nnodes * self.nspace_c) - \
-                self.dt * sp.kron(sp.eye(self.nsteps), sp.kron(Qdc, Ac)) - sp.kron(E, sp.kron(N, sp.eye(self.nspace_c)))
+            self.Tcf = np.array(np.kron(np.eye(self.nsteps), np.kron(TcfQ, TcfA)))
+            self.Tfc = np.array(np.kron(np.eye(self.nsteps), np.kron(TfcQ, TfcA)))
+
+            self.Pc = np.array(np.eye(self.nsteps * nnodesc * self.nspace_c) - \
+                self.dt * np.kron(np.eye(self.nsteps), np.kron(Qdc, Ac)) - \
+                np.kron(E, np.kron(Nc, np.eye(self.nspace_c))))
 
         self.u = np.zeros(self.nsteps * self.nnodes * self.nspace)
         self.res = np.zeros(self.nsteps * self.nnodes * self.nspace)
@@ -167,9 +172,9 @@ def build_propagation_matrix(self, niter):
         """
 
         # build smoother iteration matrix and preconditioner using nsweeps
-        Pinv = np.linalg.inv(self.P.todense())
+        Pinv = np.linalg.inv(self.P)
         precond_smoother = Pinv.copy()
-        iter_mat_smoother = np.eye(self.nsteps * self.nnodes * self.nspace) - precond_smoother.dot(self.C.todense())
+        iter_mat_smoother = np.eye(self.nsteps * self.nnodes * self.nspace) - precond_smoother.dot(self.C)
         for k in range(1, self.MS[0].levels[0].params.nsweeps):
             precond_smoother += np.linalg.matrix_power(iter_mat_smoother, k).dot(Pinv)
         iter_mat_smoother = np.linalg.matrix_power(iter_mat_smoother, self.MS[0].levels[0].params.nsweeps)
@@ -178,10 +183,10 @@ def build_propagation_matrix(self, niter):
 
         # add coarse-grid correction (single sweep, though!)
         if self.nlevels > 1:
-            precond_cgc = self.Tcf.todense().dot(np.linalg.inv(self.Pc.todense())).dot(self.Tfc.todense())
-            iter_mat_cgc = np.eye(self.nsteps * self.nnodes * self.nspace) - precond_cgc.dot(self.C.todense())
+            precond_cgc = self.Tcf.dot(np.linalg.inv(self.Pc)).dot(self.Tfc)
+            iter_mat_cgc = np.eye(self.nsteps * self.nnodes * self.nspace) - precond_cgc.dot(self.C)
             iter_mat = iter_mat.dot(iter_mat_cgc)
-            precond += precond_cgc - precond_smoother.dot(self.C.todense()).dot(precond_cgc)
+            precond += precond_cgc - precond_smoother.dot(self.C).dot(precond_cgc)
 
         # form span and reduce matrices and add to operator
         Tspread = np.kron(np.concatenate([[1] * (self.nsteps * self.nnodes)]), np.eye(self.nspace)).T
@@ -223,7 +228,7 @@ def restart_block(self, slots, time, u0):
             for lvl in self.MS[p].levels:
                 lvl.status.time = time[p]
                 P = lvl.prob
-                for m in range(1, self.nnodes + 1):
+                for m in range(1, lvl.sweep.coll.num_nodes + 1):
                     lvl.u[m] = P.dtype_u(init=P.init, val=0)
                     lvl.f[m] = P.dtype_f(init=P.init, val=0)
 
@@ -288,7 +293,7 @@ def pfasst(self, MS):
             for S in MS:
                 self.hooks.pre_sweep(step=S, level_number=0)
 
-            self.u += spla.spsolve(self.P, self.res)
+            self.u += np.linalg.solve(self.P, self.res)
             self.res = self.u0 - self.C.dot(self.u)
 
             MS = self.update_data(MS=MS, u=self.u, res=self.res, niter=niter, level=0, stage='POST_PRE_SWEEP')
@@ -310,7 +315,7 @@ def pfasst(self, MS):
                     for S in MS:
                         self.hooks.pre_sweep(step=S, level_number=1)
 
-                    self.u += self.Tcf.dot(spla.spsolve(self.Pc, self.Tfc.dot(self.res)))
+                    self.u += self.Tcf.dot(np.linalg.solve(self.Pc, self.Tfc.dot(self.res)))
                     self.res = self.u0 - self.C.dot(self.u)
 
                     MS = self.update_data(MS=MS, u=self.u, res=self.res, niter=niter, level=1,
@@ -324,7 +329,7 @@ def pfasst(self, MS):
                 for S in MS:
                     self.hooks.pre_sweep(step=S, level_number=0)
 
-                self.u += spla.spsolve(self.P, self.res)
+                self.u += np.linalg.solve(self.P, self.res)
                 self.res = self.u0 - self.C.dot(self.u)
 
                 MS = self.update_data(MS=MS, u=self.u, res=self.res, niter=niter, level=0, stage='POST_FINE_SWEEP')

pySDC/projects/matrixPFASST/compare_to_propagator.py

@@ -140,12 +140,12 @@ def testequation_setup():
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 0.25
-    level_params['nsweeps'] = [3]
+    level_params['nsweeps'] = [1, 1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
-    sweeper_params['num_nodes'] = [3]
+    sweeper_params['num_nodes'] = [3, 3]
     sweeper_params['QI'] = 'LU'
     sweeper_params['spread'] = True
 
@@ -155,19 +155,19 @@ def testequation_setup():
     # use single values like this...
     problem_params['lambdas'] = [[-1.0]]
     # .. or a list of values like this ...
-    problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
+    # problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
     # .. or a whole block of values like this
-    ilim_left = -11
-    ilim_right = 0
-    rlim_left = 0
-    rlim_right = 11
-    ilam = 1j * np.logspace(ilim_left, ilim_right, 11)
-    rlam = -1 * np.logspace(rlim_left, rlim_right, 11)
-    lambdas = []
-    for rl in rlam:
-        for il in ilam:
-            lambdas.append(rl + il)
-    problem_params['lambdas'] = [lambdas]
+    # ilim_left = -11
+    # ilim_right = 0
+    # rlim_left = 0
+    # rlim_right = 11
+    # ilam = 1j * np.logspace(ilim_left, ilim_right, 11)
+    # rlam = -1 * np.logspace(rlim_left, rlim_right, 11)
+    # lambdas = []
+    # for rl in rlam:
+    #     for il in ilam:
+    #         lambdas.append(rl + il)
+    # problem_params['lambdas'] = [lambdas]
     # note: PFASST will do all of those at once, but without interaction (realized via diagonal matrix).
     # The propagation matrix will be diagonal too, corresponding to the respective lambda value.