added paper scripts to pySDC

Author: pancetta

projects/AsympConv/PFASST_conv_tests.py

@@ -2,7 +2,7 @@
 import pickle
 import os
 import matplotlib
-# matplotlib.use('Agg')
+matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 
 from pySDC.implementations.problem_classes.HeatEquation_1D_FD import heat1d
@@ -17,20 +17,28 @@
 
 
 def main():
-    run_diffusion()
-    run_advection()
+    """
+    Main driver running diffusion and advection tests
+    """
+    nsweeps = 3
+    run_diffusion(nsweeps=nsweeps)
+    run_advection(nsweeps=nsweeps)
+    plot_results(nsweeps=nsweeps)
 
 
-def run_diffusion():
+def run_diffusion(nsweeps):
     """
-    A simple test program to test PFASST convergence for the periodic heat equation with random initial data
+    A simple test program to test PFASST convergence for the heat equation with random initial data
+
+    Args:
+        nsweeps: number of fine sweeps to perform
     """
 
     # initialize level parameters
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 0.25
-    level_params['nsweeps'] = [3, 1]
+    level_params['nsweeps'] = [nsweeps, 1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
@@ -114,6 +122,9 @@ def run_diffusion():
         out = '  Mean number of iterations: %4.2f' % np.mean(niters)
         print(out)
 
+        assert np.mean(niters) <= 6, "ERROR: number of iterations for diffusion tesst too high, got %s" \
+                                     % np.mean(niters)
+
         results[cfl] = np.mean(niters)
 
     file = open('data/results_conv_diffusion.pkl', 'wb')
@@ -123,16 +134,19 @@ def run_diffusion():
     assert os.path.isfile('data/results_conv_diffusion.pkl'), 'ERROR: pickle did not create file'
 
 
-def run_advection():
+def run_advection(nsweeps):
     """
-    A simple test program to test PFASST convergence for the periodic heat equation with random initial data
+    A simple test program to test PFASST convergence for the periodic advection equation
+
+    Args:
+        nsweeps: number of fine sweeps to perform
     """
 
     # initialize level parameters
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 0.25
-    level_params['nsweeps'] = [3, 1]
+    level_params['nsweeps'] = [nsweeps, 1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
@@ -218,6 +232,9 @@ def run_advection():
         out = '  Mean number of iterations: %4.2f' % np.mean(niters)
         print(out)
 
+        assert np.mean(niters) <= 10, "ERROR: number of iterations for advection tesst too high, got %s" \
+            % np.mean(niters)
+
         results[cfl] = np.mean(niters)
 
     file = open('data/results_conv_advection.pkl', 'wb')
@@ -227,7 +244,13 @@ def run_advection():
     assert os.path.isfile('data/results_conv_advection.pkl'), 'ERROR: pickle did not create file'
 
 
-def plot_results():
+def plot_results(nsweeps):
+    """
+    Plotting routine for iteration counts
+
+    Args:
+        nsweeps: number of fine sweeps used
+    """
 
     file = open('data/results_conv_diffusion.pkl', 'rb')
     results_diff = pickle.load(file)
@@ -269,18 +292,16 @@ def plot_results():
     # plot
     plt.semilogx(xvalues_diff, niter_diff, 'r-', marker='s', markersize=10, label='diffusion')
     plt.semilogx(xvalues_adv, niter_adv, 'b-', marker='o', markersize=10, label='advection')
-    # plt.plot()
 
     plt.legend(loc=1, ncol=1, numpoints=1)
 
     # plt.show()
     # save plot, beautify
-    fname = 'data/conv_test_niter.png'
+    fname = 'data/conv_test_niter_NS' + str(nsweeps) + '.png'
     plt.savefig(fname, rasterized=True, bbox_inches='tight')
 
     assert os.path.isfile(fname), 'ERROR: plotting did not create file'
 
 
 if __name__ == "__main__":
     main()
-    plot_results()

projects/AsympConv/conv_test_to0.py

@@ -0,0 +1,88 @@
+import numpy as np
+import scipy.linalg as LA
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pylab as plt
+from matplotlib import rc
+
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+
+
+def compute_and_plot_specrad(Nnodes, lam):
+    """
+    Compute and plot the spectral radius of the smoother for different step-sizes
+
+    Args:
+        Nnodes: number of collocation nodes
+        lam: test parameter representing the spatial problem
+    """
+
+    coll = CollGaussRadau_Right(Nnodes, 0, 1)
+    Qmat = coll.Qmat[1:, 1:]
+
+    # do LU decomposition of QT (St. Martin's trick)
+    QT = coll.Qmat[1:, 1:].T
+    [_, _, U] = LA.lu(QT, overwrite_a=True)
+    QDmat = U.T
+
+    Nmat = np.zeros((Nnodes, Nnodes))
+    Nmat[:, -1] = 1
+
+    Nsteps_list = [64, 256]
+    color_list = ['red', 'blue']
+    marker_list = ['s', 'o']
+
+    setup_list = zip(Nsteps_list, color_list, marker_list)
+
+    xlist = [0.1 ** i for i in range(11)]
+
+    plt.subplots(figsize=(15, 10))
+    rc('font', **{"sans-serif": ["Arial"], "size": 24})
+    rc('legend', fontsize='small')
+    rc('xtick', labelsize='small')
+    rc('ytick', labelsize='small')
+
+    for Nsteps, color, marker in setup_list:
+
+        Emat = np.zeros((Nsteps, Nsteps))
+        np.fill_diagonal(Emat[1:, :], 1)
+
+        Prho_list = []
+        predict_list = []
+        for x in xlist:
+
+            mat = np.linalg.inv(np.eye(Nnodes * Nsteps) - x * lam * np.kron(np.eye(Nsteps), QDmat)).dot(
+                x * lam * np.kron(np.eye(Nsteps), (Qmat - QDmat)) + np.kron(Emat, Nmat))
+
+            Prho_list.append(max(abs(np.linalg.eigvals(mat))))
+            predict_list.append((1 + x) ** (1.0 - 1.0 / (Nnodes * Nsteps)) * x ** (1.0 / (Nnodes * Nsteps)))
+
+            if len(predict_list) > 1:
+                print(x, predict_list[-1], Prho_list[-1], Prho_list[-2] / Prho_list[-1],
+                      predict_list[-2] / predict_list[-1])
+
+        plt.loglog(xlist, Prho_list, linestyle='-', linewidth=3, color=color, marker=marker, markersize=10,
+                   label='spectral radius, L=' + str(Nsteps))
+        plt.loglog(xlist, [item for item in predict_list], linestyle='--', linewidth=2, color=color, marker=marker,
+                   markersize=10, label='estimate, L=' + str(Nsteps))
+
+    ax = plt.gca()
+    ax.invert_xaxis()
+
+    plt.xlabel('time-step size')
+    plt.ylabel('spectral radius')
+    plt.legend(loc=3, numpoints=1)
+    plt.grid()
+    plt.ylim([1E-02, 1E01])
+
+    if type(lam) is complex:
+        fname = 'data/smoother_specrad_to0_L64+256_M' + str(Nnodes) + 'LU_imag.png'
+    else:
+        fname = 'data/smoother_specrad_to0_L64+256_M' + str(Nnodes) + 'LU_real.png'
+    plt.savefig(fname, rasterized=True, transparent=True, bbox_inches='tight')
+
+
+if __name__ == "__main__":
+    compute_and_plot_specrad(Nnodes=3, lam=-1)
+    compute_and_plot_specrad(Nnodes=3, lam=1j)

projects/AsympConv/conv_test_toinf.py

@@ -0,0 +1,92 @@
+import numpy as np
+import scipy.linalg as LA
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pylab as plt
+from matplotlib import rc
+
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+
+
+def compute_and_plot_specrad(Nnodes, lam):
+    """
+    Compute and plot the spectral radius of the smoother for different step-sizes
+
+    Args:
+        Nnodes: number of collocation nodes
+        lam: test parameter representing the spatial problem
+    """
+
+    Nsteps = 1
+
+    coll = CollGaussRadau_Right(Nnodes, 0, 1)
+    Qmat = coll.Qmat[1:, 1:]
+
+    # do LU decomposition of QT (St. Martin's trick)
+    QT = coll.Qmat[1:, 1:].T
+    [_, _, U] = LA.lu(QT, overwrite_a=True)
+    QDmat = U.T
+
+    Nmat = np.zeros((Nnodes, Nnodes))
+    Nmat[:, -1] = 1
+
+    Nsweep_list = [1, Nnodes - 1, Nnodes]
+    color_list = ['red', 'blue', 'green']
+    marker_list = ['s', 'o', 'd']
+
+    setup_list = zip(Nsweep_list, color_list, marker_list)
+
+    xlist = [10 ** i for i in range(11)]
+
+    plt.subplots(figsize=(15, 10))
+    rc('font', **{"sans-serif": ["Arial"], "size": 24})
+    rc('legend', fontsize='small')
+    rc('xtick', labelsize='small')
+    rc('ytick', labelsize='small')
+
+    for Nsweeps, color, marker in setup_list:
+
+        Emat = np.zeros((Nsteps, Nsteps))
+        np.fill_diagonal(Emat[1:, :], 1)
+
+        Prho_list = []
+        predict_list = []
+        for x in xlist:
+
+            mat = np.linalg.inv(np.eye(Nnodes * Nsteps) - x * lam * np.kron(np.eye(Nsteps), QDmat)).dot(
+                x * lam * np.kron(np.eye(Nsteps), (Qmat - QDmat)) + np.kron(Emat, Nmat))
+            mat = np.linalg.matrix_power(mat, Nsweeps)
+
+            Prho_list.append(max(abs(np.linalg.eigvals(mat))))
+            predict_list.append(1.0 / x)
+
+            if len(predict_list) > 1:
+                print(x, predict_list[-1], Prho_list[-1], Prho_list[-2] / Prho_list[-1],
+                      predict_list[-2] / predict_list[-1])
+
+        plt.loglog(xlist, Prho_list, linestyle='-', linewidth=3, color=color, marker=marker, markersize=10,
+                   label='spectral radius, L=' + str(Nsteps))
+
+    plt.loglog(xlist, [item / predict_list[0] for item in predict_list], linestyle='--', linewidth=2, color='k',
+               label='estimate')
+
+    plt.xlabel('time-step size')
+    plt.ylabel('spectral radius')
+    plt.legend(loc=3, numpoints=1)
+    plt.grid()
+    plt.ylim([1E-16, 1E00])
+
+    if type(lam) is complex:
+        fname = 'data/smoother_specrad_toinf_M' + str(Nnodes) + '_LU_imag.png'
+    else:
+        fname = 'data/smoother_specrad_toinf_M' + str(Nnodes) + '_LU_real.png'
+    plt.savefig(fname, rasterized=True, transparent=True, bbox_inches='tight')
+
+
+if __name__ == "__main__":
+
+    compute_and_plot_specrad(Nnodes=3, lam=-1)
+    compute_and_plot_specrad(Nnodes=3, lam=1j)
+    compute_and_plot_specrad(Nnodes=7, lam=-1)
+    compute_and_plot_specrad(Nnodes=7, lam=1j)