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

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 

9from pySDC.projects.FastWaveSlowWave.HookClass_boussinesq import gmres_tolerance 

10 

11 

12def main(cwd=''): 

13 """ 

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

15 

16 Args: 

17 cwd (string): current working directory 

18 """ 

19 

20 num_procs = 1 

21 

22 # setup parameters "in time" 

23 t0 = 0 

24 Tend = 3000 

25 Nsteps = 100 

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 controller_params['hook_class'] = gmres_tolerance 

46 

47 # initialize problem parameters 

48 problem_params = dict() 

49 problem_params['nvars'] = [(4, 300, 30)] 

50 problem_params['u_adv'] = 0.02 

51 problem_params['c_s'] = 0.3 

52 problem_params['Nfreq'] = 0.01 

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

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

55 problem_params['order'] = [4] 

56 problem_params['order_upw'] = [5] 

57 problem_params['gmres_maxiter'] = [500] 

58 problem_params['gmres_restart'] = [10] 

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

60 problem_params['gmres_tol_factor'] = [0.1] 

61 

62 # fill description dictionary for easy step instantiation 

63 description = dict() 

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

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

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

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

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

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

70 

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

72 dirk_order = step_params['maxiter'] 

73 

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

75 

76 # get initial values on finest level 

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

78 uinit = P.u_exact(t0) 

79 

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

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

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

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

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

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

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

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

88 

89 print("Running SplitExplicit ....") 

90 method_split = 'MIS4_4' 

91 # method_split = 'RK3' 

92 splitp = SplitExplicit(P, method_split, problem_params) 

93 u0 = uinit.flatten() 

94 usplit = np.copy(u0) 

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

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

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

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

99 

100 print("Running DIRK ....") 

101 dirkp = dirk(P, dirk_order) 

102 udirk = np.copy(u0) 

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

104 for _ in range(0, Nsteps): 

105 udirk = dirkp.timestep(udirk, dt) 

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

107 

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

109 rkimex = rk_imex(P, dirk_order) 

110 uimex = np.copy(u0) 

111 dt_imex = dt 

112 for _ in range(0, Nsteps): 

113 uimex = rkimex.timestep(uimex, dt_imex) 

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

115 

116 print("Running SDC...") 

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

118 

119 # For reference solution, increase GMRES tolerance 

120 P.gmres_tol_limit = 1e-10 

121 rkimexref = rk_imex(P, 5) 

122 uref = np.copy(u0) 

123 dt_ref = dt / 10.0 

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

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

126 uref = rkimexref.timestep(uref, dt_ref) 

127 

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

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

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

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

132 

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

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

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

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

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

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

139 

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

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

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

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

144 

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

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

147 

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

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

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

151 print( 

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

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

154 ) 

155 print(" ") 

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

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

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

159 print( 

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

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

162 ) 

163 print(" ") 

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

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

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

167 print( 

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

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

170 ) 

171 

172 

173if __name__ == "__main__": 

174 main()