Fix/change plots to legendre

examples/acoustic_1d_imex/plotconvdata.py

@@ -2,6 +2,7 @@
 from matplotlib import pyplot as plt
 from pylab import rcParams
 from matplotlib.ticker import ScalarFormatter
+from subprocess import call
 
 fs     = 8
 order  = np.array([])
@@ -54,5 +55,8 @@
 plt.gca().get_xaxis().get_major_formatter().labelOnlyBase = False
 plt.gca().get_xaxis().set_major_formatter(ScalarFormatter())
 plt.show()
-fig.savefig('sdc_fwsw_convergence.pdf',bbox_inches='tight')
+filename = 'sdc_fwsw_convergence.pdf'
+fig.savefig(filename,bbox_inches='tight')
+call(["pdfcrop", filename, filename])
+
 

examples/acoustic_1d_imex/runconvergence.py

@@ -31,11 +31,11 @@
 
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['nvars']     = [(2,250)]
+    pparams['nvars']     = [(2,400)]
     pparams['cadv']      = 0.05
     pparams['cs']        = 1.00
     pparams['order_adv'] = 5
-    pparams['waveno']    = 1
+    pparams['waveno']    = 2
 
     # This comes as read-in for the transfer operations
     tparams = {}
@@ -47,12 +47,11 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLobatto
+    description['collocation_class'] = collclass.CollGaussLegendre
     description['sweeper_class']     = imex_1st_order
+    description['do_coll_update']    = True
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
-    #description['transfer_class'] = mesh_to_mesh_1d
-    #description['transfer_params'] = tparams
 
     Nsteps = [15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]
 

examples/acoustic_1d_imex/runitererror.py

@@ -16,7 +16,7 @@
 
 from matplotlib import pyplot as plt
 from pylab import rcParams
-
+from subprocess import call
 
 if __name__ == "__main__":
 
@@ -62,7 +62,7 @@
       description['problem_params']    = pparams
       description['dtype_u']           = mesh
       description['dtype_f']           = rhs_imex_mesh
-      description['collocation_class'] = collclass.CollGaussLobatto
+      description['collocation_class'] = collclass.CollGaussLegendre
       description['sweeper_class']     = imex_1st_order
       description['level_params']      = lparams
       description['hook_class']        = plot_solution
@@ -127,6 +127,8 @@
     plt.yticks(fontsize=fs)
     plt.xticks(fontsize=fs)
     plt.show()
-    fig.savefig('sdc_fwsw_iteration.pdf',bbox_inches='tight')
+    filename = 'sdc_fwsw_iteration.pdf'
+    fig.savefig(filename,bbox_inches='tight')
+    call(["pdfcrop", filename, filename])
 
 

examples/acoustic_1d_imex/runmultiscale.py

@@ -16,6 +16,7 @@
 
 from matplotlib import pyplot as plt
 from pylab import rcParams
+from subprocess import call
 
 fs = 8
 
@@ -56,11 +57,9 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussRadau_Right
-    if sparams['maxiter']==2:
-      description['num_nodes']         = 2
-    else:
-      description['num_nodes']         = 3
+    description['collocation_class'] = collclass.CollGaussLegendre
+    # Number of nodes
+    description['num_nodes']         = 3
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
@@ -137,4 +136,7 @@
     plt.legend(loc='upper left', fontsize=fs, prop={'size':fs})
     fig.gca().grid()
     #plt.show()
-    plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')
+    filename = 'sdc-fwsw-multiscale-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf'
+    plt.gcf().savefig(filename, bbox_inches='tight')
+    call(["pdfcrop", filename, filename])
+

examples/boussinesq_2d_imex/plotgmrescounter.py

@@ -4,29 +4,34 @@
 from matplotlib import cm
 from matplotlib.ticker import LinearLocator, FormatStrFormatter
 from pylab import rcParams
+from subprocess import call
 
 from unflatten import unflatten
 
 if __name__ == "__main__":
   xx = np.load('xaxis.npy')
   uend = np.load('sdc.npy')
-  udirk2 = np.load('dirk2.npy')
-  udirk4 = np.load('dirk4.npy')
-  utrap   = np.load('trap.npy')
+  udirk = np.load('dirk.npy')
+  uref  = np.load('uref.npy')
+
+  print "Estimated discretisation error of DIRK: %5.3e" % ( np.linalg.norm(udirk.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )
+  print "Estimated discretisation error of SDC:  %5.3e" % ( np.linalg.norm(uend.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )
 
   fs = 8
   rcParams['figure.figsize'] = 5.0, 2.5
   fig = plt.figure()
 
+  plt.plot(xx[:,5], udirk[2,:,5], '--', color='r', markersize=fs-2, label='DIRK', dashes=(3,3))
   plt.plot(xx[:,5], uend[2,:,5], '-', color='b', label='SDC')
-  plt.plot(xx[:,5], udirk4[2,:,5], '--', color='g', markersize=fs-2, label='DIRK(4)', dashes=(3,3))
-  plt.plot(xx[:,5], udirk2[2,:,5], '--', color='r', markersize=fs-2, label='DIRK(2)', dashes=(3,3))
+  #plt.plot(xx[:,5], udirk2[2,:,5], '--', color='r', markersize=fs-2, label='DIRK(2)', dashes=(3,3))
   #plt.plot(xx[:,5], utrap[2,:,5], '--', color='k', markersize=fs-2, label='Trap', dashes=(3,3))
   plt.legend(loc='lower left', fontsize=fs, prop={'size':fs})
   plt.yticks(fontsize=fs)
   plt.xticks(fontsize=fs)
   plt.xlabel('x [km]', fontsize=fs, labelpad=0)
   plt.ylabel('Bouyancy', fontsize=fs, labelpad=1)
   #plt.show()
-  plt.savefig('boussinesq.pdf', bbox_inches='tight')
+  filename = 'sdc-fwsw-boussinesq.pdf'
+  plt.savefig(filename, bbox_inches='tight')
+  call(["pdfcrop", filename, filename])
 

examples/boussinesq_2d_imex/rungmrescounter.py

@@ -35,7 +35,7 @@
     lparams['restol'] = 1E-15
 
     swparams = {}
-    swparams['collocation_class'] = collclass.CollGaussLobatto
+    swparams['collocation_class'] = collclass.CollGaussLegendre
     swparams['num_nodes'] = 3
     swparams['do_LU'] = False
 
@@ -96,54 +96,32 @@
     print ("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor)
     print ("CFL number of acoustics (vertical):   %4.2f" % cfl_acoustic_ver)
 
-    dirk4 = dirk(P, 4)
-    dirk2 = dirk(P, 2)
-    trap  = trapezoidal(P)
-    bdf  = bdf2(P)
+    dirkp = dirk(P, dirk_order)
     u0 = uinit.values.flatten()
-    udirk4 = u0
-    udirk2 = u0
-    ubdf  = u0
-    utrap  = u0
+    udirk = u0
     for i in range(0,Nsteps):
-      udirk4 = dirk4.timestep(udirk4, dt)  
-      udirk2 = dirk2.timestep(udirk2, dt)
-      utrap  = trap.timestep(utrap, dt)
-      #if i==0:
-      #  ubdf_new = bdf.firsttimestep(ubdf, dt)
-      #  ubdf_m1  = ubdf
-      #else:        
-      #  ubdf_new  = bdf.timestep(ubdf, ubdf_m1, dt)
-      #ubdf_m1 = ubdf
-      #ubdf    = ubdf_new    
+      udirk = dirkp.timestep(udirk, dt)  
+
+    dirkref = dirk(P, 4)
+    uref    = u0
+    dt_ref  = dt/10.0
+    for i in range(0,10*Nsteps):
+      uref = dirkref.timestep(uref, dt_ref)
 
     # call main function to get things done...
     uend,stats = mp.run_pfasst(MS,u0=uinit,t0=t0,dt=dt,Tend=Tend)
-    udirk4 = unflatten(udirk4, 4, P.N[0], P.N[1])
-    udirk2 = unflatten(udirk2, 4, P.N[0], P.N[1])
-    print "Norm of final solution by trapezoidal rule: %5.3f" % np.linalg.norm( utrap, np.inf )    
-    utrap  = unflatten(utrap, 4, P.N[0], P.N[1])
-
+    udirk = unflatten(udirk, 4, P.N[0], P.N[1])
+    uref  = unflatten(uref,  4, P.N[0], P.N[1])
+
     np.save('xaxis', P.xx)
     np.save('sdc', uend.values)
-    np.save('dirk2', udirk2)
-    np.save('dirk4', udirk4)
-    np.save('trap', utrap)
+    np.save('dirk', udirk)
+    np.save('uref', uref)
 
-    print " #### Logging report for DIRK-4 #### "
-    print "Number of calls to implicit solver: %5i" % dirk4.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % dirk4.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(dirk4.logger.iterations)/float(dirk4.logger.solver_calls))
-
-    print " #### Logging report for DIRK-2 #### "
-    print "Number of calls to implicit solver: %5i" % dirk2.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % dirk2.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(dirk2.logger.iterations)/float(dirk2.logger.solver_calls))
-
-    #print " #### Logging report for BDF2 #### "
-    #print "Number of calls to implicit solver: %5i" % bdf.logger.solver_calls
-    #print "Total number of GMRES iterations: %5i" % bdf.logger.iterations
-    #print "Average number of iterations per call: %6.3f" % (float(bdf.logger.iterations)/float(bdf.logger.solver_calls))
+    print " #### Logging report for DIRK-%1i #### " % dirk_order
+    print "Number of calls to implicit solver: %5i" % dirkp.logger.solver_calls
+    print "Total number of GMRES iterations: %5i" % dirkp.logger.iterations
+    print "Average number of iterations per call: %6.3f" % (float(dirkp.logger.iterations)/float(dirkp.logger.solver_calls))
 
     print " #### Logging report for SDC-(%1i,%1i) #### " % (swparams['num_nodes'], sparams['maxiter'])
     print "Number of calls to implicit solver: %5i" % P.logger.solver_calls

examples/fwsw/plot_stifflimit_specrad.py

@@ -10,6 +10,7 @@
 
 from matplotlib import pyplot as plt
 from pylab import rcParams
+from subprocess import call
 
 if __name__ == "__main__":
 
@@ -20,7 +21,7 @@
     pparams['lambda_f'] = np.array([50.0*1j], dtype='complex')
     pparams['u0'] = 1.0
     swparams = {}
-    swparams['collocation_class'] = collclass.CollGaussLobatto
+    swparams['collocation_class'] = collclass.CollGaussLegendre
 
     nodes_v = np.arange(2,10)
     specrad = np.zeros((2,np.size(nodes_v)))
@@ -53,8 +54,10 @@
         # For Lobatto nodes, first column and row are all zeros, since q_1 = q_0; hence remove them
         QI = QI[1:,1:]
         Q  = Q[1:,1:]
-      # Eigenvalue of error propagation matrix in stiff limit: E = I - inv(QI)*Q
-      evals, evecs = np.linalg.eig( np.eye(nnodes-1) - np.linalg.inv(QI).dot(Q) )
+        # Eigenvalue of error propagation matrix in stiff limit: E = I - inv(QI)*Q
+        evals, evecs = np.linalg.eig( np.eye(nnodes-1) - np.linalg.inv(QI).dot(Q) )
+      else:
+        evals, evecs = np.linalg.eig( np.eye(nnodes) - np.linalg.inv(QI).dot(Q) )
       specrad[0,i] = np.linalg.norm( evals, np.inf )
 
   ### Plot result
@@ -71,5 +74,8 @@
     plt.yticks(fontsize=fs)
     plt.xticks(fontsize=fs)
     #plt.show()
-    fig.savefig('sdc_fwsw_stifflimit_specrad.pdf',bbox_inches='tight')
+    filename = 'sdc_fwsw_stifflimit_specrad.pdf'
+    fig.savefig(filename,bbox_inches='tight')
+    call(["pdfcrop", filename, filename])
+