Coverage for pySDC/implementations/problem_classes/Brusselator.py: 86%
44 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
2from mpi4py import MPI
4from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT
7class Brusselator(IMEX_Laplacian_MPIFFT):
8 r"""
9 2D periodic Brusselator reaction-diffusion system with mpi4py-fft, IMEX with diffusion implicit, reactions explicit.
11 This Brusselator from [1]_ is a reaction-diffusion equation with non-autonomous source term:
13 .. math::
14 \frac{\partial u}{\partial t} = \alpha \Delta u + 1 + u^2 v - 4.4u + f(x,y,t),
16 \frac{\partial v}{\partial t} = \alpha \Delta v + 3.4u - u^2 v
18 with the source term :math:`f(x,y,t) = 5` if :math:`(x-0.3)^2 + (y-0.6)^2 <= 0.1^2` and :math:`t >= 1.1` and 0 else.
19 We discretize in a periodic domain of length 1 and solve with an IMEX scheme based on a spectral method for the
20 Laplacian which we invert implicitly. We treat the reaction and source terms explicitly.
22 Parameters
23 ----------
24 alpha : float, optional
25 Diffusion coefficient in front of the Laplacian, the same for both components.
26 **kwargs
27 Passed on to ``IMEX_Laplacian_MPIFFT``: ``nvars`` (tuple of two int, spatial resolution), ``comm``
28 (communicator for ``mpi4py-fft``), ``useGPU`` and ``x0`` (left end of the domain :math:`[x_0, x_0 + 1]^2`).
29 ``spectral=False``, ``L=1.0`` and ``dtype='d'`` are fixed and cannot be passed.
31 References
32 ----------
33 .. [1] https://link.springer.com/book/10.1007/978-3-642-05221-7
34 """
36 def __init__(self, alpha=0.1, **kwargs):
37 """Initialization routine"""
38 super().__init__(spectral=False, L=1.0, dtype='d', alpha=alpha, **kwargs)
40 # prepare the array with two components
41 shape = (2,) + (self.init[0])
42 self.iU = 0
43 self.iV = 1
44 self.ncomp = 2 # needed for transfer class
45 self.init = (shape, self.comm, np.dtype('float'))
47 def _eval_explicit_part(self, u, t, f_expl):
48 iU, iV = self.iU, self.iV
49 x, y = self.X[0], self.X[1]
51 # evaluate time independent part
52 f_expl[iU, ...] = 1.0 + u[iU] ** 2 * u[iV] - 4.4 * u[iU]
53 f_expl[iV, ...] = 3.4 * u[iU] - u[iU] ** 2 * u[iV]
55 # add time-dependent part
56 if t >= 1.1:
57 mask = (x - 0.3) ** 2 + (y - 0.6) ** 2 <= 0.1**2
58 f_expl[iU][mask] += 5.0
59 return f_expl
61 def eval_f(self, u, t):
62 """
63 Routine to evaluate the right-hand side of the problem.
65 Parameters
66 ----------
67 u : dtype_u
68 Current values of the numerical solution.
69 t : float
70 Current time of the numerical solution is computed.
72 Returns
73 -------
74 f : dtype_f
75 The right-hand side of the problem.
76 """
77 f = self.dtype_f(self.init)
79 # evaluate Laplacian to be solved implicitly
80 for i in [self.iU, self.iV]:
81 f.impl[i, ...] = self._eval_Laplacian(u[i], f.impl[i])
83 f.expl[:] = self._eval_explicit_part(u, t, f.expl)
85 self.work_counters['rhs']()
87 return f
89 def solve_system(self, rhs, factor, u0, t):
90 """
91 Simple FFT solver for the diffusion part.
93 Parameters
94 ----------
95 rhs : dtype_f
96 Right-hand side for the linear system.
97 factor : float
98 Abbrev. for the node-to-node stepsize (or any other factor required).
99 u0 : dtype_u
100 Initial guess for the iterative solver (not used here so far).
101 t : float
102 Current time (e.g. for time-dependent BCs).
104 Returns
105 -------
106 me : dtype_u
107 Solution.
108 """
109 me = self.dtype_u(self.init)
111 for i in [self.iU, self.iV]:
112 me[i, ...] = self._invert_Laplacian(me[i], factor, rhs[i])
114 return me
116 def u_exact(self, t, u_init=None, t_init=None):
117 r"""
118 Initial conditions.
120 Parameters
121 ----------
122 t : float
123 Time of the exact solution.
124 u_init : dtype_u
125 Initial conditions for getting the exact solution.
126 t_init : float
127 The starting time.
129 Returns
130 -------
131 me : dtype_u
132 Exact solution.
133 """
135 iU, iV = self.iU, self.iV
136 x, y = self.X[0], self.X[1]
138 me = self.dtype_u(self.init, val=0.0)
140 if t == 0:
141 me[iU, ...] = 22.0 * y * (1 - y / self.L[0]) ** (3.0 / 2.0) / self.L[0]
142 me[iV, ...] = 27.0 * x * (1 - x / self.L[0]) ** (3.0 / 2.0) / self.L[0]
143 else:
145 def eval_rhs(t, u):
146 f = self.eval_f(u.reshape(self.init[0]), t)
147 return (f.impl + f.expl).flatten()
149 me[...] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
151 return me
153 def get_fig(self): # pragma: no cover
154 """
155 Get a figure suitable to plot the solution of this problem
157 Returns
158 -------
159 self.fig : matplotlib.pyplot.figure.Figure
160 """
161 import matplotlib.pyplot as plt
162 from mpl_toolkits.axes_grid1 import make_axes_locatable
164 self.fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=((8, 3)), constrained_layout=True)
165 divider = make_axes_locatable(axs[1])
166 self.cax = divider.append_axes('right', size='3%', pad=0.03)
167 return self.fig
169 def plot(self, u, t=None, fig=None): # pragma: no cover
170 r"""
171 Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
173 Parameters
174 ----------
175 u : dtype_u
176 Solution to be plotted
177 t : float
178 Time to display at the top of the figure
179 fig : matplotlib.pyplot.figure.Figure
180 Figure with the correct structure
182 Returns
183 -------
184 None
185 """
186 fig = self.get_fig() if fig is None else fig
187 axs = fig.axes
189 vmin = u.min()
190 vmax = u.max()
191 for i, label in zip([self.iU, self.iV], [r'$u$', r'$v$'], strict=True):
192 im = axs[i].pcolormesh(self.X[0], self.X[1], u[i], vmin=vmin, vmax=vmax)
193 axs[i].set_aspect(1)
194 axs[i].set_title(label)
196 if t is not None:
197 fig.suptitle(f't = {t:.2e}')
198 axs[0].set_xlabel(r'$x$')
199 axs[0].set_ylabel(r'$y$')
200 fig.colorbar(im, self.cax)