# Coverage for pySDC/implementations/problem_classes/Brusselator.py: 86%

## 44 statements

, created at 2024-09-20 17:10 +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 Two-dimensional Brusselator from [1]_.

10 This is a reaction-diffusion equation with non-autonomous source term:

12 .. math::

13 \frac{\partial u}{\partial t} = \varalpha \Delta u + 1 + u^2 v - 4.4u _ f(x,y,t),

14 \frac{\partial v}{\partial t} = \varalpha \Delta v + 3.4u - u^2 v

16 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.

17 We discretize in a periodic domain of length 1 and solve with an IMEX scheme based on a spectral method for the

18 Laplacian which we invert implicitly. We treat the reaction and source terms explicitly.

20 References

21 ----------

22 .. [1] https://link.springer.com/book/10.1007/978-3-642-05221-7

23 """

25 def __init__(self, alpha=0.1, **kwargs):

26 """Initialization routine"""

27 super().__init__(spectral=False, L=1.0, dtype='d', alpha=alpha, **kwargs)

29 # prepare the array with two components

30 shape = (2,) + (self.init[0])

31 self.iU = 0

32 self.iV = 1

33 self.ncomp = 2 # needed for transfer class

34 self.init = (shape, self.comm, np.dtype('float'))

36 def _eval_explicit_part(self, u, t, f_expl):

37 iU, iV = self.iU, self.iV

38 x, y = self.X[0], self.X[1]

40 # evaluate time independent part

41 f_expl[iU, ...] = 1.0 + u[iU] ** 2 * u[iV] - 4.4 * u[iU]

42 f_expl[iV, ...] = 3.4 * u[iU] - u[iU] ** 2 * u[iV]

44 # add time-dependent part

45 if t >= 1.1:

46 mask = (x - 0.3) ** 2 + (y - 0.6) ** 2 <= 0.1**2

47 f_expl[iU][mask] += 5.0

48 return f_expl

50 def eval_f(self, u, t):

51 """

52 Routine to evaluate the right-hand side of the problem.

54 Parameters

55 ----------

56 u : dtype_u

57 Current values of the numerical solution.

58 t : float

59 Current time of the numerical solution is computed.

61 Returns

62 -------

63 f : dtype_f

64 The right-hand side of the problem.

65 """

66 f = self.dtype_f(self.init)

68 # evaluate Laplacian to be solved implicitly

69 for i in [self.iU, self.iV]:

70 f.impl[i, ...] = self._eval_Laplacian(u[i], f.impl[i])

72 f.expl[:] = self._eval_explicit_part(u, t, f.expl)

74 self.work_counters['rhs']()

76 return f

78 def solve_system(self, rhs, factor, u0, t):

79 """

80 Simple FFT solver for the diffusion part.

82 Parameters

83 ----------

84 rhs : dtype_f

85 Right-hand side for the linear system.

86 factor : float

87 Abbrev. for the node-to-node stepsize (or any other factor required).

88 u0 : dtype_u

89 Initial guess for the iterative solver (not used here so far).

90 t : float

91 Current time (e.g. for time-dependent BCs).

93 Returns

94 -------

95 me : dtype_u

96 Solution.

97 """

98 me = self.dtype_u(self.init)

100 for i in [self.iU, self.iV]:

101 me[i, ...] = self._invert_Laplacian(me[i], factor, rhs[i])

103 return me

105 def u_exact(self, t, u_init=None, t_init=None):

106 r"""

107 Initial conditions.

109 Parameters

110 ----------

111 t : float

112 Time of the exact solution.

113 u_init : dtype_u

114 Initial conditions for getting the exact solution.

115 t_init : float

116 The starting time.

118 Returns

119 -------

120 me : dtype_u

121 Exact solution.

122 """

124 iU, iV = self.iU, self.iV

125 x, y = self.X[0], self.X[1]

127 me = self.dtype_u(self.init, val=0.0)

129 if t == 0:

130 me[iU, ...] = 22.0 * y * (1 - y / self.L[0]) ** (3.0 / 2.0) / self.L[0]

131 me[iV, ...] = 27.0 * x * (1 - x / self.L[0]) ** (3.0 / 2.0) / self.L[0]

132 else:

134 def eval_rhs(t, u):

135 f = self.eval_f(u.reshape(self.init[0]), t)

136 return (f.impl + f.expl).flatten()

138 me[...] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)

140 return me

142 def get_fig(self): # pragma: no cover

143 """

144 Get a figure suitable to plot the solution of this problem

146 Returns

147 -------

148 self.fig : matplotlib.pyplot.figure.Figure

149 """

150 import matplotlib.pyplot as plt

151 from mpl_toolkits.axes_grid1 import make_axes_locatable

153 plt.rcParams['figure.constrained_layout.use'] = True

154 self.fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=((8, 3)))

155 divider = make_axes_locatable(axs[1])

156 self.cax = divider.append_axes('right', size='3%', pad=0.03)

157 return self.fig

159 def plot(self, u, t=None, fig=None): # pragma: no cover

160 r"""

161 Plot the solution. Please supply a figure with the same structure as returned by self.get_fig.

163 Parameters

164 ----------

165 u : dtype_u

166 Solution to be plotted

167 t : float

168 Time to display at the top of the figure

169 fig : matplotlib.pyplot.figure.Figure

170 Figure with the correct structure

172 Returns

173 -------

174 None

175 """

176 fig = self.get_fig() if fig is None else fig

177 axs = fig.axes

179 vmin = u.min()

180 vmax = u.max()

181 for i, label in zip([self.iU, self.iV], [r'$u$', r'$v$']):

182 im = axs[i].pcolormesh(self.X[0], self.X[1], u[i], vmin=vmin, vmax=vmax)

183 axs[i].set_aspect(1)

184 axs[i].set_title(label)

186 if t is not None:

187 fig.suptitle(f't = {t:.2e}')

188 axs[0].set_xlabel(r'$x$')

189 axs[0].set_ylabel(r'$y$')

190 fig.colorbar(im, self.cax)