Coverage for pySDC/projects/FastWaveSlowWave/runmultiscale_acoustic.py: 99%

1import matplotlib 

2 

3matplotlib.use('Agg') 

4 

5import numpy as np 

6from matplotlib import pyplot as plt 

7from pylab import rcParams 

8 

9from pySDC.projects.FastWaveSlowWave.HookClass_acoustic import dump_energy 

10 

11from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

12from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

13from pySDC.implementations.problem_classes.acoustic_helpers.standard_integrators import bdf2, dirk, trapezoidal, rk_imex 

14 

15from pySDC.projects.FastWaveSlowWave.AcousticAdvection_1D_FD_imex_multiscale import acoustic_1d_imex_multiscale 

16 

17 

18def compute_and_plot_solutions(): 

19 """ 

20 Routine to compute and plot the solutions of SDC(2), IMEX, BDF-2 and RK for a multiscale problem 

21 """ 

22 

23 num_procs = 1 

24 

25 t0 = 0.0 

26 Tend = 3.0 

27 nsteps = 154 # 154 is value in Vater et al. 

28 dt = Tend / float(nsteps) 

29 

30 # initialize level parameters 

31 level_params = dict() 

32 level_params['restol'] = 1e-10 

33 level_params['dt'] = dt 

34 

35 # This comes as read-in for the step class 

36 step_params = dict() 

37 step_params['maxiter'] = 2 

38 

39 # This comes as read-in for the problem class 

40 problem_params = dict() 

41 problem_params['cadv'] = 0.05 

42 problem_params['cs'] = 1.0 

43 problem_params['nvars'] = [(2, 512)] 

44 problem_params['order_adv'] = 5 

45 problem_params['waveno'] = 5 

46 

47 # This comes as read-in for the sweeper class 

48 sweeper_params = dict() 

49 sweeper_params['quad_type'] = 'RADAU-RIGHT' 

50 sweeper_params['num_nodes'] = 2 

51 

52 # initialize controller parameters 

53 controller_params = dict() 

54 controller_params['logger_level'] = 30 

55 controller_params['hook_class'] = dump_energy 

56 

57 # Fill description dictionary for easy hierarchy creation 

58 description = dict() 

59 description['problem_class'] = acoustic_1d_imex_multiscale 

60 description['problem_params'] = problem_params 

61 description['sweeper_class'] = imex_1st_order 

62 description['sweeper_params'] = sweeper_params 

63 description['step_params'] = step_params 

64 description['level_params'] = level_params 

65 

66 # instantiate the controller 

67 controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) 

68 

69 # get initial values on finest level 

70 P = controller.MS[0].levels[0].prob 

71 uinit = P.u_exact(t0) 

72 

73 # call main function to get things done... 

74 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) 

75 

76 # instantiate standard integrators to be run for comparison 

77 trap = trapezoidal((P.A + P.Dx).astype('complex'), 0.5) 

78 bdf2_m = bdf2(P.A + P.Dx) 

79 dirk_m = dirk((P.A + P.Dx).astype('complex'), step_params['maxiter']) 

80 rkimex = rk_imex(P.A.astype('complex'), P.Dx.astype('complex'), step_params['maxiter']) 

81 

82 y0_tp = np.concatenate((uinit[0, :], uinit[1, :])) 

83 y0_bdf = y0_tp 

84 y0_dirk = y0_tp.astype('complex') 

85 y0_imex = y0_tp.astype('complex') 

86 

87 # Perform time steps with standard integrators 

88 for i in range(0, nsteps): 

89 # trapezoidal rule step 

90 ynew_tp = trap.timestep(y0_tp, dt) 

91 

92 # BDF-2 scheme 

93 if i == 0: 

94 ynew_bdf = bdf2_m.firsttimestep(y0_bdf, dt) 

95 ym1_bdf = y0_bdf 

96 else: 

97 ynew_bdf = bdf2_m.timestep(y0_bdf, ym1_bdf, dt) 

98 

99 # DIRK scheme 

100 ynew_dirk = dirk_m.timestep(y0_dirk, dt) 

101 

102 # IMEX scheme 

103 ynew_imex = rkimex.timestep(y0_imex, dt) 

104 

105 y0_tp = ynew_tp 

106 ym1_bdf = y0_bdf 

107 y0_bdf = ynew_bdf 

108 y0_dirk = ynew_dirk 

109 y0_imex = ynew_imex 

110 

111 # Finished running standard integrators 

112 unew_tp, pnew_tp = np.split(ynew_tp, 2) 

113 unew_bdf, pnew_bdf = np.split(ynew_bdf, 2) 

114 unew_dirk, pnew_dirk = np.split(ynew_dirk, 2) 

115 unew_imex, pnew_imex = np.split(ynew_imex, 2) 

116 

117 fs = 8 

118 

119 rcParams['figure.figsize'] = 2.5, 2.5 

120 # rcParams['pgf.rcfonts'] = False 

121 fig = plt.figure() 

122 

123 sigma_0 = 0.1 

124 # k = 7.0 * 2.0 * np.pi 

125 x_0 = 0.75 

126 # x_1 = 0.25 

127 

128 assert np.isclose(np.linalg.norm(uend[1, :], np.inf), 8.489e-01, 1e-03) 

129 assert np.isclose(np.linalg.norm(pnew_dirk, np.inf), 1.003e00, 1e-03) 

130 assert np.isclose(np.linalg.norm(pnew_imex, np.inf), 2.762e21, 1e-03) 

131 

132 print('Maximum pressure in SDC: %5.3e' % np.linalg.norm(uend[1, :], np.inf)) 

133 print('Maximum pressure in DIRK: %5.3e' % np.linalg.norm(pnew_dirk, np.inf)) 

134 print('Maximum pressure in RK-IMEX: %5.3e' % np.linalg.norm(pnew_imex, np.inf)) 

135 

136 if dirk_m.order == 2: 

137 plt.plot(P.mesh, pnew_bdf, 'd-', color='c', label='BDF-2', markevery=(50, 75)) 

138 p_slow = np.exp(-np.square(np.mod(P.mesh - problem_params['cadv'] * Tend, 1.0) - x_0) / (sigma_0 * sigma_0)) 

139 plt.plot(P.mesh, p_slow, '--', color='k', markersize=fs - 2, label='Slow mode', dashes=(10, 2)) 

140 if np.linalg.norm(pnew_imex, np.inf) <= 2: 

141 plt.plot( 

142 P.mesh, pnew_imex, '+-', color='r', label='IMEX(' + str(rkimex.order) + ')', markevery=(1, 75), mew=1.0 

143 ) 

144 plt.plot(P.mesh, uend[1, :], 'o-', color='b', label='SDC(' + str(step_params['maxiter']) + ')', markevery=(25, 75)) 

145 plt.plot(P.mesh, np.real(pnew_dirk), '-', color='g', label='DIRK(' + str(dirk_m.order) + ')') 

146 

147 plt.xlabel('x', fontsize=fs, labelpad=0) 

148 plt.ylabel('Pressure', fontsize=fs, labelpad=0) 

149 fig.gca().set_xlim([0, 1.0]) 

150 fig.gca().set_ylim([-0.5, 1.1]) 

151 fig.gca().tick_params(axis='both', labelsize=fs) 

152 plt.legend(loc='upper left', fontsize=fs, prop={'size': fs}, handlelength=3) 

153 fig.gca().grid() 

154 filename = 'data/multiscale-K' + str(step_params['maxiter']) + '-M' + str(sweeper_params['num_nodes']) + '.png' 

155 plt.gcf().savefig(filename, bbox_inches='tight') 

156 

157 

158if __name__ == "__main__": 

159 compute_and_plot_solutions()