Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py: 21%
70 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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`.
42 radius : float, optional
43 Radius of the circles.
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,
67 radius=0.25,
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(
87 'nvars', 'nu', 'eps', 'radius', 'L', 'init_type', localVars=locals(), readOnly=True
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`.
234 radius : float, optional
235 Radius of the circles.
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