removed run_standard; everything is now included in runmultiscale script

Author: Daniel Ruprecht

examples/acoustic_1d_imex/run_standard.py

@@ -12,11 +12,12 @@
 from pySDC import Log
 from pySDC.Stats import grep_stats, sort_stats
 
-# Sharpclaw imports
-#from clawpack import pyclaw
-#from clawpack import riemann
+from standard_integrators import bdf2, dirk, trapezoidal
+
 from matplotlib import pyplot as plt
+from pylab import rcParams
 
+fs = 8
 
 if __name__ == "__main__":
 
@@ -55,8 +56,11 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLobatto
-    description['num_nodes']         = 3
+    description['collocation_class'] = collclass.CollGaussRadau_Right
+    if sparams['maxiter']==2:
+      description['num_nodes']         = 2
+    else:
+      description['num_nodes']         = 3
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
@@ -73,25 +77,64 @@
     # call main function to get things done...
     uend,stats = mp.run_pfasst(MS,u0=uinit,t0=t0,dt=dt,Tend=Tend)
 
-    # compute exact solution and compare
-    uex = P.u_exact(Tend)
+    # instantiate standard integrators to be run for comparison
+    trap  = trapezoidal( P.A+P.Dx, 0.5 )
+    bdf2  = bdf2( P.A+P.Dx)
+    dirk  = dirk( P.A+P.Dx, sparams['maxiter'])
+
+    y0_tp = np.concatenate( (uinit.values[0,:], uinit.values[1,:]) )
+    y0_bdf  = y0_tp
+    y0_dirk = y0_tp
+      
+    # Perform 154 time steps with standard integrators
+    for i in range(0,154):  
+
+      # trapezoidal rule step
+      ynew_tp = trap.timestep(y0_tp, dt)
+
+      # BDF-2 scheme
+      if i==0:
+        ynew_bdf = bdf2.firsttimestep( y0_bdf, dt)
+        ym1_bdf = y0_bdf
+      else:
+        ynew_bdf = bdf2.timestep( y0_bdf, ym1_bdf, dt)
 
-    fig = plt.figure(figsize=(8,8))
+      # DIRK scheme
+      ynew_dirk = dirk.timestep(y0_dirk, dt)
+
+      y0_tp   = ynew_tp
+      ym1_bdf = y0_bdf
+      y0_bdf  = ynew_bdf
+      y0_dirk = ynew_dirk
+
+    # Finished running standard integrators
+      unew_tp, pnew_tp     = np.split(ynew_tp, 2)
+      unew_bdf, pnew_bdf   = np.split(ynew_bdf, 2)
+      unew_dirk, pnew_dirk = np.split(ynew_dirk, 2)
+
+    rcParams['figure.figsize'] = 5, 2.5
+    fig = plt.figure()
 
     sigma_0 = 0.1
     k       = 7.0*2.0*np.pi
     x_0     = 0.75
     x_1     = 0.25
 
-    #plt.plot(P.mesh, uex.values[0,:],  '+', color='b', label='u (exact)')
-    plt.plot(P.mesh, uend.values[1,:], '-', color='b', label='SDC')
+    plt.plot(P.mesh, pnew_tp,  '-', color='c', label='Trapezoidal')
+    plt.plot(P.mesh, uend.values[1,:], '-', color='b', label='SDC('+str(sparams['maxiter'])+')')
+    plt.plot(P.mesh, pnew_bdf, '-', color='r', label='BDF-2')
+    plt.plot(P.mesh, pnew_dirk, color='g', label='DIRK('+str(dirk.order)+')')
     #plt.plot(P.mesh, uex.values[1,:],  '+', color='r', label='p (exact)')
     #plt.plot(P.mesh, uend.values[1,:], '-', color='b', linewidth=2.0, label='p (SDC)')
-    p_slow = np.exp(-np.square(P.mesh-x_0-pparams['cadv']*Tend)/(sigma_0*sigma_0))
-    plt.plot(P.mesh, p_slow, '-', color='r', markersize=4, label='slow mode')
-    plt.legend(loc=2)
-    plt.xlim([0, 1])
-    plt.ylim([-0.1, 1.1])
+
+    p_slow = np.exp(-np.square( np.mod( P.mesh-pparams['cadv']*Tend, 1.0 ) -x_0 )/(sigma_0*sigma_0))
+    plt.plot(P.mesh, p_slow, '+', color='k', markersize=fs-2, label='Slow mode', markevery=10)
+    plt.xlabel('x', fontsize=fs, labelpad=0)
+    plt.ylabel('Pressure', fontsize=fs, labelpad=0)
+    fig.gca().set_xlim([0, 1.0])
+    fig.gca().set_ylim([-0.5, 1.1])
+    fig.gca().tick_params(axis='both', labelsize=fs)
+    plt.legend(loc='upper left', fontsize=fs, prop={'size':fs})
     fig.gca().grid()
-    plt.show()
-    #plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')
+    #plt.show()
+    plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')

examples/acoustic_1d_imex/runmultiscale.py

@@ -0,0 +1,148 @@
+import numpy as np
+import math
+import scipy.sparse.linalg as LA
+import scipy.sparse as sp
+
+#
+# Trapezoidal rule
+#
+class trapezoidal:
+ 
+  def __init__(self, M, alpha=0.5):
+    assert np.shape(M)[0]==np.shape(M)[1], "Matrix M must be quadratic"
+    self.Ndof = np.shape(M)[0]
+    self.M = M
+    self.alpha = alpha
+
+  def timestep(self, u0, dt):
+    M_trap   = sp.eye(self.Ndof) - self.alpha*dt*self.M
+    B_trap   = sp.eye(self.Ndof) + (1.0-self.alpha)*dt*self.M
+    b = B_trap.dot(u0)
+    return LA.spsolve(M_trap, b)
+#
+# A BDF-2 implicit two-step method
+#
+class bdf2:
+
+  def __init__(self, M):
+    assert np.shape(M)[0]==np.shape(M)[1], "Matrix M must be quadratic"
+    self.Ndof = np.shape(M)[0]
+    self.M = M
+
+  def firsttimestep(self, u0, dt):
+    b = u0
+    L = sp.eye(self.Ndof) - dt*self.M
+    return LA.spsolve(L, b)
+
+  def timestep(self, u0, um1, dt):
+    b = (4.0/3.0)*u0 - (1.0/3.0)*um1
+    L = sp.eye(self.Ndof) - (2.0/3.0)*dt*self.M
+    return LA.spsolve(L, b)
+#
+# A diagonally implicit Runge-Kutta method of order 2, 3 or 4
+#
+class dirk:
+
+  def __init__(self, M, order):
+
+    assert np.shape(M)[0]==np.shape(M)[1], "Matrix M must be quadratic"
+    self.Ndof = np.shape(M)[0]
+    self.M = M
+    self.order = order
+
+    assert self.order in [2,22,3,4], 'Order must be 2,22,3,4'
+    
+    if (self.order==2):
+      self.nstages = 1
+      self.A       = np.zeros((1,1))
+      self.A[0,0]  = 0.5
+      self.tau     = [0.5]
+      self.b       = [1.0]
+    
+    if (self.order==22):
+      self.nstages = 2
+      self.A       = np.zeros((2,2))
+      self.A[0,0]  = 1.0/3.0
+      self.A[1,0]  = 1.0/2.0
+      self.A[1,1]  = 1.0/2.0
+      
+      self.tau     = np.zeros(2)
+      self.tau[0]  = 1.0/3.0
+      self.tau[1]  = 1.0
+
+      self.b       = np.zeros(2)
+      self.b[0]    = 3.0/4.0
+      self.b[1]    = 1.0/4.0
+    
+    
+    if (self.order==3):
+      self.nstages = 2 
+      self.A       = np.zeros((2,2))
+      self.A[0,0]  = 0.5 + 1.0/( 2.0*math.sqrt(3.0) )
+      self.A[1,0] = -1.0/math.sqrt(3.0)
+      self.A[1,1] = self.A[0,0]
+      
+      self.tau    = np.zeros(2)
+      self.tau[0] = 0.5 + 1.0/( 2.0*math.sqrt(3.0) )
+      self.tau[1] = 0.5 - 1.0/( 2.0*math.sqrt(3.0) )
+      
+      self.b     = np.zeros(2)
+      self.b[0]  = 0.5
+      self.b[1]  = 0.5
+      
+    if (self.order==4):
+      self.nstages = 3
+      alpha = 2.0*math.cos(math.pi/18.0)/math.sqrt(3.0)
+      
+      self.A      = np.zeros((3,3))
+      self.A[0,0] = (1.0 + alpha)/2.0
+      self.A[1,0] = -alpha/2.0
+      self.A[1,1] = self.A[0,0]
+      self.A[2,0] = (1.0 + alpha)
+      self.A[2,1] =  -(1.0 + 2.0*alpha)
+      self.A[2,2] = self.A[0,0]
+      
+      self.tau    = np.zeros(3)
+      self.tau[0] = (1.0 + alpha)/2.0
+      self.tau[1] = 1.0/2.0
+      self.tau[2] = (1.0 - alpha)/2.0
+      
+      self.b      = np.zeros(3)
+      self.b[0]   = 1.0/(6.0*alpha*alpha)
+      self.b[1]   = 1.0 - 1.0/(3.0*alpha*alpha)
+      self.b[2]   = 1.0/(6.0*alpha*alpha)
+       
+    self.stages  = np.zeros((self.nstages,self.Ndof))
+
+  def timestep(self, u0, dt):
+    
+      uend           = u0
+      for i in range(0,self.nstages):  
+        
+        b = u0
+        
+        # Compute right hand side for this stage's implicit step
+        for j in range(0,i):
+          b = b + self.A[i,j]*dt*self.f(self.stages[j,:])
+        
+        # Implicit solve for current stage    
+        self.stages[i,:] = self.f_solve( b, dt*self.A[i,i] )
+        
+        # Add contribution of current stage to final value
+        uend = uend + self.b[i]*dt*self.f(self.stages[i,:])
+        
+      return uend
+      
+  # 
+  # Returns f(u) = c*u
+  #  
+  def f(self,u):
+    return self.M.dot(u)
+    
+  
+  #
+  # Solves (Id - alpha*c)*u = b for u
+  #  
+  def f_solve(self, b, alpha):
+    L = sp.eye(self.Ndof) - alpha*self.M
+    return LA.spsolve(L, b)