Merge pull request #50 from danielru/fwsw_boussinesq

Boussinesq

Author: pancetta

examples/boussinesq_2d_imex/HookClass.py

@@ -20,8 +20,11 @@ def __init__(self):
 
         # add figure object for further use
         #        self.fig = plt.figure(figsize=(18,6))
-        self.fig = plt.figure(figsize=(15,5))
-        self.counter = 0
+        #self.fig = plt.figure(figsize=(15,5))
+        #self.counter = 0
+
+    def dump_sweep(self,status):
+        pass
 
     def dump_step(self,status):
         """

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 "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])
@@ -197,8 +197,3 @@ def u_exact(self,t):
         me.values[2,:,:] = dtheta*np.sin( np.pi*self.zz/H )/( 1.0 + np.square(self.xx - x_c)/(a*a))
         me.values[3,:,:] = 0.0*self.xx
         return me
-        
-    def report_log(self):
-      print "Number of calls to implicit solver: %5i" % self.logger.solver_calls
-      print "Total number of iterations: %5i" % self.logger.iterations
-      print "Average number of iterations per call: %6.3f" % (float(self.logger.iterations)/float(self.logger.solver_calls))

examples/boussinesq_2d_imex/rungmrescounter.py

@@ -0,0 +1,132 @@
+
+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 pylab import rcParams
+
+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-15
+    
+    swparams = {}
+    swparams['collocation_class'] = collclass.CollGaussLobatto
+    swparams['num_nodes'] = 3
+    swparams['do_LU'] = False
+
+    sparams = {}
+    sparams['maxiter'] = 3
+
+    dirk_order = 2
+
+    # 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,300,20)]
+    #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'] = [500]
+    pparams['gmres_restart'] = [10]
+    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)
+
+    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)
+
+    dirk = dirk(P, dirk_order)
+    u0 = uinit.values.flatten()
+
+    for i in range(0,Nsteps):
+      u0 = dirk.timestep(u0, dt)  
+
+    # 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))
+    
+    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))

examples/boussinesq_2d_imex/standard_integrators.py

@@ -0,0 +1,122 @@
+import numpy as np
+import math
+import scipy.sparse.linalg as LA
+import scipy.sparse as sp
+from ProblemClass import Callback, logging, boussinesq_2d_imex
+
+#
+# A diagonally implicit Runge-Kutta method of order 2, 3 or 4
+#
+class dirk:
+
+  def __init__(self, problem, order):
+
+    assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object"
+    self.Ndof = np.shape(problem.M)[0]
+    self.order = order
+    self.logger = logging()
+    self.problem = problem
+
+    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.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 "DIRK: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
+      self.logger.add(cb.getcounter())    
+    return sol