Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD.py: 92%
213 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2import scipy.sparse as sp
3from scipy.sparse.linalg import cg
5from pySDC.core.Errors import ParameterError, ProblemError
6from pySDC.core.Problem import ptype, WorkCounter
7from pySDC.helpers import problem_helper
8from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
11# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
14# noinspection PyUnusedLocal
15class allencahn_fullyimplicit(ptype):
16 r"""
17 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
19 .. math::
20 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
22 for constant parameter :math:`\nu`. Initial condition are circles of the form
24 .. math::
25 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
27 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
28 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
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.
48 order : int, optional
49 Order of the finite difference matrix.
51 Attributes
52 ----------
53 A : scipy.spdiags
54 Second-order FD discretization of the 2D laplace operator.
55 dx : float
56 Distance between two spatial nodes (same for both directions).
57 xvalues : np.1darray
58 Spatial grid points, here both dimensions have the same grid points.
59 newton_itercount : int
60 Number of iterations of Newton solver.
61 lin_itercount
62 Number of iterations of linear solver.
63 newton_ncalls : int
64 Number of calls of Newton solver.
65 lin_ncalls : int
66 Number of calls of linear solver.
67 """
69 dtype_u = mesh
70 dtype_f = mesh
72 def __init__(
73 self,
74 nvars=(128, 128),
75 nu=2,
76 eps=0.04,
77 newton_maxiter=200,
78 newton_tol=1e-12,
79 lin_tol=1e-8,
80 lin_maxiter=100,
81 inexact_linear_ratio=None,
82 radius=0.25,
83 order=2,
84 ):
85 """Initialization routine"""
86 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
87 if len(nvars) != 2:
88 raise ProblemError('this is a 2d example, got %s' % nvars)
89 if nvars[0] != nvars[1]:
90 raise ProblemError('need a square domain, got %s' % nvars)
91 if nvars[0] % 2 != 0:
92 raise ProblemError('the setup requires nvars = 2^p per dimension')
94 # invoke super init, passing number of dofs, dtype_u and dtype_f
95 super().__init__((nvars, None, np.dtype('float64')))
96 self._makeAttributeAndRegister(
97 'nvars',
98 'nu',
99 'eps',
100 'radius',
101 'order',
102 localVars=locals(),
103 readOnly=True,
104 )
105 self._makeAttributeAndRegister(
106 'newton_maxiter',
107 'newton_tol',
108 'lin_tol',
109 'lin_maxiter',
110 'inexact_linear_ratio',
111 localVars=locals(),
112 readOnly=False,
113 )
115 # compute dx and get discretization matrix A
116 self.dx = 1.0 / self.nvars[0]
117 self.A, _ = problem_helper.get_finite_difference_matrix(
118 derivative=2,
119 order=self.order,
120 stencil_type='center',
121 dx=self.dx,
122 size=self.nvars[0],
123 dim=2,
124 bc='periodic',
125 )
126 self.xvalues = np.array([i * self.dx - 0.5 for i in range(self.nvars[0])])
128 self.newton_itercount = 0
129 self.lin_itercount = 0
130 self.newton_ncalls = 0
131 self.lin_ncalls = 0
133 self.work_counters['newton'] = WorkCounter()
134 self.work_counters['rhs'] = WorkCounter()
135 self.work_counters['linear'] = WorkCounter()
137 @staticmethod
138 def __get_A(N, dx):
139 """
140 Helper function to assemble FD matrix A in sparse format.
142 Parameters
143 ----------
144 N : list
145 Number of degrees of freedom.
146 dx : float
147 Distance between two spatial nodes.
149 Returns
150 -------
151 A : scipy.sparse.csc_matrix
152 Matrix in CSC format.
153 """
155 stencil = [1, -2, 1]
156 zero_pos = 2
158 dstencil = np.concatenate((stencil, np.delete(stencil, zero_pos - 1)))
159 offsets = np.concatenate(
160 (
161 [N[0] - i - 1 for i in reversed(range(zero_pos - 1))],
162 [i - zero_pos + 1 for i in range(zero_pos - 1, len(stencil))],
163 )
164 )
165 doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N[0]))
167 A = sp.diags(dstencil, doffsets, shape=(N[0], N[0]), format='csc')
168 A = sp.kron(A, sp.eye(N[0])) + sp.kron(sp.eye(N[1]), A)
169 A *= 1.0 / (dx**2)
170 return A
172 # noinspection PyTypeChecker
173 def solve_system(self, rhs, factor, u0, t):
174 """
175 Simple Newton solver.
177 Parameters
178 ----------
179 rhs : dtype_f
180 Right-hand side for the nonlinear system
181 factor : float
182 Abbrev. for the node-to-node stepsize (or any other factor required).
183 u0 : dtype_u
184 Initial guess for the iterative solver.
185 t : float
186 Current time (required here for the BC).
188 Returns
189 -------
190 me : dtype_u
191 The solution as mesh.
192 """
194 u = self.dtype_u(u0).flatten()
195 z = self.dtype_u(self.init, val=0.0).flatten()
196 nu = self.nu
197 eps2 = self.eps**2
199 Id = sp.eye(self.nvars[0] * self.nvars[1])
201 # start newton iteration
202 n = 0
203 res = 99
204 while n < self.newton_maxiter:
205 # form the function g with g(u) = 0
206 g = u - factor * (self.A.dot(u) + 1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
208 # if g is close to 0, then we are done
209 res = np.linalg.norm(g, np.inf)
211 # do inexactness in the linear solver
212 if self.inexact_linear_ratio:
213 self.lin_tol = res * self.inexact_linear_ratio
215 if res < self.newton_tol:
216 break
218 # assemble dg
219 dg = Id - factor * (self.A + 1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
221 # newton update: u1 = u0 - g/dg
222 # u -= spsolve(dg, g)
223 u -= cg(
224 dg, g, x0=z, tol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
225 )[0]
226 # increase iteration count
227 n += 1
228 # print(n, res)
230 self.work_counters['newton']()
232 # if n == self.newton_maxiter:
233 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
235 me = self.dtype_u(self.init)
236 me[:] = u.reshape(self.nvars)
238 self.newton_ncalls += 1
239 self.newton_itercount += n
241 return me
243 def eval_f(self, u, t):
244 """
245 Routine to evaluate the right-hand side of the problem.
247 Parameters
248 ----------
249 u : dtype_u
250 Current values of the numerical solution.
251 t : float
252 Current time of the numerical solution is computed (not used here).
254 Returns
255 -------
256 f : dtype_f
257 The right-hand side of the problem.
258 """
259 f = self.dtype_f(self.init)
260 v = u.flatten()
261 f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
263 self.work_counters['rhs']()
264 return f
266 def u_exact(self, t, u_init=None, t_init=None):
267 r"""
268 Routine to compute the exact solution at time :math:`t`.
270 Parameters
271 ----------
272 t : float
273 Time of the exact solution.
275 Returns
276 -------
277 me : dtype_u
278 The exact solution.
279 """
280 me = self.dtype_u(self.init, val=0.0)
281 if t > 0:
283 def eval_rhs(t, u):
284 return self.eval_f(u.reshape(self.init[0]), t).flatten()
286 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
288 else:
289 for i in range(self.nvars[0]):
290 for j in range(self.nvars[1]):
291 r2 = self.xvalues[i] ** 2 + self.xvalues[j] ** 2
292 me[i, j] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))
294 return me
297# noinspection PyUnusedLocal
298class allencahn_semiimplicit(allencahn_fullyimplicit):
299 r"""
300 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
302 .. math::
303 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
305 for constant parameter :math:`\nu`. Initial condition are circles of the form
307 .. math::
308 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
310 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
311 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
312 method, and the system containing the rest of the right-hand side is only evaluated at each time.
313 """
315 dtype_f = imex_mesh
317 def eval_f(self, u, t):
318 """
319 Routine to evaluate the right-hand side of the problem.
321 Parameters
322 ----------
323 u : dtype_u
324 Current values of the numerical solution.
325 t : float
326 Current time of the numerical solution is computed (not used here).
328 Returns
329 -------
330 f : dtype_f
331 The right-hand side of the problem.
332 """
333 f = self.dtype_f(self.init)
334 v = u.flatten()
335 f.impl[:] = self.A.dot(v).reshape(self.nvars)
336 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
338 self.work_counters['rhs']()
339 return f
341 def solve_system(self, rhs, factor, u0, t):
342 r"""
343 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
345 Parameters
346 ----------
347 rhs : dtype_f
348 Right-hand side for the linear system.
349 factor : float
350 Abbrev. for the local stepsize (or any other factor required).
351 u0 : dtype_u
352 Initial guess for the iterative solver.
353 t : float
354 Current time (e.g. for time-dependent BCs).
356 Returns
357 -------
358 me : dtype_u
359 The solution as mesh.
360 """
362 class context:
363 num_iter = 0
365 def callback(xk):
366 context.num_iter += 1
367 self.work_counters['linear']()
368 return context.num_iter
370 me = self.dtype_u(self.init)
372 Id = sp.eye(self.nvars[0] * self.nvars[1])
374 me[:] = cg(
375 Id - factor * self.A,
376 rhs.flatten(),
377 x0=u0.flatten(),
378 tol=self.lin_tol,
379 maxiter=self.lin_maxiter,
380 atol=0,
381 callback=callback,
382 )[0].reshape(self.nvars)
384 self.lin_ncalls += 1
385 self.lin_itercount += context.num_iter
387 return me
389 def u_exact(self, t, u_init=None, t_init=None):
390 """
391 Routine to compute the exact solution at time t.
393 Parameters
394 ----------
395 t : float
396 Time of the exact solution.
398 Returns
399 -------
400 me : dtype_u
401 The exact solution.
402 """
403 me = self.dtype_u(self.init, val=0.0)
404 if t > 0:
406 def eval_rhs(t, u):
407 f = self.eval_f(u.reshape(self.init[0]), t)
408 return (f.impl + f.expl).flatten()
410 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
411 else:
412 me[:] = super().u_exact(t, u_init, t_init)
413 return me
416# noinspection PyUnusedLocal
417class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
418 r"""
419 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
421 .. math::
422 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
424 for constant parameter :math:`\nu`. Initial condition are circles of the form
426 .. math::
427 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
429 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
430 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
431 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
432 is only evaluated at each time.
433 """
435 dtype_f = imex_mesh
437 def eval_f(self, u, t):
438 """
439 Routine to evaluate the right-hand side of the problem.
441 Parameters
442 ----------
443 u : dtype_u
444 Current values of the numerical solution.
445 t : float
446 Current time of the numerical solution is computed.
448 Returns
449 -------
450 f : dtype_f
451 The right-hand side of the problem.
452 """
453 f = self.dtype_f(self.init)
454 v = u.flatten()
455 f.impl[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
456 f.expl[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
458 return f
460 def solve_system(self, rhs, factor, u0, t):
461 """
462 Simple Newton solver.
464 Parameters
465 ----------
466 rhs : dtype_f
467 Right-hand side for the nonlinear system.
468 factor : float
469 Abbrev. for the node-to-node stepsize (or any other factor required).
470 u0 : dtype_u
471 Initial guess for the iterative solver.
472 t : float
473 Current time (required here for the BC).
475 Returns
476 -------
477 me : dtype_u
478 The solution as mesh.
479 """
481 u = self.dtype_u(u0).flatten()
482 z = self.dtype_u(self.init, val=0.0).flatten()
483 nu = self.nu
484 eps2 = self.eps**2
486 Id = sp.eye(self.nvars[0] * self.nvars[1])
488 # start newton iteration
489 n = 0
490 res = 99
491 while n < self.newton_maxiter:
492 # form the function g with g(u) = 0
493 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
495 # if g is close to 0, then we are done
496 # res = np.linalg.norm(g, np.inf)
497 res = np.linalg.norm(g, np.inf)
499 if res < self.newton_tol:
500 break
502 # assemble dg
503 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
505 # newton update: u1 = u0 - g/dg
506 # u -= spsolve(dg, g)
507 u -= cg(dg, g, x0=z, tol=self.lin_tol, atol=0)[0]
508 # increase iteration count
509 n += 1
510 # print(n, res)
512 # if n == self.newton_maxiter:
513 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
515 me = self.dtype_u(self.init)
516 me[:] = u.reshape(self.nvars)
518 self.newton_ncalls += 1
519 self.newton_itercount += n
521 return me
524# noinspection PyUnusedLocal
525class allencahn_multiimplicit(allencahn_fullyimplicit):
526 r"""
527 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
529 .. math::
530 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
532 for constant parameter :math:`\nu`. Initial condition are circles of the form
534 .. math::
535 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
537 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
538 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
539 method, and the system containing the rest of the right-hand side will be solved by Newton's method.
540 """
542 dtype_f = comp2_mesh
544 def eval_f(self, u, t):
545 """
546 Routine to evaluate the right-hand side of the problem.
548 Parameters
549 ----------
550 u : dtype_u
551 Current values of the numerical solution.
552 t : float
553 Current time of the numerical solution is computed.
555 Returns
556 -------
557 f : dtype_f
558 The right-hand side of the problem.
559 """
560 f = self.dtype_f(self.init)
561 v = u.flatten()
562 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
563 f.comp2[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
565 return f
567 def solve_system_1(self, rhs, factor, u0, t):
568 r"""
569 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
571 Parameters
572 ----------
573 rhs : dtype_f
574 Right-hand side for the linear system.
575 factor : float
576 Abbrev. for the local stepsize (or any other factor required).
577 u0 : dtype_u
578 Initial guess for the iterative solver.
579 t : float
580 Current time (e.g. for time-dependent BCs).
582 Returns
583 -------
584 me : dtype_u
585 The solution as mesh.
586 """
588 class context:
589 num_iter = 0
591 def callback(xk):
592 context.num_iter += 1
593 return context.num_iter
595 me = self.dtype_u(self.init)
597 Id = sp.eye(self.nvars[0] * self.nvars[1])
599 me[:] = cg(
600 Id - factor * self.A,
601 rhs.flatten(),
602 x0=u0.flatten(),
603 tol=self.lin_tol,
604 maxiter=self.lin_maxiter,
605 atol=0,
606 callback=callback,
607 )[0].reshape(self.nvars)
609 self.lin_ncalls += 1
610 self.lin_itercount += context.num_iter
612 return me
614 def solve_system_2(self, rhs, factor, u0, t):
615 """
616 Simple Newton solver.
618 Parameters
619 ----------
620 rhs : dtype_f
621 Right-hand side for the nonlinear system.
622 factor : float
623 Abbrev. for the node-to-node stepsize (or any other factor required).
624 u0 : dtype_u
625 Initial guess for the iterative solver.
626 t : float
627 Current time (required here for the BC).
629 Returns
630 -------
631 me : dtype_u
632 The solution as mesh.
633 """
635 u = self.dtype_u(u0).flatten()
636 z = self.dtype_u(self.init, val=0.0).flatten()
637 nu = self.nu
638 eps2 = self.eps**2
640 Id = sp.eye(self.nvars[0] * self.nvars[1])
642 # start newton iteration
643 n = 0
644 res = 99
645 while n < self.newton_maxiter:
646 # form the function g with g(u) = 0
647 g = u - factor * (1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
649 # if g is close to 0, then we are done
650 res = np.linalg.norm(g, np.inf)
652 if res < self.newton_tol:
653 break
655 # assemble dg
656 dg = Id - factor * (1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
658 # newton update: u1 = u0 - g/dg
659 # u -= spsolve(dg, g)
660 u -= cg(dg, g, x0=z, tol=self.lin_tol, atol=0)[0]
661 # increase iteration count
662 n += 1
663 # print(n, res)
665 # if n == self.newton_maxiter:
666 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
668 me = self.dtype_u(self.init)
669 me[:] = u.reshape(self.nvars)
671 self.newton_ncalls += 1
672 self.newton_itercount += n
674 return me
677# noinspection PyUnusedLocal
678class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
679 r"""
680 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
682 .. math::
683 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
685 for constant parameter :math:`\nu`. The initial condition has the form of circles
687 .. math::
688 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
690 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
691 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}`
692 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
693 is solved by a linear solver provided by a ``SciPy`` routine.
694 """
696 dtype_f = comp2_mesh
698 def eval_f(self, u, t):
699 """
700 Routine to evaluate the right-hand side of the problem.
702 Parameters
703 ----------
704 u : dtype_u
705 Current values of the numerical solution.
706 t : float
707 Current time of the numerical solution is computed.
709 Returns
710 -------
711 f : dtype_f
712 The right-hand side of the problem.
713 """
714 f = self.dtype_f(self.init)
715 v = u.flatten()
716 f.comp1[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
717 f.comp2[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
719 return f
721 def solve_system_1(self, rhs, factor, u0, t):
722 """
723 Simple Newton solver.
725 Parameters
726 ----------
727 rhs : dtype_f
728 Right-hand side for the nonlinear system.
729 factor : float
730 Abbrev. for the node-to-node stepsize (or any other factor required).
731 u0 : dtype_u
732 Initial guess for the iterative solver.
733 t : float
734 Current time (required here for the BC).
736 Returns
737 ------
738 me : dtype_u
739 The solution as mesh.
740 """
742 u = self.dtype_u(u0).flatten()
743 z = self.dtype_u(self.init, val=0.0).flatten()
744 nu = self.nu
745 eps2 = self.eps**2
747 Id = sp.eye(self.nvars[0] * self.nvars[1])
749 # start newton iteration
750 n = 0
751 res = 99
752 while n < self.newton_maxiter:
753 # form the function g with g(u) = 0
754 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
756 # if g is close to 0, then we are done
757 res = np.linalg.norm(g, np.inf)
759 if res < self.newton_tol:
760 break
762 # assemble dg
763 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
765 # newton update: u1 = u0 - g/dg
766 # u -= spsolve(dg, g)
767 u -= cg(
768 dg,
769 g,
770 x0=z,
771 tol=self.lin_tol,
772 atol=0,
773 )[0]
774 # increase iteration count
775 n += 1
776 # print(n, res)
778 # if n == self.newton_maxiter:
779 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
781 me = self.dtype_u(self.init)
782 me[:] = u.reshape(self.nvars)
784 self.newton_ncalls += 1
785 self.newton_itercount += n
787 return me
789 def solve_system_2(self, rhs, factor, u0, t):
790 r"""
791 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
793 Parameters
794 ----------
795 rhs : dtype_f
796 Right-hand side for the linear system.
797 factor : float
798 Abbrev. for the local stepsize (or any other factor required).
799 u0 : dtype_u
800 Initial guess for the iterative solver.
801 t : float
802 Current time (e.g. for time-dependent BCs).
804 Returns
805 -------
806 me : dtype_u
807 The solution as mesh.
808 """
810 me = self.dtype_u(self.init)
812 me[:] = (1.0 / (1.0 - factor * 1.0 / self.eps**2) * rhs).reshape(self.nvars)
813 return me