routine for plotting stability regions from stabilty function; requires beautification

Author: Daniel Ruprecht

examples/fwsw/playground.py

@@ -5,7 +5,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from examples.SWFW.ProblemClass import swfw_scalar
+from examples.fwsw.ProblemClass import swfw_scalar
 from pySDC.datatype_classes.complex_mesh import mesh, rhs_imex_mesh
 from pySDC.sweeper_classes.imex_1st_order import imex_1st_order
 import pySDC.PFASST_stepwise as mp

examples/fwsw/plot_stability.py

@@ -5,19 +5,28 @@
 
 from pySDC.datatype_classes.complex_mesh import mesh, rhs_imex_mesh
 from pySDC.sweeper_classes.imex_1st_order import imex_1st_order as imex
-from examples.SWFW.ProblemClass import swfw_scalar 
+from examples.fwsw.ProblemClass import swfw_scalar 
 import numpy as np
 
+import matplotlib.pyplot as plt
+
 if __name__ == "__main__":
 
+    N_s = 100
+    N_f = 400
+
+    lambda_s = 1j*np.linspace(0.0, 2.0, N_s)
+    lambda_f = 1j*np.linspace(0.0, 8.0, N_f)
+
     pparams = {}
-    pparams['lambda_s'] = np.array([-0.0], dtype='complex')
-    pparams['lambda_f'] = np.array([-1.0], dtype='complex')
+    # the following are not used in the computation
+    pparams['lambda_s'] = np.array([0.0])
+    pparams['lambda_f'] = np.array([0.0])
     pparams['u0'] = 1.0
     swparams = {}
     swparams['collocation_class'] = collclass.CollGaussLobatto
-    swparams['num_nodes'] = 9
-    K = 1
+    swparams['num_nodes'] = 2
+    K = 2
 
     #
     # ...this is functionality copied from test_imexsweeper. Ideally, it should be available in one place.
@@ -32,17 +41,37 @@
     nnodes  = step.levels[0].sweep.coll.num_nodes
     level   = step.levels[0]
     problem = level.prob
+  
     QE = level.sweep.QE[1:,1:]
     QI = level.sweep.QI[1:,1:]
     Q  = level.sweep.coll.Qmat[1:,1:]
 
-    dt = step.status.dt
-    LHS = np.eye(nnodes) - step.status.dt*( problem.lambda_f[0]*QI + problem.lambda_s[0]*QE )
-    RHS = step.status.dt*( (problem.lambda_f[0]+problem.lambda_s[0])*Q - (problem.lambda_f[0]*QI + problem.lambda_s[0]*QE) )
+    stab = np.zeros((N_f, N_s), dtype='complex')
+
+    for i in range(0,N_s):
+      for j in range(0,N_f):
+        lambda_fast = lambda_f[j]
+        lambda_slow = lambda_s[i]
+        LHS = np.eye(nnodes) - step.status.dt*( lambda_fast*QI + lambda_slow*QE )
+        RHS = step.status.dt*( (lambda_fast+lambda_slow)*Q - (lambda_fast*QI + lambda_slow*QE) )
+
+        Pinv = np.linalg.inv(LHS)
+        Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
+        for k in range(0,K):
+          Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),k).dot(Pinv)
+          stab_fh = 1.0 + (lambda_fast + lambda_slow)*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
+        stab[j,i] = stab_fh
 
-    Pinv = np.linalg.inv(LHS)
-    Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
-    for i in range(0,K):
-      Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),i).dot(Pinv)
-    stab_fh = 1.0 + (pparams['lambda_s'][0] + pparams['lambda_f'][0])*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
-    print abs(stab_fh)
+    ###
+    fig  = plt.figure(figsize=(12,12))
+    #pcol = plt.pcolor(lambda_s.imag, lambda_f.imag, np.absolute(stab), vmin=0.99, vmax=2.01)
+    #pcol.set_edgecolor('face')
+    levels = np.array([0.95, 0.99, 1.01])
+    CS = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), levels, colors='k')
+    plt.clabel(CS, fontsize=9)
+    plt.plot([0, 4], [0, 4], color='k', linewidth=2.5)
+    plt.gca().set_xticks([0.0, 1.0, 2.0, 3.0])
+    plt.xlim([np.min(lambda_s.imag), np.max(lambda_s.imag)])
+    plt.xlabel('$\Delta t \lambda_{slow}$', fontsize=18, labelpad=0.0)
+    plt.ylabel('$\Delta t \lambda_{fast}$', fontsize=18)
+    plt.show()