Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD_gpu.py: 0%
130 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import numpy as np
2import cupy as cp
3import cupyx.scipy.sparse as csp
4from cupyx.scipy.sparse.linalg import cg # , spsolve, gmres, minres
6from pySDC.core.errors import ProblemError
7from pySDC.core.problem import Problem
8from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh, comp2_cupy_mesh
10# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
13class allencahn_fullyimplicit(Problem): # pragma: no cover
14 r"""
15 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
17 .. math::
18 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
20 for constant parameter :math:`\nu`. Initial condition are circles of the form
22 .. math::
23 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
25 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
26 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
28 This class is especially developed for solving it on GPUs using ``CuPy``.
30 Parameters
31 ----------
32 nvars : tuple of int, optional
33 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``.
34 nu : float, optional
35 Problem parameter :math:`\nu`.
36 eps : float, optional
37 Scaling parameter :math:`\varepsilon`.
38 newton_maxiter : int, optional
39 Maximum number of iterations for the Newton solver.
40 newton_tol : float, optional
41 Tolerance for Newton's method to terminate.
42 lin_tol : float, optional
43 Tolerance for linear solver to terminate.
44 lin_maxiter : int, optional
45 Maximum number of iterations for the linear solver.
46 radius : float, optional
47 Radius of the circles.
49 Attributes
50 ----------
51 A : scipy.spdiags
52 Second-order FD discretization of the 2D laplace operator.
53 dx : float
54 Distance between two spatial nodes (same for both directions).
55 xvalues : np.1darray
56 Spatial grid points, here both dimensions have the same grid points.
57 newton_itercount : int
58 Number of iterations of Newton solver.
59 lin_itercount
60 Number of iterations of linear solver.
61 newton_ncalls : int
62 Number of calls of Newton solver.
63 lin_ncalls : int
64 Number of calls of linear solver.
65 """
67 dtype_u = cupy_mesh
68 dtype_f = cupy_mesh
70 def __init__(
71 self,
72 nvars=(128, 128),
73 nu=2,
74 eps=0.04,
75 newton_maxiter=1e-12,
76 newton_tol=100,
77 lin_tol=1e-8,
78 lin_maxiter=100,
79 radius=0.25,
80 ):
81 """Initialization routine"""
82 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
83 if len(nvars) != 2:
84 raise ProblemError('this is a 2d example, got %s' % nvars)
85 if nvars[0] != nvars[1]:
86 raise ProblemError('need a square domain, got %s' % nvars)
87 if nvars[0] % 2 != 0:
88 raise ProblemError('the setup requires nvars = 2^p per dimension')
90 # invoke super init, passing number of dofs, dtype_u and dtype_f
91 super().__init__((nvars, None, cp.dtype('float64')))
92 self._makeAttributeAndRegister(
93 'nvars',
94 'nu',
95 'eps',
96 'newton_maxiter',
97 'newton_tol',
98 'lin_tol',
99 'lin_maxiter',
100 'radius',
101 localVars=locals(),
102 readOnly=True,
103 )
105 # compute dx and get discretization matrix A
106 self.dx = 1.0 / self.nvars[0]
107 self.A = self.__get_A(self.nvars, self.dx)
108 self.xvalues = cp.array([i * self.dx - 0.5 for i in range(self.nvars[0])])
110 self.newton_itercount = 0
111 self.lin_itercount = 0
112 self.newton_ncalls = 0
113 self.lin_ncalls = 0
115 @staticmethod
116 def __get_A(N, dx):
117 """
118 Helper function to assemble FD matrix A in sparse format.
120 Parameters
121 ----------
122 N : list
123 Number of degrees of freedom.
124 dx : float
125 Distance between two spatial nodes.
127 Returns
128 -------
129 A : cupyx.scipy.sparse.csr_matrix
130 CuPy-matrix A in CSR format.
131 """
133 stencil = cp.asarray([-2, 1])
134 A = stencil[0] * csp.eye(N[0], format='csr')
135 for i in range(1, len(stencil)):
136 A += stencil[i] * csp.eye(N[0], k=-i, format='csr')
137 A += stencil[i] * csp.eye(N[0], k=+i, format='csr')
138 A += stencil[i] * csp.eye(N[0], k=N[0] - i, format='csr')
139 A += stencil[i] * csp.eye(N[0], k=-N[0] + i, format='csr')
140 A = csp.kron(A, csp.eye(N[0])) + csp.kron(csp.eye(N[1]), A)
141 A *= 1.0 / (dx**2)
142 return A
144 # noinspection PyTypeChecker
145 def solve_system(self, rhs, factor, u0, t):
146 """
147 Simple Newton solver.
149 Parameters
150 ----------
151 rhs : dtype_f
152 Right-hand side for the nonlinear system.
153 factor : float
154 Abbrev. for the node-to-node stepsize (or any other factor required).
155 u0 : dtype_u
156 Initial guess for the iterative solver.
157 t : float
158 Current time (required here for the BC).
160 Returns
161 -------
162 u : dtype_u
163 The solution as mesh.
164 """
166 u = self.dtype_u(u0).flatten()
167 z = self.dtype_u(self.init, val=0.0).flatten()
168 nu = self.nu
169 eps2 = self.eps**2
171 Id = csp.eye(self.nvars[0] * self.nvars[1])
173 # start newton iteration
174 n = 0
175 res = 99
176 while n < self.newton_maxiter:
177 # form the function g with g(u) = 0
178 g = u - factor * (self.A.dot(u) + 1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
180 # if g is close to 0, then we are done
181 res = cp.linalg.norm(g, np.inf)
183 if res < self.newton_tol:
184 break
186 # assemble dg
187 dg = Id - factor * (self.A + 1.0 / eps2 * csp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
189 # newton update: u1 = u0 - g/dg
190 # u -= spsolve(dg, g)
191 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0]
192 # increase iteration count
193 n += 1
194 # print(n, res)
196 # if n == self.newton_maxiter:
197 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
199 me = self.dtype_u(self.init)
200 me[:] = u.reshape(self.nvars)
202 self.newton_ncalls += 1
203 self.newton_itercount += n
205 return me
207 def eval_f(self, u, t):
208 """
209 Routine to evaluate the right-hand side of the problem.
211 Parameters
212 ----------
213 u : dtype_u
214 Current values of the numerical solution.
215 t : float
216 Current time of the numerical solution is computed (not used here).
218 Returns
219 -------
220 f : dtype_f
221 The right-hand side of the problem.
222 """
223 f = self.dtype_f(self.init)
224 v = u.flatten()
225 f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
227 return f
229 def u_exact(self, t):
230 r"""
231 Routine to compute the exact solution at time :math:`t`.
233 Parameters
234 ----------
235 t : float
236 Time of the exact solution.
238 Returns
239 -------
240 me : dtype_u
241 The exact solution.
242 """
244 assert t == 0, 'ERROR: u_exact only valid for t=0'
245 me = self.dtype_u(self.init, val=0.0)
246 mx, my = cp.meshgrid(self.xvalues, self.xvalues)
247 me[:] = cp.tanh((self.radius - cp.sqrt(mx**2 + my**2)) / (cp.sqrt(2) * self.eps))
248 # print(type(me))
249 return me
252# noinspection PyUnusedLocal
253class allencahn_semiimplicit(allencahn_fullyimplicit):
254 r"""
255 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
257 .. math::
258 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
260 for constant parameter :math:`\nu`. Initial condition are circles of the form
262 .. math::
263 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
265 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
266 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the
267 conjugate gradients method, and the system containing the rest of the right-hand side is only evaluated at each time.
269 This class is especially developed for solving it on GPUs using ``CuPy``.
270 """
272 dtype_f = imex_cupy_mesh
274 def eval_f(self, u, t):
275 """
276 Routine to evaluate the right-hand side of the problem.
278 Parameters
279 ----------
280 u : dtype_u
281 Current values of the numerical solution.
282 t : float
283 Current time of the numerical solution is computed.
285 Returns
286 -------
287 f : dtype_f
288 The right-hand side of the problem.
289 """
290 f = self.dtype_f(self.init)
291 v = u.flatten()
292 f.impl[:] = self.A.dot(v).reshape(self.nvars)
293 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
295 return f
297 def solve_system(self, rhs, factor, u0, t):
298 r"""
299 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
301 Parameters
302 ----------
303 rhs : dtype_f
304 Right-hand side for the linear system.
305 factor : float
306 Abbrev. for the local stepsize (or any other factor required).
307 u0 : dtype_u
308 Initial guess for the iterative solver.
309 t : float
310 Current time (e.g. for time-dependent BCs).
312 Returns
313 -------
314 me : dtype_u
315 The solution as mesh.
316 """
318 class context:
319 num_iter = 0
321 def callback(xk):
322 context.num_iter += 1
323 return context.num_iter
325 me = self.dtype_u(self.init)
327 Id = csp.eye(self.nvars[0] * self.nvars[1])
329 me[:] = cg(
330 Id - factor * self.A,
331 rhs.flatten(),
332 x0=u0.flatten(),
333 tol=self.lin_tol,
334 maxiter=self.lin_maxiter,
335 callback=callback,
336 )[0].reshape(self.nvars)
338 self.lin_ncalls += 1
339 self.lin_itercount += context.num_iter
341 return me
344# noinspection PyUnusedLocal
345class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
346 r"""
347 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
349 .. math::
350 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
352 for constant parameter :math:`\nu`. Initial condition are circles of the form
354 .. math::
355 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
357 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
358 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
359 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
360 is only evaluated at each time.
362 This class is especially developed for solving it on GPUs using ``CuPy``.
363 """
365 dtype_f = imex_cupy_mesh
367 def eval_f(self, u, t):
368 """
369 Routine to evaluate the right-hand side of the problem.
371 Parameters
372 ----------
373 u : dtype_u
374 Current values of the numerical solution.
375 t : float
376 Current time of the numerical solution is computed.
378 Returns
379 -------
380 f : dtype_f
381 The right-hand side of the problem.
382 """
383 f = self.dtype_f(self.init)
384 v = u.flatten()
385 f.impl[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
386 f.expl[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
388 return f
390 def solve_system(self, rhs, factor, u0, t):
391 """
392 Simple Newton solver.
394 Parameters
395 ----------
396 rhs : dtype_f
397 Right-hand side for the nonlinear system.
398 factor : float
399 Abbrev. for the node-to-node stepsize (or any other factor required).
400 u0 : dtype_u
401 Initial guess for the iterative solver.
402 t : float
403 Current time (required here for the BC).
405 Returns
406 -------
407 me : dtype_u
408 The solution as mesh.
409 """
411 u = self.dtype_u(u0).flatten()
412 z = self.dtype_u(self.init, val=0.0).flatten()
413 nu = self.nu
414 eps2 = self.eps**2
416 Id = csp.eye(self.nvars[0] * self.nvars[1])
418 # start newton iteration
419 n = 0
420 res = 99
421 while n < self.newton_maxiter:
422 # form the function g with g(u) = 0
423 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
425 # if g is close to 0, then we are done
426 res = cp.linalg.norm(g, np.inf)
428 if res < self.newton_tol:
429 break
431 # assemble dg
432 dg = Id - factor * (self.A - 1.0 / eps2 * csp.diags(((nu + 1) * u**nu), offsets=0))
434 # newton update: u1 = u0 - g/dg
435 # u -= spsolve(dg, g)
436 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0]
437 # increase iteration count
438 n += 1
439 # print(n, res)
441 # if n == self.newton_maxiter:
442 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
444 me = self.dtype_u(self.init)
445 me[:] = u.reshape(self.nvars)
447 self.newton_ncalls += 1
448 self.newton_itercount += n
450 return me
453# noinspection PyUnusedLocal
454class allencahn_multiimplicit(allencahn_fullyimplicit):
455 r"""
456 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
458 .. math::
459 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
461 for constant parameter :math:`\nu`. Initial condition are circles of the form
463 .. math::
464 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
466 for :math::`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
467 treated in *multi-implicit* fashion, i.e., the nonlinear system containing the second order term is solved by the
468 conjugate gradients method, and the system containing the rest of the right-hand side will be solved by Newton's method.
470 This class is especially developed for solving it on GPUs using ``CuPy``.
471 """
473 dtype_f = comp2_cupy_mesh
475 def eval_f(self, u, t):
476 """
477 Routine to evaluate the right-hand side of the problem.
479 Parameters
480 ----------
481 u : dtype_u
482 Current values of the numerical solution.
483 t : float
484 Current time of the numerical solution is computed (not used here).
486 Returns
487 -------
488 f : dtype_f
489 The right-hand side of the problem.
490 """
491 f = self.dtype_f(self.init)
492 v = u.flatten()
493 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
494 f.comp2[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
496 return f
498 def solve_system_1(self, rhs, factor, u0, t):
499 r"""
500 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
502 Parameters
503 ----------
504 rhs : dtype_f
505 Right-hand side for the linear system.
506 factor : float
507 Abbrev. for the local stepsize (or any other factor required).
508 u0 : dtype_u
509 Initial guess for the iterative solver.
510 t : float
511 Current time (e.g. for time-dependent BCs).
513 Returns
514 -------
515 me : dtype_u
516 The solution as mesh.
517 """
519 class context:
520 num_iter = 0
522 def callback(xk):
523 context.num_iter += 1
524 return context.num_iter
526 me = self.dtype_u(self.init)
528 Id = csp.eye(self.nvars[0] * self.nvars[1])
530 me[:] = cg(
531 Id - factor * self.A,
532 rhs.flatten(),
533 x0=u0.flatten(),
534 tol=self.lin_tol,
535 maxiter=self.lin_maxiter,
536 callback=callback,
537 )[0].reshape(self.nvars)
539 self.lin_ncalls += 1
540 self.lin_itercount += context.num_iter
542 return me
544 def solve_system_2(self, rhs, factor, u0, t):
545 """
546 Simple Newton solver.
548 Parameters
549 ----------
550 rhs : dtype_f
551 Right-hand side for the nonlinear system
552 factor : float
553 Abbrev. for the node-to-node stepsize (or any other factor required).
554 u0 : dtype_u
555 Initial guess for the iterative solver.
556 t : float
557 Current time (required here for the BC).
559 Returns
560 -------
561 me : dtype_u
562 The solution as mesh.
563 """
565 u = self.dtype_u(u0).flatten()
566 z = self.dtype_u(self.init, val=0.0).flatten()
567 nu = self.nu
568 eps2 = self.eps**2
570 Id = csp.eye(self.nvars[0] * self.nvars[1])
572 # start newton iteration
573 n = 0
574 res = 99
575 while n < self.newton_maxiter:
576 # form the function g with g(u) = 0
577 g = u - factor * (1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
579 # if g is close to 0, then we are done
580 res = cp.linalg.norm(g, np.inf)
582 if res < self.newton_tol:
583 break
585 # assemble dg
586 dg = Id - factor * (1.0 / eps2 * csp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
588 # newton update: u1 = u0 - g/dg
589 # u -= spsolve(dg, g)
590 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0]
591 # increase iteration count
592 n += 1
593 # print(n, res)
595 # if n == self.newton_maxiter:
596 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
598 me = self.dtype_u(self.init)
599 me[:] = u.reshape(self.nvars)
601 self.newton_ncalls += 1
602 self.newton_itercount += n
604 return me
607# noinspection PyUnusedLocal
608class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
609 r"""
610 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
612 .. math::
613 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
615 for constant parameter :math:`\nu`. The initial condition has the form of circles
617 .. math::
618 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
620 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
621 is used here to get another kind of *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
622 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
623 is solved by a linear solver provided by a ``SciPy`` routine.
625 This class is especially developed for solving it on GPUs using ``CuPy``.
626 """
628 dtype_f = comp2_cupy_mesh
630 def eval_f(self, u, t):
631 """
632 Routine to evaluate the right-hand side of the problem.
634 Parameters
635 ----------
636 u : dtype_u
637 Current values of the numerical solution.
638 t : float
639 Current time of the numerical solution is computed (not used here).
641 Returns
642 -------
643 f : dtype_f
644 The right-hand side of the problem.
645 """
646 f = self.dtype_f(self.init)
647 v = u.flatten()
648 f.comp1[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
649 f.comp2[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
651 return f
653 def solve_system_1(self, rhs, factor, u0, t):
654 """
655 Simple Newton solver.
657 Parameters
658 ----------
659 rhs : dtype_f
660 Right-hand side for the nonlinear system.
661 factor : float
662 Abbrev. for the node-to-node stepsize (or any other factor required).
663 u0 : dtype_u
664 Initial guess for the iterative solver.
665 t : float
666 Current time (required here for the BC).
668 Returns
669 -------
670 me : dtype_u
671 The solution as mesh.
672 """
674 u = self.dtype_u(u0).flatten()
675 z = self.dtype_u(self.init, val=0.0).flatten()
676 nu = self.nu
677 eps2 = self.eps**2
679 Id = csp.eye(self.nvars[0] * self.nvars[1])
681 # start newton iteration
682 n = 0
683 res = 99
684 while n < self.newton_maxiter:
685 # form the function g with g(u) = 0
686 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
688 # if g is close to 0, then we are done
689 res = cp.linalg.norm(g, np.inf)
691 if res < self.newton_tol:
692 break
694 # assemble dg
695 dg = Id - factor * (self.A - 1.0 / eps2 * csp.diags(((nu + 1) * u**nu), offsets=0))
697 # newton update: u1 = u0 - g/dg
698 # u -= spsolve(dg, g)
699 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0]
700 # increase iteration count
701 n += 1
702 # print(n, res)
704 # if n == self.newton_maxiter:
705 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
707 me = self.dtype_u(self.init)
708 me[:] = u.reshape(self.nvars)
710 self.newton_ncalls += 1
711 self.newton_itercount += n
713 return me
715 def solve_system_2(self, rhs, factor, u0, t):
716 r"""
717 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
719 Parameters
720 ----------
721 rhs : dtype_f
722 Right-hand side for the linear system.
723 factor : float
724 Abbrev. for the local stepsize (or any other factor required).
725 u0 : dtype_u
726 Initial guess for the iterative solver.
727 t : float
728 Current time (e.g. for time-dependent BCs).
730 Returns
731 -------
732 me : dtype_u
733 The solution as mesh.
734 """
736 me = self.dtype_u(self.init)
738 me[:] = (1.0 / (1.0 - factor * 1.0 / self.eps**2) * rhs).reshape(self.nvars)
739 return me