Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD.py: 96%
208 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +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 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
17 phases at :math:`u = 0` and :math:`u = 1`
19 .. math::
20 \frac{\partial u}{\partial t} = \Delta u
21 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
23 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
24 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
26 Initial condition are circles of the form
28 .. math::
29 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
30 {\sqrt{2}\varepsilon}\right)\right)
32 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
33 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
35 Parameters
36 ----------
37 nvars : tuple of int, optional
38 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``.
39 nu : int, optional
40 Exponent of the double well; :math:`\nu = 2` is the standard Allen-Cahn nonlinearity.
41 eps : float, optional
42 Scaling parameter :math:`\varepsilon`.
43 newton_maxiter : int, optional
44 Maximum number of iterations for the Newton solver.
45 newton_tol : float, optional
46 Tolerance for Newton's method to terminate.
47 lin_tol : float, optional
48 Tolerance for linear solver to terminate.
49 lin_maxiter : int, optional
50 Maximum number of iterations for the linear solver.
51 radius : float, optional
52 Radius of the circles.
53 order : int, optional
54 Order of the finite difference matrix.
55 useGPU : bool, optional
56 Run on the GPU with CuPy instead of on the CPU with NumPy.
58 Attributes
59 ----------
60 A : scipy.spdiags
61 Second-order FD discretization of the 2D laplace operator.
62 dx : float
63 Distance between two spatial nodes (same for both directions).
64 xvalues : np.1darray
65 Spatial grid points, here both dimensions have the same grid points.
66 newton_ncalls : int
67 Number of calls of the Newton solver. The iterations themselves are counted in
68 ``work_counters['newton']``, and the linear ones in ``work_counters['linear']``.
69 lin_ncalls : int
70 Number of calls of the linear solver.
71 """
73 dtype_u = mesh
74 dtype_f = mesh
76 xp = np
77 xsp = sp
78 linalg = spla
80 def setup_GPU(self):
81 """
82 Switch the array, sparse and solver modules and the datatypes over to CuPy.
84 This changes the class, not the instance, as everything else in pySDC that does this
85 does: once one instance of a class runs on the GPU, they all do.
86 """
87 import cupy as cp
88 import cupyx.scipy.sparse as csp
89 import cupyx.scipy.sparse.linalg as cspla
90 from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh, comp2_cupy_mesh
92 self.xp = cp
93 self.xsp = csp
94 self.linalg = cspla
95 self.dtype_u = cupy_mesh
96 # .get, not [], because this runs once per instance and the class keeps what it is given
97 GPU_versions = {mesh: cupy_mesh, imex_mesh: imex_cupy_mesh, comp2_mesh: comp2_cupy_mesh}
98 self.dtype_f = GPU_versions.get(self.dtype_f, self.dtype_f)
100 def __init__(
101 self,
102 nvars=(128, 128),
103 nu=2,
104 eps=0.04,
105 newton_maxiter=200,
106 newton_tol=1e-12,
107 lin_tol=1e-8,
108 lin_maxiter=100,
109 inexact_linear_ratio=None,
110 radius=0.25,
111 order=2,
112 useGPU=False,
113 ):
114 """Initialization routine"""
115 if useGPU:
116 self.setup_GPU()
118 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
119 if len(nvars) != 2:
120 raise ProblemError('this is a 2d example, got %s' % nvars)
121 if nvars[0] != nvars[1]:
122 raise ProblemError('need a square domain, got %s' % nvars)
123 if nvars[0] % 2 != 0:
124 raise ProblemError('the setup requires nvars = 2^p per dimension')
126 # invoke super init, passing number of dofs, dtype_u and dtype_f
127 super().__init__((nvars, None, np.dtype('float64')))
128 self._makeAttributeAndRegister(
129 'nvars',
130 'nu',
131 'eps',
132 'radius',
133 'order',
134 'useGPU',
135 localVars=locals(),
136 readOnly=True,
137 )
138 self._makeAttributeAndRegister(
139 'newton_maxiter',
140 'newton_tol',
141 'lin_tol',
142 'lin_maxiter',
143 'inexact_linear_ratio',
144 localVars=locals(),
145 readOnly=False,
146 )
148 # compute dx and get discretization matrix A
149 self.dx = 1.0 / self.nvars[0]
150 self.A, _ = problem_helper.get_finite_difference_matrix(
151 derivative=2,
152 order=self.order,
153 stencil_type='center',
154 dx=self.dx,
155 size=self.nvars[0],
156 dim=2,
157 bc='periodic',
158 cupy=self.useGPU,
159 )
160 self.xvalues = self.xp.arange(self.nvars[0]) * self.dx - 0.5
162 self.newton_ncalls = 0
163 self.lin_ncalls = 0
165 self.work_counters['newton'] = WorkCounter()
166 self.work_counters['rhs'] = WorkCounter()
167 self.work_counters['linear'] = WorkCounter()
169 def reaction(self, u):
170 r"""
171 The reaction term, :math:`\frac{1}{2\varepsilon^2}(2u - 1)\left(1 - (2u - 1)^\nu\right)`.
173 The wells sit at :math:`u = 0` and :math:`u = 1`, so the double well is symmetric about
174 :math:`2u - 1`; writing the term in that variable is what lets :math:`\nu` keep the meaning
175 it has always had here. For the default :math:`\nu = 2` this is
176 :math:`-\frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
177 """
178 v = 2.0 * u - 1.0
179 return 0.5 / self.eps**2 * v * (1.0 - v**self.nu)
181 def reaction_prime(self, u):
182 """Derivative of :meth:`reaction`, ready to go on the diagonal of a Jacobian."""
183 v = 2.0 * u - 1.0
184 return 1.0 / self.eps**2 * (1.0 - (self.nu + 1.0) * v**self.nu)
186 def reaction_cubic(self, u):
187 r"""The stiff part of :meth:`reaction`, :math:`-\frac{1}{2\varepsilon^2}(2u - 1)^{\nu + 1}`."""
188 return -0.5 / self.eps**2 * (2.0 * u - 1.0) ** (self.nu + 1)
190 def reaction_cubic_prime(self, u):
191 """Derivative of :meth:`reaction_cubic`."""
192 return -(self.nu + 1.0) / self.eps**2 * (2.0 * u - 1.0) ** self.nu
194 def reaction_linear(self, u):
195 r"""The rest of :meth:`reaction`, :math:`\frac{1}{2\varepsilon^2}(2u - 1)`, so the two sum back to it."""
196 return 0.5 / self.eps**2 * (2.0 * u - 1.0)
198 # noinspection PyTypeChecker
199 def solve_system(self, rhs, factor, u0, t):
200 """
201 Simple Newton solver.
203 Parameters
204 ----------
205 rhs : dtype_f
206 Right-hand side for the nonlinear system
207 factor : float
208 Abbrev. for the node-to-node stepsize (or any other factor required).
209 u0 : dtype_u
210 Initial guess for the iterative solver.
211 t : float
212 Current time (required here for the BC).
214 Returns
215 -------
216 me : dtype_u
217 The solution as mesh.
218 """
220 u = self.dtype_u(u0).flatten()
221 z = self.dtype_u(self.init, val=0.0).flatten()
223 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
225 # start newton iteration
226 n = 0
227 res = 99
228 while n < self.newton_maxiter:
229 # form the function g with g(u) = 0
230 g = u - factor * (self.A.dot(u) + self.reaction(u)) - rhs.flatten()
232 # if g is close to 0, then we are done
233 res = self.xp.linalg.norm(g, self.xp.inf)
235 # do inexactness in the linear solver
236 if self.inexact_linear_ratio:
237 self.lin_tol = res * self.inexact_linear_ratio
239 if res < self.newton_tol:
240 break
242 # assemble dg
243 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_prime(u), offsets=0))
245 # newton update: u1 = u0 - g/dg
246 # u -= spsolve(dg, g)
247 u -= self.linalg.cg(
248 dg, g, x0=z, rtol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
249 )[0]
250 # increase iteration count
251 n += 1
252 # print(n, res)
254 self.work_counters['newton']()
256 # if n == self.newton_maxiter:
257 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
259 me = self.dtype_u(self.init)
260 me[:] = u.reshape(self.nvars)
262 self.newton_ncalls += 1
264 return me
266 def eval_f(self, u, t):
267 """
268 Routine to evaluate the right-hand side of the problem.
270 Parameters
271 ----------
272 u : dtype_u
273 Current values of the numerical solution.
274 t : float
275 Current time of the numerical solution is computed (not used here).
277 Returns
278 -------
279 f : dtype_f
280 The right-hand side of the problem.
281 """
282 f = self.dtype_f(self.init)
283 v = u.flatten()
284 f[:] = (self.A.dot(v) + self.reaction(v)).reshape(self.nvars)
286 self.work_counters['rhs']()
287 return f
289 def u_exact(self, t, u_init=None, t_init=None):
290 r"""
291 Routine to compute the exact solution at time :math:`t`.
293 Parameters
294 ----------
295 t : float
296 Time of the exact solution.
298 Returns
299 -------
300 me : dtype_u
301 The exact solution.
302 """
303 me = self.dtype_u(self.init, val=0.0)
304 if t > 0:
306 def eval_rhs(t, u):
307 return self.eval_f(u.reshape(self.init[0]), t).flatten()
309 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
311 else:
312 X, Y = self.xp.meshgrid(self.xvalues, self.xvalues)
313 r2 = X**2 + Y**2
314 me[:] = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)))
316 return me
319# noinspection PyUnusedLocal
320class allencahn_semiimplicit(allencahn_fullyimplicit):
321 r"""
322 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
323 phases at :math:`u = 0` and :math:`u = 1`
325 .. math::
326 \frac{\partial u}{\partial t} = \Delta u
327 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
329 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
330 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
332 Initial condition are circles of the form
334 .. math::
335 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
336 {\sqrt{2}\varepsilon}\right)\right)
338 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
339 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
340 method, and the system containing the rest of the right-hand side is only evaluated at each time.
341 """
343 dtype_f = imex_mesh
345 def eval_f(self, u, t):
346 """
347 Routine to evaluate the right-hand side of the problem.
349 Parameters
350 ----------
351 u : dtype_u
352 Current values of the numerical solution.
353 t : float
354 Current time of the numerical solution is computed (not used here).
356 Returns
357 -------
358 f : dtype_f
359 The right-hand side of the problem.
360 """
361 f = self.dtype_f(self.init)
362 v = u.flatten()
363 f.impl[:] = self.A.dot(v).reshape(self.nvars)
364 f.expl[:] = self.reaction(v).reshape(self.nvars)
366 self.work_counters['rhs']()
367 return f
369 def solve_system(self, rhs, factor, u0, t):
370 r"""
371 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
373 Parameters
374 ----------
375 rhs : dtype_f
376 Right-hand side for the linear system.
377 factor : float
378 Abbrev. for the local stepsize (or any other factor required).
379 u0 : dtype_u
380 Initial guess for the iterative solver.
381 t : float
382 Current time (e.g. for time-dependent BCs).
384 Returns
385 -------
386 me : dtype_u
387 The solution as mesh.
388 """
390 me = self.dtype_u(self.init)
392 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
394 me[:] = self.linalg.cg(
395 Id - factor * self.A,
396 rhs.flatten(),
397 x0=u0.flatten(),
398 rtol=self.lin_tol,
399 maxiter=self.lin_maxiter,
400 atol=0,
401 callback=self.work_counters['linear'],
402 )[0].reshape(self.nvars)
404 self.lin_ncalls += 1
406 return me
408 def u_exact(self, t, u_init=None, t_init=None):
409 """
410 Routine to compute the exact solution at time t.
412 Parameters
413 ----------
414 t : float
415 Time of the exact solution.
417 Returns
418 -------
419 me : dtype_u
420 The exact solution.
421 """
422 me = self.dtype_u(self.init, val=0.0)
423 if t > 0:
425 def eval_rhs(t, u):
426 f = self.eval_f(u.reshape(self.init[0]), t)
427 return (f.impl + f.expl).flatten()
429 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
430 else:
431 me[:] = super().u_exact(t, u_init, t_init)
432 return me
435# noinspection PyUnusedLocal
436class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
437 r"""
438 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions, with the two
439 phases at :math:`u = 0` and :math:`u = 1`
441 .. math::
442 \frac{\partial u}{\partial t} = \Delta u
443 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
445 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
446 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
448 Initial condition are circles of the form
450 .. math::
451 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
452 {\sqrt{2}\varepsilon}\right)\right)
454 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
455 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{2\varepsilon^2}(2u - 1)^3`
456 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{2\varepsilon^2}(2u - 1)`
457 is only evaluated at each time.
458 """
460 dtype_f = imex_mesh
462 def eval_f(self, u, t):
463 """
464 Routine to evaluate the right-hand side of the problem.
466 Parameters
467 ----------
468 u : dtype_u
469 Current values of the numerical solution.
470 t : float
471 Current time of the numerical solution is computed.
473 Returns
474 -------
475 f : dtype_f
476 The right-hand side of the problem.
477 """
478 f = self.dtype_f(self.init)
479 v = u.flatten()
480 f.impl[:] = (self.A.dot(v) + self.reaction_cubic(v)).reshape(self.nvars)
481 f.expl[:] = self.reaction_linear(v).reshape(self.nvars)
483 self.work_counters['rhs']()
484 return f
486 def solve_system(self, rhs, factor, u0, t):
487 """
488 Simple Newton solver.
490 Parameters
491 ----------
492 rhs : dtype_f
493 Right-hand side for the nonlinear system.
494 factor : float
495 Abbrev. for the node-to-node stepsize (or any other factor required).
496 u0 : dtype_u
497 Initial guess for the iterative solver.
498 t : float
499 Current time (required here for the BC).
501 Returns
502 -------
503 me : dtype_u
504 The solution as mesh.
505 """
507 u = self.dtype_u(u0).flatten()
508 z = self.dtype_u(self.init, val=0.0).flatten()
510 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
512 # start newton iteration
513 n = 0
514 res = 99
515 while n < self.newton_maxiter:
516 # form the function g with g(u) = 0
517 g = u - factor * (self.A.dot(u) + self.reaction_cubic(u)) - rhs.flatten()
519 # if g is close to 0, then we are done
520 res = self.xp.linalg.norm(g, self.xp.inf)
522 if res < self.newton_tol:
523 break
525 # assemble dg
526 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_cubic_prime(u), offsets=0))
528 # newton update: u1 = u0 - g/dg
529 # u -= spsolve(dg, g)
530 u -= self.linalg.cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
531 # increase iteration count
532 n += 1
533 # print(n, res)
535 self.work_counters['newton']()
537 # if n == self.newton_maxiter:
538 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
540 me = self.dtype_u(self.init)
541 me[:] = u.reshape(self.nvars)
543 self.newton_ncalls += 1
545 return me
548# noinspection PyUnusedLocal
549class allencahn_multiimplicit(allencahn_fullyimplicit):
550 r"""
551 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions, with the two
552 phases at :math:`u = 0` and :math:`u = 1`
554 .. math::
555 \frac{\partial u}{\partial t} = \Delta u
556 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
558 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
559 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
561 Initial condition are circles of the form
563 .. math::
564 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
565 {\sqrt{2}\varepsilon}\right)\right)
567 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
568 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
569 method, and the system containing the rest of the right-hand side will be solved by Newton's method.
570 """
572 dtype_f = comp2_mesh
574 def eval_f(self, u, t):
575 """
576 Routine to evaluate the right-hand side of the problem.
578 Parameters
579 ----------
580 u : dtype_u
581 Current values of the numerical solution.
582 t : float
583 Current time of the numerical solution is computed.
585 Returns
586 -------
587 f : dtype_f
588 The right-hand side of the problem.
589 """
590 f = self.dtype_f(self.init)
591 v = u.flatten()
592 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
593 f.comp2[:] = self.reaction(v).reshape(self.nvars)
595 self.work_counters['rhs']()
596 return f
598 def solve_system_1(self, rhs, factor, u0, t):
599 r"""
600 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
602 Parameters
603 ----------
604 rhs : dtype_f
605 Right-hand side for the linear system.
606 factor : float
607 Abbrev. for the local stepsize (or any other factor required).
608 u0 : dtype_u
609 Initial guess for the iterative solver.
610 t : float
611 Current time (e.g. for time-dependent BCs).
613 Returns
614 -------
615 me : dtype_u
616 The solution as mesh.
617 """
619 me = self.dtype_u(self.init)
621 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
623 me[:] = self.linalg.cg(
624 Id - factor * self.A,
625 rhs.flatten(),
626 x0=u0.flatten(),
627 rtol=self.lin_tol,
628 maxiter=self.lin_maxiter,
629 atol=0,
630 callback=self.work_counters['linear'],
631 )[0].reshape(self.nvars)
633 self.lin_ncalls += 1
635 return me
637 def solve_system_2(self, rhs, factor, u0, t):
638 """
639 Simple Newton solver.
641 Parameters
642 ----------
643 rhs : dtype_f
644 Right-hand side for the nonlinear system.
645 factor : float
646 Abbrev. for the node-to-node stepsize (or any other factor required).
647 u0 : dtype_u
648 Initial guess for the iterative solver.
649 t : float
650 Current time (required here for the BC).
652 Returns
653 -------
654 me : dtype_u
655 The solution as mesh.
656 """
658 u = self.dtype_u(u0).flatten()
659 z = self.dtype_u(self.init, val=0.0).flatten()
661 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
663 # start newton iteration
664 n = 0
665 res = 99
666 while n < self.newton_maxiter:
667 # form the function g with g(u) = 0
668 g = u - factor * self.reaction(u) - rhs.flatten()
670 # if g is close to 0, then we are done
671 res = self.xp.linalg.norm(g, self.xp.inf)
673 if res < self.newton_tol:
674 break
676 # assemble dg
677 dg = Id - factor * self.xsp.diags(self.reaction_prime(u), offsets=0)
679 # newton update: u1 = u0 - g/dg
680 # u -= spsolve(dg, g)
681 u -= self.linalg.cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
682 # increase iteration count
683 n += 1
684 # print(n, res)
686 self.work_counters['newton']()
688 # if n == self.newton_maxiter:
689 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
691 me = self.dtype_u(self.init)
692 me[:] = u.reshape(self.nvars)
694 self.newton_ncalls += 1
696 return me
699# noinspection PyUnusedLocal
700class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
701 r"""
702 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions, with the two
703 phases at :math:`u = 0` and :math:`u = 1`
705 .. math::
706 \frac{\partial u}{\partial t} = \Delta u
707 + \frac{1}{2\varepsilon^2} (2u - 1)\left(1 - (2u - 1)^\nu\right)
709 for a constant :math:`\nu`, which at the default :math:`\nu = 2` is the usual
710 :math:`\Delta u - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u)`.
712 The initial condition has the form of circles
714 .. math::
715 u({\bf x}, 0) = \frac{1}{2}\left(1 + \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}
716 {\sqrt{2}\varepsilon}\right)\right)
718 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
719 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`
720 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{2\varepsilon^2}(2u - 1)`
721 is solved by a linear solver provided by a ``SciPy`` routine.
722 """
724 dtype_f = comp2_mesh
726 def eval_f(self, u, t):
727 """
728 Routine to evaluate the right-hand side of the problem.
730 Parameters
731 ----------
732 u : dtype_u
733 Current values of the numerical solution.
734 t : float
735 Current time of the numerical solution is computed.
737 Returns
738 -------
739 f : dtype_f
740 The right-hand side of the problem.
741 """
742 f = self.dtype_f(self.init)
743 v = u.flatten()
744 f.comp1[:] = (self.A.dot(v) + self.reaction_cubic(v)).reshape(self.nvars)
745 f.comp2[:] = self.reaction_linear(v).reshape(self.nvars)
747 self.work_counters['rhs']()
748 return f
750 def solve_system_1(self, rhs, factor, u0, t):
751 """
752 Simple Newton solver.
754 Parameters
755 ----------
756 rhs : dtype_f
757 Right-hand side for the nonlinear system.
758 factor : float
759 Abbrev. for the node-to-node stepsize (or any other factor required).
760 u0 : dtype_u
761 Initial guess for the iterative solver.
762 t : float
763 Current time (required here for the BC).
765 Returns
766 ------
767 me : dtype_u
768 The solution as mesh.
769 """
771 u = self.dtype_u(u0).flatten()
772 z = self.dtype_u(self.init, val=0.0).flatten()
774 Id = self.xsp.eye(self.nvars[0] * self.nvars[1])
776 # start newton iteration
777 n = 0
778 res = 99
779 while n < self.newton_maxiter:
780 # form the function g with g(u) = 0
781 g = u - factor * (self.A.dot(u) + self.reaction_cubic(u)) - rhs.flatten()
783 # if g is close to 0, then we are done
784 res = self.xp.linalg.norm(g, self.xp.inf)
786 if res < self.newton_tol:
787 break
789 # assemble dg
790 dg = Id - factor * (self.A + self.xsp.diags(self.reaction_cubic_prime(u), offsets=0))
792 # newton update: u1 = u0 - g/dg
793 # u -= spsolve(dg, g)
794 u -= self.linalg.cg(
795 dg,
796 g,
797 x0=z,
798 rtol=self.lin_tol,
799 atol=0,
800 )[0]
801 # increase iteration count
802 n += 1
803 # print(n, res)
805 self.work_counters['newton']()
807 # if n == self.newton_maxiter:
808 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
810 me = self.dtype_u(self.init)
811 me[:] = u.reshape(self.nvars)
813 self.newton_ncalls += 1
815 return me
817 def solve_system_2(self, rhs, factor, u0, t):
818 r"""
819 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
821 Parameters
822 ----------
823 rhs : dtype_f
824 Right-hand side for the linear system.
825 factor : float
826 Abbrev. for the local stepsize (or any other factor required).
827 u0 : dtype_u
828 Initial guess for the iterative solver.
829 t : float
830 Current time (e.g. for time-dependent BCs).
832 Returns
833 -------
834 me : dtype_u
835 The solution as mesh.
836 """
838 me = self.dtype_u(self.init)
840 me[:] = ((rhs - 0.5 * factor / self.eps**2) / (1.0 - factor / self.eps**2)).reshape(self.nvars)
841 return me