many things, incl. the test equation to check with Thibaut

Author: pancetta

pySDC/playgrounds/Gander/HookClass_error_output.py

@@ -0,0 +1,32 @@
+from pySDC.core.Hooks import hooks
+
+
+class error_output(hooks):
+    """
+    Hook class to add output of error
+    """
+
+    def post_step(self, step, level_number):
+        """
+        Default routine called after each step
+        Args:
+            step: the current step
+            level_number: the current level number
+        """
+
+        super(error_output, self).post_step(step, level_number)
+
+        # some abbreviations
+        L = step.levels[level_number]
+        P = L.prob
+
+        L.sweep.compute_end_point()
+
+        # compute and save errors
+        upde = P.u_exact(step.time + step.dt)
+        pde_err = abs(upde - L.uend)
+
+        self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=L.level_index, iter=step.status.iter,
+                          sweep=L.status.sweep, type='error_after_step', value=pde_err)
+        self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=L.level_index, iter=step.status.iter,
+                          sweep=L.status.sweep, type='residual_after_step', value=L.status.residual)

pySDC/playgrounds/Gander/__init__.py

@@ -0,0 +1,132 @@
+import numpy as np
+
+from pySDC.helpers.stats_helper import filter_stats, sort_stats
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.playgrounds.Gander.HookClass_error_output import error_output
+
+
+def testequation_setup(prec_type=None, maxiter=None):
+    """
+    Setup routine for the test equation
+
+    Args:
+        par (float): parameter for controlling stiffness
+    """
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 0.0
+    level_params['dt'] = 0.25
+    level_params['nsweeps'] = [1]
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussRadau_Right
+    sweeper_params['num_nodes'] = [3]
+    sweeper_params['QI'] = prec_type
+    sweeper_params['initial_guess'] = 'spread'
+
+    # 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]]
+    problem_params['lambdas'] = [[-1000]]
+    # 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'] = maxiter
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 30
+    controller_params['hook_class'] = error_output
+
+    # 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['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
+
+    return description, controller_params
+
+
+def compare_preconditioners(f=None, list_of_k=None):
+
+    # set time parameters
+    t0 = 0.0
+    Tend = 1.0
+
+    for k in list_of_k:
+
+        description_IE, controller_params_IE = testequation_setup(prec_type='IE', maxiter=k)
+        description_LU, controller_params_LU = testequation_setup(prec_type='LU', maxiter=k)
+
+        out = f'\nWorking with maxiter = {k}'
+        f.write(out + '\n')
+        print(out)
+
+        # instantiate controller
+        controller_IE = controller_nonMPI(num_procs=1, controller_params=controller_params_IE, description=description_IE)
+        controller_LU = controller_nonMPI(num_procs=1, controller_params=controller_params_LU, description=description_LU)
+
+        # get initial values on finest level
+        P = controller_IE.MS[0].levels[0].prob
+        uinit = P.u_exact(t0)
+        uex = P.u_exact(Tend)
+
+        # this is where the iteration is happening
+        uend_IE, stats_IE = controller_IE.run(u0=uinit, t0=t0, Tend=Tend)
+        uend_LU, stats_LU = controller_LU.run(u0=uinit, t0=t0, Tend=Tend)
+
+        diff = abs(uend_IE - uend_LU)
+
+        err_IE = abs(uend_IE - uex)
+        err_LU = abs(uend_LU - uex)
+
+        out = '  Error (IE/LU) vs. exact solution: %6.4e -- %6.4e' % (err_IE, err_LU)
+        f.write(out + '\n')
+        print(out)
+        out = '  Difference between both results: %6.4e' % diff
+        f.write(out + '\n')
+        print(out)
+
+        # filter statistics by type
+        filtered_stats_IE = filter_stats(stats_IE, type='error_after_step')
+        filtered_stats_LU = filter_stats(stats_LU, type='error_after_step')
+
+        # convert filtered statistics to list
+        errors_IE = sort_stats(filtered_stats_IE, sortby='time')
+        errors_LU = sort_stats(filtered_stats_LU, sortby='time')
+        print(errors_IE)
+        print(errors_LU)
+
+    # out = '  Iteration counts for matrix-based version: %s' % iter_counts_mat
+    # f.write(out + '\n')
+    # print(out)
+    # out = '  Iteration counts for matrix-free version: %s' % iter_counts
+    # f.write(out + '\n')
+    # print(out)
+    #
+    # assert iter_counts == iter_counts_mat, \
+    #     'ERROR: number of iterations differ between matrix-based and matrix-free controller'
+
+
+def main():
+
+    f = open('comparison_IE_vs_LU.txt', 'w')
+    compare_preconditioners(f=f, list_of_k=[1, 2, 3, 4])
+    f.close()
+
+
+if __name__ == "__main__":
+    main()

pySDC/playgrounds/Gander/testequation.py

@@ -0,0 +1,154 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pySDC.implementations.collocation_classes.gauss_lobatto import \
+    CollGaussLobatto
+from pySDC.implementations.collocation_classes.gauss_radau_left import \
+    CollGaussRadau_Left
+from pySDC.implementations.collocation_classes.gauss_radau_right import \
+    CollGaussRadau_Right
+from pySDC.implementations.collocation_classes.equidistant import \
+    Equidistant
+
+collDict = {'LOBATTO': CollGaussLobatto,
+            'RADAU_LEFT': CollGaussRadau_Left,
+            'RADAU_RIGHT': CollGaussRadau_Right,
+            'EQUID': Equidistant}
+
+# Problem parameters
+tBeg = 0
+tEnd = 1
+lam = -1.0+1j  # \lambda coefficient from the space operator
+numType = np.array(lam).dtype
+u0 = 1  # initial solution
+
+# Discretization parameters
+L = 8  # number of time steps
+M = 5  # number of nodes
+sweepType = 'BE'
+nodesType = 'EQUID'
+
+# Time interval decomposition
+times = np.linspace(tBeg, tEnd, num=L)
+dt = times[1]-times[0]
+
+# Generate nodes, deltas and Q using pySDC
+coll = collDict[nodesType](M, 0, 1)
+nodes = coll._getNodes
+deltas = coll._gen_deltas
+Q = coll._gen_Qmatrix[1:, 1:]
+Q = Q*lam*dt
+
+# Generate Q_\Delta
+QDelta = np.zeros((M, M))
+if sweepType in ['BE', 'FE']:
+    offset = 1 if sweepType == 'FE' else 0
+    for i in range(offset, M):
+        QDelta[i:, :M-i] += np.diag(deltas[:M-i])
+else:
+    raise ValueError(f'sweepType={sweepType}')
+QDelta = QDelta*lam*dt
+
+print(QDelta.real)
+exit()
+
+# Generate other matrices
+H = np.zeros((M, M), dtype=numType)
+H[:, -1] = 1
+I = np.identity(M)
+
+# Generate operators
+ImQDelta = I - QDelta
+ImQ = I - Q
+
+
+def fineSweep(u, uStar, f):
+    u += np.linalg.solve(ImQDelta, H.dot(uStar) + f - ImQ.dot(u))
+
+
+def coarseSweep(u, uStar, f):
+    u += np.linalg.solve(ImQDelta, H.dot(uStar) + f)
+
+
+# Variables
+f = np.zeros((L, M), dtype=numType)
+f[0] = u0
+u = np.zeros((L, M), dtype=numType)
+t = np.array([t + dt*nodes for t in times])
+
+# Exact solution
+uExact = np.exp(t*lam)
+
+
+nSweep = 3  # also number of iteration for Parareal-SDC and PFASST
+
+# SDC solution
+uSDC = u.copy()
+# uSDC[:] = u0
+uStar = u[0]*0
+for l in range(L):
+    for n in range(nSweep):
+        fineSweep(uSDC[l], uStar, f[l])
+    uStar = uSDC[l]
+
+# Parareal-SDC solution <=> pipelined SDC <=> RIDC ?
+uPar = u.copy()
+# uPar[:] = u0
+for n in range(nSweep):
+    uStar = u[0]*0
+    for l in range(L):
+        fineSweep(uPar[l], uStar, f[l])
+        uStar = uPar[l]
+
+# PFASST solution
+uPFA = u.copy()
+if True:
+    # -- initialization with coarse sweep
+    # (equivalent to fine sweep, since uPFA[l] = 0)
+    uStar = u[0]*0
+    for l in range(L):
+        coarseSweep(uPFA[l], uStar, f[l])
+        uStar = uPFA[l]
+# -- PFASST iteration
+for n in range(nSweep):
+    uStar = u[0]*0
+    for l in range(L):
+        uStarPrev = uPFA[l].copy()
+        fineSweep(uPFA[l], uStar, f[l])
+        uStar = uStarPrev
+
+# Plot solution and error
+if True:
+    plt.figure('Solution (amplitude)')
+    for sol, lbl, sym in [[uExact, 'Exact', 'o-'],
+                          [uSDC, 'SDC', 's-'],
+                          [uPar, 'Parareal', '>-'],
+                          [uPFA, 'PFASST', 'p-']]:
+        plt.plot(t.ravel(), sol.ravel().real, sym, label=lbl)
+    plt.legend()
+    plt.grid(True)
+    plt.xlabel('Time')
+
+    plt.figure('Solution (phase)')
+    for sol, lbl, sym in [[uExact, 'Exact', 'o-'],
+                          [uSDC, 'SDC', 's-'],
+                          [uPar, 'Parareal', '>-'],
+                          [uPFA, 'PFASST', 'p-']]:
+        plt.plot(t.ravel(), sol.ravel().imag, sym, label=lbl)
+    plt.legend()
+    plt.grid(True)
+    plt.xlabel('Time')
+
+eSDC = np.abs(uExact.ravel()-uSDC.ravel())
+ePar = np.abs(uExact.ravel()-uPar.ravel())
+ePFA = np.abs(uExact.ravel()-uPFA.ravel())
+
+plt.figure('Error')
+plt.semilogy(t.ravel(), eSDC, 's-', label=f'SDC {nSweep} sweep')
+plt.semilogy(t.ravel(), ePar, '>-', label=f'Parareal {nSweep} iter')
+plt.semilogy(t.ravel(), ePFA, 'p-', label=f'PFASST {nSweep} iter')
+plt.legend()
+plt.grid(True)
+plt.xlabel('Time')
+
+plt.show()

pySDC/playgrounds/Gander/thibaut_algorithms.py

@@ -0,0 +1,24 @@
+import numpy as np
+import scipy.optimize as opt
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+
+
+m = 5
+coll = CollGaussRadau_Right(num_nodes=m, tleft=0.0, tright=1.0)
+Q = coll.Qmat[1:, 1:]
+
+def func(x):
+    return max(abs(np.linalg.eigvals(np.eye(m) - np.diag(x).dot(Q))))
+
+
+def func2(x0):
+    # x0 = 1E-01 * np.ones(m)
+    d = opt.minimize(func, x0, method='Nelder-Mead', tol=1E-08)
+    print(d.fun, np.diag(np.linalg.inv(np.diag(d.x))))
+    return d.fun
+
+
+y0 = 1.0 / np.asarray([0.2818591930905709, 0.2011358490453793, 0.06274536689514164, 0.11790265267514095, 0.1571629578515223])
+d2 = opt.minimize(func2, y0, method='Nelder-Mead', tol=1E-08)
+print(d2)
+# print(np.linalg.inv(np.diag(d.x)))