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

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 Example implementing the :math:`N`-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [0, 1]^2` 

11 

12 .. math:: 

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

14 - 6 d_w u (1 - u) 

15 

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

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

18 

19 .. math:: 

20 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{(x_i-0.5)^2 + (y_j-0.5)^2}}{\sqrt{2}\varepsilon}\right), 

21 

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

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

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

25 

26 Parameters 

27 ---------- 

28 nvars : List of int tuples, optional 

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

30 eps : float, optional 

31 Scaling parameter :math:`\varepsilon`. 

32 radius : float, optional 

33 Radius of the circles. 

34 spectral : bool, optional 

35 Indicates if spectral initial condition is used. 

36 dw : float, optional 

37 Driving force. 

38 L : float, optional 

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

40 init_type : str, optional 

41 Initialises type of initial state. 

42 comm : bool, optional 

43 Communicator for parallelization. 

44 

45 Attributes 

46 ---------- 

47 fft : fft object 

48 Object for FFT. 

49 X : np.ogrid 

50 Grid coordinates in real space. 

51 K2 : np.1darray 

52 Laplace operator in spectral space. 

53 dx : float 

54 Mesh width in x direction. 

55 dy : float 

56 Mesh width in y direction. 

57 

58 References 

59 ---------- 

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

61 """ 

62 

63 def __init__( 

64 self, 

65 eps=0.04, 

66 radius=0.25, 

67 dw=0.0, 

68 init_type='circle', 

69 **kwargs, 

70 ): 

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

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

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

74 

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

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

77 return f_expl 

78 

79 def eval_f(self, u, t): 

80 """ 

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

82 

83 Parameters 

84 ---------- 

85 u : dtype_u 

86 Current values of the numerical solution. 

87 t : float 

88 Current time of the numerical solution is computed. 

89 

90 Returns 

91 ------- 

92 f : dtype_f 

93 The right-hand side of the problem. 

94 """ 

95 

96 f = self.dtype_f(self.init) 

97 

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

99 

100 if self.spectral: 

101 f.impl = -self.K2 * u 

102 

103 if self.eps > 0: 

104 tmp = self.fft.backward(u) 

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

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

107 

108 else: 

109 

110 if self.eps > 0: 

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

112 

113 self.work_counters['rhs']() 

114 return f 

115 

116 def u_exact(self, t): 

117 r""" 

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

119 

120 Parameters 

121 ---------- 

122 t : float 

123 Time of the exact solution. 

124 

125 Returns 

126 ------- 

127 me : dtype_u 

128 The exact solution. 

129 """ 

130 

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

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

133 if self.init_type == 'circle': 

134 r2 = (self.X[0] - 0.5) ** 2 + (self.X[1] - 0.5) ** 2 

135 if self.spectral: 

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

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

138 else: 

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

140 elif self.init_type == 'circle_rand': 

141 ndim = len(me.shape) 

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

143 # get random radii for circles/spheres 

144 self.xp.random.seed(1) 

145 lbound = 3.0 * self.eps 

146 ubound = 0.5 - self.eps 

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

148 # distribute circles/spheres 

149 tmp = newDistArray(self.fft, False) 

150 if ndim == 2: 

151 for i in range(0, L): 

152 for j in range(0, L): 

153 # build radius 

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

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

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

157 # normalize to [0,1] 

158 tmp *= 0.5 

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

160 if self.spectral: 

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

162 else: 

163 self.xp.copyto(me, tmp) 

164 else: 

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

166 

167 return me 

168 

169 

170class allencahn_imex_timeforcing(allencahn_imex): 

171 r""" 

172 Example implementing the :math:`N`-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [0, 1]^2` 

173 using time-dependent forcing 

174 

175 .. math:: 

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

177 - 6 d_w u (1 - u) 

178 

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

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

181 

182 .. math:: 

183 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{(x_i-0.5)^2 + (y_j-0.5)^2}}{\sqrt{2}\varepsilon}\right), 

184 

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

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

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

188 """ 

189 

190 def eval_f(self, u, t): 

191 """ 

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

193 

194 Parameters 

195 ---------- 

196 u : dtype_u 

197 Current values of the numerical solution. 

198 t : float 

199 Current time of the numerical solution is computed. 

200 

201 Returns 

202 ------- 

203 f : dtype_f 

204 The right-hand side of the problem. 

205 """ 

206 

207 f = self.dtype_f(self.init) 

208 

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

210 

211 if self.spectral: 

212 

213 tmp = newDistArray(self.fft, False) 

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

215 

216 if self.eps > 0: 

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

218 else: 

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

220 

221 # build sum over RHS without driving force 

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

223 if self.comm is not None: 

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

225 else: 

226 Rt_global = Rt_local 

227 

228 # build sum over driving force term 

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

230 if self.comm is not None: 

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

232 else: 

233 Ht_global = Rt_local 

234 

235 # add/substract time-dependent driving force 

236 if Ht_global != 0.0: 

237 dw = Rt_global / Ht_global 

238 else: 

239 dw = 0.0 

240 

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

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

243 

244 else: 

245 

246 if self.eps > 0: 

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

248 

249 # build sum over RHS without driving force 

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

251 if self.comm is not None: 

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

253 else: 

254 Rt_global = Rt_local 

255 

256 # build sum over driving force term 

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

258 if self.comm is not None: 

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

260 else: 

261 Ht_global = Rt_local 

262 

263 # add/substract time-dependent driving force 

264 if Ht_global != 0.0: 

265 dw = Rt_global / Ht_global 

266 else: 

267 dw = 0.0 

268 

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

270 

271 self.work_counters['rhs']() 

272 return f