Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py: 21%
72 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
3from pySDC.core.Errors import ProblemError
4from pySDC.core.Problem import ptype
5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
8# noinspection PyUnusedLocal
9class allencahn2d_imex(ptype):
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 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 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 v = u.flatten()
124 tmp = self.lap * np.fft.rfft2(u)
125 f.impl[:] = np.fft.irfft2(tmp)
126 if self.eps > 0:
127 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
128 return f
130 def solve_system(self, rhs, factor, u0, t):
131 """
132 Simple FFT solver for the diffusion part.
134 Parameters
135 ----------
136 rhs : dtype_f
137 Right-hand side for the linear system.
138 factor : float
139 Abbrev. for the node-to-node stepsize (or any other factor required).
140 u0 : dtype_u
141 Initial guess for the iterative solver (not used here so far).
142 t : float
143 Current time (e.g. for time-dependent BCs).
145 Returns
146 -------
147 me : dtype_u
148 The solution as mesh.
149 """
151 me = self.dtype_u(self.init)
153 tmp = np.fft.rfft2(rhs) / (1.0 - factor * self.lap)
154 me[:] = np.fft.irfft2(tmp)
156 return me
158 def u_exact(self, t, u_init=None, t_init=None):
159 r"""
160 Routine to compute the exact solution at time :math:`t`.
162 Parameters
163 ----------
164 t : float
165 Time of the exact solution.
166 u_init : pySDC.implementations.problem_classes.allencahn2d_imex.dtype_u
167 Initial conditions for getting the exact solution.
168 t_init : float
169 The starting time.
171 Returns
172 -------
173 me : dtype_u
174 The exact solution.
175 """
177 me = self.dtype_u(self.init, val=0.0)
179 if t == 0:
180 if self.init_type == 'circle':
181 xv, yv = np.meshgrid(self.xvalues, self.xvalues, indexing='ij')
182 me[:, :] = np.tanh((self.radius - np.sqrt(xv**2 + yv**2)) / (np.sqrt(2) * self.eps))
183 elif self.init_type == 'checkerboard':
184 xv, yv = np.meshgrid(self.xvalues, self.xvalues)
185 me[:, :] = np.sin(2.0 * np.pi * xv) * np.sin(2.0 * np.pi * yv)
186 elif self.init_type == 'random':
187 me[:, :] = np.random.uniform(-1, 1, self.init)
188 else:
189 raise NotImplementedError('type of initial value not implemented, got %s' % self.init_type)
190 else:
192 def eval_rhs(t, u):
193 f = self.eval_f(u.reshape(self.init[0]), t)
194 return (f.impl + f.expl).flatten()
196 me[:, :] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
198 return me
201class allencahn2d_imex_stab(allencahn2d_imex):
202 r"""
203 This implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
204 with stabilized splitting
206 .. math::
207 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) + \frac{2}{\varepsilon^2}u
209 on a spatial domain :math:`[-\frac{L}{2}, \frac{L}{2}]^2`, and constant parameter :math:`\nu`. Different initial conditions
210 can be used here, for example, circles of the form
212 .. math::
213 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right),
215 or *checker-board*
217 .. math::
218 u({\bf x}, 0) = \sin(2 \pi x_i) \sin(2 \pi y_j),
220 or uniform distributed random numbers in :math:`[-1, 1]` for :math:`i, j=0,..,N-1`, where :math:`N` is the number of
221 spatial grid points. For time-stepping, the problem is treated *semi-implicitly*, i.e., the diffusion part is solved with
222 Fast-Fourier Tranform (FFT) and the nonlinear term is treated explicitly.
224 An exact solution is not known, but instead the numerical solution can be compared via a generated reference solution computed
225 by a ``SciPy`` routine.
227 Parameters
228 ----------
229 nvars : List of int tuples, optional
230 Number of unknowns in the problem, e.g. ``nvars=[(128, 128), (128, 128)]``.
231 nu : float, optional
232 Problem parameter :math:`\nu`.
233 eps : float, optional
234 Scaling parameter :math:`\varepsilon`.
235 radius : float, optional
236 Radius of the circles.
237 L : float, optional
238 Denotes the period of the function to be approximated for the Fourier transform.
239 init_type : str, optional
240 Indicates which type of initial condition is used.
242 Attributes
243 ----------
244 xvalues : np.1darray
245 Grid points in space.
246 dx : float
247 Mesh width.
248 lap : np.1darray
249 Spectral operator for Laplacian.
250 """
252 def __init__(self, nvars=None, nu=2, eps=0.04, radius=0.25, L=1.0, init_type='circle'):
253 """Initialization routine"""
255 if nvars is None:
256 nvars = [(256, 256), (64, 64)]
258 super().__init__(nvars, nu, eps, radius, L, init_type)
259 self.lap -= 2.0 / self.eps**2
261 def eval_f(self, u, t):
262 """
263 Routine to evaluate the right-hand side of the problem.
265 Parameters
266 ----------
267 u : dtype_u
268 Current values of the numerical solution.
269 t : float
270 Current time of the numerical solution is computed.
272 Returns
273 -------
274 f : dtype_f
275 The right-hand side of the problem.
276 """
278 f = self.dtype_f(self.init)
279 v = u.flatten()
280 tmp = self.lap * np.fft.rfft2(u)
281 f.impl[:] = np.fft.irfft2(tmp)
282 if self.eps > 0:
283 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu) + 2.0 / self.eps**2 * v).reshape(self.nvars)
284 return f
286 def solve_system(self, rhs, factor, u0, t):
287 """
288 Simple FFT solver for the diffusion part.
290 Parameters
291 ----------
292 rhs : dtype_f
293 Right-hand side for the linear system.
294 factor : float
295 Abbrev. for the node-to-node stepsize (or any other factor required).
296 u0 : dtype_u
297 Initial guess for the iterative solver (not used here so far).
298 t : float
299 Current time (e.g. for time-dependent BCs).
301 Returns
302 -------
303 me : dtype_u
304 The solution as mesh.
305 """
307 me = self.dtype_u(self.init)
309 tmp = np.fft.rfft2(rhs) / (1.0 - factor * self.lap)
310 me[:] = np.fft.irfft2(tmp)
312 return me