Coverage for pySDC/tutorial/step_8/B_multistep_SDC.py: 100%

1# --- 

2# jupyter: 

3# jupytext: 

4# formats: py:percent 

5# kernelspec: 

6# display_name: Python 3 

7# name: python3 

8# --- 

9 

10# %% [markdown] 

11# # Part B: Multi-step SDC 

12# 

13# What happens if we want parallel time steps, but only a single level? The result is **multi-step SDC** 

14# (MSSDC): after each sweep, a step sends its result forward, and the next step picks it up in its next iteration, 

15# not in the current one. This is like doing only the smoother of a multigrid method. Parallelization is dead 

16# simple and needs no coarsening. Yet, without the stabilization of a coarse level, the iteration counts grow 

17# significantly when more steps are computed in parallel. The information can also be picked up immediately, but 

18# then the steps wait for each other and the method is not parallel any more. 

19# 

20# The controller parameter `mssdc_jac` chooses between the two: `True` for the "parallel", Jacobi-like variant, 

21# `False` for the "serial", Gauss-Seidel-like one. We compare both with PFASST, with 8 parallel steps. 

22 

23# %% 

24import os 

25from pathlib import Path 

26 

27import matplotlib.pyplot as plt 

28 

29from pySDC.helpers.stats_helper import get_sorted 

30from pySDC.helpers.visualization_tools import show_residual_across_simulation 

31from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

32from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

33from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

34from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh 

35 

36# initialize level parameters 

37level_params = {'restol': 5e-10, 'dt': 0.125} 

38 

39# initialize sweeper parameters 

40sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': [3], 'QI': 'LU'} 

41 

42# initialize problem parameters 

43problem_params = { 

44 'nu': 0.1, # diffusion coefficient 

45 'freq': 2, # frequency for the test value 

46 'bc': 'dirichlet-zero', # boundary conditions 

47} 

48 

49# initialize step parameters 

50step_params = {'maxiter': 50} 

51 

52# initialize space transfer parameters 

53space_transfer_params = {'rorder': 2, 'iorder': 6} 

54 

55# initialize controller parameters 

56controller_params = {'logger_level': 40} 

57 

58# fill description dictionary for easy step instantiation 

59description = { 

60 'problem_class': heatNd_unforced, 

61 'sweeper_class': generic_implicit, 

62 'sweeper_params': sweeper_params, 

63 'level_params': level_params, 

64 'step_params': step_params, 

65 'space_transfer_class': mesh_to_mesh, 

66 'space_transfer_params': space_transfer_params, 

67} 

68 

69# set up parameters for PFASST run 

70problem_params['nvars'] = [63, 31] 

71description['problem_params'] = problem_params.copy() 

72description_pfasst = description.copy() 

73 

74# set up parameters for MSSDC run 

75problem_params['nvars'] = [63] 

76description['problem_params'] = problem_params.copy() 

77description_mssdc = description.copy() 

78 

79controller_params['mssdc_jac'] = True 

80controller_params_jac = controller_params.copy() 

81controller_params['mssdc_jac'] = False 

82controller_params_gs = controller_params.copy() 

83 

84# %% [markdown] 

85# The multi-step SDC descriptions are copies of the PFASST one, down to the transfer class, which a single level 

86# does not need. pySDC warns about that, and the `logger_level` of 40 silences these (safe) warnings. 

87# 

88# ## Three runs 

89 

90# %% 

91# set time parameters 

92t0 = 0.0 

93Tend = 1.0 

94 

95# set up list of parallel time-steps to run PFASST/MSSDC with 

96num_proc = 8 

97 

98# instantiate controllers 

99controller_mssdc_jac = controller_nonMPI( 

100 num_procs=num_proc, controller_params=controller_params_jac, description=description_mssdc 

101) 

102controller_mssdc_gs = controller_nonMPI( 

103 num_procs=num_proc, controller_params=controller_params_gs, description=description_mssdc 

104) 

105controller_pfasst = controller_nonMPI( 

106 num_procs=num_proc, controller_params=controller_params, description=description_pfasst 

107) 

108 

109# get initial values on finest level 

110P = controller_mssdc_jac.MS[0].levels[0].prob 

111uinit = P.u_exact(t0) 

112 

113# call main functions to get things done... 

114uend_pfasst, stats_pfasst = controller_pfasst.run(u0=uinit, t0=t0, Tend=Tend) 

115uend_mssdc_jac, stats_mssdc_jac = controller_mssdc_jac.run(u0=uinit, t0=t0, Tend=Tend) 

116uend_mssdc_gs, stats_mssdc_gs = controller_mssdc_gs.run(u0=uinit, t0=t0, Tend=Tend) 

117 

118# compute exact solution and compare for both runs 

119uex = P.u_exact(Tend) 

120err_mssdc_jac = abs(uex - uend_mssdc_jac) 

121err_mssdc_gs = abs(uex - uend_mssdc_gs) 

122err_pfasst = abs(uex - uend_pfasst) 

123diff_jac = abs(uend_mssdc_jac - uend_pfasst) 

124diff_gs = abs(uend_mssdc_gs - uend_pfasst) 

125diff_jac_gs = abs(uend_mssdc_gs - uend_mssdc_jac) 

126 

127print('Error PFASST: %12.8e' % err_pfasst) 

128print('Error parallel MSSDC: %12.8e' % err_mssdc_jac) 

129print('Error serial MSSDC: %12.8e' % err_mssdc_gs) 

130print('Diff PFASST vs. parallel MSSDC: %12.8e' % diff_jac) 

131print('Diff PFASST vs. serial MSSDC: %12.8e' % diff_gs) 

132print('Diff parallel vs. serial MSSDC: %12.8e' % diff_jac_gs) 

133 

134# %% [markdown] 

135# All three arrive at the same solution. The difference is in how long it takes: 

136 

137# %% 

138# convert filtered statistics to list of iterations count, sorted by process 

139iter_counts_pfasst = get_sorted(stats_pfasst, type='niter', sortby='time') 

140iter_counts_mssdc_jac = get_sorted(stats_mssdc_jac, type='niter', sortby='time') 

141iter_counts_mssdc_gs = get_sorted(stats_mssdc_gs, type='niter', sortby='time') 

142 

143# compute and print statistics 

144for item_pfasst, item_mssdc_jac, item_mssdc_gs in zip( 

145 iter_counts_pfasst, iter_counts_mssdc_jac, iter_counts_mssdc_gs, strict=True 

146): 

147 print( 

148 'Number of iterations for time %4.2f (PFASST/parMSSDC/serMSSDC): %2i / %2i / %2i' 

149 % (item_pfasst[0], item_pfasst[1], item_mssdc_jac[1], item_mssdc_gs[1]) 

150 ) 

151 

152# %% tags=["hide-input"] 

153fig, ax = plt.subplots(figsize=(6, 3.5)) 

154for name, counts in [ 

155 ('PFASST', iter_counts_pfasst), 

156 ('parallel MSSDC', iter_counts_mssdc_jac), 

157 ('serial MSSDC', iter_counts_mssdc_gs), 

158]: 

159 ax.plot(range(len(counts)), [n for _, n in counts], 'o-', label=name) 

160ax.set_xlabel('process (time step)') 

161ax.set_ylabel('iterations') 

162ax.legend(frameon=False) 

163ax.grid(alpha=0.3) 

164fig.tight_layout() 

165 

166# %% [markdown] 

167# Parallel multi-step SDC needs more and more iterations along the block, up to more than twice as many as PFASST. 

168# The serial variant needs fewer, but its steps wait for each other. The residuals, with the helper 

169# `show_residual_across_simulation` from `pySDC.helpers.visualization_tools`, for the parallel and then the serial 

170# variant, over the processes (rows) and iterations (columns): 

171 

172# %% 

173Path("data").mkdir(parents=True, exist_ok=True) 

174# call helper routine to produce residual plot 

175fig_jac = show_residual_across_simulation(stats_mssdc_jac, 'data/step_8_residuals_mssdc_jac.png') 

176fig_gs = show_residual_across_simulation(stats_mssdc_gs, 'data/step_8_residuals_mssdc_gs.png') 

177 

178# %% [markdown] 

179# The checks the tests run: 

180 

181# %% 

182assert os.path.isfile('data/step_8_residuals_mssdc_jac.png') 

183assert os.path.isfile('data/step_8_residuals_mssdc_gs.png') 

184assert ( 

185 diff_jac < 3.1e-10 

186), f"ERROR: difference between PFASST and parallel MSSDC controller is too large, got {diff_jac}" 

187assert diff_gs < 3.1e-10, f"ERROR: difference between PFASST and serial MSSDC controller is too large, got {diff_gs}" 

188assert ( 

189 diff_jac_gs < 3.1e-10 

190), f"ERROR: difference between parallel and serial MSSDC controller is too large, got {diff_jac_gs}"