fixed parameter for convergence plot

Author: Daniel Ruprecht

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/plotconvdata.py

@@ -49,9 +49,9 @@
 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-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.ylim([1e-5, 1e0])
+plt.yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0],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()

examples/acoustic_1d_imex/runconvergence.py

@@ -33,7 +33,7 @@
     pparams['cadv']      = 0.1
     pparams['cs']        = 1.00
     pparams['order_adv'] = 5
-    pparams['waveno']    = 0
+    pparams['waveno']    = 5
 
     # Fill description dictionary for easy hierarchy creation
     description = {}
@@ -49,8 +49,8 @@
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
 
-    nsteps = np.zeros((3,7))
-    nsteps[0,:] = [10, 15, 20, 25, 30, 35, 40]
+    nsteps = np.zeros((3,9))
+    nsteps[0,:] = [20, 30, 40, 50, 60, 70, 80, 90, 100]
     nsteps[1,:] = nsteps[0,:]
     nsteps[2,:] = nsteps[0,:]
 
@@ -73,15 +73,16 @@
         description['num_nodes'] = 3
       elif order==5:
         description['num_nodes'] = 3
+      
       sparams['maxiter'] = order
 
       for ii in range(0,np.shape(nsteps)[1]):
 
         ns = nsteps[order-3,ii]
         if ((order==3) or (order==4)):
-          pparams['nvars']     = [(2,4*ns)]
+          pparams['nvars']     = [(2,5*ns)]
         elif order==5:
-          pparams['navrs'] = [(2,3*ns)]
+          pparams['nvars'] = [(2,5*ns)]
 
         # quickly generate block of steps
         MS = mp.generate_steps(num_procs,sparams,description)

examples/acoustic_1d_imex/runitererror.py

@@ -29,7 +29,7 @@
 
     # This comes as read-in for the level class
     lparams = {}
-    lparams['restol'] = 1E-14
+    lparams['restol'] = 1E-10
 
     sparams = {}
 
@@ -62,7 +62,9 @@
       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['level_params']      = lparams
       description['hook_class']        = plot_solution
@@ -97,7 +99,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 +129,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

@@ -38,7 +38,7 @@
     t0   = 0.0
     Tend = 3.0
     nsteps = 154 # 154 is value in Vater et al.
-    nsteps = 2*154    
+    #nsteps = 2*154
     dt = Tend/float(nsteps)
 
     # This comes as read-in for the problem class

examples/ruprecht_speck_2016.README

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