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

1import numpy as np 

2from mpi4py import MPI 

3 

4from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT 

5 

6 

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: 

11 

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 

15 

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. 

19 

20 References 

21 ---------- 

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

23 """ 

24 

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

26 """Initialization routine""" 

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

28 

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')) 

35 

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] 

39 

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] 

43 

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 

49 

50 def eval_f(self, u, t): 

51 """ 

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

53 

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. 

60 

61 Returns 

62 ------- 

63 f : dtype_f 

64 The right-hand side of the problem. 

65 """ 

66 f = self.dtype_f(self.init) 

67 

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]) 

71 

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

73 

74 self.work_counters['rhs']() 

75 

76 return f 

77 

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

79 """ 

80 Simple FFT solver for the diffusion part. 

81 

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

92 

93 Returns 

94 ------- 

95 me : dtype_u 

96 Solution. 

97 """ 

98 me = self.dtype_u(self.init) 

99 

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

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

102 

103 return me 

104 

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

106 r""" 

107 Initial conditions. 

108 

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. 

117 

118 Returns 

119 ------- 

120 me : dtype_u 

121 Exact solution. 

122 """ 

123 

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

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

126 

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

128 

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: 

133 

134 def eval_rhs(t, u): 

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

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

137 

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

139 

140 return me 

141 

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

143 """ 

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

145 

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 

152 

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 

158 

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

162 

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 

171 

172 Returns 

173 ------- 

174 None 

175 """ 

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

177 axs = fig.axes 

178 

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) 

185 

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)