changed run routines to legendre nodes; automatically cropping generated pdfs

Author: Daniel Ruprecht

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/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