Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FFT_gpu.py: 0%

1import cupy as cp 

2import numpy as np 

3 

4from pySDC.core.Errors import ProblemError 

5from pySDC.core.Problem import ptype 

6from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh 

7 

8 

9class allencahn2d_imex(ptype): # pragma: no cover 

10 r""" 

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

12 

13 .. math:: 

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

15 

16 on a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`, and constant parameter :math:`\nu`. Different initial conditions 

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

18 

19 .. math:: 

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

21 

22 or *checker-board* 

23 

24 .. math:: 

25 u({\bf x}, 0) = \sin(2 \pi x_i) \sin(2 \pi y_j), 

26 

27 or uniform distributed random numbers in :math:`[-1, 1]` for :math:`i, j=0,..,N-1`, where :math:`N` is the number of 

28 spatial grid points. For time-stepping, the problem is treated *semi-implicitly*, i.e., the diffusion part is solved with 

29 Fast-Fourier Tranform (FFT) and the nonlinear term is treated explicitly. 

30 

31 An exact solution is not known, but instead the numerical solution can be compared via a generated reference solution computed 

32 by a ``SciPy`` routine. 

33 

34 This class is especially developed for solving it on GPUs using ``CuPy``. 

35 

36 Parameters 

37 ---------- 

38 nvars : List of int tuples, optional 

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

40 nu : float, optional 

41 Problem parameter :math:`\nu`. 

42 eps : float, optional 

43 Scaling parameter :math:`\varepsilon`. 

44 radius : float, optional 

45 Radius of the circles. 

46 L : float, optional 

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

48 init_type : str, optional 

49 Indicates which type of initial condition is used. 

50 

51 Attributes 

52 ---------- 

53 xvalues : cp.1darray 

54 Grid points in space. 

55 dx : float 

56 Cupy mesh width. 

57 lap : cp.1darray 

58 Spectral operator for Laplacian. 

59 """ 

60 

61 dtype_u = cupy_mesh 

62 dtype_f = imex_cupy_mesh 

63 

64 def __init__(self, nvars=None, nu=2, eps=0.04, radius=0.25, L=1.0, init_type='circle'): 

65 """Initialization routine""" 

66 

67 if nvars is None: 

68 nvars = [(256, 256), (64, 64)] 

69 

70 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on 

71 if len(nvars) != 2: 

72 raise ProblemError('this is a 2d example, got %s' % nvars) 

73 if nvars[0] != nvars[1]: 

74 raise ProblemError('need a square domain, got %s' % nvars) 

75 if nvars[0] % 2 != 0: 

76 raise ProblemError('the setup requires nvars = 2^p per dimension') 

77 

78 # invoke super init, passing number of dofs, dtype_u and dtype_f 

79 super().__init__(init=(nvars, None, cp.dtype('float64'))) 

80 self._makeAttributeAndRegister( 

81 'nvars', 'nu', 'eps', 'radius', 'L', 'init_type', localVars=locals(), readOnly=True 

82 ) 

83 

84 self.dx = self.L / self.nvars[0] # could be useful for hooks, too. 

85 self.xvalues = cp.array([i * self.dx - self.L / 2.0 for i in range(self.nvars[0])]) 

86 

87 kx = cp.zeros(self.init[0][0]) 

88 ky = cp.zeros(self.init[0][1] // 2 + 1) 

89 

90 kx[: int(self.init[0][0] / 2) + 1] = 2 * np.pi / self.L * cp.arange(0, int(self.init[0][0] / 2) + 1) 

91 kx[int(self.init[0][0] / 2) + 1 :] = ( 

92 2 * np.pi / self.L * cp.arange(int(self.init[0][0] / 2) + 1 - self.init[0][0], 0) 

93 ) 

94 ky[:] = 2 * np.pi / self.L * cp.arange(0, self.init[0][1] // 2 + 1) 

95 

96 xv, yv = cp.meshgrid(kx, ky, indexing='ij') 

97 self.lap = -(xv**2) - yv**2 

98 

99 def eval_f(self, u, t): 

100 """ 

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

102 

103 Parameters 

104 ---------- 

105 u : dtype_u 

106 Current values of the numerical solution. 

107 t : float 

108 Current time of the numerical solution is computed. 

109 

110 Returns 

111 ------- 

112 f : dtype_f 

113 The right-hand side of the problem. 

114 """ 

115 

116 f = self.dtype_f(self.init) 

117 v = u.flatten() 

118 tmp = self.lap * cp.fft.rfft2(u) 

119 f.impl[:] = cp.fft.irfft2(tmp) 

120 if self.eps > 0: 

121 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars) 

122 return f 

123 

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

125 """ 

126 Simple FFT solver for the diffusion part. 

127 

128 Parameters 

129 ---------- 

130 rhs : dtype_f 

131 Right-hand side for the linear system. 

132 factor : float 

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

134 u0 : dtype_u 

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

136 t : float 

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

138 

139 Parameters 

140 ---------- 

141 me : dtype_u 

142 The solution as mesh. 

143 """ 

144 

145 me = self.dtype_u(self.init) 

146 tmp = cp.fft.rfft2(rhs) / (1.0 - factor * self.lap) 

147 me[:] = cp.fft.irfft2(tmp) 

148 

149 return me 

150 

151 def u_exact(self, t): 

152 r""" 

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

154 

155 Parameters 

156 ---------- 

157 t : float 

158 Time of the exact solution. 

159 

160 Returns 

161 ------- 

162 me : dtype_u 

163 The exact solution. 

164 """ 

165 

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

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

168 if self.init_type == 'circle': 

169 xv, yv = cp.meshgrid(self.xvalues, self.xvalues, indexing='ij') 

170 me[:, :] = cp.tanh((self.radius - cp.sqrt(xv**2 + yv**2)) / (cp.sqrt(2) * self.eps)) 

171 elif self.init_type == 'checkerboard': 

172 xv, yv = cp.meshgrid(self.xvalues, self.xvalues) 

173 me[:, :] = cp.sin(2.0 * np.pi * xv) * cp.sin(2.0 * np.pi * yv) 

174 elif self.init_type == 'random': 

175 me[:, :] = cp.random.uniform(-1, 1, self.init) 

176 else: 

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

178 

179 return me 

180 

181 

182class allencahn2d_imex_stab(allencahn2d_imex): 

183 r""" 

184 This implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

185 with stabilized splitting 

186 

187 .. math:: 

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

189 

190 on a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`, and constant parameter :math:`\nu`. Different initial conditions 

191 can be used here, for example, circles of the form 

192 

193 .. math:: 

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

195 

196 or *checker-board* 

197 

198 .. math:: 

199 u({\bf x}, 0) = \sin(2 \pi x_i) \sin(2 \pi y_j), 

200 

201 or uniform distributed random numbers in :math:`[-1, 1]` for :math:`i, j=0,..,N-1`, where :math:`N` is the number of 

202 spatial grid points. For time-stepping, the problem is treated *semi-implicitly*, i.e., the diffusion part is solved with 

203 Fast-Fourier Tranform (FFT) and the nonlinear term is treated explicitly. 

204 

205 An exact solution is not known, but instead the numerical solution can be compared via a generated reference solution computed 

206 by a ``SciPy`` routine. 

207 

208 This class is especially developed for solving it on GPUs using ``CuPy``. 

209 

210 Parameters 

211 ---------- 

212 nvars : List of int tuples, optional 

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

214 nu : float, optional 

215 Problem parameter :math:`\nu`. 

216 eps : float, optional 

217 Scaling parameter :math:`\varepsilon`. 

218 radius : float, optional 

219 Radius of the circles. 

220 L : float, optional 

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

222 init_type : str, optional 

223 Indicates which type of initial condition is used. 

224 

225 Attributes 

226 ---------- 

227 xvalues : cp.1darray 

228 Grid points in space. 

229 dx : float 

230 Cupy mesh width. 

231 lap : cp.1darray 

232 Spectral operator for Laplacian. 

233 """ 

234 

235 def __init__(self, nvars=None, nu=2, eps=0.04, radius=0.25, L=1.0, init_type='circle'): 

236 """Initialization routine""" 

237 

238 if nvars is None: 

239 nvars = [(256, 256), (64, 64)] 

240 

241 super().__init__(nvars, nu, eps, radius, L, init_type) 

242 self.lap -= 2.0 / self.eps**2 

243 

244 def eval_f(self, u, t): 

245 """ 

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

247 

248 Parameters 

249 ---------- 

250 u : dtype_u 

251 Current values of the numerical solution. 

252 t : float 

253 Current time of the numerical solution is computed. 

254 

255 Returns 

256 ------- 

257 f : dtype_f 

258 The right-hand side of the problem. 

259 """ 

260 

261 f = self.dtype_f(self.init) 

262 v = u.flatten() 

263 tmp = self.lap * cp.fft.rfft2(u) 

264 f.impl[:] = cp.fft.irfft2(tmp) 

265 if self.eps > 0: 

266 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu) + 2.0 / self.eps**2 * v).reshape(self.nvars) 

267 return f 

268 

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

270 """ 

271 Simple FFT solver for the diffusion part. 

272 

273 Parameters 

274 ---------- 

275 rhs : dtype_f 

276 Right-hand side for the linear system. 

277 factor : float 

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

279 u0 : dtype_u 

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

281 t : float 

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

283 

284 Returns 

285 ------- 

286 me : dtype_u 

287 The solution as mesh. 

288 """ 

289 

290 me = self.dtype_u(self.init) 

291 

292 tmp = cp.fft.rfft2(rhs) / (1.0 - factor * self.lap) 

293 me[:] = cp.fft.irfft2(tmp) 

294 

295 return me