|
17 | 17 | from mpl_toolkits.mplot3d import Axes3D |
18 | 18 | from matplotlib import cm |
19 | 19 | from matplotlib.ticker import LinearLocator, FormatStrFormatter |
| 20 | +from pylab import rcParams |
| 21 | + |
20 | 22 | from unflatten import unflatten |
21 | 23 |
|
22 | 24 | from standard_integrators import dirk |
|
30 | 32 |
|
31 | 33 | # This comes as read-in for the level class |
32 | 34 | lparams = {} |
33 | | - lparams['restol'] = 1E-8 |
| 35 | + lparams['restol'] = 1E-15 |
34 | 36 |
|
35 | 37 | swparams = {} |
36 | 38 | swparams['collocation_class'] = collclass.CollGaussLobatto |
37 | 39 | swparams['num_nodes'] = 3 |
38 | | - swparams['do_LU'] = True |
| 40 | + swparams['do_LU'] = False |
39 | 41 |
|
40 | 42 | sparams = {} |
41 | | - sparams['maxiter'] = 4 |
| 43 | + sparams['maxiter'] = 3 |
42 | 44 |
|
43 | | - dirk_order = 4 |
| 45 | + dirk_order = 2 |
44 | 46 |
|
45 | 47 | # setup parameters "in time" |
46 | 48 | t0 = 0 |
|
52 | 54 |
|
53 | 55 | # This comes as read-in for the problem class |
54 | 56 | pparams = {} |
55 | | - #pparams['nvars'] = [(4,450,30)] |
56 | | - pparams['nvars'] = [(4,150,10)] |
| 57 | + pparams['nvars'] = [(4,300,20)] |
| 58 | + #pparams['nvars'] = [(4,150,10)] |
57 | 59 | pparams['u_adv'] = 0.02 |
58 | 60 | pparams['c_s'] = 0.3 |
59 | 61 | pparams['Nfreq'] = 0.01 |
60 | 62 | pparams['x_bounds'] = [(-150.0, 150.0)] |
61 | 63 | pparams['z_bounds'] = [( 0.0, 10.0)] |
62 | 64 | pparams['order'] = [4] # [fine_level, coarse_level] |
63 | 65 | pparams['order_upw'] = [5] |
64 | | - pparams['gmres_maxiter'] = [50] |
| 66 | + pparams['gmres_maxiter'] = [500] |
65 | 67 | pparams['gmres_restart'] = [10] |
66 | 68 | pparams['gmres_tol'] = [1e-6] |
67 | 69 |
|
|
87 | 89 | P = MS[0].levels[0].prob |
88 | 90 | uinit = P.u_exact(t0) |
89 | 91 |
|
90 | | - dirk = dirk(P, dirk_order) |
91 | | - u0 = uinit.values.flatten() |
92 | | - |
93 | | - for i in range(0,Nsteps): |
94 | | - u0 = dirk.timestep(u0, dt) |
95 | | - |
96 | 92 | cfl_advection = pparams['u_adv']*dt/P.h[0] |
97 | 93 | cfl_acoustic_hor = pparams['c_s']*dt/P.h[0] |
98 | 94 | cfl_acoustic_ver = pparams['c_s']*dt/P.h[1] |
99 | 95 | print ("CFL number of advection: %4.2f" % cfl_advection) |
100 | 96 | print ("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor) |
101 | 97 | print ("CFL number of acoustics (vertical): %4.2f" % cfl_acoustic_ver) |
102 | 98 |
|
| 99 | + dirk = dirk(P, dirk_order) |
| 100 | + u0 = uinit.values.flatten() |
| 101 | + |
| 102 | + for i in range(0,Nsteps): |
| 103 | + u0 = dirk.timestep(u0, dt) |
| 104 | + |
103 | 105 | # call main function to get things done... |
104 | 106 | uend,stats = mp.run_pfasst(MS,u0=uinit,t0=t0,dt=dt,Tend=Tend) |
105 | 107 |
|
106 | 108 | u0 = unflatten(u0, 4, P.N[0], P.N[1]) |
| 109 | + |
| 110 | + fs = 8 |
| 111 | + rcParams['figure.figsize'] = 5.0, 2.5 |
107 | 112 | fig = plt.figure() |
108 | 113 |
|
109 | | - plt.plot(P.xx[:,5], uend.values[2,:,5], color='r', marker='+', markevery=3) |
110 | | - plt.plot(P.xx[:,5], u0[2,:,5], color='b', markevery=5) |
111 | | - plt.show() |
| 114 | + plt.plot(P.xx[:,5], uend.values[2,:,5], '-', color='b', label='SDC') |
| 115 | + plt.plot(P.xx[:,5], u0[2,:,5], '+', color='g', markevery=5, markersize=fs-2, label='DIRK') |
| 116 | + plt.legend(loc='lower left', fontsize=fs, prop={'size':fs}) |
| 117 | + plt.yticks(fontsize=fs) |
| 118 | + plt.xticks(fontsize=fs) |
| 119 | + plt.xlabel('x', fontsize=fs, labelpad=0) |
| 120 | + plt.ylabel('Bouyancy', fontsize=fs, labelpad=1) |
| 121 | + #plt.show() |
| 122 | + plt.savefig('boussinesq.pdf', bbox_inches='tight') |
112 | 123 |
|
113 | 124 | print " #### Logging report for DIRK #### " |
114 | 125 | print "Number of calls to implicit solver: %5i" % dirk.logger.solver_calls |
|
0 commit comments