added functions for standard integrators

Daniel Ruprecht · Daniel Ruprecht · commit e93845d167a2 · 2015-12-04T11:51:32.000Z
diff --git a/examples/boussinesq_2d_imex/rungmrescounter.py b/examples/boussinesq_2d_imex/rungmrescounter.py
@@ -0,0 +1,121 @@
+
+from pySDC import CollocationClasses as collclass
+
+import numpy as np
+
+from ProblemClass import boussinesq_2d_imex
+from examples.boussinesq_2d_imex.TransferClass import mesh_to_mesh_2d
+from examples.boussinesq_2d_imex.HookClass import plot_solution
+
+from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
+from pySDC.sweeper_classes.imex_1st_order import imex_1st_order
+import pySDC.PFASST_stepwise as mp
+from pySDC import Log
+from pySDC.Stats import grep_stats, sort_stats
+
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+from unflatten import unflatten
+
+from standard_integrators import dirk
+
+if __name__ == "__main__":
+
+    # set global logger (remove this if you do not want the output at all)
+    logger = Log.setup_custom_logger('root')
+
+    num_procs = 1
+
+    # This comes as read-in for the level class
+    lparams = {}
+    lparams['restol'] = 1E-8
+    
+    swparams = {}
+    swparams['collocation_class'] = collclass.CollGaussLobatto
+    swparams['num_nodes'] = 3
+    swparams['do_LU'] = True
+
+    sparams = {}
+    sparams['maxiter'] = 4
+
+    dirk_order = 4
+
+    # setup parameters "in time"
+    t0     = 0
+    Tend   = 3000
+    Nsteps =  500
+    #Tend   = 30
+    #Nsteps =  5
+    dt = Tend/float(Nsteps)
+
+    # This comes as read-in for the problem class
+    pparams = {}
+    pparams['nvars']    = [(4,450,30)]
+    #pparams['nvars']    = [(4,150,10)]
+    pparams['u_adv']    = 0.02
+    pparams['c_s']      = 0.3
+    pparams['Nfreq']    = 0.01
+    pparams['x_bounds'] = [(-150.0, 150.0)]
+    pparams['z_bounds'] = [(   0.0,  10.0)]
+    pparams['order']    = [4] # [fine_level, coarse_level]
+    pparams['order_upw'] = [5]
+    pparams['gmres_maxiter'] = [50]
+    pparams['gmres_restart'] = [20]
+    pparams['gmres_tol']     = [1e-6]
+
+    # This comes as read-in for the transfer operations
+    tparams = {}
+    tparams['finter'] = False
+
+    # Fill description dictionary for easy hierarchy creation
+    description = {}
+    description['problem_class']     = boussinesq_2d_imex
+    description['problem_params']    = pparams
+    description['dtype_u']           = mesh
+    description['dtype_f']           = rhs_imex_mesh
+    description['sweeper_params']    = swparams
+    description['sweeper_class']     = imex_1st_order
+    description['level_params']      = lparams
+    description['hook_class']        = plot_solution
+
+    # quickly generate block of steps
+    MS = mp.generate_steps(num_procs,sparams,description)
+
+    # get initial values on finest level
+    P = MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+
+    dirk = dirk(P.D_upwind + P.M, dirk_order, pparams['gmres_maxiter'], pparams['gmres_restart'], pparams['gmres_tol'])
+    u0 = uinit.values.flatten()
+
+    for i in range(0,Nsteps):
+      u0 = dirk.timestep(u0, dt)
+    
+    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 ("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)
+
+    # 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])
+    fig = plt.figure()
+
+    plt.plot(P.xx[:,5], uend.values[2,:,5], color='r', marker='+', markevery=3)
+    plt.plot(P.xx[:,5], u0[2,:,5], color='b', markevery=5)
+    plt.show()
+
+    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))
+    
+    print " #### Logging report for SDC #### "
+    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))
diff --git a/examples/boussinesq_2d_imex/standard_integrators.py b/examples/boussinesq_2d_imex/standard_integrators.py
@@ -0,0 +1,160 @@
+import numpy as np
+import math
+import scipy.sparse.linalg as LA
+import scipy.sparse as sp
+from ProblemClass import Callback, logging
+#
+# 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, gmres_maxiter, gmres_restart, gmres_tol):
+
+    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
+    self.logger = logging()
+    self.gmres_maxiter = gmres_maxiter
+    self.gmres_restart = gmres_restart
+    self.gmres_tol = gmres_tol
+
+    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    
+        #if i==0:
+        self.stages[i,:] = self.f_solve( b, dt*self.A[i,i] , u0 )
+        #else:
+        #  self.stages[i,:] = self.f_solve( b, dt*self.A[i,i] , self.stages[i-1,:] )
+        
+        # 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, u0):
+    cb = Callback()
+    L = sp.eye(self.Ndof) - alpha*self.M
+    sol, info = LA.gmres( L, b, x0=u0, tol=self.gmres_tol, restart=self.gmres_restart, maxiter=self.gmres_maxiter, callback=cb)
+    if alpha!=0.0:
+      print "DIRK: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+      self.logger.add(cb.getcounter())    
+    return sol
