diff --git a/examples/acoustic_1d_imex/plotconvdata.py b/examples/acoustic_1d_imex/plotconvdata.py index d5665af939..67757f98b4 100644 --- a/examples/acoustic_1d_imex/plotconvdata.py +++ b/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]) + diff --git a/examples/acoustic_1d_imex/runconvergence.py b/examples/acoustic_1d_imex/runconvergence.py index ebe3697326..064b70df0c 100644 --- a/examples/acoustic_1d_imex/runconvergence.py +++ b/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] diff --git a/examples/acoustic_1d_imex/runitererror.py b/examples/acoustic_1d_imex/runitererror.py index 5c72bcbe9e..b06811a22e 100644 --- a/examples/acoustic_1d_imex/runitererror.py +++ b/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]) diff --git a/examples/acoustic_1d_imex/runmultiscale.py b/examples/acoustic_1d_imex/runmultiscale.py index 1ed81ec921..e1d0fc5ddb 100644 --- a/examples/acoustic_1d_imex/runmultiscale.py +++ b/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]) + diff --git a/examples/boussinesq_2d_imex/plotgmrescounter.py b/examples/boussinesq_2d_imex/plotgmrescounter.py index c436929acf..4dcac27f72 100644 --- a/examples/boussinesq_2d_imex/plotgmrescounter.py +++ b/examples/boussinesq_2d_imex/plotgmrescounter.py @@ -4,23 +4,26 @@ 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) @@ -28,5 +31,7 @@ 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]) diff --git a/examples/boussinesq_2d_imex/rungmrescounter.py b/examples/boussinesq_2d_imex/rungmrescounter.py index 17cfdf2438..c74138ff1f 100644 --- a/examples/boussinesq_2d_imex/rungmrescounter.py +++ b/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 diff --git a/examples/fwsw/plot_stifflimit_specrad.py b/examples/fwsw/plot_stifflimit_specrad.py index 17347b644d..67605071e3 100644 --- a/examples/fwsw/plot_stifflimit_specrad.py +++ b/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]) +