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

## 24 statements

, created at 2024-09-09 14:59 +0000

1import cupy as cp

2import numpy as np

4from pySDC.core.errors import ProblemError

5from pySDC.core.problem import Problem

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

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

10 r"""

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

13 .. math::

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

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

19 .. math::

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

22 or *checker-board*

24 .. math::

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

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.

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.

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

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.

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.

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 """

61 dtype_u = cupy_mesh

62 dtype_f = imex_cupy_mesh

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

65 """Initialization routine"""

67 if nvars is None:

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

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

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(

82 )

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

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

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

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)

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

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

99 def eval_f(self, u, t):

100 """

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

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.

110 Returns

111 -------

112 f : dtype_f

113 The right-hand side of the problem.

114 """

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

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

125 """

126 Simple FFT solver for the diffusion part.

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

139 Parameters

140 ----------

141 me : dtype_u

142 The solution as mesh.

143 """

145 me = self.dtype_u(self.init)

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

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

149 return me

151 def u_exact(self, t):

152 r"""

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

155 Parameters

156 ----------

157 t : float

158 Time of the exact solution.

160 Returns

161 -------

162 me : dtype_u

163 The exact solution.

164 """

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)

179 return me

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

187 .. math::

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

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

193 .. math::

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

196 or *checker-board*

198 .. math::

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

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.

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.

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

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.

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.

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 """

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

236 """Initialization routine"""

238 if nvars is None:

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

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

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

244 def eval_f(self, u, t):

245 """

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

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.

255 Returns

256 -------

257 f : dtype_f

258 The right-hand side of the problem.

259 """

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

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

270 """

271 Simple FFT solver for the diffusion part.

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

284 Returns

285 -------

286 me : dtype_u

287 The solution as mesh.

288 """

290 me = self.dtype_u(self.init)

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

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

295 return me