Coverage for pySDC/projects/TOMS/pySDC_with_PETSc.py: 0%

1import sys 

2 

3import numpy as np 

4from mpi4py import MPI 

5 

6from pySDC.helpers.stats_helper import get_sorted 

7 

8from pySDC.implementations.controller_classes.controller_MPI import controller_MPI 

9from pySDC.implementations.problem_classes.HeatEquation_2D_PETSc_forced import heat2d_petsc_forced 

10from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

11from pySDC.implementations.transfer_classes.TransferPETScDMDA import mesh_to_mesh_petsc_dmda 

12 

13 

14def main(): 

15 """ 

16 Program to demonstrate usage of PETSc data structures and spatial parallelization, 

17 combined with parallelization in time. 

18 """ 

19 # set MPI communicator 

20 comm = MPI.COMM_WORLD 

21 

22 world_rank = comm.Get_rank() 

23 world_size = comm.Get_size() 

24 

25 # split world communicator to create space-communicators 

26 if len(sys.argv) >= 2: 

27 color = int(world_rank / int(sys.argv[1])) 

28 else: 

29 color = int(world_rank / 1) 

30 space_comm = comm.Split(color=color) 

31 space_size = space_comm.Get_size() 

32 space_rank = space_comm.Get_rank() 

33 

34 # split world communicator to create time-communicators 

35 if len(sys.argv) >= 2: 

36 color = int(world_rank % int(sys.argv[1])) 

37 else: 

38 color = int(world_rank / world_size) 

39 time_comm = comm.Split(color=color) 

40 time_size = time_comm.Get_size() 

41 time_rank = time_comm.Get_rank() 

42 

43 print( 

44 "IDs (world, space, time): %i / %i -- %i / %i -- %i / %i" 

45 % (world_rank, world_size, space_rank, space_size, time_rank, time_size) 

46 ) 

47 

48 # initialize level parameters 

49 level_params = dict() 

50 level_params['restol'] = 1e-08 

51 level_params['dt'] = 0.125 

52 level_params['nsweeps'] = [3, 1] 

53 

54 # initialize sweeper parameters 

55 sweeper_params = dict() 

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

57 sweeper_params['num_nodes'] = [5] 

58 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part 

59 sweeper_params['initial_guess'] = 'zero' 

60 

61 # initialize problem parameters 

62 problem_params = dict() 

63 problem_params['nu'] = 1.0 # diffusion coefficient 

64 problem_params['freq'] = 2 # frequency for the test value 

65 problem_params['cnvars'] = [(129, 129)] # number of degrees of freedom on coarse level 

66 problem_params['refine'] = [1, 0] # number of refinements 

67 problem_params['comm'] = space_comm # pass space-communicator to problem class 

68 problem_params['sol_tol'] = 1e-10 # set tolerance to PETSc' linear solver 

69 

70 # initialize step parameters 

71 step_params = dict() 

72 step_params['maxiter'] = 50 

73 

74 # initialize space transfer parameters 

75 space_transfer_params = dict() 

76 space_transfer_params['rorder'] = 2 

77 space_transfer_params['iorder'] = 2 

78 space_transfer_params['periodic'] = False 

79 

80 # initialize controller parameters 

81 controller_params = dict() 

82 controller_params['logger_level'] = 20 if space_rank == 0 else 99 # set level depending on rank 

83 controller_params['dump_setup'] = False 

84 

85 # fill description dictionary for easy step instantiation 

86 description = dict() 

87 description['problem_class'] = heat2d_petsc_forced # pass problem class 

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

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

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

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

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

93 description['space_transfer_class'] = mesh_to_mesh_petsc_dmda # pass spatial transfer class 

94 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer 

95 

96 # set time parameters 

97 t0 = 0.0 

98 Tend = 3.0 

99 

100 # instantiate controller 

101 controller = controller_MPI(controller_params=controller_params, description=description, comm=time_comm) 

102 

103 # get initial values on finest level 

104 P = controller.S.levels[0].prob 

105 uinit = P.u_exact(t0) 

106 

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

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

109 

110 # compute exact solution and compare 

111 uex = P.u_exact(Tend) 

112 err = abs(uex - uend) 

113 

114 # filter statistics by type (number of iterations) 

115 iter_counts = get_sorted(stats, type='niter', sortby='time') 

116 

117 niters = np.array([item[1] for item in iter_counts]) 

118 

119 # limit output to space-rank 0 (as before when setting the logger level) 

120 if space_rank == 0: 

121 out = 'This is time-rank %i...' % time_rank 

122 print(out) 

123 

124 # compute and print statistics 

125 for item in iter_counts: 

126 out = 'Number of iterations for time %4.2f: %2i' % item 

127 print(out) 

128 

129 out = ' Mean number of iterations: %4.2f' % np.mean(niters) 

130 print(out) 

131 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters) 

132 print(out) 

133 out = ' Position of max/min number of iterations: %2i -- %2i' % ( 

134 int(np.argmax(niters)), 

135 int(np.argmin(niters)), 

136 ) 

137 print(out) 

138 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters))) 

139 print(out) 

140 

141 print(' Iteration count linear solver: %i' % P.ksp_itercount) 

142 print(' Mean Iteration count per call: %4.2f' % (P.ksp_itercount / max(P.ksp_ncalls, 1))) 

143 

144 timing = get_sorted(stats, type='timing_run', sortby='time') 

145 

146 out = 'Time to solution: %6.4f sec.' % timing[0][1] 

147 print(out) 

148 out = 'Error vs. PDE solution: %6.4e' % err 

149 print(out) 

150 

151 space_comm.Free() 

152 time_comm.Free() 

153 

154 

155if __name__ == "__main__": 

156 main()