added dirk method to run_standard for comparison against fwsw-SDC

Author: Daniel Ruprecht

examples/acoustic_1d_imex/dirk.py

@@ -0,0 +1,110 @@
+import numpy as np
+import math
+import scipy.sparse.linalg as LA
+import scipy.sparse as sp
+
+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)

examples/acoustic_1d_imex/run_standard.py

@@ -4,6 +4,7 @@
 import os
 from matplotlib import pyplot as plt
 from buildWave1DMatrix import getWave1DMatrix, getWave1DAdvectionMatrix
+from dirk import dirk
 
 sigma_0 = 0.1
 k       = 7.0*2*np.pi
@@ -30,6 +31,8 @@ def u(x,t, multiscale):
 Id, A  = getWave1DMatrix(nvars, dx, ['periodic','periodic'], ['periodic','periodic'])
 A      = -cs*A
 
+dirk     = dirk( (A+Dx), 4)
+
 M_ieuler = Id - dt*(A + Dx)
 
 M_bdf    = Id - (2.0/3.0)*dt*(A + Dx)
@@ -38,11 +41,11 @@ def u(x,t, multiscale):
 M_trap   = Id - alpha*dt*(A+Dx)
 B_trap   = Id + (1-alpha)*dt*(A+Dx)
 
-u0, p0 = u(mesh, 0.0, 1.0)
-y0_ie  = np.concatenate( (u0, p0) )
-y0_tp  = y0_ie
-
-y0_bdf = y0_ie
+u0, p0  = u(mesh, 0.0, 1.0)
+y0_ie   = np.concatenate( (u0, p0) )
+y0_tp   = y0_ie
+y0_bdf  = y0_ie
+y0_dirk = y0_ie
 
 fig = plt.figure(figsize=(8,8))
 
@@ -62,15 +65,20 @@ def u(x,t, multiscale):
     b_bdf    = (4.0/3.0)*y0_bdf - (1.0/3.0)*ym1_bdf
     ynew_bdf = LA.spsolve(M_bdf, b_bdf)
 
-  unew_ie, pnew_ie = np.split(ynew_ie, 2)
-  unew_tp, pnew_tp = np.split(ynew_tp, 2)
+  # DIRK scheme
+  ynew_dirk = dirk.timestep(y0_dirk, dt)
+
+  unew_ie, pnew_ie   = np.split(ynew_ie, 2)
+  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)
   uex, pex = u(mesh, float(i+1)*dt, 0.0)
 
   if i==Nsteps-1:
     fig.gca().clear()
     plt.plot(mesh, pnew_bdf, 'b', label='BDF-2')
-    plt.plot(mesh, pnew_tp, 'r', label='Trapezoidal')
+    #plt.plot(mesh, pnew_tp, 'r', label='Trapezoidal')
+    plt.plot(mesh, pnew_dirk, 'r', label='DIRK')
     plt.plot(mesh, pex, 'k', label='Slow Mode')
     fig.gca().set_xlim([0, 1.0])
     fig.gca().set_ylim([-0.5, 1.1])
@@ -89,6 +97,8 @@ def u(x,t, multiscale):
   y0_tp   = ynew_tp
   ym1_bdf = y0_bdf
   y0_bdf  = ynew_bdf
+  y0_dirk = ynew_dirk
+
 plt.show()
 #lt.gcf().savefig('final.pdf', bbox_inches='tight')
 #os.system('ffmpeg -r 25 -i images/standard-%03d.jpeg -vcodec libx264 -crf 25 movie.avi')