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:`2`-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 Transform (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, **kwargs): 

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] - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)) + 1 

157 else: 

158 raise NotImplementedError 

159 # normalize to [0,1] 

160 tmp *= 0.5 

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

162 if self.spectral: 

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

164 else: 

165 self.xp.copyto(me, tmp) 

166 else: 

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

168 

169 return me 

170 

171 

172class allencahn_imex_timeforcing(allencahn_imex): 

173 r""" 

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

175 using time-dependent forcing 

176 

177 .. math:: 

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

179 - 6 d_w u (1 - u) 

180 

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

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

183 

184 .. math:: 

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

186 

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

188 *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 

189 side will be treated explicitly. 

190 """ 

191 

192 def eval_f(self, u, t): 

193 """ 

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

195 

196 Parameters 

197 ---------- 

198 u : dtype_u 

199 Current values of the numerical solution. 

200 t : float 

201 Current time of the numerical solution is computed. 

202 

203 Returns 

204 ------- 

205 f : dtype_f 

206 The right-hand side of the problem. 

207 """ 

208 

209 f = self.dtype_f(self.init) 

210 

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

212 

213 if self.spectral: 

214 

215 tmp = newDistArray(self.fft, False) 

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

217 

218 if self.eps > 0: 

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

220 else: 

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

222 

223 # build sum over RHS without driving force 

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

225 if self.comm is not None: 

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

227 else: 

228 Rt_global = Rt_local 

229 

230 # build sum over driving force term 

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

232 if self.comm is not None: 

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

234 else: 

235 Ht_global = Ht_local 

236 

237 # add/substract time-dependent driving force 

238 if Ht_global != 0.0: 

239 dw = Rt_global / Ht_global 

240 else: 

241 dw = 0.0 

242 

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

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

245 

246 else: 

247 

248 if self.eps > 0: 

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

250 

251 # build sum over RHS without driving force 

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

253 if self.comm is not None: 

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

255 else: 

256 Rt_global = Rt_local 

257 

258 # build sum over driving force term 

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

260 if self.comm is not None: 

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

262 else: 

263 Ht_global = Ht_local 

264 

265 # add/substract time-dependent driving force 

266 if Ht_global != 0.0: 

267 dw = Rt_global / Ht_global 

268 else: 

269 dw = 0.0 

270 

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

272 

273 self.work_counters['rhs']() 

274 return f