Coverage for pySDC / projects / RayleighBenard / RBC3D_configs.py: 84%

232 statements  

« prev     ^ index     » next       coverage.py v7.13.4, created at 2026-02-12 11:13 +0000

1from pySDC.projects.GPU.configs.base_config import Config 

2 

3 

4def get_config(args): 

5 name = args['config'] 

6 if name == 'RBC3D': 

7 return RayleighBenard3DRegular(args) 

8 elif name in globals().keys(): 

9 return globals()[name](args) 

10 else: 

11 raise NotImplementedError(f'There is no configuration called {name!r}!') 

12 

13 

14class RayleighBenard3DRegular(Config): 

15 sweeper_type = 'IMEX' 

16 Tend = 50 

17 gamma = 1 

18 res_ratio = 1 

19 dealiasing = 3.0 / 2.0 

20 

21 def get_file_name(self): 

22 res = self.args['res'] 

23 return f'{self.base_path}/data/{type(self).__name__}-res{res}.pySDC' 

24 

25 def get_LogToFile(self, *args, **kwargs): 

26 if self.comms[1].rank > 0: 

27 return None 

28 import numpy as np 

29 from pySDC.implementations.hooks.log_solution import LogToFile 

30 

31 LogToFile.filename = self.get_file_name() 

32 LogToFile.time_increment = 5e-1 

33 # LogToFile.allow_overwriting = True 

34 

35 return LogToFile 

36 

37 def get_controller_params(self, *args, **kwargs): 

38 from pySDC.implementations.hooks.log_step_size import LogStepSize 

39 

40 controller_params = super().get_controller_params(*args, **kwargs) 

41 controller_params['hook_class'] += [LogStepSize] 

42 return controller_params 

43 

44 def get_description(self, *args, MPIsweeper=False, res=-1, **kwargs): 

45 from pySDC.implementations.problem_classes.RayleighBenard3D import ( 

46 RayleighBenard3D, 

47 ) 

48 from pySDC.implementations.problem_classes.generic_spectral import ( 

49 compute_residual_DAE, 

50 compute_residual_DAE_MPI, 

51 ) 

52 from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeSlopeLimiter 

53 from pySDC.implementations.convergence_controller_classes.crash import StopAtNan 

54 

55 desc = super().get_description(*args, MPIsweeper=MPIsweeper, **kwargs) 

56 

57 if MPIsweeper: 

58 desc['sweeper_class'].compute_residual = compute_residual_DAE_MPI 

59 else: 

60 desc['sweeper_class'].compute_residual = compute_residual_DAE 

61 

62 desc['level_params']['dt'] = 0.01 

63 desc['level_params']['restol'] = 1e-7 

64 

65 desc['convergence_controllers'][StepSizeSlopeLimiter] = {'dt_rel_min_slope': 0.1} 

66 desc['convergence_controllers'][StopAtNan] = {} 

67 

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

69 desc['sweeper_params']['num_nodes'] = 2 

70 desc['sweeper_params']['QI'] = 'MIN-SR-S' 

71 desc['sweeper_params']['QE'] = 'PIC' 

72 

73 res = 64 if res == -1 else res 

74 desc['problem_params']['Rayleigh'] = 1e8 

75 desc['problem_params']['nx'] = self.res_ratio * res 

76 desc['problem_params']['ny'] = self.res_ratio * res 

77 desc['problem_params']['nz'] = res 

78 desc['problem_params']['Lx'] = self.gamma 

79 desc['problem_params']['Ly'] = self.gamma 

80 desc['problem_params']['Lz'] = 1 

81 desc['problem_params']['heterogeneous'] = True 

82 desc['problem_params']['dealiasing'] = self.dealiasing 

83 

84 desc['step_params']['maxiter'] = 3 

85 

86 desc['problem_class'] = RayleighBenard3D 

87 

88 return desc 

89 

90 def get_initial_condition(self, P, *args, restart_idx=0, **kwargs): 

91 

92 if restart_idx == 0: 

93 u0 = P.u_exact(t=0, seed=P.comm.rank, noise_level=1e-3) 

94 u0_with_pressure = P.solve_system(u0, 1e-9, u0) 

95 P.cached_factorizations.pop(1e-9) 

96 return u0_with_pressure, 0 

97 else: 

98 from pySDC.helpers.fieldsIO import FieldsIO 

99 

100 P.setUpFieldsIO() 

101 outfile = FieldsIO.fromFile(self.get_file_name()) 

102 

103 t0, solution = outfile.readField(restart_idx) 

104 solution = solution[: P.spectral.ncomponents, ...] 

105 

106 u0 = P.u_init 

107 

108 if P.spectral_space: 

109 u0[...] = P.transform(solution) 

110 else: 

111 u0[...] = solution 

112 

113 return u0, t0 

114 

115 def prepare_caches(self, prob): 

116 """ 

117 Cache the fft objects, which are expensive to create on GPU because graphs have to be initialized. 

118 """ 

119 prob.eval_f(prob.u_init) 

120 

121 

122class RBC3Dverification(RayleighBenard3DRegular): 

123 converged = 0 

124 dt = 1e-2 

125 ic_config = { 

126 'config': None, 

127 'res': -1, 

128 'dt': -1, 

129 } 

130 res = None 

131 Ra = None 

132 Tend = 100 

133 res_ratio = 4 

134 gamma = 4 

135 

136 def get_file_name(self): 

137 res = self.args['res'] 

138 dt = self.args['dt'] 

139 return f'{self.base_path}/data/{type(self).__name__}-res{res}-dt{dt:.0e}.pySDC' 

140 

141 def get_description(self, *args, res=-1, dt=-1, **kwargs): 

142 desc = super().get_description(*args, **kwargs) 

143 desc['level_params']['nsweeps'] = 4 

144 desc['level_params']['restol'] = -1 

145 desc['step_params']['maxiter'] = 1 

146 desc['sweeper_params']['QI'] = 'MIN-SR-S' 

147 desc['sweeper_params']['skip_residual_computation'] = ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE') 

148 desc['sweeper_params']['num_nodes'] = 4 

149 Ra = int(type(self).__name__[-3]) * 10 ** int(type(self).__name__[-1]) 

150 desc['problem_params']['Rayleigh'] = Ra 

151 desc['problem_params']['Prandtl'] = 0.7 

152 

153 _res = self.res if res == -1 else res 

154 desc['problem_params']['nx'] = _res * self.res_ratio 

155 desc['problem_params']['ny'] = _res * self.res_ratio 

156 desc['problem_params']['nz'] = _res 

157 

158 _dt = self.dt if dt == -1 else dt 

159 desc['level_params']['dt'] = _dt 

160 

161 desc['problem_params']['Lx'] = float(self.gamma) 

162 desc['problem_params']['Ly'] = float(self.gamma) 

163 desc['problem_params']['Lz'] = 1.0 

164 return desc 

165 

166 def get_initial_condition(self, P, *args, restart_idx=0, **kwargs): 

167 if self.ic_config['config'] is None or restart_idx != 0: 

168 return super().get_initial_condition(P, *args, restart_idx=restart_idx, **kwargs) 

169 

170 # read initial conditions 

171 from pySDC.helpers.fieldsIO import FieldsIO 

172 

173 ic_config = self.ic_config['config']( 

174 args={**self.args, 'res': self.ic_config['res'], 'dt': self.ic_config['dt']} 

175 ) 

176 desc = ic_config.get_description(res=self.ic_config['res'], dt=self.ic_config['dt']) 

177 ic_nx = desc['problem_params']['nx'] 

178 ic_ny = desc['problem_params']['ny'] 

179 ic_nz = desc['problem_params']['nz'] 

180 

181 _P = type(P)(nx=ic_nx, ny=ic_ny, nz=ic_nz, comm=P.comm, useGPU=P.useGPU) 

182 _P.setUpFieldsIO() 

183 filename = ic_config.get_file_name() 

184 ic_file = FieldsIO.fromFile(filename) 

185 t0, ics = ic_file.readField(-1) 

186 ics = ics[: P.spectral.ncomponents, ...] 

187 P.logger.info(f'Loaded initial conditions from {filename!r} at t={t0}.') 

188 

189 # interpolate the initial conditions using padded transforms 

190 padding = (P.nx / ic_nx, P.ny / ic_ny, P.nz / ic_nz) 

191 P.logger.info(f'Interpolating initial conditions from {ic_nx}x{ic_ny}x{ic_nz} to {P.nx}x{P.ny}x{P.nz}') 

192 

193 ics = _P.xp.array(ics) 

194 _ics_hat = _P.transform(ics) 

195 ics_interpolated = _P.itransform(_ics_hat, padding=padding) 

196 

197 self.get_LogToFile() 

198 

199 P.setUpFieldsIO() 

200 if P.spectral_space: 

201 u0_hat = P.u_init_forward 

202 u0_hat[...] = P.transform(ics_interpolated) 

203 return u0_hat, 0 

204 else: 

205 return ics_interpolated, 0 

206 

207 

208class RBC3DM2K3(RBC3Dverification): 

209 

210 def get_description(self, *args, **kwargs): 

211 desc = super().get_description(*args, **kwargs) 

212 desc['level_params']['nsweeps'] = 3 

213 desc['sweeper_params']['num_nodes'] = 2 

214 return desc 

215 

216 

217class RBC3DM2K2(RBC3Dverification): 

218 

219 def get_description(self, *args, **kwargs): 

220 desc = super().get_description(*args, **kwargs) 

221 desc['level_params']['nsweeps'] = 2 

222 desc['sweeper_params']['num_nodes'] = 2 

223 return desc 

224 

225 

226class RBC3DM3K4(RBC3Dverification): 

227 

228 def get_description(self, *args, **kwargs): 

229 desc = super().get_description(*args, **kwargs) 

230 desc['level_params']['nsweeps'] = 4 

231 desc['sweeper_params']['num_nodes'] = 3 

232 return desc 

233 

234 

235class RBC3DM4K4(RBC3Dverification): 

236 

237 def get_description(self, *args, **kwargs): 

238 desc = super().get_description(*args, **kwargs) 

239 desc['level_params']['nsweeps'] = 4 

240 desc['sweeper_params']['num_nodes'] = 4 

241 return desc 

242 

243 

244class RBC3DverificationRK(RBC3Dverification): 

245 

246 def get_description(self, *args, res=-1, dt=-1, **kwargs): 

247 from pySDC.implementations.sweeper_classes.Runge_Kutta import ARK3 

248 

249 desc = super().get_description(*args, res=res, dt=dt, **kwargs) 

250 desc['level_params']['nsweeps'] = 1 

251 desc['level_params']['restol'] = -1 

252 desc['step_params']['maxiter'] = 1 

253 desc['sweeper_params']['skip_residual_computation'] = ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE') 

254 desc['sweeper_params'].pop('QI') 

255 desc['sweeper_params'].pop('num_nodes') 

256 desc['sweeper_class'] = ARK3 

257 return desc 

258 

259 

260class RBC3DverificationEuler(RBC3DverificationRK): 

261 

262 def get_description(self, *args, res=-1, dt=-1, **kwargs): 

263 from pySDC.implementations.sweeper_classes.Runge_Kutta import IMEXEulerStifflyAccurate 

264 

265 desc = super().get_description(*args, res=res, dt=dt, **kwargs) 

266 desc['sweeper_class'] = IMEXEulerStifflyAccurate 

267 return desc 

268 

269 

270# --- Ra 1e5 --- 

271class RBC3DG4R4SDC22Ra1e5(RBC3DM2K2): 

272 Tend = 200 

273 dt = 4e-2 

274 res = 32 

275 converged = 50 

276 

277 

278class RBC3DG4R4SDC23Ra1e5(RBC3DM2K3): 

279 Tend = 200 

280 dt = 6e-2 

281 res = 32 

282 converged = 50 

283 

284 

285class RBC3DG4R4SDC34Ra1e5(RBC3DM3K4): 

286 Tend = 200 

287 dt = 6e-2 

288 res = 32 

289 converged = 50 

290 

291 

292class RBC3DG4R4SDC44Ra1e5(RBC3DM4K4): 

293 Tend = 200 

294 dt = 6e-2 

295 res = 32 

296 converged = 50 

297 

298 

299class RBC3DG4R4RKRa1e5(RBC3DverificationRK): 

300 Tend = 200 

301 dt = 8e-2 

302 res = 32 

303 converged = 50 

304 

305 

306class RBC3DG4R4EulerRa1e5(RBC3DverificationEuler): 

307 Tend = 200 

308 dt = 8e-2 

309 res = 32 

310 converged = 50 

311 

312 

313# --- Ra 1e6 --- 

314class RBC3DG4R4SDC44Ra1e6(RBC3DM4K4): 

315 Tend = 75 

316 dt = 2e-2 

317 res = 64 

318 # converged = 22 

319 ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02} 

320 

321 

322class RBC3DG4R4SDC23Ra1e6(RBC3DM2K3): 

323 Tend = 75 

324 dt = 1e-2 

325 res = 64 

326 converged = 22 

327 ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02} 

328 

329 

330class RBC3DG4R4RKRa1e6(RBC3DverificationRK): 

331 Tend = 75 

332 dt = 1e-2 

333 res = 64 

334 ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02} 

335 # converged = 22 

336 

337 

338class RBC3DG4R4EulerRa1e6(RBC3DverificationEuler): 

339 Tend = 75 

340 dt = 1e-2 

341 res = 64 

342 ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02} 

343 # converged = 22