Coverage for pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py: 100%

1from pathlib import Path 

2import numpy as np 

3 

4from pySDC.helpers.stats_helper import get_sorted 

5 

6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

7from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import ( 

8 fenics_heat_mass, 

9 fenics_heat, 

10 fenics_heat_mass_timebc, 

11) 

12from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass, imex_1st_order 

13from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics 

14 

15 

16def setup(t0=None, ml=None): 

17 """ 

18 Helper routine to set up parameters 

19 

20 Args: 

21 t0 (float): initial time 

22 ml (bool): use single or multiple levels 

23 

24 Returns: 

25 description and controller_params parameter dictionaries 

26 """ 

27 

28 # initialize level parameters 

29 level_params = dict() 

30 level_params['restol'] = 5e-10 

31 level_params['dt'] = 0.2 

32 

33 # initialize step parameters 

34 step_params = dict() 

35 step_params['maxiter'] = 20 

36 

37 # initialize sweeper parameters 

38 sweeper_params = dict() 

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

40 if ml: 

41 # Note that coarsening in the nodes actually HELPS MLSDC to converge (M=1 is exact on the coarse level) 

42 sweeper_params['num_nodes'] = [3, 1] 

43 else: 

44 sweeper_params['num_nodes'] = [3] 

45 

46 problem_params = dict() 

47 problem_params['nu'] = 0.1 

48 problem_params['t0'] = t0 # ugly, but necessary to set up this ProblemClass 

49 problem_params['c_nvars'] = [128] 

50 problem_params['family'] = 'CG' 

51 problem_params['c'] = 1.0 

52 if ml: 

53 # We can do rather aggressive coarsening here. As long as we have 1 node on the coarse level, all is "well" (ie. 

54 # MLSDC does not take more iterations than SDC, but also not less). If we just coarsen in the refinement (and 

55 # not in the nodes and order, the mass inverse approach is way better, ie. halves the number of iterations! 

56 problem_params['order'] = [4, 1] 

57 problem_params['refinements'] = [1, 0] 

58 else: 

59 problem_params['order'] = [4] 

60 problem_params['refinements'] = [1] 

61 

62 # initialize controller parameters 

63 controller_params = dict() 

64 controller_params['logger_level'] = 30 

65 

66 base_transfer_params = dict() 

67 base_transfer_params['finter'] = True 

68 

69 # Fill description dictionary for easy hierarchy creation 

70 description = dict() 

71 description['problem_class'] = None 

72 description['problem_params'] = problem_params 

73 description['sweeper_class'] = None 

74 description['sweeper_params'] = sweeper_params 

75 description['level_params'] = level_params 

76 description['step_params'] = step_params 

77 description['space_transfer_class'] = mesh_to_mesh_fenics 

78 description['base_transfer_params'] = base_transfer_params 

79 

80 return description, controller_params 

81 

82 

83def run_variants(variant=None, ml=None, num_procs=None): 

84 """ 

85 Main routine to run the different implementations of the heat equation with FEniCS 

86 

87 Args: 

88 variant (str): specifies the variant 

89 ml (bool): use single or multiple levels 

90 num_procs (int): number of processors in time 

91 """ 

92 Tend = 1.0 

93 t0 = 0.0 

94 

95 description, controller_params = setup(t0=t0, ml=ml) 

96 

97 if variant == 'mass': 

98 # Note that we need to reduce the tolerance for the residual here, since otherwise the error will be too high 

99 description['level_params']['restol'] /= 500 

100 description['problem_class'] = fenics_heat_mass 

101 description['sweeper_class'] = imex_1st_order_mass 

102 elif variant == 'mass_inv': 

103 description['problem_class'] = fenics_heat 

104 description['sweeper_class'] = imex_1st_order 

105 elif variant == 'mass_timebc': 

106 # Can increase the tolerance here, errors are higher anyway 

107 description['level_params']['restol'] *= 20 

108 description['problem_class'] = fenics_heat_mass_timebc 

109 description['sweeper_class'] = imex_1st_order_mass 

110 else: 

111 raise NotImplementedError('Variant %s is not implemented' % variant) 

112 

113 # quickly generate block of steps 

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

115 

116 # get initial values on finest level 

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

118 uinit = P.u_exact(0.0) 

119 

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

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

122 

123 # compute exact solution and compare 

124 uex = P.u_exact(Tend) 

125 err = abs(uex - uend) / abs(uex) 

126 

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

128 f = open('data/step_7_A_out.txt', 'a') 

129 

130 out = f'Variant {variant} with ml={ml} and num_procs={num_procs} -- error at time {Tend}: {err}' 

131 f.write(out + '\n') 

132 print(out) 

133 

134 # filter statistics by type (number of iterations) 

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

136 

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

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

139 f.write(out + '\n') 

140 print(out) 

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

142 f.write(out + '\n') 

143 print(out) 

144 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters))) 

145 f.write(out + '\n') 

146 print(out) 

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

148 f.write(out + '\n') 

149 print(out) 

150 

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

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

153 f.write(out + '\n') 

154 print(out) 

155 

156 if num_procs == 1: 

157 assert np.mean(niters) <= 6.0, 'Mean number of iterations is too high, got %s' % np.mean(niters) 

158 if variant == 'mass' or variant == 'mass_inv': 

159 assert err <= 1.14e-08, 'Error is too high, got %s' % err 

160 else: 

161 assert err <= 3.25e-07, 'Error is too high, got %s' % err 

162 else: 

163 assert np.mean(niters) <= 11.6, 'Mean number of iterations is too high, got %s' % np.mean(niters) 

164 assert err <= 1.14e-08, 'Error is too high, got %s' % err 

165 

166 f.write('\n') 

167 print() 

168 f.close() 

169 

170 

171def main(): 

172 run_variants(variant='mass_inv', ml=False, num_procs=1) 

173 run_variants(variant='mass', ml=False, num_procs=1) 

174 run_variants(variant='mass_timebc', ml=False, num_procs=1) 

175 run_variants(variant='mass_inv', ml=True, num_procs=1) 

176 run_variants(variant='mass', ml=True, num_procs=1) 

177 run_variants(variant='mass_timebc', ml=True, num_procs=1) 

178 run_variants(variant='mass_inv', ml=True, num_procs=5) 

179 

180 # WARNING: all other variants do NOT work, either because of FEniCS restrictions (weak forms with different meshes 

181 # will not work together) or because of inconsistent use of the mass matrix (locality condition for the tau 

182 # correction is not satisfied, mass matrix does not permute with restriction). 

183 # run_pfasst_variants(variant='mass', ml=True, num_procs=5) 

184 

185 

186if __name__ == "__main__": 

187 main()