Coverage for pySDC/projects/DAE/run/accuracy_check_MPI.py: 100%

1import numpy as np 

2from mpi4py import MPI 

3 

4from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

5from pySDC.projects.DAE.misc.hooksDAE import ( 

6 LogGlobalErrorPostStepDifferentialVariable, 

7 LogGlobalErrorPostStepAlgebraicVariable, 

8) 

9from pySDC.helpers.stats_helper import get_sorted 

10 

11 

12def run(dt, num_nodes, use_MPI, semi_implicit, residual_type, index_case, initial_guess='spread', comm=None): 

13 r""" 

14 Prepares the controller with all the description needed. Here, the function decides to choose the correct sweeper 

15 for the test. 

16 

17 Parameters 

18 ---------- 

19 dt : float 

20 Time step size chosen for simulation. 

21 num_nodes : int 

22 Number of collocation nodes. 

23 use_MPI : bool 

24 If True, the MPI sweeper classes are used. 

25 semi_implicit : bool 

26 Modules are loaded either for fully-implicit case or semi-implicit case. 

27 residual_type : str 

28 Choose how to compute the residual. 

29 index_case : int 

30 Denotes the index case of a DAE to be tested here, can be either ``1`` or ``2``. 

31 initial_guess : str, optional 

32 Type of initial guess for simulation. 

33 comm : mpi4py.MPI.COMM_WORLD 

34 Communicator. 

35 """ 

36 

37 if not semi_implicit: 

38 if use_MPI: 

39 from pySDC.projects.DAE.sweepers.fullyImplicitDAEMPI import FullyImplicitDAEMPI as sweeper 

40 

41 else: 

42 from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE as sweeper 

43 

44 else: 

45 if use_MPI: 

46 from pySDC.projects.DAE.sweepers.semiImplicitDAEMPI import SemiImplicitDAEMPI as sweeper 

47 

48 else: 

49 from pySDC.projects.DAE.sweepers.semiImplicitDAE import SemiImplicitDAE as sweeper 

50 

51 if index_case == 1: 

52 from pySDC.projects.DAE.problems.discontinuousTestDAE import DiscontinuousTestDAE as problem 

53 

54 t0 = 1.0 

55 Tend = 1.5 

56 

57 elif index_case == 2: 

58 from pySDC.projects.DAE.problems.simpleDAE import SimpleDAE as problem 

59 

60 t0 = 0.0 

61 Tend = 0.4 

62 

63 else: 

64 raise NotImplementedError(f"DAE case of index {index_case} is not implemented!") 

65 

66 # initialize level parameters 

67 level_params = { 

68 'restol': 1e-12, 

69 'residual_type': residual_type, 

70 'dt': dt, 

71 } 

72 

73 # initialize problem parameters 

74 problem_params = { 

75 'newton_tol': 1e-6, 

76 } 

77 

78 # initialize sweeper parameters 

79 sweeper_params = { 

80 'quad_type': 'RADAU-RIGHT', 

81 'num_nodes': num_nodes, 

82 'QI': 'MIN-SR-S', # use a diagonal Q_Delta here! 

83 'initial_guess': initial_guess, 

84 } 

85 

86 # check if number of processes requested matches with number of nodes 

87 if comm is not None: 

88 sweeper_params.update({'comm': comm}) 

89 assert ( 

90 sweeper_params['num_nodes'] == comm.Get_size() 

91 ), f"Number of nodes does not match with number of processes! Expected {sweeper_params['num_nodes']}, got {comm.Get_size()}!" 

92 

93 # initialize step parameters 

94 step_params = { 

95 'maxiter': 20, 

96 } 

97 

98 # initialize controller parameters 

99 controller_params = { 

100 'logger_level': 30, 

101 'hook_class': [LogGlobalErrorPostStepDifferentialVariable, LogGlobalErrorPostStepAlgebraicVariable], 

102 } 

103 

104 # fill description dictionary for easy step instantiation 

105 description = { 

106 'problem_class': problem, 

107 'problem_params': problem_params, 

108 'sweeper_class': sweeper, 

109 'sweeper_params': sweeper_params, 

110 'level_params': level_params, 

111 'step_params': step_params, 

112 } 

113 

114 # instantiate controller 

115 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) 

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

117 

118 uinit = P.u_exact(t0) 

119 

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

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

122 controller.MS[0].levels[0].sweep.compute_end_point() 

123 

124 residual = controller.MS[0].levels[0].status.residual 

125 

126 return uend, residual, stats 

127 

128 

129def check_order(comm): 

130 num_nodes = comm.Get_size() 

131 use_MPI = True 

132 residual_type = 'full_abs' 

133 for semi_implicit in [False, True]: 

134 for index_case in [1, 2]: 

135 dt_list = np.logspace(-1.7, -1.0, num=5) 

136 

137 errorsDiff, errorsAlg = np.zeros(len(dt_list)), np.zeros(len(dt_list)) 

138 for i, dt in enumerate(dt_list): 

139 _, _, stats = run( 

140 dt=dt, 

141 num_nodes=num_nodes, 

142 use_MPI=use_MPI, 

143 semi_implicit=semi_implicit, 

144 residual_type=residual_type, 

145 index_case=index_case, 

146 comm=comm, 

147 ) 

148 

149 errorsDiff[i] = max( 

150 np.array( 

151 get_sorted(stats, type='e_global_differential_post_step', sortby='time', recomputed=False) 

152 )[:, 1] 

153 ) 

154 errorsAlg[i] = max( 

155 np.array(get_sorted(stats, type='e_global_algebraic_post_step', sortby='time', recomputed=False))[ 

156 :, 1 

157 ] 

158 ) 

159 

160 # only process with index 0 should plot 

161 if comm.Get_rank() == 0: 

162 orderDiff = np.mean( 

163 [ 

164 np.log(errorsDiff[i] / errorsDiff[i - 1]) / np.log(dt_list[i] / dt_list[i - 1]) 

165 for i in range(1, len(dt_list)) 

166 ] 

167 ) 

168 orderAlg = np.mean( 

169 [ 

170 np.log(errorsAlg[i] / errorsAlg[i - 1]) / np.log(dt_list[i] / dt_list[i - 1]) 

171 for i in range(1, len(dt_list)) 

172 ] 

173 ) 

174 

175 refOrderDiff = 2 * comm.Get_size() - 1 

176 refOrderAlg = 2 * comm.Get_size() - 1 if index_case == 1 else comm.Get_size() 

177 assert np.isclose( 

178 orderDiff, refOrderDiff, atol=1e0 

179 ), f"Expected order {refOrderDiff} in differential variable, got {orderDiff}" 

180 assert np.isclose( 

181 orderAlg, refOrderAlg, atol=1e0 

182 ), f"Expected order {refOrderAlg} in algebraic variable, got {orderAlg}" 

183 

184 

185if __name__ == "__main__": 

186 comm = MPI.COMM_WORLD 

187 check_order(comm)