Merge pull request #83 from danielru/fix/plots

Fix/plots

Author: pancetta

.gitignore

@@ -2,6 +2,11 @@
 /doc/latex
 /doc/html
 
+*.npz
+*.npy
+*.pdf
+*.png
+
 # Created by https://www.gitignore.io
 
 ### Python ###

examples/acoustic_1d_imex/buildFDMatrix.py

@@ -54,6 +54,10 @@ def getMatrix(N, dx, bc_left, bc_right, order):
     stencil = [1.0, -8.0, 0.0, 8.0, -1.0]
     range   = [ -2,   -1,   0,   1,    2]
     coeff   = 1.0/12.0
+  elif order==6:
+    stencil = [-1.0, 9.0, -45.0, 0.0, 45.0, -9.0, 1.0]
+    range   =[ -3, -2, -1, 0, 1, 2, 3]
+    coeff   = 1.0/60.0
 
   A       = sp.diags(stencil, range, shape=(N,N))
   A       = sp.lil_matrix(A)
@@ -73,13 +77,28 @@ def getMatrix(N, dx, bc_left, bc_right, order):
       A[0,N-1] = stencil[1]
       A[1,N-1] = stencil[0]
 
+    elif order==6:
+      A[0,N-3] = stencil[0]
+      A[0,N-2] = stencil[1]
+      A[0,N-1] = stencil[2]
+      A[1,N-2] = stencil[0]
+      A[1,N-1] = stencil[1]
+      A[2,N-1] = stencil[0]
+      
   if bc_right in ['periodic']:
     if order==2:
       A[N-1,0] = stencil[2]
     elif order==4:
       A[N-2,0] = stencil[4]
       A[N-1,0] = stencil[3]
       A[N-1,1] = stencil[4]
+    elif order==6:
+      A[N-3,0] = stencil[6]
+      A[N-2,0] = stencil[5]
+      A[N-2,1] = stencil[6]
+      A[N-1,0] = stencil[4]
+      A[N-1,1] = stencil[5]
+      A[N-1,2] = stencil[6]
 
   #
   # Neumann boundary conditions

examples/acoustic_1d_imex/buildWave1DMatrix.py

@@ -4,7 +4,7 @@
 import scipy.sparse as sp
 from buildFDMatrix import getMatrix, getHorizontalDx, getBCLeft, getBCRight
 
-wave_order = 4
+wave_order = 6
 
 def getWave1DMatrix(N, dx, bc_left, bc_right):
 

examples/acoustic_1d_imex/plot_dispersion.py

@@ -39,11 +39,22 @@ def findomega(stab_fh):
     # the following are not used in the computation
     pparams['lambda_s'] = np.array([0.0])
     pparams['lambda_f'] = np.array([0.0])
+    
     pparams['u0'] = 1.0
     swparams = {}
-    swparams['collocation_class'] = collclass.CollGaussLegendre
+    
+    ### SET TYPE OF QUADRATURE NODES ###
+    #swparams['collocation_class'] = collclass.CollGaussLobatto
+    #swparams['collocation_class'] = collclass.CollGaussLegendre
+    swparams['collocation_class'] = collclass.CollGaussRadau_Right
+    
+    ### SET NUMBER OF NODES ###
     swparams['num_nodes'] = 3
-    K = 4
+    
+    ### SET NUMBER OF ITERATIONS ###
+    K = 5
+    
+    ### ORDER OF DIRK/IMEX IS EQUAL TO NUMBER OF ITERATIONS AND THUS ORDER OF SDC ###
     dirk_order = K
 
     c_speed = 1.0
@@ -66,7 +77,7 @@ def findomega(stab_fh):
     QE = level.sweep.QE[1:,1:]
     QI = level.sweep.QI[1:,1:]
     Q  = level.sweep.coll.Qmat[1:,1:]
-    Nsamples = 30
+    Nsamples = 15
     k_vec = np.linspace(0, np.pi, Nsamples+1, endpoint=False)
     k_vec = k_vec[1:]
     phase = np.zeros((3,Nsamples))
@@ -100,7 +111,7 @@ def findomega(stab_fh):
       # For testing, insert exact stability function exp(-dt*i*k*(Cs+Uadv)
       #stab_fh = la.expm(Cs+Uadv)
 
-      dirkts = dirk(Cs+Uadv, np.min([4,dirk_order]))
+      dirkts = dirk(Cs+Uadv, dirk_order)
       stab_fh1 = dirkts.timestep(y1, 1.0)
       stab_fh2 = dirkts.timestep(y2, 1.0)
       stab_dirk = np.column_stack((stab_fh1, stab_fh2))
@@ -128,7 +139,7 @@ def findomega(stab_fh):
     fig  = plt.figure()
     plt.plot(k_vec, (U_speed+c_speed)+np.zeros(np.size(k_vec)), '--', color='k', linewidth=1.5, label='Exact')
     plt.plot(k_vec, phase[1,:], '-', color='g', linewidth=1.5, label='DIRK('+str(dirkts.order)+')')
-    plt.plot(k_vec, phase[2,:], '-+', color='r', linewidth=1.5, label='RK-IMEX('+str(rkimex.order)+')', markevery=5, mew=1.0)
+    plt.plot(k_vec, phase[2,:], '-+', color='r', linewidth=1.5, label='IMEX('+str(rkimex.order)+')', markevery=5, mew=1.0)
     plt.plot(k_vec, phase[0,:], '-o', color='b', linewidth=1.5, label='SDC('+str(K)+')', markevery=5, markersize=fs/2)
     plt.xlabel('Wave number', fontsize=fs, labelpad=0.25)
     plt.ylabel('Phase speed', fontsize=fs, labelpad=0.5)
@@ -138,14 +149,14 @@ def findomega(stab_fh):
     plt.legend(loc='lower left', fontsize=fs, prop={'size':fs-2})
     plt.xticks([0, 1, 2, 3], fontsize=fs)
     #plt.show()
-    filename = 'sdc-fwsw-disprel-phase-K'+str(K)+'-M'+str(swparams['num_nodes'])+'.pdf'
+    filename = 'phase-K'+str(K)+'-M'+str(swparams['num_nodes'])+'.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
     call(["pdfcrop", filename, filename])
 
     fig  = plt.figure()
     plt.plot(k_vec, 1.0+np.zeros(np.size(k_vec)), '--', color='k', linewidth=1.5, label='Exact')
     plt.plot(k_vec, amp_factor[1,:], '-', color='g', linewidth=1.5, label='DIRK('+str(dirkts.order)+')')
-    plt.plot(k_vec, amp_factor[2,:], '-+', color='r', linewidth=1.5, label='RK-IMEX('+str(rkimex.order)+')', markevery=5, mew=1.0)
+    plt.plot(k_vec, amp_factor[2,:], '-+', color='r', linewidth=1.5, label='IMEX('+str(rkimex.order)+')', markevery=5, mew=1.0)
     plt.plot(k_vec, amp_factor[0,:], '-o', color='b', linewidth=1.5, label='SDC('+str(K)+')', markevery=5, markersize=fs/2)
     plt.xlabel('Wave number', fontsize=fs, labelpad=0.25)
     plt.ylabel('Amplification factor', fontsize=fs, labelpad=0.5)
@@ -156,7 +167,7 @@ def findomega(stab_fh):
     plt.gca().set_ylim([0.0, 1.1])
     plt.xticks([0, 1, 2, 3], fontsize=fs)
     #plt.show()
-    filename = 'sdc-fwsw-disprel-ampfac-K'+str(K)+'-M'+str(swparams['num_nodes'])+'.pdf'
+    filename = 'ampfactor-K'+str(K)+'-M'+str(swparams['num_nodes'])+'.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
     call(["pdfcrop", filename, filename])
 

examples/acoustic_1d_imex/plotconvdata.py

@@ -15,7 +15,7 @@
   if not line: break
   items  = str.split(line, "    ", 3)
   order  = np.append(order,  int(items[0]))
-  nsteps = np.append(nsteps, int(items[1]))
+  nsteps = np.append(nsteps, int(float(items[1])))
   error  = np.append(error,  float(items[2]))
 
 assert np.size(order)==np.size(nsteps), 'Found different number of entries in order and nsteps'
@@ -51,11 +51,11 @@
 plt.xlim([0.9*np.min(nsteps_plot), 1.1*np.max(nsteps_plot)])
 plt.ylim([1e-5, 1e0])
 plt.yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0],fontsize=fs)
-plt.xticks([30, 40, 60, 80], fontsize=fs)
+plt.xticks([20, 30, 40, 60, 80, 100], fontsize=fs)
 plt.gca().get_xaxis().get_major_formatter().labelOnlyBase = False
 plt.gca().get_xaxis().set_major_formatter(ScalarFormatter())
-plt.show()
-filename = 'sdc-fwsw-convergence.pdf'
+#plt.show()
+filename = 'convergence.pdf'
 fig.savefig(filename,bbox_inches='tight')
 call(["pdfcrop", filename, filename])
 

examples/acoustic_1d_imex/runconvergence.py

@@ -28,10 +28,8 @@
 
     sparams = {}
 
-
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['nvars']     = [(2,300)]
     pparams['cadv']      = 0.1
     pparams['cs']        = 1.00
     pparams['order_adv'] = 5
@@ -43,41 +41,57 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLegendre
+    
+    ### SET TYPE OF QUADRATURE NODES ###
+    #description['collocation_class'] = collclass.CollGaussLobatto
+    #description['collocation_class'] = collclass.CollGaussLegendre
+    description['collocation_class'] = collclass.CollGaussRadau_Right
+    
     description['sweeper_class']     = imex_1st_order
     description['do_coll_update']    = True
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
 
-    Nsteps = [30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]
+    nsteps = np.zeros((3,9))
+    nsteps[0,:] = [20, 30, 40, 50, 60, 70, 80, 90, 100]
+    nsteps[1,:] = nsteps[0,:]
+    nsteps[2,:] = nsteps[0,:]
 
-    for order in [2, 3, 4]:
+    for order in [3, 4, 5]:
 
-      error  = np.zeros(np.size(Nsteps))
+      error  = np.zeros(np.shape(nsteps)[1])
 
       # setup parameters "in time"
       t0   = 0
       Tend = 1.0
 
-      if order==2:
+      if order==3:
         file = open('conv-data.txt', 'w')
       else:
         file = open('conv-data.txt', 'a')
 
-      if order==2:
-        description['num_nodes'] = 2
-      elif order==3:
+      ### SET NUMBER OF NODES DEPENDING ON REQUESTED ORDER ###
+      if order==3:
         description['num_nodes'] = 3
       elif order==4:
         description['num_nodes'] = 3
-
+      elif order==5:
+        description['num_nodes'] = 3
+      
       sparams['maxiter'] = order
 
-       # quickly generate block of steps
-      MS = mp.generate_steps(num_procs,sparams,description)
-      for ii in range(0,np.size(Nsteps)):
+      for ii in range(0,np.shape(nsteps)[1]):
 
-        dt = Tend/float(Nsteps[ii])
+        ns = nsteps[order-3,ii]
+        if ((order==3) or (order==4)):
+          pparams['nvars']     = [(2,5*ns)]
+        elif order==5:
+          pparams['nvars'] = [(2,5*ns)]
+          
+        # quickly generate block of steps
+        MS = mp.generate_steps(num_procs,sparams,description)
+        
+        dt = Tend/float(ns)
 
         # get initial values on finest level
         P = MS[0].levels[0].prob
@@ -94,11 +108,11 @@
         uex = P.u_exact(Tend)
 
         error[ii] = np.linalg.norm(uex.values-uend.values,np.inf)/np.linalg.norm(uex.values,np.inf)
-        file.write(str(order)+"    "+str(Nsteps[ii])+"    "+str(error[ii])+"\n")
+        file.write(str(order)+"    "+str(ns)+"    "+str(error[ii])+"\n")
 
       file.close()
 
-    if np.size(Nsteps)==1:
+    if np.shape(nsteps)[1]==1:
       fig = plt.figure(figsize=(8,8))
 
       plt.plot(P.mesh, uex.values[0,:],  '+', color='b', label='u (exact)')
@@ -110,5 +124,5 @@
       plt.ylim([-1.0, 1.0])
       plt.show()
     else:
-      for ii in range(0,np.size(Nsteps)):
-        print('error for nsteps= %s: %s' % (Nsteps[ii], error[ii]))
+      for ii in range(0,np.shape(nsteps)[1]):
+        print('error for nsteps= %s: %s' % (nsteps[order-3,ii], error[ii]))

examples/acoustic_1d_imex/runitererror.py

@@ -29,12 +29,14 @@
 
     # This comes as read-in for the level class
     lparams = {}
-    lparams['restol'] = 1E-14
+    lparams['restol'] = 1E-10
 
     sparams = {}
 
     cs_v = [0.5, 1.0, 1.5, 5.0]
     sparams['maxiter'] = 15
+    
+    ### SET NUMBER OF NODES ###
     nodes_v = [3]
 
     residual = np.zeros((np.size(cs_v), np.size(nodes_v), sparams['maxiter'])) 
@@ -62,7 +64,12 @@
       description['problem_params']    = pparams
       description['dtype_u']           = mesh
       description['dtype_f']           = rhs_imex_mesh
-      description['collocation_class'] = collclass.CollGaussLegendre
+      
+      ### SELECT TYPE OF QUADRATURE NODES ###
+      #description['collocation_class'] = collclass.CollGaussLobatto
+      #description['collocation_class'] = collclass.CollGaussLegendre
+      description['collocation_class'] = collclass.CollGaussRadau_Right
+      
       description['sweeper_class']     = imex_1st_order
       description['level_params']      = lparams
       description['hook_class']        = plot_solution
@@ -97,7 +104,7 @@
 
         # Compute convergence rates
         for iter in range(0,sparams['maxiter']-1):
-            if residual[cs_ind, nodes_ind, iter]<lparams['restol']:
+            if residual[cs_ind, nodes_ind, iter]< lparams['restol']:
               lastiter[cs_ind,nodes_ind] = iter
             else:
               convrate[cs_ind, nodes_ind, iter] = residual[cs_ind, nodes_ind, iter+1]/residual[cs_ind, nodes_ind, iter]
@@ -127,7 +134,7 @@
     plt.yticks(fontsize=fs)
     plt.xticks(fontsize=fs)
     plt.show()
-    filename = 'sdc-fwsw-iteration.pdf'
+    filename = 'iteration.pdf'
     fig.savefig(filename,bbox_inches='tight')
     call(["pdfcrop", filename, filename])
 

examples/acoustic_1d_imex/runmultiscale.py

@@ -32,13 +32,14 @@
     lparams['restol'] = 1E-10
 
     sparams = {}
-    sparams['maxiter'] = 5
+    
+    ### SET NUMBER OF ITERATIONS ###
+    sparams['maxiter'] = 2
 
     # setup parameters "in time"
     t0   = 0.0
     Tend = 3.0
     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
@@ -59,9 +60,15 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLegendre
-    # Number of nodes
-    description['num_nodes']         = 3
+    
+    ### SET TYPE OF QUADRATURE NODES ###
+    #description['collocation_class'] = collclass.CollGaussLobatto
+    #description['collocation_class'] = collclass.CollGaussLegendre
+    description['collocation_class'] = collclass.CollGaussRadau_Right
+    
+    ### SET NUMBER OF NODES ###
+    description['num_nodes']         = 2
+    
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
@@ -152,7 +159,7 @@
     plt.legend(loc='upper left', fontsize=fs, prop={'size':fs}, handlelength=3)
     fig.gca().grid()
     #plt.show()
-    filename = 'sdc-fwsw-multiscale-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf'
+    filename = 'multiscale-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
     call(["pdfcrop", filename, filename])
 
@@ -168,6 +175,6 @@
     fig.gca().tick_params(axis='both', labelsize=fs)
     fig.gca().grid()
     #plt.show()
-    filename = 'sdc-fwsw-multiscale-initial.pdf'
+    filename = 'multiscale-initial.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
     call(["pdfcrop", filename, filename])