Coverage for pySDC/projects/FastWaveSlowWave/rungmrescounter_boussinesq.py: 100%

1import numpy as np 

2 

3 

4from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

5from pySDC.implementations.problem_classes.Boussinesq_2D_FD_imex import boussinesq_2d_imex 

6from pySDC.implementations.problem_classes.boussinesq_helpers.standard_integrators import SplitExplicit, dirk, rk_imex 

7from pySDC.implementations.problem_classes.boussinesq_helpers.unflatten import unflatten 

8from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

9 

10 

11def main(cwd='', nvars=(4, 300, 30), Tend=3000, Nsteps=100): 

12 """ 

13 Example running/comparing SDC and different standard integrators for the 2D Boussinesq equation 

14 

15 Args: 

16 cwd (string): current working directory 

17 nvars (tuple): number of degrees of freedom, reduced by the test to keep it cheap 

18 Tend (float): end time 

19 Nsteps (int): number of time steps 

20 """ 

21 

22 num_procs = 1 

23 

24 # setup parameters "in time" 

25 t0 = 0 

26 dt = Tend / float(Nsteps) 

27 

28 # initialize level parameters 

29 level_params = dict() 

30 level_params['restol'] = 1e-15 

31 level_params['dt'] = dt 

32 

33 # initialize step parameters 

34 step_params = dict() 

35 step_params['maxiter'] = 4 

36 

37 # initialize sweeper parameters 

38 sweeper_params = dict() 

39 sweeper_params['quad_type'] = 'GAUSS' 

40 sweeper_params['num_nodes'] = 3 

41 

42 # initialize controller parameters 

43 controller_params = dict() 

44 controller_params['logger_level'] = 20 

45 

46 # initialize problem parameters 

47 problem_params = dict() 

48 problem_params['nvars'] = [nvars] 

49 problem_params['u_adv'] = 0.02 

50 problem_params['c_s'] = 0.3 

51 problem_params['Nfreq'] = 0.01 

52 problem_params['x_bounds'] = [(-150.0, 150.0)] 

53 problem_params['z_bounds'] = [(0.0, 10.0)] 

54 problem_params['order'] = [4] 

55 problem_params['order_upw'] = [5] 

56 problem_params['gmres_maxiter'] = [500] 

57 problem_params['gmres_restart'] = [10] 

58 problem_params['gmres_tol_limit'] = [1e-05] 

59 

60 # fill description dictionary for easy step instantiation 

61 description = dict() 

62 description['problem_class'] = boussinesq_2d_imex # pass problem class 

63 description['problem_params'] = problem_params # pass problem parameters 

64 description['sweeper_class'] = imex_1st_order # pass sweeper (see part B) 

65 description['sweeper_params'] = sweeper_params # pass sweeper parameters 

66 description['level_params'] = level_params # pass level parameters 

67 description['step_params'] = step_params # pass step parameters 

68 

69 # ORDER OF DIRK/IMEX EQUAL TO NUMBER OF SDC ITERATIONS AND THUS SDC ORDER 

70 dirk_order = step_params['maxiter'] 

71 

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

73 

74 # get initial values on finest level 

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

76 uinit = P.u_exact(t0) 

77 

78 cfl_advection = P.u_adv * dt / P.h[0] 

79 cfl_acoustic_hor = P.c_s * dt / P.h[0] 

80 cfl_acoustic_ver = P.c_s * dt / P.h[1] 

81 print("Horizontal resolution: %4.2f" % P.h[0]) 

82 print("Vertical resolution: %4.2f" % P.h[1]) 

83 print("CFL number of advection: %4.2f" % cfl_advection) 

84 print("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor) 

85 print("CFL number of acoustics (vertical): %4.2f" % cfl_acoustic_ver) 

86 

87 print("Running SplitExplicit ....") 

88 method_split = 'MIS4_4' 

89 # method_split = 'RK3' 

90 splitp = SplitExplicit(P, method_split, problem_params) 

91 u0 = uinit.flatten() 

92 usplit = np.copy(u0) 

93 print(np.linalg.norm(usplit)) 

94 for _ in range(0, 2 * Nsteps): 

95 usplit = splitp.timestep(usplit, dt / 2) 

96 print(np.linalg.norm(usplit)) 

97 

98 print("Running DIRK ....") 

99 dirkp = dirk(P, dirk_order) 

100 udirk = np.copy(u0) 

101 print(np.linalg.norm(udirk)) 

102 for _ in range(0, Nsteps): 

103 udirk = dirkp.timestep(udirk, dt) 

104 print(np.linalg.norm(udirk)) 

105 

106 print("Running RK-IMEX ....") 

107 rkimex = rk_imex(P, dirk_order) 

108 uimex = np.copy(u0) 

109 dt_imex = dt 

110 for _ in range(0, Nsteps): 

111 uimex = rkimex.timestep(uimex, dt_imex) 

112 print(np.linalg.norm(uimex)) 

113 

114 print("Running SDC...") 

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

116 

117 # For reference solution, increase GMRES tolerance 

118 P.gmres_tol_limit = 1e-10 

119 rkimexref = rk_imex(P, 5) 

120 uref = np.copy(u0) 

121 dt_ref = dt / 10.0 

122 print("Running RK-IMEX reference....") 

123 for _ in range(0, 10 * Nsteps): 

124 uref = rkimexref.timestep(uref, dt_ref) 

125 

126 usplit = unflatten(usplit, 4, P.N[0], P.N[1]) 

127 udirk = unflatten(udirk, 4, P.N[0], P.N[1]) 

128 uimex = unflatten(uimex, 4, P.N[0], P.N[1]) 

129 uref = unflatten(uref, 4, P.N[0], P.N[1]) 

130 

131 np.save(cwd + 'data/xaxis', P.xx) 

132 np.save(cwd + 'data/sdc', uend) 

133 np.save(cwd + 'data/dirk', udirk) 

134 np.save(cwd + 'data/rkimex', uimex) 

135 np.save(cwd + 'data/split', usplit) 

136 np.save(cwd + 'data/uref', uref) 

137 

138 print("diff split ", np.linalg.norm(uref - usplit)) 

139 print("diff dirk ", np.linalg.norm(uref - udirk)) 

140 print("diff rkimex ", np.linalg.norm(uref - uimex)) 

141 print("diff sdc ", np.linalg.norm(uref - uend)) 

142 

143 print(" #### Logging report for Split #### ") 

144 print("Total number of matrix multiplications: %5i" % splitp.logger.nsmall) 

145 

146 print(" #### Logging report for DIRK-%1i #### " % dirkp.order) 

147 print("Number of calls to implicit solver: %5i" % dirkp.logger.solver_calls) 

148 print("Total number of GMRES iterations: %5i" % dirkp.logger.iterations) 

149 print( 

150 "Average number of iterations per call: %6.3f" 

151 % (float(dirkp.logger.iterations) / float(dirkp.logger.solver_calls)) 

152 ) 

153 print(" ") 

154 print(" #### Logging report for RK-IMEX-%1i #### " % rkimex.order) 

155 print("Number of calls to implicit solver: %5i" % rkimex.logger.solver_calls) 

156 print("Total number of GMRES iterations: %5i" % rkimex.logger.iterations) 

157 print( 

158 "Average number of iterations per call: %6.3f" 

159 % (float(rkimex.logger.iterations) / float(rkimex.logger.solver_calls)) 

160 ) 

161 print(" ") 

162 print(" #### Logging report for SDC-(%1i,%1i) #### " % (sweeper_params['num_nodes'], step_params['maxiter'])) 

163 print("Number of calls to implicit solver: %5i" % P.gmres_logger.solver_calls) 

164 print("Total number of GMRES iterations: %5i" % P.gmres_logger.iterations) 

165 print( 

166 "Average number of iterations per call: %6.3f" 

167 % (float(P.gmres_logger.iterations) / float(P.gmres_logger.solver_calls)) 

168 ) 

169 

170 

171if __name__ == "__main__": 

172 main()