fixed propagator (now for real), added test, bugfixes

Author: pancetta

pySDC/implementations/transfer_classes/TransferMesh_NoCoarse.py

@@ -3,6 +3,7 @@
 import scipy.sparse as sp
 
 from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh
+from pySDC.implementations.datatype_classes.complex_mesh import mesh as cmesh
 from pySDC.core.SpaceTransfer import space_transfer
 from pySDC.core.Errors import TransferError
 
@@ -43,6 +44,8 @@ def restrict(self, F):
         """
         if isinstance(F, mesh):
             G = mesh(F)
+        elif isinstance(F, cmesh):
+            G = cmesh(F)
         elif isinstance(F, rhs_imex_mesh):
             G = rhs_imex_mesh(F)
         else:
@@ -58,6 +61,8 @@ def prolong(self, G):
         """
         if isinstance(G, mesh):
             F = mesh(G)
+        elif isinstance(G, cmesh):
+            F = cmesh(G)
         elif isinstance(G, rhs_imex_mesh):
             F = rhs_imex_mesh(G)
         else:

pySDC/projects/matrixPFASST/allinclusive_matrix_nonMPI.py

@@ -85,7 +85,6 @@ def __init__(self, num_procs, controller_params, description):
 
             Ac = prob_c.A.todense()
             Qdc = self.MS[0].levels[1].sweep.QI[1:, 1:]
-
             nnodesc = self.MS[0].levels[1].sweep.coll.num_nodes
             Nc = np.zeros((nnodesc, nnodesc))
             Nc[:, -1] = 1
@@ -181,15 +180,16 @@ def build_propagation_matrix(self, niter):
         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)
-        iter_mat = iter_mat_smoother
-        precond = precond_smoother
 
         # add coarse-grid correction (single sweep, though!)
         if self.nlevels > 1:
             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).dot(precond_cgc)
+            iter_mat = iter_mat_smoother.dot(iter_mat_cgc)
+            precond = precond_smoother + precond_cgc - precond_smoother.dot(self.C).dot(precond_cgc)
+        else:
+            iter_mat = iter_mat_smoother
+            precond = precond_smoother
 
         # form span and reduce matrices and add to operator
         Tspread = np.kron(np.concatenate([[1] * (self.nsteps * self.nnodes)]), np.eye(self.nspace)).T

pySDC/projects/matrixPFASST/compare_to_matrixbased.py

@@ -1,21 +1,26 @@
+import numpy as np
+
 from pySDC.implementations.problem_classes.HeatEquation_1D_FD import heat1d
 from pySDC.implementations.problem_classes.AdvectionEquation_1D_FD import advection1d
+from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d
 from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.complex_mesh import mesh as cmesh
 from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
 from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
+from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nocoarse
 from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
 from pySDC.projects.matrixPFASST.allinclusive_matrix_nonMPI import allinclusive_matrix_nonMPI
 
 from pySDC.helpers.stats_helper import filter_stats, sort_stats
 
 
-def diffusion_setup(mu=0.0):
+def diffusion_setup(par=0.0):
     """
     Setup routine for advection test
 
     Args:
-        mu (float): parameter for controlling stiffness
+        par (float): parameter for controlling stiffness
     """
     # initialize level parameters
     level_params = dict()
@@ -32,7 +37,7 @@ def diffusion_setup(mu=0.0):
 
     # initialize problem parameters
     problem_params = dict()
-    problem_params['nu'] = mu  # diffusion coefficient
+    problem_params['nu'] = par  # diffusion coefficient
     problem_params['freq'] = 4  # frequency for the test value
     problem_params['nvars'] = [127, 63]  # number of degrees of freedom for each level
 
@@ -66,12 +71,12 @@ def diffusion_setup(mu=0.0):
     return description, controller_params
 
 
-def advection_setup(mu=0.0):
+def advection_setup(par=0.0):
     """
     Setup routine for advection test
 
     Args:
-        mu (float): parameter for controlling stiffness
+        par (float): parameter for controlling stiffness
     """
     # initialize level parameters
     level_params = dict()
@@ -88,7 +93,7 @@ def advection_setup(mu=0.0):
 
     # initialize problem parameters
     problem_params = dict()
-    problem_params['c'] = mu
+    problem_params['c'] = par
     problem_params['freq'] = 4  # frequency for the test value
     problem_params['nvars'] = [128, 64]  # number of degrees of freedom for each level
     problem_params['order'] = 2
@@ -125,13 +130,80 @@ def advection_setup(mu=0.0):
     return description, controller_params
 
 
-def compare_controllers(type=None, mu=0.0, f=None):
+def testequation_setup():
+    """
+    Setup routine for the test equation
+
+    Args:
+        par (float): parameter for controlling stiffness
+    """
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-08
+    level_params['dt'] = 0.25
+    level_params['nsweeps'] = [3, 1]
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussRadau_Right
+    sweeper_params['num_nodes'] = [3, 2]
+    sweeper_params['QI'] = 'LU'
+    sweeper_params['spread'] = True
+
+    # initialize problem parameters
+    problem_params = dict()
+    problem_params['u0'] = 1.0  # initial value (for all instances)
+    # 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]]
+    # .. 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]
+    # 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.
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 50
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 30
+    controller_params['predict'] = False
+
+    # fill description dictionary for easy step instantiation
+    description = dict()
+    description['problem_class'] = testequation0d  # pass problem class
+    description['problem_params'] = problem_params  # pass problem parameters
+    description['dtype_u'] = cmesh  # pass data type for u
+    description['dtype_f'] = cmesh  # pass data type for f
+    description['sweeper_class'] = generic_implicit  # pass sweeper
+    description['sweeper_params'] = sweeper_params  # pass sweeper parameters
+    description['level_params'] = level_params  # pass level parameters
+    description['step_params'] = step_params  # pass step parameters
+    description['space_transfer_class'] = mesh_to_mesh_nocoarse  # pass spatial transfer class
+    description['space_transfer_params'] = dict()  # pass paramters for spatial transfer
+
+    return description, controller_params
+
+
+def compare_controllers(type=None, par=0.0, f=None):
     """
     A simple test program to compare PFASST runs with matrix-based and matrix-free controllers
 
     Args:
         type (str): setup type
-        mu (float) parameter for controlling stiffness
+        par (float) parameter for controlling stiffness
         f: file handler
     """
 
@@ -140,13 +212,15 @@ def compare_controllers(type=None, mu=0.0, f=None):
     Tend = 1.0
 
     if type == 'diffusion':
-        description, controller_params = diffusion_setup(mu)
+        description, controller_params = diffusion_setup(par)
     elif type == 'advection':
-        description, controller_params = advection_setup(mu)
+        description, controller_params = advection_setup(par)
+    elif type == 'testequation':
+        description, controller_params = testequation_setup()
     else:
         raise ValueError('No valis setup type provided, aborting..')
 
-    out = '\nWorking with %s setup and parameter %3.1e..' % (type, mu)
+    out = '\nWorking with %s setup and parameter %3.1e..' % (type, par)
     f.write(out + '\n')
     print(out)
 
@@ -178,7 +252,7 @@ def compare_controllers(type=None, mu=0.0, f=None):
     f.write(out + '\n')
     print(out)
 
-    assert diff < 2.0E-15, 'ERROR: difference between matrix-based and matrix-free result is too large, got %s' % diff
+    assert diff < 2.3E-15, 'ERROR: difference between matrix-based and matrix-free result is too large, got %s' % diff
 
     # filter statistics by type (number of iterations)
     filtered_stats_mat = filter_stats(stats_mat, type='niter')
@@ -201,12 +275,13 @@ def compare_controllers(type=None, mu=0.0, f=None):
 
 def main():
 
-    mu_list = [1E-02, 1.0, 1E+02]
+    par_list = [1E-02, 1.0, 1E+02]
 
-    f = open('comparison_matrix_vs_nomat_detail.txt', 'a')
-    for mu in mu_list:
-        compare_controllers(type='diffusion', mu=mu, f=f)
-        compare_controllers(type='advection', mu=mu, f=f)
+    f = open('comparison_matrix_vs_nomat_detail.txt', 'w')
+    for par in par_list:
+        compare_controllers(type='diffusion', par=par, f=f)
+        compare_controllers(type='advection', par=par, f=f)
+    compare_controllers(type='testequation', par=0.0, f=f)
     f.close()
 
 

pySDC/projects/matrixPFASST/compare_to_propagator.py

@@ -38,7 +38,7 @@ def diffusion_setup(par=0.0):
     problem_params = dict()
     problem_params['nu'] = par  # diffusion coefficient
     problem_params['freq'] = 4  # frequency for the test value
-    problem_params['nvars'] = [127, 63]  # number of degrees of freedom for each level
+    problem_params['nvars'] = [127]  # number of degrees of freedom for each level
 
     # initialize step parameters
     step_params = dict()
@@ -140,34 +140,34 @@ def testequation_setup():
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 0.25
-    level_params['nsweeps'] = [1, 1]
+    level_params['nsweeps'] = [3, 1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
-    sweeper_params['num_nodes'] = [3, 3]
+    sweeper_params['num_nodes'] = [3, 2]
     sweeper_params['QI'] = 'LU'
     sweeper_params['spread'] = True
 
     # initialize problem parameters
     problem_params = dict()
     problem_params['u0'] = 1.0  # initial value (for all instances)
     # use single values like this...
-    problem_params['lambdas'] = [[-1.0]]
+    # problem_params['lambdas'] = [[-1.0]]
     # .. or a list of values like this ...
     # 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.
 
@@ -253,24 +253,27 @@ def compare_controllers(type=None, par=0.0, f=None):
     prop = controller.build_propagation_matrix(niter=niter)
 
     err_prop_ex = np.linalg.norm(prop.dot(uinit.values) - uex.values)
-    out = '  Difference between propagation and exact solution: %6.4e' % err_prop_ex
-    f.write(out + '\n')
-    print(out)
     err_mat_ex = np.linalg.norm(uend_mat.values - uex.values)
-    out = '  Difference between matrix-PFASST and exact solution: %6.4e' % err_mat_ex
+    out = '  Error (mat/prop) vs. exact solution: %6.4e -- %6.4e' % (err_mat_ex, err_prop_ex)
     f.write(out + '\n')
     print(out)
     err_mat_prop = np.linalg.norm(prop.dot(uinit.values) - uend_mat.values)
     out = '  Difference between matrix-PFASST and propagator: %6.4e' % err_mat_prop
     f.write(out + '\n')
     print(out)
 
+    assert err_mat_prop < 2.0E-14, \
+        'ERROR: difference between matrix-based and propagator result is too large, got %s' % err_mat_prop
+
 
 def main():
 
+    par_list = [1E-02, 1.0, 1E+02]
+
     f = open('comparison_matrix_vs_propagator_detail.txt', 'w')
-    # compare_controllers(type='diffusion', par=1E-00, f=f)
-    # compare_controllers(type='advection', par=1E-00, f=f)
+    for par in par_list:
+        compare_controllers(type='diffusion', par=par, f=f)
+        compare_controllers(type='advection', par=par, f=f)
     compare_controllers(type='testequation', par=0.0, f=f)
     f.close()
 

pySDC/tests/test_projects/test_matrixPFASST/test_controllers.py

@@ -0,0 +1,10 @@
+from pySDC.projects.matrixPFASST.compare_to_matrixbased import main as A
+from pySDC.projects.matrixPFASST.compare_to_propagator import main as B
+
+def test_matrixbased():
+    A()
+
+def test_propagator():
+    B()
+
+