# Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py: 21%

## 70 statements

, created at 2024-09-20 17:10 +0000

1import numpy as np

3from pySDC.core.errors import ProblemError

4from pySDC.core.problem import Problem

5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh

8# noinspection PyUnusedLocal

9class allencahn2d_imex(Problem):

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 by

29 Fast-Fourier Transform (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 Parameters

35 ----------

36 nvars : List of int tuples, optional

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

38 nu : float, optional

39 Problem parameter :math:\nu.

40 eps : float, optional

41 Scaling parameter :math:\varepsilon.

44 L : float, optional

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

46 init_type : str, optional

47 Indicates which type of initial condition is used.

49 Attributes

50 ----------

51 xvalues : np.1darray

52 Grid points in space.

53 dx : float

54 Mesh width.

55 lap : np.1darray

56 Spectral operator for Laplacian.

57 """

59 dtype_u = mesh

60 dtype_f = imex_mesh

62 def __init__(

63 self,

64 nvars=None,

65 nu=2,

66 eps=0.04,

68 L=1.0,

69 init_type='circle',

70 ):

71 """Initialization routine"""

73 if nvars is None:

74 nvars = (128, 128)

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

77 if len(nvars) != 2:

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

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

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

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

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

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

85 super().__init__(init=(nvars, None, np.dtype('float64')))

86 self._makeAttributeAndRegister(

88 )

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

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

93 kx = np.zeros(self.init[0][0])

94 ky = np.zeros(self.init[0][1] // 2 + 1)

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

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

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

99 )

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

102 xv, yv = np.meshgrid(kx, ky, indexing='ij')

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

105 def eval_f(self, u, t):

106 """

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

109 Parameters

110 ----------

111 u : dtype_u

112 Current values of the numerical solution.

113 t : float

114 Current time of the numerical solution is computed.

116 Returns

117 -------

118 f : dtype_f

119 The right-hand side of the problem.

120 """

122 f = self.dtype_f(self.init)

123 tmp = self.lap * np.fft.rfft2(u)

124 f.impl[:] = np.fft.irfft2(tmp)

125 if self.eps > 0:

126 f.expl[:] = 1.0 / self.eps**2 * u * (1.0 - u**self.nu)

127 return f

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

130 """

131 Simple FFT solver for the diffusion part.

133 Parameters

134 ----------

135 rhs : dtype_f

136 Right-hand side for the linear system.

137 factor : float

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

139 u0 : dtype_u

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

141 t : float

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

144 Returns

145 -------

146 me : dtype_u

147 The solution as mesh.

148 """

150 me = self.dtype_u(self.init)

152 tmp = np.fft.rfft2(rhs) / (1.0 - factor * self.lap)

153 me[:] = np.fft.irfft2(tmp)

155 return me

157 def u_exact(self, t, u_init=None, t_init=None):

158 r"""

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

161 Parameters

162 ----------

163 t : float

164 Time of the exact solution.

165 u_init : pySDC.implementations.problem_classes.allencahn2d_imex.dtype_u

166 Initial conditions for getting the exact solution.

167 t_init : float

168 The starting time.

170 Returns

171 -------

172 me : dtype_u

173 The exact solution.

174 """

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

178 if t == 0:

179 if self.init_type == 'circle':

180 xv, yv = np.meshgrid(self.xvalues, self.xvalues, indexing='ij')

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

182 elif self.init_type == 'checkerboard':

183 xv, yv = np.meshgrid(self.xvalues, self.xvalues)

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

185 elif self.init_type == 'random':

186 me[:, :] = np.random.uniform(-1, 1, self.init)

187 else:

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

189 else:

191 def eval_rhs(t, u):

192 f = self.eval_f(u.reshape(self.init[0]), t)

193 return (f.impl + f.expl).flatten()

195 me[:, :] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)

197 return me

200class allencahn2d_imex_stab(allencahn2d_imex):

201 r"""

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

203 with stabilized splitting

205 .. math::

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

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

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

211 .. math::

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

214 or *checker-board*

216 .. math::

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

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

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

221 Fast-Fourier Transform (FFT) and the nonlinear term is treated explicitly.

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

224 by a SciPy routine.

226 Parameters

227 ----------

228 nvars : List of int tuples, optional

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

230 nu : float, optional

231 Problem parameter :math:\nu.

232 eps : float, optional

233 Scaling parameter :math:\varepsilon.

236 L : float, optional

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

238 init_type : str, optional

239 Indicates which type of initial condition is used.

241 Attributes

242 ----------

243 xvalues : np.1darray

244 Grid points in space.

245 dx : float

246 Mesh width.

247 lap : np.1darray

248 Spectral operator for Laplacian.

249 """

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

252 """Initialization routine"""

254 if nvars is None:

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

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

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

260 def eval_f(self, u, t):

261 """

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

264 Parameters

265 ----------

266 u : dtype_u

267 Current values of the numerical solution.

268 t : float

269 Current time of the numerical solution is computed.

271 Returns

272 -------

273 f : dtype_f

274 The right-hand side of the problem.

275 """

277 f = self.dtype_f(self.init)

278 tmp = self.lap * np.fft.rfft2(u)

279 f.impl[:] = np.fft.irfft2(tmp)

280 if self.eps > 0:

281 f.expl[:] = 1.0 / self.eps**2 * u * (1.0 - u**self.nu) + 2.0 / self.eps**2 * u

282 return f

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

285 """

286 Simple FFT solver for the diffusion part.

288 Parameters

289 ----------

290 rhs : dtype_f

291 Right-hand side for the linear system.

292 factor : float

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

294 u0 : dtype_u

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

296 t : float

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

299 Returns

300 -------

301 me : dtype_u

302 The solution as mesh.

303 """

305 me = self.dtype_u(self.init)

307 tmp = np.fft.rfft2(rhs) / (1.0 - factor * self.lap)

308 me[:] = np.fft.irfft2(tmp)

310 return me