Coverage for pySDC/projects/Resilience/RBC.py: 66%
113 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1# script to run a Rayleigh-Benard Convection problem
2from pySDC.implementations.problem_classes.generic_spectral import compute_residual_DAE
3from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
4from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
5from pySDC.projects.Resilience.hook import hook_collection, LogData
6from pySDC.projects.Resilience.strategies import merge_descriptions
7from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient
8from pySDC.projects.Resilience.reachTendExactly import ReachTendExactly
10from pySDC.core.errors import ConvergenceError
12import numpy as np
14PROBLEM_PARAMS = {'Rayleigh': 3.2e5, 'nx': 256, 'nz': 128, 'max_cached_factorizations': 30, 'Lx': 8, 'Lz': 4, 'z0': 0}
17def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None):
18 import pickle
19 import os
21 path = f'data/stats/RBC-u_init-{t:.8f}.pickle'
22 if os.path.exists(path) and not recompute and t_init is None:
23 with open(path, 'rb') as file:
24 data = pickle.load(file)
25 else:
26 from pySDC.helpers.stats_helper import get_sorted
27 from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
29 convergence_controllers = {
30 Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5},
31 }
32 desc = {'convergence_controllers': convergence_controllers}
34 u0 = self._u_exact(0) if u_init is None else u_init
35 t0 = 0 if t_init is None else t_init
37 if t == t0:
38 return u0
39 else:
40 u0 = u0 if _t0 is None else self.u_exact(_t0)
41 _t0 = t0 if _t0 is None else _t0
43 stats, _, _ = run_RBC(Tend=t, u0=u0, t0=_t0, custom_description=desc)
45 u = get_sorted(stats, type='u', recomputed=False)[-1]
46 data = u[1]
48 if t0 == 0:
49 with open(path, 'wb') as file:
50 pickle.dump(data, file)
52 return data
55class RayleighBenardReference(RayleighBenard):
56 """Rayleigh-Benard whose solution at t > 0 is a reference solution computed with `run_RBC` and cached on disk"""
58 _u_exact = RayleighBenard.u_exact
59 u_exact = u_exact
62class imex_1st_order_efficient_DAE(imex_1st_order_efficient):
63 compute_residual = compute_residual_DAE
66def run_RBC(
67 custom_description=None,
68 num_procs=1,
69 Tend=21.0,
70 hook_class=LogData,
71 fault_stuff=None,
72 custom_controller_params=None,
73 u0=None,
74 t0=20.0,
75 use_MPI=False,
76 step_size_rounding=False,
77 **kwargs,
78):
79 """
80 Args:
81 custom_description (dict): Overwrite presets
82 num_procs (int): Number of steps for MSSDC
83 Tend (float): Time to integrate to
84 hook_class (pySDC.Hook): A hook to store data
85 fault_stuff (dict): A dictionary with information on how to add faults
86 custom_controller_params (dict): Overwrite presets
87 u0 (dtype_u): Initial value
88 t0 (float): Starting time
89 use_MPI (bool): Whether or not to use MPI
91 Returns:
92 dict: The stats object
93 controller: The controller
94 bool: If the code crashed
95 """
96 level_params = {}
97 level_params['dt'] = 1e-3
98 level_params['restol'] = -1
100 sweeper_params = {}
101 sweeper_params['quad_type'] = 'RADAU-RIGHT'
102 sweeper_params['num_nodes'] = 3
103 sweeper_params['QI'] = 'LU'
104 sweeper_params['QE'] = 'PIC'
105 sweeper_params['initial_guess'] = 'copy'
107 from mpi4py import MPI
109 problem_params = {'comm': MPI.COMM_SELF, **PROBLEM_PARAMS}
111 step_params = {}
112 step_params['maxiter'] = 5
114 convergence_controllers = {}
115 convergence_controllers[ReachTendExactly] = {'Tend': Tend}
116 from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding
118 if step_size_rounding:
119 convergence_controllers[StepSizeRounding] = {}
121 controller_params = {}
122 controller_params['logger_level'] = 15
123 controller_params['hook_class'] = hook_collection + (hook_class if type(hook_class) == list else [hook_class])
124 controller_params['mssdc_jac'] = False
126 if custom_controller_params is not None:
127 controller_params = {**controller_params, **custom_controller_params}
129 description = {}
130 description['problem_class'] = RayleighBenardReference
131 description['problem_params'] = problem_params
132 description['sweeper_class'] = imex_1st_order_efficient_DAE
133 description['sweeper_params'] = sweeper_params
134 description['level_params'] = level_params
135 description['step_params'] = step_params
136 description['convergence_controllers'] = convergence_controllers
138 if custom_description is not None:
139 description = merge_descriptions(description, custom_description)
141 controller_args = {
142 'controller_params': controller_params,
143 'description': description,
144 }
145 if use_MPI:
146 from mpi4py import MPI
147 from pySDC.implementations.controller_classes.controller_MPI import controller_MPI
149 comm = kwargs.get('comm', MPI.COMM_WORLD)
150 controller = controller_MPI(**controller_args, comm=comm)
151 P = controller.S.levels[0].prob
152 else:
153 controller = controller_nonMPI(**controller_args, num_procs=num_procs)
154 P = controller.MS[0].levels[0].prob
156 t0 = 0.0 if t0 is None else t0
157 uinit = P.u_exact(t0) if u0 is None else u0
159 if fault_stuff is not None:
160 from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults
162 prepare_controller_for_faults(controller, fault_stuff)
164 crash = False
165 try:
166 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
167 except ConvergenceError as e:
168 print(f'Warning: Premature termination!: {e}')
169 stats = controller.return_stats()
170 crash = True
171 return stats, controller, crash
174def generate_data_for_fault_stats(Tend):
175 prob = RayleighBenardReference(**PROBLEM_PARAMS)
176 _ts = np.linspace(0, Tend, Tend * 10 + 1, dtype=float)
177 for i in range(len(_ts) - 1):
178 print(f'Generating reference solution from {_ts[i]:.4e} to {_ts[i+1]:.4e}')
179 prob.u_exact(_ts[i + 1], _t0=_ts[i], recompute=False)
182def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False): # pragma: no cover
183 from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
184 from pySDC.implementations.hooks.log_work import LogSDCIterations, LogWork
185 from pySDC.implementations.convergence_controller_classes.crash import StopAtNan
186 from pySDC.helpers.stats_helper import get_sorted
187 import pickle
188 import os
190 sweeper_params = {'num_nodes': num_nodes}
191 level_params = {'e_tol': e_tol, 'restol': restol}
192 step_params = {'maxiter': 99}
193 convergence_controllers = {StopAtNan: {}}
195 errors = []
196 dts = []
197 for i in range(steps):
198 dts += [dt / 2**i]
200 path = f'data/stats/RBC-u_order-{t:.8f}-{dts[-1]:.8e}-{num_nodes}-{e_tol:.2e}-{restol:.2e}.pickle'
202 if os.path.exists(path) and not recompute:
203 with open(path, 'rb') as file:
204 stats = pickle.load(file)
205 else:
207 level_params['dt'] = dts[-1]
209 desc = {
210 'sweeper_params': sweeper_params,
211 'level_params': level_params,
212 'step_params': step_params,
213 'convergence_controllers': convergence_controllers,
214 }
216 stats, _, _ = run_RBC(
217 Tend=t + dt,
218 t0=t,
219 custom_description=desc,
220 hook_class=[LogGlobalErrorPostRun, LogSDCIterations, LogWork],
221 )
223 with open(path, 'wb') as file:
224 pickle.dump(stats, file)
226 e = get_sorted(stats, type='e_global_post_run')
227 # k = get_sorted(stats, type='k')
229 if len(e) > 0:
230 errors += [e[-1][1]]
231 else:
232 errors += [np.nan]
234 errors = np.array(errors)
235 dts = np.array(dts)
236 mask = np.isfinite(errors)
237 max_error = np.nanmax(errors)
239 errors = errors[mask]
240 dts = dts[mask]
241 ax.loglog(dts, errors, label=f'{num_nodes} nodes', marker='x')
242 ax.loglog(
243 dts, [max_error * (me / dts[0]) ** (2 * num_nodes - 1) for me in dts], ls='--', label=f'order {2*num_nodes-1}'
244 )
245 ax.set_xlabel(r'$\Delta t$')
246 ax.set_ylabel(r'global error')
247 ax.legend(frameon=False)
250def check_order(t=14, dt=1e-1, steps=6):
251 prob = RayleighBenardReference(**PROBLEM_PARAMS)
252 _ts = [0, t, t + dt]
253 for i in range(len(_ts) - 1):
254 prob.u_exact(_ts[i + 1], _t0=_ts[i])
256 import matplotlib.pyplot as plt
258 fig, ax = plt.subplots(1, 1)
259 for num_nodes in [1, 2, 3, 4]:
260 plot_order(t=t, dt=dt, steps=steps, num_nodes=num_nodes, ax=ax, restol=1e-9)
261 ax.set_title(f't={t:.2f}, dt={dt:.2f}')
262 plt.show()
265def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False): # pragma: no cover
266 import matplotlib.pyplot as plt
267 import pickle
268 import os
269 from pySDC.helpers.stats_helper import get_sorted
271 path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}.pickle'
272 if os.path.exists(path) and not recompute:
273 with open(path, 'rb') as file:
274 stats = pickle.load(file)
275 else:
276 from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
278 convergence_controllers = {
279 Adaptivity: {'e_tol': e_tol, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5},
280 }
281 desc = {'convergence_controllers': convergence_controllers}
283 stats, _, _ = run_RBC(Tend=Tend, t0=t0, custom_description=desc)
285 with open(path, 'wb') as file:
286 pickle.dump(stats, file)
288 fig, ax = plt.subplots(1, 1)
290 dt = get_sorted(stats, type='dt', recomputed=False)
291 ax.plot([me[0] for me in dt], [me[1] for me in dt])
292 ax.set_ylabel(r'$\Delta t$')
293 ax.set_xlabel(r'$t$')
294 ax.set_yscale('log')
295 plt.show()
298def plot_factorizations_over_time(t0=0, Tend=50, e_tol=1e-3, recompute=False, adaptivity_mode='dt'): # pragma: no cover
299 import matplotlib.pyplot as plt
300 import pickle
301 import os
302 from pySDC.helpers.stats_helper import get_sorted
303 from pySDC.implementations.hooks.log_work import LogWork
304 from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity, AdaptivityPolynomialError
305 from pySDC.helpers.plot_helper import figsize_by_journal, setup_mpl
306 from pySDC.projects.Resilience.paper_plots import savefig
308 setup_mpl()
310 fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal('TUHH_thesis', 1.0, 0.4))
312 if adaptivity_mode == 'dt':
313 adaptivity = Adaptivity
314 elif adaptivity_mode == 'dt_k':
315 adaptivity = AdaptivityPolynomialError
317 dt_controllers = {
318 'basic': {
319 adaptivity: {
320 'e_tol': e_tol,
321 }
322 },
323 'min slope': {adaptivity: {'e_tol': e_tol, 'beta': 0.5, 'dt_rel_min_slope': 2}},
324 'fixed': {},
325 # 'rounding': {adaptivity: {'e_tol': e_tol, 'beta': 0.5, 'dt_rel_min_slope': 2}, StepSizeRounding: {}},
326 }
328 for name, params in dt_controllers.items():
329 if adaptivity_mode == 'dt':
330 path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}-{name}.pickle'
331 elif adaptivity_mode == 'dt_k':
332 path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}-{name}-dtk.pickle'
334 if os.path.exists(path) and not recompute:
335 with open(path, 'rb') as file:
336 stats = pickle.load(file)
337 else:
339 convergence_controllers = {
340 **params,
341 }
342 desc = {'convergence_controllers': convergence_controllers}
344 if name == 'fixed':
345 if adaptivity_mode == 'dt':
346 desc['level_params'] = {'dt': 2e-2}
347 elif adaptivity_mode == 'dt_k':
348 desc['level_params'] = {'dt': 2e-3}
349 elif adaptivity_mode == 'dt_k':
350 desc['level_params'] = {'restol': 1e-7}
352 stats, _, _ = run_RBC(
353 Tend=Tend, t0=t0, custom_description=desc, hook_class=LogWork, step_size_rounding=False
354 )
356 with open(path, 'wb') as file:
357 pickle.dump(stats, file)
359 factorizations = get_sorted(stats, type='work_factorizations')
360 rhs_evals = get_sorted(stats, type='work_rhs')
361 axs[0].plot([me[0] for me in factorizations], np.cumsum([me[1] for me in factorizations]), label=name)
362 axs[1].plot([me[0] for me in rhs_evals], np.cumsum([me[1] for me in rhs_evals]), label=name)
364 axs[0].set_ylabel(r'matrix factorizations')
365 axs[1].set_ylabel(r'right hand side evaluations')
366 axs[0].set_xlabel(r'$t$')
367 axs[1].set_xlabel(r'$t$')
368 axs[0].set_yscale('log')
369 axs[1].set_yscale('log')
370 axs[0].legend(frameon=False)
371 savefig(fig, f'RBC_step_size_controller_{adaptivity_mode}')
374if __name__ == '__main__':
375 # plot_step_size(0, 30)
376 generate_data_for_fault_stats(Tend=30)
377 # plot_factorizations_over_time(e_tol=1e-3, adaptivity_mode='dt')
378 # plot_factorizations_over_time(recompute=False, e_tol=1e-5, adaptivity_mode='dt_k')
379 # check_order(t=20, dt=1., steps=7)
380 # stats, _, _ = run_RBC()