Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py: 88%
85 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
3from pySDC.core.errors import ProblemError
4from pySDC.core.problem import Problem, WorkCounter
5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
8# noinspection PyUnusedLocal
9class allencahn2d_imex(Problem):
10 r"""
11 2D periodic Allen-Cahn equation with FFTs, IMEX with the Laplacian implicit and the reaction explicit.
13 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
14 phases at :math:`u = 0` and :math:`u = 1`
16 .. math::
17 \frac{\partial u}{\partial t} = \Delta u
18 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
20 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
21 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
23 On a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`. Different initial conditions
24 can be used, for example, circles of the form
26 .. math::
27 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
28 {\sqrt{2}\varepsilon}\right)\right),
30 or *checker-board*
32 .. math::
33 u({\bf x}, 0) = \frac{1}{2}\left(1 + \sin(2 \pi x_i) \sin(2 \pi y_j)\right),
35 or uniform distributed random numbers in :math:`[0, 1]` for :math:`i, j=0,..,N-1`, where :math:`N` is the number of
36 spatial grid points. For time-stepping, the problem is treated *semi-implicitly*, i.e., the diffusion part is solved by
37 Fast-Fourier Transform (FFT) and the nonlinear term is treated explicitly.
39 An exact solution is not known, but instead the numerical solution can be compared via a generated reference solution computed
40 by a ``SciPy`` routine.
42 Parameters
43 ----------
44 nvars : tuple of int, optional
45 Number of unknowns in each spatial direction, e.g. ``nvars=(128, 128)``.
46 nu : int, optional
47 Exponent of the double well; :math:`\nu = 2` is the standard Allen-Cahn nonlinearity.
48 eps : float, optional
49 Scaling parameter :math:`\varepsilon`.
50 radius : float, optional
51 Radius of the circles.
52 L : float, optional
53 Denotes the period of the function to be approximated for the Fourier transform.
54 init_type : str, optional
55 Indicates which type of initial condition is used.
56 useGPU : bool, optional
57 Run on the GPU with CuPy instead of on the CPU with NumPy.
59 Attributes
60 ----------
61 xvalues : np.1darray
62 Grid points in space.
63 dx : float
64 Mesh width.
65 lap : np.1darray
66 Spectral operator for Laplacian.
67 work_counters : WorkCounter
68 Counts the right-hand side evaluations.
69 """
71 dtype_u = mesh
72 dtype_f = imex_mesh
74 xp = np
76 def setup_GPU(self):
77 """
78 Switch the array module and the datatypes over to CuPy.
80 This changes the class, not the instance, as everything else in pySDC that does this
81 does: once one instance of a class runs on the GPU, they all do.
82 """
83 import cupy as cp
84 from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh
86 self.xp = cp
87 self.dtype_u = cupy_mesh
88 self.dtype_f = imex_cupy_mesh
90 def __init__(
91 self,
92 nvars=None,
93 nu=2,
94 eps=0.04,
95 radius=0.25,
96 L=1.0,
97 init_type='circle',
98 useGPU=False,
99 ):
100 """Initialization routine"""
102 if useGPU:
103 self.setup_GPU()
105 if nvars is None:
106 nvars = (128, 128)
108 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
109 if len(nvars) != 2:
110 raise ProblemError('this is a 2d example, got %s' % nvars)
111 if nvars[0] != nvars[1]:
112 raise ProblemError('need a square domain, got %s' % nvars)
113 if nvars[0] % 2 != 0:
114 raise ProblemError('the setup requires nvars = 2^p per dimension')
116 # invoke super init, passing number of dofs, dtype_u and dtype_f
117 super().__init__(init=(nvars, None, np.dtype('float64')))
118 self._makeAttributeAndRegister(
119 'nvars', 'nu', 'eps', 'radius', 'L', 'init_type', 'useGPU', localVars=locals(), readOnly=True
120 )
122 self.dx = self.L / self.nvars[0] # could be useful for hooks, too.
123 self.xvalues = self.xp.arange(self.nvars[0]) * self.dx - self.L / 2.0
125 kx = self.xp.zeros(self.init[0][0])
126 ky = self.xp.zeros(self.init[0][1] // 2 + 1)
128 kx[: int(self.init[0][0] / 2) + 1] = 2 * np.pi / self.L * self.xp.arange(0, int(self.init[0][0] / 2) + 1)
129 kx[int(self.init[0][0] / 2) + 1 :] = (
130 2 * np.pi / self.L * self.xp.arange(int(self.init[0][0] / 2) + 1 - self.init[0][0], 0)
131 )
132 ky[:] = 2 * np.pi / self.L * self.xp.arange(0, self.init[0][1] // 2 + 1)
134 xv, yv = self.xp.meshgrid(kx, ky, indexing='ij')
135 self.lap = -(xv**2) - yv**2
137 self.work_counters['rhs'] = WorkCounter()
139 def reaction(self, u):
140 r"""
141 The reaction term, :math:`\frac{1}{2\varepsilon^2}(2u - 1)\left(1 - (2u - 1)^\nu\right)`.
143 The wells sit at :math:`u = 0` and :math:`u = 1`, so the double well is symmetric about
144 :math:`2u - 1`; writing the term in that variable is what lets :math:`\nu` keep the meaning
145 it has always had here. For the default :math:`\nu = 2` this is
146 :math:`-\frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
147 """
148 v = 2.0 * u - 1.0
149 return 0.5 / self.eps**2 * v * (1.0 - v**self.nu)
151 def eval_f(self, u, t):
152 """
153 Routine to evaluate the right-hand side of the problem.
155 Parameters
156 ----------
157 u : dtype_u
158 Current values of the numerical solution.
159 t : float
160 Current time of the numerical solution is computed.
162 Returns
163 -------
164 f : dtype_f
165 The right-hand side of the problem.
166 """
168 f = self.dtype_f(self.init)
169 tmp = self.lap * self.xp.fft.rfft2(u)
170 f.impl[:] = self.xp.fft.irfft2(tmp)
171 if self.eps > 0:
172 f.expl[:] = self.reaction(u)
174 self.work_counters['rhs']()
175 return f
177 def solve_system(self, rhs, factor, u0, t):
178 """
179 Simple FFT solver for the diffusion part.
181 Parameters
182 ----------
183 rhs : dtype_f
184 Right-hand side for the linear system.
185 factor : float
186 Abbrev. for the node-to-node stepsize (or any other factor required).
187 u0 : dtype_u
188 Initial guess for the iterative solver (not used here so far).
189 t : float
190 Current time (e.g. for time-dependent BCs).
192 Returns
193 -------
194 me : dtype_u
195 The solution as mesh.
196 """
198 me = self.dtype_u(self.init)
200 tmp = self.xp.fft.rfft2(rhs) / (1.0 - factor * self.lap)
201 me[:] = self.xp.fft.irfft2(tmp)
203 return me
205 def u_exact(self, t, u_init=None, t_init=None):
206 r"""
207 Routine to compute the exact solution at time :math:`t`.
209 Parameters
210 ----------
211 t : float
212 Time of the exact solution.
213 u_init : pySDC.implementations.problem_classes.allencahn2d_imex.dtype_u
214 Initial conditions for getting the exact solution.
215 t_init : float
216 The starting time.
218 Returns
219 -------
220 me : dtype_u
221 The exact solution.
222 """
224 me = self.dtype_u(self.init, val=0.0)
226 if t == 0:
227 if self.init_type == 'circle':
228 xv, yv = self.xp.meshgrid(self.xvalues, self.xvalues, indexing='ij')
229 me[:, :] = 0.5 * (
230 1.0 + self.xp.tanh((self.radius - self.xp.sqrt(xv**2 + yv**2)) / (np.sqrt(2) * self.eps))
231 )
232 elif self.init_type == 'checkerboard':
233 xv, yv = self.xp.meshgrid(self.xvalues, self.xvalues)
234 me[:, :] = 0.5 * (1.0 + self.xp.sin(2.0 * np.pi * xv) * self.xp.sin(2.0 * np.pi * yv))
235 elif self.init_type == 'random':
236 me[:, :] = self.xp.random.uniform(0, 1, self.init)
237 else:
238 raise NotImplementedError('type of initial value not implemented, got %s' % self.init_type)
239 else:
241 def eval_rhs(t, u):
242 f = self.eval_f(u.reshape(self.init[0]), t)
243 return (f.impl + f.expl).flatten()
245 me[:, :] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
247 return me
250class allencahn2d_imex_stab(allencahn2d_imex):
251 r"""
252 2D periodic Allen-Cahn equation with FFTs, IMEX with a stabilizing linear term shifted into the implicit part.
254 This implements the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
255 phases at :math:`u = 0` and :math:`u = 1`
257 .. math::
258 \frac{\partial u}{\partial t} = \Delta u
259 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
261 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
262 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
264 On a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`. Different initial conditions
265 can be used here, for example, circles of the form
267 .. math::
268 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
269 {\sqrt{2}\varepsilon}\right)\right),
271 or *checker-board*
273 .. math::
274 u({\bf x}, 0) = \frac{1}{2}\left(1 + \sin(2 \pi x_i) \sin(2 \pi y_j)\right),
276 or uniform distributed random numbers in :math:`[0, 1]` for :math:`i, j=0,..,N-1`, where :math:`N` is the number of
277 spatial grid points. For time-stepping, the problem is treated *semi-implicitly*, i.e., the diffusion part is solved with
278 Fast-Fourier Transform (FFT) and the nonlinear term is treated explicitly.
280 An exact solution is not known, but instead the numerical solution can be compared via a generated reference solution computed
281 by a ``SciPy`` routine.
283 Same equation as :class:`allencahn2d_imex`, split differently: :math:`-\frac{2}{\varepsilon^2}`
284 is folded into the implicit operator and :math:`\frac{2}{\varepsilon^2}u` added back to the
285 explicit part, which leaves the sum unchanged but flattens the explicit part at both wells,
286 since the reaction has derivative :math:`-\frac{2}{\varepsilon^2}` at each of them.
289 Parameters
290 ----------
291 nvars : tuple of int, optional
292 Number of unknowns in each spatial direction, e.g. ``nvars=(128, 128)``.
293 nu : int, optional
294 Exponent of the double well; :math:`\nu = 2` is the standard Allen-Cahn nonlinearity.
295 eps : float, optional
296 Scaling parameter :math:`\varepsilon`.
297 radius : float, optional
298 Radius of the circles.
299 L : float, optional
300 Denotes the period of the function to be approximated for the Fourier transform.
301 init_type : str, optional
302 Indicates which type of initial condition is used.
303 useGPU : bool, optional
304 Run on the GPU with CuPy instead of on the CPU with NumPy.
306 Attributes
307 ----------
308 xvalues : np.1darray
309 Grid points in space.
310 dx : float
311 Mesh width.
312 lap : np.1darray
313 Spectral operator for Laplacian.
314 work_counters : WorkCounter
315 Counts the right-hand side evaluations.
316 """
318 def __init__(self, nvars=None, nu=2, eps=0.04, radius=0.25, L=1.0, init_type='circle', useGPU=False):
319 """Initialization routine"""
321 if nvars is None:
322 nvars = (256, 256)
324 super().__init__(nvars, nu, eps, radius, L, init_type, useGPU)
325 self.lap -= 2.0 / self.eps**2
327 def eval_f(self, u, t):
328 """
329 Routine to evaluate the right-hand side of the problem.
331 Parameters
332 ----------
333 u : dtype_u
334 Current values of the numerical solution.
335 t : float
336 Current time of the numerical solution is computed.
338 Returns
339 -------
340 f : dtype_f
341 The right-hand side of the problem.
342 """
344 f = self.dtype_f(self.init)
345 tmp = self.lap * self.xp.fft.rfft2(u)
346 f.impl[:] = self.xp.fft.irfft2(tmp)
347 if self.eps > 0:
348 f.expl[:] = self.reaction(u) + 2.0 / self.eps**2 * u
350 self.work_counters['rhs']()
351 return f
353 def solve_system(self, rhs, factor, u0, t):
354 """
355 Simple FFT solver for the diffusion part.
357 Parameters
358 ----------
359 rhs : dtype_f
360 Right-hand side for the linear system.
361 factor : float
362 Abbrev. for the node-to-node stepsize (or any other factor required).
363 u0 : dtype_u
364 Initial guess for the iterative solver (not used here so far).
365 t : float
366 Current time (e.g. for time-dependent BCs).
368 Returns
369 -------
370 me : dtype_u
371 The solution as mesh.
372 """
374 me = self.dtype_u(self.init)
376 tmp = self.xp.fft.rfft2(rhs) / (1.0 - factor * self.lap)
377 me[:] = self.xp.fft.irfft2(tmp)
379 return me