Merge branch 'feature/plot_stability' of https://github.com/danielru/pySDC into danielru-feature/plot_stability

pancetta · pancetta · commit 8085ed83ecb0 · 2015-12-21T15:00:05.000+01:00
diff --git a/examples/acoustic_1d_imex/runconvergence.py b/examples/acoustic_1d_imex/runconvergence.py
@@ -54,7 +54,7 @@
     #description['transfer_class'] = mesh_to_mesh_1d
     #description['transfer_params'] = tparams
     
-    Nsteps = [20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]
+    Nsteps = [15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]
     
     for order in [2, 3, 4]:
 
@@ -87,6 +87,10 @@
         # get initial values on finest level
         P = MS[0].levels[0].prob
         uinit = P.u_exact(t0)
+        if ii==0:
+          print "Time step: %4.2f" % dt
+          print "Fast CFL number: %4.2f" % (pparams['cs']*dt/P.dx) 
+          print "Slow CFL number: %4.2f" % (pparams['cadv']*dt/P.dx) 
 
         # call main function to get things done...
         uend,stats = mp.run_pfasst(MS, u0=uinit, t0=t0, dt=dt, Tend=Tend)
diff --git a/examples/acoustic_1d_imex/runitererror.py b/examples/acoustic_1d_imex/runitererror.py
@@ -83,6 +83,9 @@
         # get initial values on finest level
         P = MS[0].levels[0].prob
         uinit = P.u_exact(t0)
+        
+        print "Fast CFL number: %4.2f" % (pparams['cs']*dt/P.dx) 
+        print "Slow CFL number: %4.2f" % (pparams['cadv']*dt/P.dx) 
 
         # call main function to get things done...
         uend,stats = mp.run_pfasst(MS, u0=uinit, t0=t0, dt=dt, Tend=Tend)
@@ -113,12 +116,12 @@
     for ii in range(0,np.size(cs_v)):
       x = np.arange(1,lastiter[ii,0])
       y = convrate[ii, 0, 0:lastiter[ii,0]-1]
-      plt.plot(x, y, shape[ii], markersize=fs-2, color=color[ii], label=r'$c_{s}$=%4.2f' % cs_v[ii])       
+      plt.plot(x, y, shape[ii], markersize=fs-2, color=color[ii], label=r'$C_{\rm fast}$=%4.2f' % (cs_v[ii]*dt/P.dx))       
       #plt.plot(x, 0.0*y+avg_convrate[ii,0], '--', color=color[ii])
 
     plt.legend(loc='upper right', fontsize=fs, prop={'size':fs})
     plt.xlabel('Iteration', fontsize=fs)
-    plt.ylabel('Convergence rate', fontsize=fs, labelpad=2)
+    plt.ylabel(r'$|| r^{k+1} ||_{\infty}/|| r^k ||_{\infty}$', fontsize=fs, labelpad=2)
     plt.xlim([0, sparams['maxiter']])
     plt.ylim([0, 0.8])
     plt.yticks(fontsize=fs)
diff --git a/examples/boussinesq_2d_imex/ProblemClass.py b/examples/boussinesq_2d_imex/ProblemClass.py
@@ -106,7 +106,7 @@ def solve_system(self,rhs,factor,u0,t):
         sol, info = LA.gmres( self.Id - factor*self.M, b, x0=u0.values.flatten(), tol=self.gmres_tol, restart=self.gmres_restart, maxiter=self.gmres_maxiter, callback=cb)
         # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
         if factor!=0.0:
-          print "SDC: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+          #print "SDC: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
           self.logger.add(cb.getcounter())
         me        = mesh(self.nvars)
         me.values = unflatten(sol, 4, self.N[0], self.N[1])
diff --git a/examples/boussinesq_2d_imex/plotgmrescounter.py b/examples/boussinesq_2d_imex/plotgmrescounter.py
@@ -0,0 +1,32 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+from pylab import rcParams
+
+from unflatten import unflatten
+
+if __name__ == "__main__":
+  xx = np.load('xaxis.npy')
+  uend = np.load('sdc.npy')
+  udirk2 = np.load('dirk2.npy')
+  udirk4 = np.load('dirk4.npy')
+  utrap   = np.load('trap.npy')
+
+  fs = 8
+  rcParams['figure.figsize'] = 5.0, 2.5
+  fig = plt.figure()
+
+  plt.plot(xx[:,5], uend[2,:,5], '-', color='b', label='SDC')
+  plt.plot(xx[:,5], udirk4[2,:,5], '--', color='g', markersize=fs-2, label='DIRK(4)', dashes=(3,3))
+  plt.plot(xx[:,5], udirk2[2,:,5], '--', color='r', markersize=fs-2, label='DIRK(2)', dashes=(3,3))
+  #plt.plot(xx[:,5], utrap[2,:,5], '--', color='k', markersize=fs-2, label='Trap', dashes=(3,3))
+  plt.legend(loc='lower left', fontsize=fs, prop={'size':fs})
+  plt.yticks(fontsize=fs)
+  plt.xticks(fontsize=fs)
+  plt.xlabel('x [km]', fontsize=fs, labelpad=0)
+  plt.ylabel('Bouyancy', fontsize=fs, labelpad=1)
+  #plt.show()
+  plt.savefig('boussinesq.pdf', bbox_inches='tight')
+
diff --git a/examples/boussinesq_2d_imex/rungmrescounter.py b/examples/boussinesq_2d_imex/rungmrescounter.py
@@ -21,7 +21,7 @@
 
 from unflatten import unflatten
 
-from standard_integrators import dirk
+from standard_integrators import dirk, bdf2, trapezoidal
 
 if __name__ == "__main__":
 
@@ -40,16 +40,14 @@
     swparams['do_LU'] = False
 
     sparams = {}
-    sparams['maxiter'] = 3
+    sparams['maxiter'] = 4
 
-    dirk_order = 2
+    dirk_order = 4
 
     # setup parameters "in time"
     t0     = 0
     Tend   = 3000
-    Nsteps =  500
-    #Tend   = 30
-    #Nsteps =  5
+    Nsteps =  100
     dt = Tend/float(Nsteps)
 
     # This comes as read-in for the problem class
@@ -92,41 +90,63 @@
     cfl_advection    = pparams['u_adv']*dt/P.h[0]
     cfl_acoustic_hor = pparams['c_s']*dt/P.h[0]
     cfl_acoustic_ver = pparams['c_s']*dt/P.h[1]
+    print "Horizontal resolution: %4.2f" % P.h[0]
+    print "Vertical resolution:   %4.2f" % P.h[1]
     print ("CFL number of advection: %4.2f" % cfl_advection)
     print ("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor)
     print ("CFL number of acoustics (vertical):   %4.2f" % cfl_acoustic_ver)
 
-    dirk = dirk(P, dirk_order)
+    dirk4 = dirk(P, 4)
+    dirk2 = dirk(P, 2)
+    trap  = trapezoidal(P)
+    bdf  = bdf2(P)
     u0 = uinit.values.flatten()
-
+    udirk4 = u0
+    udirk2 = u0
+    ubdf  = u0
+    utrap  = u0
     for i in range(0,Nsteps):
-      u0 = dirk.timestep(u0, dt)  
-
+      udirk4 = dirk4.timestep(udirk4, dt)  
+      udirk2 = dirk2.timestep(udirk2, dt)
+      utrap  = trap.timestep(utrap, dt)
+      #if i==0:
+      #  ubdf_new = bdf.firsttimestep(ubdf, dt)
+      #  ubdf_m1  = ubdf
+      #else:        
+      #  ubdf_new  = bdf.timestep(ubdf, ubdf_m1, dt)
+      #ubdf_m1 = ubdf
+      #ubdf    = ubdf_new    
+  
     # call main function to get things done...
     uend,stats = mp.run_pfasst(MS,u0=uinit,t0=t0,dt=dt,Tend=Tend)
-
-    u0 = unflatten(u0, 4, P.N[0], P.N[1])
-
-    fs = 8
-    rcParams['figure.figsize'] = 5.0, 2.5
-    fig = plt.figure()
-
-    plt.plot(P.xx[:,5], uend.values[2,:,5], '-', color='b', label='SDC')
-    plt.plot(P.xx[:,5], u0[2,:,5], '+', color='g', markevery=5, markersize=fs-2, label='DIRK')
-    plt.legend(loc='lower left', fontsize=fs, prop={'size':fs})
-    plt.yticks(fontsize=fs)
-    plt.xticks(fontsize=fs)
-    plt.xlabel('x', fontsize=fs, labelpad=0)
-    plt.ylabel('Bouyancy', fontsize=fs, labelpad=1)
-    #plt.show()
-    plt.savefig('boussinesq.pdf', bbox_inches='tight')
-
-    print " #### Logging report for DIRK #### "
-    print "Number of calls to implicit solver: %5i" % dirk.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % dirk.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(dirk.logger.iterations)/float(dirk.logger.solver_calls))
+    udirk4 = unflatten(udirk4, 4, P.N[0], P.N[1])
+    udirk2 = unflatten(udirk2, 4, P.N[0], P.N[1])
+    print "Norm of final solution by trapezoidal rule: %5.3f" % np.linalg.norm( utrap, np.inf )    
+    utrap  = unflatten(utrap, 4, P.N[0], P.N[1])
+  
+    np.save('xaxis', P.xx)
+    np.save('sdc', uend.values)
+    np.save('dirk2', udirk2)
+    np.save('dirk4', udirk4)
+    np.save('trap', utrap)
     
-    print " #### Logging report for SDC #### "
+    print " #### Logging report for DIRK-4 #### "
+    print "Number of calls to implicit solver: %5i" % dirk4.logger.solver_calls
+    print "Total number of GMRES iterations: %5i" % dirk4.logger.iterations
+    print "Average number of iterations per call: %6.3f" % (float(dirk4.logger.iterations)/float(dirk4.logger.solver_calls))
+
+    print " #### Logging report for DIRK-2 #### "
+    print "Number of calls to implicit solver: %5i" % dirk2.logger.solver_calls
+    print "Total number of GMRES iterations: %5i" % dirk2.logger.iterations
+    print "Average number of iterations per call: %6.3f" % (float(dirk2.logger.iterations)/float(dirk2.logger.solver_calls))
+
+    #print " #### Logging report for BDF2 #### "
+    #print "Number of calls to implicit solver: %5i" % bdf.logger.solver_calls
+    #print "Total number of GMRES iterations: %5i" % bdf.logger.iterations
+    #print "Average number of iterations per call: %6.3f" % (float(bdf.logger.iterations)/float(bdf.logger.solver_calls))
+
+    print " #### Logging report for SDC-(%1i,%1i) #### " % (swparams['num_nodes'], sparams['maxiter'])
     print "Number of calls to implicit solver: %5i" % P.logger.solver_calls
     print "Total number of GMRES iterations: %5i" % P.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(P.logger.iterations)/float(P.logger.solver_calls))
+    print "Average number of iterations per call: %6.3f" % (float(P.logger.iterations)/float(P.logger.solver_calls))  
+
diff --git a/examples/boussinesq_2d_imex/standard_integrators.py b/examples/boussinesq_2d_imex/standard_integrators.py
@@ -4,6 +4,79 @@
 import scipy.sparse as sp
 from ProblemClass import Callback, logging, boussinesq_2d_imex
 
+#
+# Trapezoidal rule
+#
+class trapezoidal:
+ 
+  def __init__(self, problem, alpha=0.5):
+    assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object"
+    self.Ndof = np.shape(problem.M)[0]
+    self.order = 2
+    self.logger = logging()
+    self.problem = problem
+    self.alpha = alpha
+
+  def timestep(self, u0, dt):
+    B_trap   = sp.eye(self.Ndof) + self.alpha*dt*(self.problem.D_upwind + self.problem.M)
+    b         = B_trap.dot(u0)
+    return self.f_solve(b, alpha=(1.0-self.alpha)*dt, u0 = u0)
+
+  # 
+  # Returns f(u) = c*u
+  #  
+  def f(self,u):
+    return self.problem.D_upwind.dot(u)+self.problem.M.dot(u)
+    
+  
+  #
+  # Solves (Id - alpha*c)*u = b for u
+  #  
+  def f_solve(self, b, alpha, u0):
+    cb = Callback()
+    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    if alpha!=0.0:
+      #print "BDF-2: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+      self.logger.add(cb.getcounter())    
+    return sol
+
+#
+# A BDF-2 implicit two-step method
+#
+class bdf2:
+
+  def __init__(self, problem):
+    assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object"
+    self.Ndof = np.shape(problem.M)[0]
+    self.order = 2
+    self.logger = logging()
+    self.problem = problem
+
+  def firsttimestep(self, u0, dt):
+    return self.f_solve(b = u0, alpha = dt, u0 = u0)
+
+  def timestep(self, u0, um1, dt):
+    b = (4.0/3.0)*u0 - (1.0/3.0)*um1
+    return self.f_solve(b = b, alpha = (2.0/3.0)*dt, u0 = u0)
+
+  # 
+  # Returns f(u) = c*u
+  #  
+  def f(self,u):
+    return self.problem.D_upwind.dot(u)+self.problem.M.dot(u)
+    
+  
+  #
+  # Solves (Id - alpha*c)*u = b for u
+  #  
+  def f_solve(self, b, alpha, u0):
+    cb = Callback()
+    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    if alpha!=0.0:
+      #print "BDF-2: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+      self.logger.add(cb.getcounter())    
+    return sol
+
 #
 # A diagonally implicit Runge-Kutta method of order 2, 3 or 4
 #
@@ -117,6 +190,6 @@ def f_solve(self, b, alpha, u0):
     cb = Callback()
     sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
     if alpha!=0.0:
-      print "DIRK: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+      #print "DIRK-%1i: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( self.order, cb.getcounter(), cb.getresidual() )
       self.logger.add(cb.getcounter())    
     return sol
diff --git a/examples/fwsw/ProblemClass.py b/examples/fwsw/ProblemClass.py
diff --git a/examples/fwsw/__init__.py b/examples/fwsw/__init__.py
diff --git a/examples/fwsw/playground.py b/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
diff --git a/examples/fwsw/plot_stability.py b/examples/fwsw/plot_stability.py
@@ -0,0 +1,89 @@
+from pySDC import Level as lvl
+from pySDC import Hooks as hookclass
+from pySDC import CollocationClasses as collclass
+from pySDC import Step as stepclass
+
+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.fwsw.ProblemClass import swfw_scalar 
+import numpy as np
+
+from pylab import rcParams
+import matplotlib.pyplot as plt
+from subprocess import call
+
+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(2.0, 8.0, N_f)
+
+    pparams = {}
+    # 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'] = 4
+    K = 5
+    
+    #
+    # ...this is functionality copied from test_imexsweeper. Ideally, it should be available in one place.
+    #
+    step = stepclass.step(params={})
+    L = lvl.level(problem_class=swfw_scalar, problem_params=pparams, dtype_u=mesh, dtype_f=rhs_imex_mesh, sweeper_class=imex, sweeper_params=swparams, level_params={}, hook_class=hookclass.hooks, id="stability")
+    step.register_level(L)
+    step.status.dt   = 1.0
+    step.status.time = 0.0
+    u0 = step.levels[0].prob.u_exact(step.status.time)
+    step.init_step(u0)
+    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:]
+
+    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
+
+    ###
+    rcParams['figure.figsize'] = 2.5, 2.5
+    fs = 8
+    fig  = plt.figure()
+    #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, 1.05])
+#    levels = np.array([1.0])
+    CS1 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), levels, colors='k',linestyles='dashed')
+    CS2 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), [1.0], colors='k')
+    plt.clabel(CS1, fontsize=fs-2)
+    plt.clabel(CS2, fontsize=fs-2)
+    #plt.plot([0, 2], [0, 2], color='k', linewidth=1)
+    plt.gca().set_xticks([0.0, 1.0, 2.0, 3.0])
+    plt.gca().tick_params(axis='both', which='both', labelsize=fs)
+    plt.xlim([np.min(lambda_s.imag), np.max(lambda_s.imag)])
+    plt.xlabel('$\Delta t \lambda_{slow}$', fontsize=fs, labelpad=2.0)
+    plt.ylabel('$\Delta t \lambda_{fast}$', fontsize=fs)
+    plt.title(r'$M=%1i$, $K=%1i$' % (swparams['num_nodes'],K), fontsize=fs)
+    filename = 'sdc-fwsw-stability-K'+str(K)+'-M'+str(swparams['num_nodes'])+'.pdf'
+    fig.savefig(filename, bbox_inches='tight')
+    call(["pdfcrop", filename, filename])
+
diff --git a/examples/fwsw/plot_stifflimit_specrad.py b/examples/fwsw/plot_stifflimit_specrad.py
