Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD.py: 96%
212 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
2import scipy.sparse as sp
3import scipy.sparse.linalg as spla
5from pySDC.core.errors import ParameterError, ProblemError
6from pySDC.core.problem import Problem, WorkCounter
7from pySDC.helpers import problem_helper
8from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
10# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
13# noinspection PyUnusedLocal
14class allencahn_fullyimplicit(Problem):
15 r"""
16 2D periodic Allen-Cahn equation with finite differences, fully implicit with Newton.
18 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
19 phases at :math:`u = 0` and :math:`u = 1`
21 .. math::
22 \frac{\partial u}{\partial t} = \Delta u
23 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
25 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
26 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
28 Initial condition are circles of the form
30 .. math::
31 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
32 {\sqrt{2}\varepsilon}\right)\right)
34 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
35 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
37 Parameters
38 ----------
39 nvars : tuple of int, optional
40 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``.
41 nu : int, optional
42 Exponent of the double well; :math:`\nu = 2` is the standard Allen-Cahn nonlinearity.
43 eps : float, optional
44 Scaling parameter :math:`\varepsilon`.
45 newton_maxiter : int, optional
46 Maximum number of iterations for the Newton solver.
47 newton_tol : float, optional
48 Tolerance for Newton's method to terminate.
49 lin_tol : float, optional
50 Tolerance for linear solver to terminate.
51 lin_maxiter : int, optional
52 Maximum number of iterations for the linear solver.
53 inexact_linear_ratio : float, optional
54 Ratio of tolerance of linear solver to the Newton residual, overrides ``lin_tol`` if set.
55 radius : float, optional
56 Radius of the circles.
57 order : int, optional
58 Order of the finite difference matrix.
59 useGPU : bool, optional
60 Run on the GPU with CuPy instead of on the CPU with NumPy.
62 Attributes
63 ----------
64 A : scipy.spdiags
65 Second-order FD discretization of the 2D laplace operator.
66 dx : float
67 Distance between two spatial nodes (same for both directions).
68 xvalues : np.1darray
69 Spatial grid points, here both dimensions have the same grid points.
70 newton_ncalls : int
71 Number of calls of the Newton solver. The iterations themselves are counted in
72 ``work_counters['newton']``, and the linear ones in ``work_counters['linear']``.
73 lin_ncalls : int
74 Number of calls of the linear solver.
75 """
77 dtype_u = mesh
78 dtype_f = mesh
80 xp = np
81 xsp = sp
82 linalg = spla
84 def setup_GPU(self):
85 """
86 Switch the array, sparse and solver modules and the datatypes over to CuPy.
88 This changes the class, not the instance, as everything else in pySDC that does this
89 does: once one instance of a class runs on the GPU, they all do.
90 """
91 import cupy as cp
92 import cupyx.scipy.sparse as csp
93 import cupyx.scipy.sparse.linalg as cspla
94 from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh, comp2_cupy_mesh
96 self.xp = cp
97 self.xsp = csp
98 self.linalg = cspla
99 self.dtype_u = cupy_mesh
100 # .get, not [], because this runs once per instance and the class keeps what it is given
101 GPU_versions = {mesh: cupy_mesh, imex_mesh: imex_cupy_mesh, comp2_mesh: comp2_cupy_mesh}
102 self.dtype_f = GPU_versions.get(self.dtype_f, self.dtype_f)
104 def __init__(
105 self,
106 nvars=(128, 128),
107 nu=2,
108 eps=0.04,
109 newton_maxiter=200,
110 newton_tol=1e-12,
111 lin_tol=1e-8,
112 lin_maxiter=100,
113 inexact_linear_ratio=None,
114 radius=0.25,
115 order=2,
116 useGPU=False,
117 ):
118 """Initialization routine"""
119 if useGPU:
120 self.setup_GPU()
122 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
123 if len(nvars) != 2:
124 raise ProblemError('this is a 2d example, got %s' % nvars)
125 if nvars[0] != nvars[1]:
126 raise ProblemError('need a square domain, got %s' % nvars)
127 if nvars[0] % 2 != 0:
128 raise ProblemError('the setup requires nvars = 2^p per dimension')
130 # invoke super init, passing number of dofs, dtype_u and dtype_f
131 super().__init__((nvars, None, np.dtype('float64')))
132 self._makeAttributeAndRegister(
133 'nvars',
134 'nu',
135 'eps',
136 'radius',
137 'order',
138 'useGPU',
139 localVars=locals(),
140 readOnly=True,
141 )
142 self._makeAttributeAndRegister(
143 'newton_maxiter',
144 'newton_tol',
145 'lin_tol',
146 'lin_maxiter',
147 'inexact_linear_ratio',
148 localVars=locals(),
149 readOnly=False,
150 )
152 # compute dx and get discretization matrix A
153 self.dx = 1.0 / self.nvars[0]
154 self.A, _ = problem_helper.get_finite_difference_matrix(
155 derivative=2,
156 order=self.order,
157 stencil_type='center',
158 dx=self.dx,
159 size=self.nvars[0],
160 dim=2,
161 bc='periodic',
162 cupy=self.useGPU,
163 )
164 self.xvalues = self.xp.arange(self.nvars[0]) * self.dx - 0.5
166 self.newton_ncalls = 0
167 self.lin_ncalls = 0
169 self.work_counters['newton'] = WorkCounter()
170 self.work_counters['rhs'] = WorkCounter()
171 self.work_counters['linear'] = WorkCounter()
173 def reaction(self, u):
174 r"""
175 The reaction term, :math:`\frac{1}{2\varepsilon^2}(2u - 1)\left(1 - (2u - 1)^\nu\right)`.
177 The wells sit at :math:`u = 0` and :math:`u = 1`, so the double well is symmetric about
178 :math:`2u - 1`; writing the term in that variable is what lets :math:`\nu` keep the meaning
179 it has always had here. For the default :math:`\nu = 2` this is
180 :math:`-\frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
181 """
182 v = 2.0 * u - 1.0
183 return 0.5 / self.eps**2 * v * (1.0 - v**self.nu)
185 def reaction_prime(self, u):
186 """Derivative of :meth:`reaction`, ready to go on the diagonal of a Jacobian."""
187 v = 2.0 * u - 1.0
188 return 1.0 / self.eps**2 * (1.0 - (self.nu + 1.0) * v**self.nu)
190 def reaction_cubic(self, u):
191 r"""The stiff part of :meth:`reaction`, :math:`-\frac{1}{2\varepsilon^2}(2u - 1)^{\nu + 1}`."""
192 return -0.5 / self.eps**2 * (2.0 * u - 1.0) ** (self.nu + 1)
194 def reaction_cubic_prime(self, u):
195 """Derivative of :meth:`reaction_cubic`."""
196 return -(self.nu + 1.0) / self.eps**2 * (2.0 * u - 1.0) ** self.nu
198 def reaction_linear(self, u):
199 r"""The rest of :meth:`reaction`, :math:`\frac{1}{2\varepsilon^2}(2u - 1)`, so the two sum back to it."""
200 return 0.5 / self.eps**2 * (2.0 * u - 1.0)
202 # noinspection PyTypeChecker
203 def solve_system(self, rhs, factor, u0, t):
204 """
205 Simple Newton solver.
207 Parameters
208 ----------
209 rhs : dtype_f
210 Right-hand side for the nonlinear system
211 factor : float
212 Abbrev. for the node-to-node stepsize (or any other factor required).
213 u0 : dtype_u
214 Initial guess for the iterative solver.
215 t : float
216 Current time (required here for the BC).
218 Returns
219 -------
220 me : dtype_u
221 The solution as mesh.
222 """
224 # plain arrays, not the datatype: its `__abs__` is a norm, which CuPy's `linalg.norm` trips over
225 u = u0.view(self.xp.ndarray).flatten()
226 b = rhs.view(self.xp.ndarray).flatten()
227 z = self.xp.zeros_like(u)
229 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
231 # start newton iteration
232 n = 0
233 res = 99
234 while n < self.newton_maxiter:
235 # form the function g with g(u) = 0
236 g = u - factor * (self.A.dot(u) + self.reaction(u)) - b
238 # if g is close to 0, then we are done
239 res = self.xp.linalg.norm(g, self.xp.inf)
241 # do inexactness in the linear solver
242 if self.inexact_linear_ratio:
243 self.lin_tol = res * self.inexact_linear_ratio
245 if res < self.newton_tol:
246 break
248 # assemble dg
249 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_prime(u), offsets=0))
251 # newton update: u1 = u0 - g/dg
252 # u -= spsolve(dg, g)
253 u -= self.linalg.cg(
254 dg, g, x0=z, rtol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
255 )[0]
256 # increase iteration count
257 n += 1
258 # print(n, res)
260 self.work_counters['newton']()
262 # if n == self.newton_maxiter:
263 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
265 me = self.dtype_u(self.init)
266 me[:] = u.reshape(self.nvars)
268 self.newton_ncalls += 1
270 return me
272 def eval_f(self, u, t):
273 """
274 Routine to evaluate the right-hand side of the problem.
276 Parameters
277 ----------
278 u : dtype_u
279 Current values of the numerical solution.
280 t : float
281 Current time of the numerical solution is computed (not used here).
283 Returns
284 -------
285 f : dtype_f
286 The right-hand side of the problem.
287 """
288 f = self.dtype_f(self.init)
289 v = u.flatten()
290 f[:] = (self.A.dot(v) + self.reaction(v)).reshape(self.nvars)
292 self.work_counters['rhs']()
293 return f
295 def u_exact(self, t, u_init=None, t_init=None):
296 r"""
297 Routine to compute the exact solution at time :math:`t`.
299 Parameters
300 ----------
301 t : float
302 Time of the exact solution.
303 u_init : dtype_u, optional
304 Initial conditions for getting the exact solution.
305 t_init : float, optional
306 The starting time.
308 Returns
309 -------
310 me : dtype_u
311 The exact solution.
312 """
313 me = self.dtype_u(self.init, val=0.0)
314 if t > 0:
316 def eval_rhs(t, u):
317 return self.eval_f(u.reshape(self.init[0]), t).flatten()
319 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
321 else:
322 X, Y = self.xp.meshgrid(self.xvalues, self.xvalues)
323 r2 = X**2 + Y**2
324 me[:] = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)))
326 return me
329# noinspection PyUnusedLocal
330class allencahn_semiimplicit(allencahn_fullyimplicit):
331 r"""
332 2D periodic Allen-Cahn equation with finite differences, IMEX with Laplacian implicit, reaction explicit.
334 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
335 phases at :math:`u = 0` and :math:`u = 1`
337 .. math::
338 \frac{\partial u}{\partial t} = \Delta u
339 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
341 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
342 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
344 Initial condition are circles of the form
346 .. math::
347 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
348 {\sqrt{2}\varepsilon}\right)\right)
350 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
351 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
352 method, and the system containing the rest of the right-hand side is only evaluated at each time.
353 """
355 dtype_f = imex_mesh
357 def eval_f(self, u, t):
358 """
359 Routine to evaluate the right-hand side of the problem.
361 Parameters
362 ----------
363 u : dtype_u
364 Current values of the numerical solution.
365 t : float
366 Current time of the numerical solution is computed (not used here).
368 Returns
369 -------
370 f : dtype_f
371 The right-hand side of the problem.
372 """
373 f = self.dtype_f(self.init)
374 v = u.flatten()
375 f.impl[:] = self.A.dot(v).reshape(self.nvars)
376 f.expl[:] = self.reaction(v).reshape(self.nvars)
378 self.work_counters['rhs']()
379 return f
381 def solve_system(self, rhs, factor, u0, t):
382 r"""
383 Pointwise solver for :math:`\vec{u} - factor \cdot \frac{1}{2\varepsilon^2}(2\vec{u} - 1) = \vec{rhs}`.
385 Parameters
386 ----------
387 rhs : dtype_f
388 Right-hand side for the linear system.
389 factor : float
390 Abbrev. for the local stepsize (or any other factor required).
391 u0 : dtype_u
392 Initial guess for the iterative solver.
393 t : float
394 Current time (e.g. for time-dependent BCs).
396 Returns
397 -------
398 me : dtype_u
399 The solution as mesh.
400 """
402 me = self.dtype_u(self.init)
404 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
406 me[:] = self.linalg.cg(
407 Id - factor * self.A,
408 rhs.flatten(),
409 x0=u0.flatten(),
410 rtol=self.lin_tol,
411 maxiter=self.lin_maxiter,
412 atol=0,
413 callback=self.work_counters['linear'],
414 )[0].reshape(self.nvars)
416 self.lin_ncalls += 1
418 return me
420 def u_exact(self, t, u_init=None, t_init=None):
421 """
422 Routine to compute the exact solution at time t.
424 Parameters
425 ----------
426 t : float
427 Time of the exact solution.
428 u_init : dtype_u, optional
429 Initial conditions for getting the exact solution.
430 t_init : float, optional
431 The starting time.
433 Returns
434 -------
435 me : dtype_u
436 The exact solution.
437 """
438 me = self.dtype_u(self.init, val=0.0)
439 if t > 0:
441 def eval_rhs(t, u):
442 f = self.eval_f(u.reshape(self.init[0]), t)
443 return (f.impl + f.expl).flatten()
445 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
446 else:
447 me[:] = super().u_exact(t, u_init, t_init)
448 return me
451# noinspection PyUnusedLocal
452class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
453 r"""
454 2D periodic Allen-Cahn equation with finite differences, IMEX with Laplacian and cubic term implicit.
456 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions, with the two
457 phases at :math:`u = 0` and :math:`u = 1`
459 .. math::
460 \frac{\partial u}{\partial t} = \Delta u
461 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
463 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
464 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
466 Initial condition are circles of the form
468 .. math::
469 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
470 {\sqrt{2}\varepsilon}\right)\right)
472 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
473 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{2\varepsilon^2}(2u - 1)^3`
474 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{2\varepsilon^2}(2u - 1)`
475 is only evaluated at each time.
476 """
478 dtype_f = imex_mesh
480 def eval_f(self, u, t):
481 """
482 Routine to evaluate the right-hand side of the problem.
484 Parameters
485 ----------
486 u : dtype_u
487 Current values of the numerical solution.
488 t : float
489 Current time of the numerical solution is computed.
491 Returns
492 -------
493 f : dtype_f
494 The right-hand side of the problem.
495 """
496 f = self.dtype_f(self.init)
497 v = u.flatten()
498 f.impl[:] = (self.A.dot(v) + self.reaction_cubic(v)).reshape(self.nvars)
499 f.expl[:] = self.reaction_linear(v).reshape(self.nvars)
501 self.work_counters['rhs']()
502 return f
504 def solve_system(self, rhs, factor, u0, t):
505 """
506 Simple Newton solver.
508 Parameters
509 ----------
510 rhs : dtype_f
511 Right-hand side for the nonlinear system.
512 factor : float
513 Abbrev. for the node-to-node stepsize (or any other factor required).
514 u0 : dtype_u
515 Initial guess for the iterative solver.
516 t : float
517 Current time (required here for the BC).
519 Returns
520 -------
521 me : dtype_u
522 The solution as mesh.
523 """
525 # plain arrays, not the datatype: its `__abs__` is a norm, which CuPy's `linalg.norm` trips over
526 u = u0.view(self.xp.ndarray).flatten()
527 b = rhs.view(self.xp.ndarray).flatten()
528 z = self.xp.zeros_like(u)
530 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
532 # start newton iteration
533 n = 0
534 res = 99
535 while n < self.newton_maxiter:
536 # form the function g with g(u) = 0
537 g = u - factor * (self.A.dot(u) + self.reaction_cubic(u)) - b
539 # if g is close to 0, then we are done
540 res = self.xp.linalg.norm(g, self.xp.inf)
542 if res < self.newton_tol:
543 break
545 # assemble dg
546 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_cubic_prime(u), offsets=0))
548 # newton update: u1 = u0 - g/dg
549 # u -= spsolve(dg, g)
550 u -= self.linalg.cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
551 # increase iteration count
552 n += 1
553 # print(n, res)
555 self.work_counters['newton']()
557 # if n == self.newton_maxiter:
558 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
560 me = self.dtype_u(self.init)
561 me[:] = u.reshape(self.nvars)
563 self.newton_ncalls += 1
565 return me
568# noinspection PyUnusedLocal
569class allencahn_multiimplicit(allencahn_fullyimplicit):
570 r"""
571 2D periodic Allen-Cahn equation with finite differences, multi-implicit: Laplacian (CG), reaction (Newton).
573 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
574 phases at :math:`u = 0` and :math:`u = 1`
576 .. math::
577 \frac{\partial u}{\partial t} = \Delta u
578 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
580 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
581 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
583 Initial condition are circles of the form
585 .. math::
586 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
587 {\sqrt{2}\varepsilon}\right)\right)
589 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
590 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
591 method, and the system containing the rest of the right-hand side will be solved by Newton's method.
592 """
594 dtype_f = comp2_mesh
596 def eval_f(self, u, t):
597 """
598 Routine to evaluate the right-hand side of the problem.
600 Parameters
601 ----------
602 u : dtype_u
603 Current values of the numerical solution.
604 t : float
605 Current time of the numerical solution is computed.
607 Returns
608 -------
609 f : dtype_f
610 The right-hand side of the problem.
611 """
612 f = self.dtype_f(self.init)
613 v = u.flatten()
614 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
615 f.comp2[:] = self.reaction(v).reshape(self.nvars)
617 self.work_counters['rhs']()
618 return f
620 def solve_system_1(self, rhs, factor, u0, t):
621 r"""
622 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
624 Parameters
625 ----------
626 rhs : dtype_f
627 Right-hand side for the linear system.
628 factor : float
629 Abbrev. for the local stepsize (or any other factor required).
630 u0 : dtype_u
631 Initial guess for the iterative solver.
632 t : float
633 Current time (e.g. for time-dependent BCs).
635 Returns
636 -------
637 me : dtype_u
638 The solution as mesh.
639 """
641 me = self.dtype_u(self.init)
643 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
645 me[:] = self.linalg.cg(
646 Id - factor * self.A,
647 rhs.flatten(),
648 x0=u0.flatten(),
649 rtol=self.lin_tol,
650 maxiter=self.lin_maxiter,
651 atol=0,
652 callback=self.work_counters['linear'],
653 )[0].reshape(self.nvars)
655 self.lin_ncalls += 1
657 return me
659 def solve_system_2(self, rhs, factor, u0, t):
660 """
661 Simple Newton solver.
663 Parameters
664 ----------
665 rhs : dtype_f
666 Right-hand side for the nonlinear system.
667 factor : float
668 Abbrev. for the node-to-node stepsize (or any other factor required).
669 u0 : dtype_u
670 Initial guess for the iterative solver.
671 t : float
672 Current time (required here for the BC).
674 Returns
675 -------
676 me : dtype_u
677 The solution as mesh.
678 """
680 # plain arrays, not the datatype: its `__abs__` is a norm, which CuPy's `linalg.norm` trips over
681 u = u0.view(self.xp.ndarray).flatten()
682 b = rhs.view(self.xp.ndarray).flatten()
683 z = self.xp.zeros_like(u)
685 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
687 # start newton iteration
688 n = 0
689 res = 99
690 while n < self.newton_maxiter:
691 # form the function g with g(u) = 0
692 g = u - factor * self.reaction(u) - b
694 # if g is close to 0, then we are done
695 res = self.xp.linalg.norm(g, self.xp.inf)
697 if res < self.newton_tol:
698 break
700 # assemble dg
701 dg = Id - factor * self.xsp.diags(self.reaction_prime(u), offsets=0)
703 # newton update: u1 = u0 - g/dg
704 # u -= spsolve(dg, g)
705 u -= self.linalg.cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
706 # increase iteration count
707 n += 1
708 # print(n, res)
710 self.work_counters['newton']()
712 # if n == self.newton_maxiter:
713 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
715 me = self.dtype_u(self.init)
716 me[:] = u.reshape(self.nvars)
718 self.newton_ncalls += 1
720 return me
723# noinspection PyUnusedLocal
724class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
725 r"""
726 2D periodic Allen-Cahn with finite differences, multi-implicit: Laplacian plus cubic term, linear term.
728 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions, with the two
729 phases at :math:`u = 0` and :math:`u = 1`
731 .. math::
732 \frac{\partial u}{\partial t} = \Delta u
733 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
735 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
736 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
738 The initial condition has the form of circles
740 .. math::
741 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
742 {\sqrt{2}\varepsilon}\right)\right)
744 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
745 is used here to get another kind of *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{2\varepsilon^2}(2u - 1)^3`
746 is handled implicitly and the nonlinear system including this part will be solved by Newton. The linear reaction
747 term :math:`\frac{1}{2\varepsilon^2}(2u - 1)` acts pointwise, so its implicit system is solved by a pointwise
748 division.
749 """
751 dtype_f = comp2_mesh
753 def eval_f(self, u, t):
754 """
755 Routine to evaluate the right-hand side of the problem.
757 Parameters
758 ----------
759 u : dtype_u
760 Current values of the numerical solution.
761 t : float
762 Current time of the numerical solution is computed.
764 Returns
765 -------
766 f : dtype_f
767 The right-hand side of the problem.
768 """
769 f = self.dtype_f(self.init)
770 v = u.flatten()
771 f.comp1[:] = (self.A.dot(v) + self.reaction_cubic(v)).reshape(self.nvars)
772 f.comp2[:] = self.reaction_linear(v).reshape(self.nvars)
774 self.work_counters['rhs']()
775 return f
777 def solve_system_1(self, rhs, factor, u0, t):
778 """
779 Simple Newton solver.
781 Parameters
782 ----------
783 rhs : dtype_f
784 Right-hand side for the nonlinear system.
785 factor : float
786 Abbrev. for the node-to-node stepsize (or any other factor required).
787 u0 : dtype_u
788 Initial guess for the iterative solver.
789 t : float
790 Current time (required here for the BC).
792 Returns
793 ------
794 me : dtype_u
795 The solution as mesh.
796 """
798 # plain arrays, not the datatype: its `__abs__` is a norm, which CuPy's `linalg.norm` trips over
799 u = u0.view(self.xp.ndarray).flatten()
800 b = rhs.view(self.xp.ndarray).flatten()
801 z = self.xp.zeros_like(u)
803 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
805 # start newton iteration
806 n = 0
807 res = 99
808 while n < self.newton_maxiter:
809 # form the function g with g(u) = 0
810 g = u - factor * (self.A.dot(u) + self.reaction_cubic(u)) - b
812 # if g is close to 0, then we are done
813 res = self.xp.linalg.norm(g, self.xp.inf)
815 if res < self.newton_tol:
816 break
818 # assemble dg
819 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_cubic_prime(u), offsets=0))
821 # newton update: u1 = u0 - g/dg
822 # u -= spsolve(dg, g)
823 u -= self.linalg.cg(
824 dg,
825 g,
826 x0=z,
827 rtol=self.lin_tol,
828 atol=0,
829 )[0]
830 # increase iteration count
831 n += 1
832 # print(n, res)
834 self.work_counters['newton']()
836 # if n == self.newton_maxiter:
837 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
839 me = self.dtype_u(self.init)
840 me[:] = u.reshape(self.nvars)
842 self.newton_ncalls += 1
844 return me
846 def solve_system_2(self, rhs, factor, u0, t):
847 r"""
848 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
850 Parameters
851 ----------
852 rhs : dtype_f
853 Right-hand side for the linear system.
854 factor : float
855 Abbrev. for the local stepsize (or any other factor required).
856 u0 : dtype_u
857 Initial guess for the iterative solver.
858 t : float
859 Current time (e.g. for time-dependent BCs).
861 Returns
862 -------
863 me : dtype_u
864 The solution as mesh.
865 """
867 me = self.dtype_u(self.init)
869 me[:] = ((rhs - 0.5 * factor / self.eps**2) / (1.0 - factor / self.eps**2)).reshape(self.nvars)
870 return me