Coverage for pySDC/projects/AllenCahn_Bayreuth/run_temp_forcing_verification.py: 99%

1from argparse import ArgumentParser 

2import json 

3import numpy as np 

4from mpi4py import MPI 

5from mpi4py_fft import newDistArray 

6 

7from pySDC.helpers.stats_helper import get_sorted 

8 

9from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

10from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

11from pySDC.implementations.problem_classes.AllenCahn_Temp_MPIFFT import allencahn_temp_imex 

12from pySDC.implementations.transfer_classes.TransferMesh_MPIFFT import fft_to_fft 

13 

14 

15def run_simulation(name='', spectral=None, nprocs_time=None, nprocs_space=None, dt=None, cwd='.'): 

16 """ 

17 A test program to do PFASST runs for the AC equation with temperature-based forcing 

18 

19 (slightly inefficient, but will run for a few seconds only) 

20 

21 Args: 

22 name (str): name of the run, will be used to distinguish different setups 

23 spectral (bool): run in real or spectral space 

24 nprocs_time (int): number of processors in time 

25 nprocs_space (int): number of processors in space (None if serial) 

26 dt (float): time-step size 

27 cwd (str): current working directory 

28 """ 

29 

30 # set MPI communicator 

31 comm = MPI.COMM_WORLD 

32 

33 world_rank = comm.Get_rank() 

34 world_size = comm.Get_size() 

35 

36 # split world communicator to create space-communicators 

37 if nprocs_space is not None: 

38 color = int(world_rank / nprocs_space) 

39 else: 

40 color = int(world_rank / 1) 

41 space_comm = comm.Split(color=color) 

42 space_rank = space_comm.Get_rank() 

43 space_size = space_comm.Get_size() 

44 

45 assert world_size == space_size, 'This script cannot run parallel-in-time with MPI, only spatial parallelism' 

46 

47 # initialize level parameters 

48 level_params = dict() 

49 level_params['restol'] = 1e-12 

50 level_params['dt'] = dt 

51 level_params['nsweeps'] = [1] 

52 

53 # initialize sweeper parameters 

54 sweeper_params = dict() 

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

56 sweeper_params['num_nodes'] = [3] 

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

58 sweeper_params['initial_guess'] = 'spread' 

59 

60 # initialize problem parameters 

61 problem_params = dict() 

62 problem_params['L'] = 1.0 

63 problem_params['nvars'] = [(128, 128), (32, 32)] 

64 problem_params['eps'] = [0.04] 

65 problem_params['radius'] = 0.25 

66 problem_params['TM'] = 1.0 

67 problem_params['D'] = 0.1 

68 problem_params['dw'] = [21.0] 

69 problem_params['comm'] = space_comm 

70 problem_params['init_type'] = 'circle' 

71 problem_params['spectral'] = spectral 

72 

73 # initialize step parameters 

74 step_params = dict() 

75 step_params['maxiter'] = 50 

76 

77 # initialize controller parameters 

78 controller_params = dict() 

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

80 controller_params['predict_type'] = 'pfasst_burnin' 

81 

82 # fill description dictionary for easy step instantiation 

83 description = dict() 

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

85 description['sweeper_class'] = imex_1st_order 

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

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

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

89 description['space_transfer_class'] = fft_to_fft 

90 description['problem_class'] = allencahn_temp_imex 

91 

92 # set time parameters 

93 t0 = 0.0 

94 Tend = 1 * 0.001 

95 

96 if space_rank == 0: 

97 out = f'---------> Running {name} with spectral={spectral} and {space_size} process(es) in space...' 

98 print(out) 

99 

100 # instantiate controller 

101 controller = controller_nonMPI(num_procs=nprocs_time, controller_params=controller_params, description=description) 

102 

103 # get initial values on finest level 

104 P = controller.MS[0].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 if space_rank == 0: 

111 # convert filtered statistics of iterations count, sorted by time 

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

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

114 out = f'Mean number of iterations: {niters:.4f}' 

115 print(out) 

116 

117 # get setup time 

118 timing = get_sorted(stats, type='timing_setup', sortby='time') 

119 out = f'Setup time: {timing[0][1]:.4f} sec.' 

120 print(out) 

121 

122 # get running time 

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

124 out = f'Time to solution: {timing[0][1]:.4f} sec.' 

125 print(out) 

126 

127 refname = f'{cwd}/data/AC-reference-tempforce_00001000' 

128 with open(f'{refname}.json', 'r') as fp: 

129 obj = json.load(fp) 

130 

131 array = np.fromfile(f'{refname}.dat', dtype=obj['datatype']) 

132 array = array.reshape(obj['shape'], order='C') 

133 

134 if spectral: 

135 ureal = newDistArray(P.fft, False) 

136 ureal = P.fft.backward(uend[..., 0], ureal) 

137 Treal = newDistArray(P.fft, False) 

138 Treal = P.fft.backward(uend[..., 1], Treal) 

139 err = max(np.amax(abs(ureal - array[..., 0])), np.amax(abs(Treal - array[..., 1]))) 

140 else: 

141 err = abs(array - uend) 

142 

143 out = '...Done <---------\n' 

144 print(out) 

145 

146 return err 

147 

148 

149def main(nprocs_space=None, cwd='.'): 

150 """ 

151 Little helper routine to run the whole thing 

152 

153 Args: 

154 nprocs_space (int): number of processors in space (None if serial) 

155 cwd (str): current working directory 

156 """ 

157 name = 'AC-test-tempforce' 

158 

159 nsteps = [2**i for i in range(4)] 

160 

161 errors = [1] 

162 orders = [] 

163 for n in nsteps: 

164 err = run_simulation(name=name, spectral=False, nprocs_time=n, nprocs_space=nprocs_space, dt=1e-03 / n, cwd=cwd) 

165 errors.append(err) 

166 orders.append(np.log(errors[-1] / errors[-2]) / np.log(0.5)) 

167 print(f'Error: {errors[-1]:6.4e}') 

168 print(f'Order of accuracy: {orders[-1]:4.2f}\n') 

169 

170 assert errors[2 + 1] < 8e-10, f'Errors are too high, got {errors[2 + 1]}' 

171 assert np.isclose(orders[3], 5.3, rtol=2e-02), f'Order of accuracy is not within tolerance, got {orders[3]}' 

172 

173 print() 

174 

175 errors = [1] 

176 orders = [] 

177 for n in nsteps: 

178 err = run_simulation(name=name, spectral=True, nprocs_time=n, nprocs_space=nprocs_space, dt=1e-03 / n, cwd=cwd) 

179 errors.append(err) 

180 orders.append(np.log(errors[-1] / errors[-2]) / np.log(0.5)) 

181 print(f'Error: {errors[-1]:6.4e}') 

182 print(f'Order of accuracy: {orders[-1]:4.2f}\n') 

183 

184 assert errors[2 + 1] < 8e-10, f'Errors are too high, got {errors[2 + 1]}' 

185 assert np.isclose( 

186 orders[1].view(np.ndarray), 4.6, rtol=7e-02 

187 ), f'Order of accuracy is not within tolerance, got {orders[1]}' 

188 

189 

190if __name__ == "__main__": 

191 # Add parser to get number of processors in space (have to do this here to enable automatic testing) 

192 parser = ArgumentParser() 

193 parser.add_argument("-n", "--nprocs_space", help='Specifies the number of processors in space', type=int) 

194 args = parser.parse_args() 

195 

196 main(nprocs_space=args.nprocs_space)