some fixes in convergence plots

Author: Daniel Ruprecht

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'
@@ -49,13 +49,13 @@
 plt.xlabel('Number of time steps', fontsize=fs)
 plt.ylabel('Relative error', fontsize=fs, labelpad=2)
 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.ylim([1e-7, 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.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,56 +28,65 @@
 
     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
-    pparams['waveno']    = 5
+    pparams['waveno']    = 0
 
     # Fill description dictionary for easy hierarchy creation
     description = {}
     description['problem_class']     = acoustic_1d_imex
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLegendre
+    #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,7))
+    nsteps[0,:] = [10, 15, 20, 25, 30, 35, 40]
+    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:
+      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,4*ns)]
+        elif order==5:
+          pparams['navrs'] = [(2,3*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 +103,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 +119,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/ruprecht_speck_2016.README

@@ -1,7 +1,7 @@
 Figure 1: fwsw/plot_stifflimit_specrad.py
 Figure 2: fwsw/plot_stability.py
 Figure 3: acoustic_1d_imex/plot_dispersion.py 
-Figure 4:
+Figure 4: acoustic_1d_imex/runconvergence.py and then plotconvergence.py
 Figure 5:
 Figure 6: boussinesq_2d_imex/rungmrescounter.py and then plotgmrescounter.py for sparams['maxiter']=4 and dirk_order=4