Coverage for pySDC/projects/Resilience/AC.py: 74%
98 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1# script to run an Allen-Cahn problem
2from pySDC.implementations.problem_classes.AllenCahn_2D_FFT import allencahn2d_imex
3from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
4from pySDC.core.hooks import Hooks
5from pySDC.projects.Resilience.hook import hook_collection, LogData
6from pySDC.projects.Resilience.strategies import merge_descriptions
7import matplotlib.pyplot as plt
8import numpy as np
10from pySDC.core.errors import ConvergenceError
13class allencahn_imex_timeforcing_adaptivity(allencahn2d_imex):
14 r"""
15 Add more source terms to `allencahn_imex_timeforcing` such that the time-scale changes and we can benefit from adaptivity.
16 """
18 def __init__(self, time_freq=2.0, time_dep_strength=1e-2, *args, **kwargs):
19 super().__init__(*args, **kwargs)
20 self._makeAttributeAndRegister('time_freq', 'time_dep_strength', localVars=locals(), readOnly=True)
22 def eval_f(self, u, t):
23 f = super().eval_f(u, t)
24 time_mod = self.get_time_dep_fac(self.time_freq, self.time_dep_strength, t)
26 if self.eps > 0:
27 f.expl = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u)
29 # build sum over RHS without driving force
30 Rt = float(np.sum(f.impl + f.expl))
32 # build sum over driving force term
33 Ht = float(np.sum(6.0 * u * (1.0 - u)))
35 # add/subtract time-dependent driving force
36 if Ht != 0.0:
37 dw = Rt / Ht * time_mod
38 else:
39 dw = 0.0
41 f.expl -= 6.0 * dw * u * (1.0 - u)
43 return f
45 @staticmethod
46 def get_time_dep_fac(time_freq, time_dep_strength, t):
47 return 1 - time_dep_strength * np.sin(time_freq * 2 * np.pi / 0.032 * t)
50def run_AC(
51 custom_description=None,
52 num_procs=1,
53 Tend=1e-2,
54 hook_class=LogData,
55 fault_stuff=None,
56 custom_controller_params=None,
57 imex=False,
58 u0=None,
59 t0=None,
60 use_MPI=False,
61 live_plot=False,
62 FFT=True,
63 time_forcing=True,
64 **kwargs,
65):
66 """
67 Args:
68 custom_description (dict): Overwrite presets
69 num_procs (int): Number of steps for MSSDC
70 Tend (float): Time to integrate to
71 hook_class (pySDC.Hook): A hook to store data
72 fault_stuff (dict): A dictionary with information on how to add faults
73 custom_controller_params (dict): Overwrite presets
74 imex (bool): Solve the problem IMEX or fully implicit
75 u0 (dtype_u): Initial value
76 t0 (float): Starting time
77 use_MPI (bool): Whether or not to use MPI
79 Returns:
80 dict: The stats object
81 controller: The controller
82 bool: If the code crashed
83 """
84 if custom_description is not None:
85 problem_params = custom_description.get('problem_params', {})
86 if 'imex' in problem_params.keys():
87 imex = problem_params['imex']
88 problem_params.pop('imex', None)
89 if 'FFT' in problem_params.keys():
90 FFT = problem_params['FFT']
91 problem_params.pop('FFT', None)
93 # import problem and sweeper class
94 if time_forcing:
95 problem_class = allencahn_imex_timeforcing_adaptivity
96 from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient as sweeper_class
97 elif FFT:
98 from pySDC.implementations.problem_classes.AllenCahn_2D_FFT import allencahn2d_imex as problem_class
99 from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient as sweeper_class
100 elif imex:
101 from pySDC.implementations.problem_classes.AllenCahn_2D_FD import allencahn_semiimplicit as problem_class
102 from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient as sweeper_class
103 else:
104 from pySDC.implementations.problem_classes.AllenCahn_2D_FD import allencahn_fullyimplicit as problem_class
105 from pySDC.projects.Resilience.sweepers import generic_implicit_efficient as sweeper_class
107 level_params = {}
108 level_params['dt'] = 1e-4
109 level_params['restol'] = 1e-8
111 sweeper_params = {}
112 sweeper_params['quad_type'] = 'RADAU-RIGHT'
113 sweeper_params['num_nodes'] = 3
114 sweeper_params['QI'] = 'LU'
115 sweeper_params['QE'] = 'PIC'
117 # problem params
118 fd_params = {
119 'newton_tol': 1e-9,
120 'order': 2,
121 }
122 problem_params = {
123 'nvars': (128, 128),
124 'init_type': 'circle',
125 }
126 if not FFT:
127 problem_params = {**problem_params, **fd_params}
129 step_params = {}
130 step_params['maxiter'] = 5
132 controller_params = {}
133 controller_params['logger_level'] = 30
134 controller_params['hook_class'] = (
135 hook_collection + (hook_class if type(hook_class) == list else [hook_class]) + ([LivePlot] if live_plot else [])
136 )
137 controller_params['mssdc_jac'] = False
139 if custom_controller_params is not None:
140 controller_params = {**controller_params, **custom_controller_params}
142 description = {}
143 description['problem_class'] = problem_class
144 description['problem_params'] = problem_params
145 description['sweeper_class'] = sweeper_class
146 description['sweeper_params'] = sweeper_params
147 description['level_params'] = level_params
148 description['step_params'] = step_params
150 if custom_description is not None:
151 description = merge_descriptions(description, custom_description)
153 t0 = 0.0 if t0 is None else t0
155 controller_args = {
156 'controller_params': controller_params,
157 'description': description,
158 }
159 if use_MPI:
160 from mpi4py import MPI
161 from pySDC.implementations.controller_classes.controller_MPI import controller_MPI
163 comm = kwargs.get('comm', MPI.COMM_WORLD)
164 controller = controller_MPI(**controller_args, comm=comm)
165 P = controller.S.levels[0].prob
166 else:
167 controller = controller_nonMPI(**controller_args, num_procs=num_procs)
168 P = controller.MS[0].levels[0].prob
170 uinit = P.u_exact(t0) if u0 is None else u0
172 if fault_stuff is not None:
173 from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults
175 prepare_controller_for_faults(controller, fault_stuff)
177 crash = False
178 try:
179 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
180 except ConvergenceError as e:
181 print(f'Warning: Premature termination!: {e}')
182 stats = controller.return_stats()
183 crash = True
184 return stats, controller, crash
187def plot_solution(stats): # pragma: no cover
188 import matplotlib.pyplot as plt
189 from pySDC.helpers.stats_helper import get_sorted
191 fig, ax = plt.subplots(1, 1)
193 u = get_sorted(stats, type='u', recomputed=False)
194 for me in u: # pun intended
195 ax.imshow(me[1], vmin=0.0, vmax=1.0)
196 ax.set_title(f't={me[0]:.2e}')
197 plt.pause(1e-1)
199 plt.show()
202class LivePlot(Hooks): # pragma: no cover
203 def __init__(self):
204 super().__init__()
205 self.fig, self.axs = plt.subplots(1, 3, figsize=(12, 4))
206 self.radius = []
207 self.exact_radius = []
208 self.t = []
209 self.dt = []
211 def post_step(self, step, level_number):
212 super().post_step(step, level_number)
213 L = step.levels[level_number]
214 self.t += [step.time + step.dt]
216 # plot solution
217 self.axs[0].cla()
218 if len(L.uend.shape) > 1:
219 self.axs[0].imshow(L.uend, vmin=0.0, vmax=1.0)
221 # plot radius
222 self.axs[1].cla()
223 radius = np.sqrt(np.count_nonzero(L.uend > 0.5) / np.pi) * L.prob.dx
224 exact_radius = np.sqrt(max(L.prob.radius**2 - 2.0 * (L.time + L.dt), 0))
226 self.radius += [radius]
227 self.exact_radius += [exact_radius]
228 self.axs[1].plot(self.t, self.exact_radius, label='exact')
229 self.axs[1].plot(self.t, self.radius, label='numerical')
230 self.axs[1].set_ylim([0, 0.26])
231 self.axs[1].set_xlim([0, 0.03])
232 self.axs[1].legend(frameon=False)
233 self.axs[1].set_title(r'Radius')
234 else:
235 self.axs[0].plot(L.prob.xvalues, L.prob.u_exact(t=L.time + L.dt), label='exact')
236 self.axs[0].plot(L.prob.xvalues, L.uend, label='numerical')
237 self.axs[0].set_title(f't = {step.time + step.dt:.2e}')
239 # plot step size
240 self.axs[2].cla()
241 self.dt += [step.dt]
242 self.axs[2].plot(self.t, self.dt)
243 self.axs[2].set_yscale('log')
244 self.axs[2].axhline(step.levels[level_number].prob.eps**2, label=r'$\epsilon^2$', color='black', ls='--')
245 self.axs[2].legend(frameon=False)
246 self.axs[2].set_xlim([0, 0.03])
247 self.axs[2].set_title(r'$\Delta t$')
249 if step.status.restart:
250 for me in [self.radius, self.exact_radius, self.t, self.dt]:
251 try:
252 me.pop(-1)
253 except (TypeError, IndexError):
254 pass
256 plt.pause(1e-9)
259if __name__ == '__main__':
261 stats, _, _ = run_AC()
262 plot_solution(stats)