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

1# --- 

2# jupyter: 

3# jupytext: 

4# formats: py:percent 

5# kernelspec: 

6# display_name: Python 3 

7# name: python3 

8# language_info: 

9# name: python 

10# --- 

11 

12# %% [markdown] 

13# # Part A: pySDC and FEniCS 

14# 

15# In this example, pySDC is coupled with the [FEniCS framework](https://fenicsproject.org/) for finite elements in 

16# space. This implies significant changes to the algorithm, depending on whether or not the mass matrix should be 

17# inverted. SDC, MLSDC and PFASST can be used without changes when the right-hand side of the ODE is defined with the 

18# inverse of the mass matrix. Otherwise, the mass matrix has to be used, e.g. in the tau-correction. This example tests 

19# different variants of this methodology for SDC, MLSDC and PFASST. 

20# 

21# ## The setup 

22# 

23# The forced heat equation in 1D, with continuous Lagrange elements of order 4, and for MLSDC a second level on a 

24# mesh half as fine, with the same element order and the same collocation nodes: coarsening in the mesh is where the 

25# multilevel gain is, while coarsening in the nodes or in the element order throws it away. The problem and sweeper 

26# classes are set per variant below. 

27 

28# %% 

29from pathlib import Path 

30import numpy as np 

31 

32from pySDC.helpers.stats_helper import get_sorted 

33 

34from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

35from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import ( 

36 fenics_heat_mass, 

37 fenics_heat, 

38 fenics_heat_mass_timebc, 

39) 

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

41from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics 

42from pySDC.implementations.transfer_classes.BaseTransfer_mass import base_transfer_mass 

43 

44 

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

46 """ 

47 Helper routine to set up parameters 

48 

49 Args: 

50 t0 (float): initial time 

51 ml (bool): use single or multiple levels 

52 

53 Returns: 

54 description and controller_params parameter dictionaries 

55 """ 

56 

57 # initialize level parameters 

58 level_params = dict() 

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

60 level_params['dt'] = 0.2 

61 

62 # initialize step parameters 

63 step_params = dict() 

64 step_params['maxiter'] = 20 

65 

66 # initialize sweeper parameters 

67 sweeper_params = dict() 

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

69 if ml: 

70 # Keep the collocation nodes on both levels. Coarsening them does not help: with M=1 the coarse 

71 # level is asymptotically inert, so MLSDC takes neither more nor fewer iterations than SDC. 

72 sweeper_params['num_nodes'] = [3, 3] 

73 else: 

74 sweeper_params['num_nodes'] = [3] 

75 

76 problem_params = dict() 

77 problem_params['nu'] = 0.1 

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

79 problem_params['c_nvars'] = [128] 

80 problem_params['family'] = 'CG' 

81 problem_params['c'] = 1.0 

82 if ml: 

83 # Coarsen in the mesh only. This is where the multilevel gain actually is: it halves the number of 

84 # iterations (6 -> 3 here, and 11.6 -> 3.8 for PFASST below), for both the mass-inverse and the mass 

85 # formulation. Coarsening in the nodes and the element order instead throws that away. 

86 problem_params['order'] = [4, 4] 

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

88 else: 

89 problem_params['order'] = [4] 

90 problem_params['refinements'] = [1] 

91 

92 # initialize controller parameters 

93 controller_params = dict() 

94 controller_params['logger_level'] = 30 

95 

96 base_transfer_params = dict() 

97 base_transfer_params['finter'] = True 

98 

99 # Fill description dictionary for easy hierarchy creation 

100 description = dict() 

101 description['problem_class'] = None 

102 description['problem_params'] = problem_params 

103 description['sweeper_class'] = None 

104 description['sweeper_params'] = sweeper_params 

105 description['level_params'] = level_params 

106 description['step_params'] = step_params 

107 description['space_transfer_class'] = mesh_to_mesh_fenics 

108 description['base_transfer_params'] = base_transfer_params 

109 

110 return description, controller_params 

111 

112 

113# %% [markdown] 

114# ## The variants 

115# 

116# - `'mass_inv'`: the right-hand side includes the inverse of the mass matrix, with the problem class `fenics_heat`, 

117# so that the standard IMEX sweeper works unchanged. 

118# - `'mass'`: the mass matrix stays on the left, with `fenics_heat_mass` and the sweeper `imex_1st_order_mass`, which 

119# applies it to the initial value and in the residual. Between the levels, `base_transfer_mass` restricts the 

120# quantities that carry the mass matrix, the tau-correction and the initial value, as load vectors. The right-hand 

121# side is one as well, so it is re-evaluated on the fine level instead of being interpolated (`finter=False`). 

122# - `'mass_timebc'`: as `'mass'`, but with time-dependent boundary conditions, `fenics_heat_mass_timebc`. 

123# 

124# Each run prints the error, statistics of the iterations and the time to solution, and appends them to 

125# `data/step_7_A_out.txt`. 

126 

127 

128# %% 

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

130 """ 

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

132 

133 Args: 

134 variant (str): specifies the variant 

135 ml (bool): use single or multiple levels 

136 num_procs (int): number of processors in time 

137 """ 

138 Tend = 1.0 

139 t0 = 0.0 

140 

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

142 

143 if variant == 'mass': 

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

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

146 description['problem_class'] = fenics_heat_mass 

147 description['sweeper_class'] = imex_1st_order_mass 

148 description['base_transfer_class'] = base_transfer_mass 

149 # prolong_f is not available for the mass formulation: f is a load vector there, so it cannot be 

150 # interpolated. base_transfer_mass falls back to prolong, which re-evaluates f on the fine level. 

151 description['base_transfer_params']['finter'] = False 

152 elif variant == 'mass_inv': 

153 description['problem_class'] = fenics_heat 

154 description['sweeper_class'] = imex_1st_order 

155 elif variant == 'mass_timebc': 

156 # Trades accuracy for iterations: converged this runs to 1.7e-07 in 9.4 iterations, and 

157 # stopping 20x earlier costs about a factor two in error to get back to 6. 

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

159 description['problem_class'] = fenics_heat_mass_timebc 

160 description['sweeper_class'] = imex_1st_order_mass 

161 description['base_transfer_class'] = base_transfer_mass 

162 description['base_transfer_params']['finter'] = False 

163 else: 

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

165 

166 # quickly generate block of steps 

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

168 

169 # get initial values on finest level 

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

171 uinit = P.u_exact(0.0) 

172 

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

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

175 

176 # compute exact solution and compare 

177 uex = P.u_exact(Tend) 

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

179 

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

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

182 

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

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

185 print(out) 

186 

187 # filter statistics by type (number of iterations) 

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

189 

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

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

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

193 print(out) 

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

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

196 print(out) 

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

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

199 print(out) 

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

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

202 print(out) 

203 

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

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

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

207 print(out) 

208 

209 # Bounds are meant to catch a regression, not to pin the current numbers: at the committed 

210 # settings the errors are 1.14e-08, or 2.8-3.2e-07 for mass_timebc, whose time-dependent 

211 # boundary data makes it a harder problem; the iteration counts are 6.00 serial, 3.00-3.20 with 

212 # a coarse level and 3.80 on five parallel steps. 

213 max_err = 2e-08 

214 if variant == 'mass_timebc': 

215 # the loosened tolerance stops the parallel run a little earlier still, at a larger error 

216 max_err = 5e-07 if num_procs == 1 else 2e-06 

217 max_niter = (5.0 if ml else 8.0) if num_procs == 1 else 6.0 

218 

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

220 assert err <= max_err, 'Error is too high, got %s' % err 

221 

222 f.write('\n') 

223 print() 

224 f.close() 

225 

226 

227# %% [markdown] 

228# SDC, MLSDC and PFASST with 5 steps in parallel, emulated in one process, each with all three variants. 

229 

230 

231# %% 

232def main(): 

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

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

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

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

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

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

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

240 run_variants(variant='mass', ml=True, num_procs=5) 

241 run_variants(variant='mass_timebc', ml=True, num_procs=5) 

242 

243 

244if __name__ == "__main__": 

245 main() 

246 

247# %% [markdown] 

248# ## Results 

249# 

250# FEniCS does not run in the browser, nor in the environment this website is built in. These are the results of our 

251# CI, which runs this part in an environment with FEniCS, in the run that built this page: 

252# 

253# :::{literalinclude} /../../data/step_7_A_out.txt 

254# :language: text 

255# ::: 

256# 

257# :::{admonition} Important things to note 

258# - Even core routines can be replaced where a method needs it: for the mass-matrix formulation, pySDC also has 

259# `base_transfer_mass`, which the mass variants here use for MLSDC and PFASST. 

260# - The project [Finite elements, the mass-matrix route](../../projects/FEM_with_FEniCS) takes the mass-matrix 

261# formulation further: nonlinear problems, discontinuous elements, three levels, and which coarsening pays. 

262# - It is also valuable to check out the data type and transfer classes required to work with FEniCS. Both can be 

263# found in the `implementations` folder. 

264# :::