Coverage for pySDC/implementations/problem_classes/AllenCahn_MPIFFT.py: 91%

1import numpy as np 

2from mpi4py import MPI 

3 

4from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT 

5from mpi4py_fft import newDistArray 

6 

7 

8class allencahn_imex(IMEX_Laplacian_MPIFFT): 

9 r""" 

10 Periodic Allen-Cahn equation with driving force and mpi4py-fft FFTs, IMEX with the Laplacian implicit. 

11 

12 Example implementing the :math:`2`-dimensional Allen-Cahn equation with periodic boundary conditions, with the two 

13 phases at :math:`u = 0` and :math:`u = 1` 

14 

15 .. math:: 

16 \frac{\partial u}{\partial t} = \Delta u - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

17 - 6 d_w u (1 - u) 

18 

19 on a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`, driving force :math:`d_w`, and :math:`N=2,3`. Different initial 

20 conditions can be used, for example, circles of the form 

21 

22 .. math:: 

23 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}} 

24 {\sqrt{2}\varepsilon}\right)\right), 

25 

26 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is treated 

27 *semi-implicitly*, i.e., the linear part is solved with Fast-Fourier Transform (FFT) and the nonlinear part in the right-hand 

28 side will be treated explicitly using ``mpi4py-fft`` [1]_ to solve them. 

29 

30 Parameters 

31 ---------- 

32 nvars : List of int tuples, optional 

33 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``. 

34 eps : float, optional 

35 Scaling parameter :math:`\varepsilon`. 

36 radius : float, optional 

37 Radius of the circles. 

38 spectral : bool, optional 

39 Indicates if spectral initial condition is used. 

40 dw : float, optional 

41 Driving force. 

42 L : float, optional 

43 Denotes the period of the function to be approximated for the Fourier transform. 

44 init_type : str, optional 

45 Initialises type of initial state. 

46 comm : MPI.Intracomm, optional 

47 Communicator for parallelization. 

48 

49 Attributes 

50 ---------- 

51 fft : fft object 

52 Object for FFT. 

53 X : np.ogrid 

54 Grid coordinates in real space. 

55 K2 : np.1darray 

56 Laplace operator in spectral space. 

57 dx : float 

58 Mesh width in x direction. 

59 dy : float 

60 Mesh width in y direction. 

61 

62 References 

63 ---------- 

64 .. [1] https://mpi4py-fft.readthedocs.io/en/latest/ 

65 """ 

66 

67 def __init__( 

68 self, 

69 eps=0.04, 

70 radius=0.25, 

71 dw=0.0, 

72 init_type='circle', 

73 **kwargs, 

74 ): 

75 """Initialization routine. ``kwargs`` are passed on to ``IMEX_Laplacian_MPIFFT``.""" 

76 kwargs['L'] = kwargs.get('L', 1.0) 

77 kwargs['x0'] = kwargs.get('x0', -kwargs['L'] / 2.0) 

78 super().__init__(alpha=1.0, dtype=np.dtype('float'), **kwargs) 

79 self._makeAttributeAndRegister('eps', 'radius', 'dw', 'init_type', localVars=locals(), readOnly=True) 

80 

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

82 f_expl[:] = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u) - 6.0 * self.dw * u * (1.0 - u) 

83 return f_expl 

84 

85 def eval_f(self, u, t): 

86 """ 

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

88 

89 Parameters 

90 ---------- 

91 u : dtype_u 

92 Current values of the numerical solution. 

93 t : float 

94 Current time of the numerical solution is computed. 

95 

96 Returns 

97 ------- 

98 f : dtype_f 

99 The right-hand side of the problem. 

100 """ 

101 

102 f = self.dtype_f(self.init) 

103 

104 f.impl[:] = self._eval_Laplacian(u, f.impl) 

105 

106 if self.spectral: 

107 if self.eps > 0: 

108 tmp = self.fft.backward(u) 

109 tmp[:] = self._eval_explicit_part(tmp, t, tmp) 

110 f.expl[:] = self.fft.forward(tmp) 

111 

112 else: 

113 

114 if self.eps > 0: 

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

116 

117 self.work_counters['rhs']() 

118 return f 

119 

120 def u_exact(self, t, **kwargs): 

121 r""" 

122 Routine to compute the exact solution at time :math:`t`. 

123 

124 Parameters 

125 ---------- 

126 t : float 

127 Time of the exact solution. 

128 

129 Returns 

130 ------- 

131 me : dtype_u 

132 The exact solution. 

133 """ 

134 

135 assert t == 0, 'ERROR: u_exact only valid for t=0' 

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

137 if self.init_type == 'circle': 

138 r2 = self.X[0] ** 2 + self.X[1] ** 2 

139 if self.spectral: 

140 tmp = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps))) 

141 me[:] = self.fft.forward(tmp) 

142 else: 

143 me[:] = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps))) 

144 elif self.init_type == 'circle_rand': 

145 ndim = len(me.shape) 

146 L = int(self.L[0]) 

147 # get random radii for circles/spheres 

148 self.xp.random.seed(1) 

149 lbound = 3.0 * self.eps 

150 ubound = 0.5 - self.eps 

151 rand_radii = (ubound - lbound) * self.xp.random.random_sample(size=tuple([L] * ndim)) + lbound 

152 # distribute circles/spheres 

153 tmp = newDistArray(self.fft, False) 

154 if ndim == 2: 

155 for i in range(0, L): 

156 for j in range(0, L): 

157 # build radius 

158 r2 = (self.X[0] + i - L / 2 + 0.5) ** 2 + (self.X[1] + j - L / 2 + 0.5) ** 2 

159 # add this blob, shifted by 1 to avoid issues with adding up negative contributions 

160 tmp += self.xp.tanh((rand_radii[i, j] - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)) + 1 

161 else: 

162 raise NotImplementedError 

163 # normalize to [0,1] 

164 tmp *= 0.5 

165 assert self.xp.all(tmp <= 1.0) 

166 if self.spectral: 

167 me[:] = self.fft.forward(tmp) 

168 else: 

169 self.xp.copyto(me, tmp) 

170 else: 

171 raise NotImplementedError('type of initial value not implemented, got %s' % self.init_type) 

172 

173 return me 

174 

175 

176class allencahn_imex_timeforcing(allencahn_imex): 

177 r""" 

178 Periodic Allen-Cahn equation with mpi4py-fft FFTs, IMEX, with a time-dependent, mass-conserving driving force. 

179 

180 Example implementing the :math:`N`-dimensional Allen-Cahn equation with periodic boundary conditions, with the two 

181 phases at :math:`u = 0` and :math:`u = 1` 

182 using time-dependent forcing 

183 

184 .. math:: 

185 \frac{\partial u}{\partial t} = \Delta u - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

186 - 6 d_w u (1 - u) 

187 

188 on a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`, driving force :math:`d_w`, and :math:`N=2,3`. Different initial 

189 conditions can be used, for example, circles of the form 

190 

191 .. math:: 

192 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}} 

193 {\sqrt{2}\varepsilon}\right)\right), 

194 

195 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is treated 

196 *semi-implicitly*, i.e., the linear part is solved with Fast-Fourier Transform (FFT) using ``mpi4py-fft`` [1]_ and the nonlinear part in the right-hand 

197 side will be treated explicitly. 

198 """ 

199 

200 def eval_f(self, u, t): 

201 """ 

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

203 

204 Parameters 

205 ---------- 

206 u : dtype_u 

207 Current values of the numerical solution. 

208 t : float 

209 Current time of the numerical solution is computed. 

210 

211 Returns 

212 ------- 

213 f : dtype_f 

214 The right-hand side of the problem. 

215 """ 

216 

217 f = self.dtype_f(self.init) 

218 

219 f.impl[:] = self._eval_Laplacian(u, f.impl) 

220 

221 if self.spectral: 

222 

223 tmp = newDistArray(self.fft, False) 

224 tmp[:] = self.fft.backward(u, tmp) 

225 

226 if self.eps > 0: 

227 tmpf = -2.0 / self.eps**2 * tmp * (1.0 - tmp) * (1.0 - 2.0 * tmp) 

228 else: 

229 tmpf = self.dtype_f(self.init, val=0.0) 

230 

231 # build sum over RHS without driving force 

232 Rt_local = float(self.xp.sum(self.fft.backward(f.impl) + tmpf)) 

233 if self.comm is not None: 

234 Rt_global = self.comm.allreduce(sendobj=Rt_local, op=MPI.SUM) 

235 else: 

236 Rt_global = Rt_local 

237 

238 # build sum over driving force term 

239 Ht_local = float(self.xp.sum(6.0 * tmp * (1.0 - tmp))) 

240 if self.comm is not None: 

241 Ht_global = self.comm.allreduce(sendobj=Ht_local, op=MPI.SUM) 

242 else: 

243 Ht_global = Ht_local 

244 

245 # add/substract time-dependent driving force 

246 if Ht_global != 0.0: 

247 dw = Rt_global / Ht_global 

248 else: 

249 dw = 0.0 

250 

251 tmpf -= 6.0 * dw * tmp * (1.0 - tmp) 

252 f.expl[:] = self.fft.forward(tmpf) 

253 

254 else: 

255 

256 if self.eps > 0: 

257 f.expl = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u) 

258 

259 # build sum over RHS without driving force 

260 Rt_local = float(self.xp.sum(f.impl + f.expl)) 

261 if self.comm is not None: 

262 Rt_global = self.comm.allreduce(sendobj=Rt_local, op=MPI.SUM) 

263 else: 

264 Rt_global = Rt_local 

265 

266 # build sum over driving force term 

267 Ht_local = float(self.xp.sum(6.0 * u * (1.0 - u))) 

268 if self.comm is not None: 

269 Ht_global = self.comm.allreduce(sendobj=Ht_local, op=MPI.SUM) 

270 else: 

271 Ht_global = Ht_local 

272 

273 # add/substract time-dependent driving force 

274 if Ht_global != 0.0: 

275 dw = Rt_global / Ht_global 

276 else: 

277 dw = 0.0 

278 

279 f.expl -= 6.0 * dw * u * (1.0 - u) 

280 

281 self.work_counters['rhs']() 

282 return f