Merge pull request #82 from danielru/add/dirk_order_five

add fifth order DIRK method to 1d acoustic IMEX and boussinesq example

Author: pancetta

examples/acoustic_1d_imex/runmultiscale.py

@@ -32,12 +32,14 @@
     lparams['restol'] = 1E-10
 
     sparams = {}
-    sparams['maxiter'] = 4
+    sparams['maxiter'] = 5
 
     # setup parameters "in time"
     t0   = 0.0
     Tend = 3.0
-    dt = Tend/float(154)
+    nsteps = 154 # 154 is value in Vater et al.
+    nsteps = 2*154    
+    dt = Tend/float(nsteps)
 
     # This comes as read-in for the problem class
     pparams = {}
@@ -87,8 +89,8 @@
     y0_dirk = y0_tp.astype('complex')
     y0_imex = y0_tp.astype('complex')
 
-    # Perform 154 time steps with standard integrators
-    for i in range(0,154):  
+    # Perform time steps with standard integrators
+    for i in range(0,nsteps):  
 
       # trapezoidal rule step
       ynew_tp = trap.timestep(y0_tp, dt)

examples/acoustic_1d_imex/standard_integrators.py

@@ -241,7 +241,7 @@ def __init__(self, M, order):
     self.M = M
     self.order = order
 
-    assert self.order in [2,22,3,4], 'Order must be 2,22,3,4'
+    assert self.order in [2,22,3,4,5], 'Order must be 2,22,3,4'
 
     if (self.order==2):
       self.nstages = 1
@@ -302,6 +302,27 @@ def __init__(self, M, order):
       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)
+
+    if (self.order==5):
+      # From A. H. Al-Rabeh, "OPTIMAL ORDER DIAGONALLY IMPLICIT RUNGE-KUTTA METHODS"
+      self.nstages = 4
+      self.A   = np.zeros((4,4))
+      self.A[1,0] = 0.1090390091
+      self.A[1,1] = 0.1090390091
+      self.A[2,0] = 0.0177359481 
+      self.A[2,1] = 0.4277445474
+      self.A[2,2] = 0.1090390091
+      self.A[3,0] = 0.1173343519
+      self.A[3,1] = 0.2044057169
+      self.A[3,2] = 0.4601819131
+      self.A[3,3] = 0.1090390091
+
+      self.tau = np.zeros(4)
+      self.b   = np.zeros(4)
+      self.b[0] = 0.0707307044
+      self.b[1] = 0.3078968440
+      self.b[2] = 0.3589454736
+      self.b[3] = 0.2624269779
 
     self.stages  = np.zeros((self.nstages,self.Ndof), dtype='complex')
 

examples/boussinesq_2d_imex/rungmrescounter.py

@@ -46,8 +46,10 @@
 
     # setup parameters "in time"
     t0     = 0
-    Tend   = 3000 
-    Nsteps = 500
+
+    Tend   = 3000   
+    Nsteps =  100
+
     dt = Tend/float(Nsteps)
 
     # This comes as read-in for the problem class
@@ -63,7 +65,7 @@
     pparams['gmres_maxiter'] = [500]
     pparams['gmres_restart'] = [10]
     pparams['gmres_tol_limit'] = [1e-5]
-    pparams['gmres_tol_factor'] = [0.1]
+    pparams['gmres_tol_factor'] = [0.05]
 
     # This comes as read-in for the transfer operations
     tparams = {}
@@ -99,6 +101,7 @@
     method_split = 'MIS4_4'
 #   method_split = 'RK3'
     splitp = SplitExplicit(P, method_split, pparams) 
+
     u0 = uinit.values.flatten()
     usplit = np.copy(u0)
 #   for i in range(0,Nsteps):
@@ -107,7 +110,7 @@
       usplit = splitp.timestep(usplit, dt/2)  
     print(np.linalg.norm(usplit))
 
-    dirkp = dirk(P, np.min([4,dirk_order]))
+    dirkp = dirk(P, dirk_order)
     udirk = np.copy(u0)
     print("Running DIRK ....")
     print(np.linalg.norm(udirk))  

examples/boussinesq_2d_imex/standard_integrators.py

@@ -423,8 +423,8 @@ def __init__(self, problem, order):
     self.logger = logging()
     self.problem = problem
 
-    assert self.order in [2,22,3,4], 'Order must be 2,22,3,4'
-   
+    assert self.order in [2,22,3,4,5], 'Order must be 2,22,3,4'
+    
     if (self.order==2):
       self.nstages = 1
       self.A       = np.zeros((1,1))
@@ -485,6 +485,26 @@ def __init__(self, problem, order):
       self.b[1]   = 1.0 - 1.0/(3.0*alpha*alpha)
       self.b[2]   = 1.0/(6.0*alpha*alpha)
 
+    if (self.order==5):
+      # From A. H. Al-Rabeh, "OPTIMAL ORDER DIAGONALLY IMPLICIT RUNGE-KUTTA METHODS"
+      self.nstages = 4
+      self.A   = np.zeros((4,4))
+      self.A[1,0] = 0.1090390091
+      self.A[1,1] = 0.1090390091
+      self.A[2,0] = 0.0177359481 
+      self.A[2,1] = 0.4277445474
+      self.A[2,2] = 0.1090390091
+      self.A[3,0] = 0.1173343519
+      self.A[3,1] = 0.2044057169
+      self.A[3,2] = 0.4601819131
+      self.A[3,3] = 0.1090390091
+      self.tau = np.zeros(4)
+      self.b   = np.zeros(4)
+      self.b[0] = 0.0707307044
+      self.b[1] = 0.3078968440
+      self.b[2] = 0.3589454736
+      self.b[3] = 0.2624269779
+
     self.stages  = np.zeros((self.nstages,self.Ndof))
 
   def timestep(self, u0, dt):