Coverage for pySDC/implementations/problem_classes/GeneralizedFisher_1D_PETSc.py: 99%
243 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
2from petsc4py import PETSc
4from pySDC.core.problem import Problem
5from pySDC.implementations.datatype_classes.petsc_vec import petsc_vec, petsc_vec_imex, petsc_vec_comp2
8class Fisher_full(object):
9 r"""
10 Helper class to generate residual and Jacobian matrix for PETSc's nonlinear solver SNES.
12 Parameters
13 ----------
14 da : DMDA object
15 Object of ``PETSc``.
16 prob : problem instance
17 Contains problem information for ``PETSc``.
18 factor : float
19 Temporal factor :math:`\Delta t Q_\Delta`.
20 dx : float
21 Grid spacing in x direction.
23 Attributes
24 ----------
25 localX : PETSc vector object
26 Local vector for ``PETSc``.
27 xs, xe : int
28 Defines the ranges for spatial domain.
29 mx : tuple
30 Get sizes for the vector containing the spatial points.
31 row : PETSc matrix stencil object
32 Row for a matrix.
33 col : PETSc matrix stencil object
34 Column for a matrix.
35 """
37 def __init__(self, da, prob, factor, dx):
38 """Initialization routine"""
39 assert da.getDim() == 1
40 self.da = da
41 self.factor = factor
42 self.dx = dx
43 self.prob = prob
44 self.localX = da.createLocalVec()
45 self.xs, self.xe = self.da.getRanges()[0]
46 self.mx = self.da.getSizes()[0]
47 self.row = PETSc.Mat.Stencil()
48 self.col = PETSc.Mat.Stencil()
50 def formFunction(self, snes, X, F):
51 r"""
52 Function to evaluate the residual for the Newton solver. This function should be equal to the RHS
53 in the solution.
55 Parameters
56 ----------
57 snes : PETSc solver object
58 Nonlinear solver object.
59 X : PETSc vector object
60 Input vector.
61 F : PETSc vector object
62 Output vector :math:`F(X)`.
64 Returns
65 -------
66 None
67 Overwrites F.
68 """
69 self.da.globalToLocal(X, self.localX)
70 x = self.da.getVecArray(self.localX)
71 f = self.da.getVecArray(F)
73 for i in range(self.xs, self.xe):
74 if i == 0 or i == self.mx - 1:
75 f[i] = x[i]
76 else:
77 u = x[i] # center
78 u_e = x[i + 1] # east
79 u_w = x[i - 1] # west
80 u_xx = (u_e - 2 * u + u_w) / self.dx**2
81 f[i] = x[i] - self.factor * (u_xx + self.prob.lambda0**2 * x[i] * (1 - x[i] ** self.prob.nu))
83 def formJacobian(self, snes, X, J, P):
84 """
85 Function to return the Jacobian matrix.
87 Parameters
88 ----------
89 snes : PETSc solver object
90 Nonlinear solver object.
91 X : PETSc vector object
92 Input vector.
93 J : PETSc matrix object
94 Current Jacobian matrix.
95 P : PETSc matrix object
96 New Jacobian matrix.
98 Returns
99 -------
100 PETSc.Mat.Structure.SAME_NONZERO_PATTERN
101 Matrix status.
102 """
103 self.da.globalToLocal(X, self.localX)
104 x = self.da.getVecArray(self.localX)
105 P.zeroEntries()
107 for i in range(self.xs, self.xe):
108 self.row.i = i
109 self.row.field = 0
110 if i == 0 or i == self.mx - 1:
111 P.setValueStencil(self.row, self.row, 1.0)
112 else:
113 diag = 1.0 - self.factor * (
114 -2.0 / self.dx**2 + self.prob.lambda0**2 * (1.0 - (self.prob.nu + 1) * x[i] ** self.prob.nu)
115 )
116 for index, value in [
117 (i - 1, -self.factor / self.dx**2),
118 (i, diag),
119 (i + 1, -self.factor / self.dx**2),
120 ]:
121 self.col.i = index
122 self.col.field = 0
123 P.setValueStencil(self.row, self.col, value)
124 P.assemble()
125 if J != P:
126 J.assemble() # matrix-free operator
127 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN
130class Fisher_reaction(object):
131 r"""
132 Helper class to generate residual and Jacobian matrix for ``PETSc``'s nonlinear solver SNES.
134 Parameters
135 ----------
136 da : DMDA object
137 Object of ``PETSc``.
138 prob : problem instance
139 Contains problem information for ``PETSc``.
140 factor : float
141 Temporal factor :math:`\Delta t Q_\Delta`.
143 Attributes
144 ----------
145 localX : PETSc vector object
146 Local vector for ``PETSc``.
147 """
149 def __init__(self, da, prob, factor):
150 """Initialization routine"""
151 assert da.getDim() == 1
152 self.da = da
153 self.prob = prob
154 self.factor = factor
155 self.localX = da.createLocalVec()
157 def formFunction(self, snes, X, F):
158 r"""
159 Function to evaluate the residual for the Newton solver. This function should be equal to the RHS
160 in the solution.
162 Parameters
163 ----------
164 snes : PETSc solver object
165 Nonlinear solver object.
166 X : PETSc vector object
167 Input vector.
168 F : PETSc vector object
169 Output vector :math:`F(X)`.
171 Returns
172 -------
173 None
174 Overwrites F.
175 """
176 self.da.globalToLocal(X, self.localX)
177 x = self.da.getVecArray(self.localX)
178 f = self.da.getVecArray(F)
179 mx = self.da.getSizes()[0]
180 xs, xe = self.da.getRanges()[0]
181 for i in range(xs, xe):
182 if i == 0 or i == mx - 1:
183 f[i] = x[i]
184 else:
185 f[i] = x[i] - self.factor * self.prob.lambda0**2 * x[i] * (1 - x[i] ** self.prob.nu)
187 def formJacobian(self, snes, X, J, P):
188 """
189 Function to return the Jacobian matrix.
191 Parameters
192 ----------
193 snes : PETSc solver object
194 Nonlinear solver object.
195 X : PETSc vector object
196 Input vector.
197 J : PETSc matrix object
198 Current Jacobian matrix.
199 P : PETSc matrix object
200 New Jacobian matrix.
202 Returns
203 -------
204 PETSc.Mat.Structure.SAME_NONZERO_PATTERN
205 Matrix status.
206 """
207 self.da.globalToLocal(X, self.localX)
208 x = self.da.getVecArray(self.localX)
209 P.zeroEntries()
210 row = PETSc.Mat.Stencil()
211 mx = self.da.getSizes()[0]
212 xs, xe = self.da.getRanges()[0]
213 for i in range(xs, xe):
214 row.i = i
215 row.field = 0
216 if i == 0 or i == mx - 1:
217 P.setValueStencil(row, row, 1.0)
218 else:
219 diag = 1.0 - self.factor * self.prob.lambda0**2 * (1.0 - (self.prob.nu + 1) * x[i] ** self.prob.nu)
220 P.setValueStencil(row, row, diag)
221 P.assemble()
222 if J != P:
223 J.assemble() # matrix-free operator
224 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN
227class petsc_fisher_multiimplicit(Problem):
228 r"""
229 1D generalized Fisher equation with PETSc finite differences, multi-implicit: diffusion (CG), reaction (SNES).
231 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
232 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
233 problem [1]_ on a finite interval, where the first and last grid points keep their initial values (Dirichlet
234 boundary conditions)
236 .. math::
237 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
239 with exact solution
241 .. math::
242 u(x, 0) = \left[
243 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
244 \right]^{-2 / \nu}
246 for :math:`x \in \mathbb{R}`, and
248 .. math::
249 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
251 .. math::
252 \lambda_1 = \frac{\lambda_0}{2} \left[
253 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
254 \right].
256 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem will be solved
257 *multi-implicitly*.
259 Parameters
260 ----------
261 nvars : int, optional
262 Spatial resolution.
263 lambda0 : float, optional
264 Problem parameter : math:`\lambda_0`.
265 nu : float, optional
266 Problem parameter :math:`\nu`.
267 interval : tuple, optional
268 Defines the spatial domain.
269 comm : PETSc.COMM_WORLD, optional
270 Communicator for PETSc.
271 lsol_tol : float, optional
272 Tolerance for linear solver to terminate.
273 nlsol_tol : float, optional
274 Tolerance for nonlinear solver to terminate.
275 lsol_maxiter : int, optional
276 Maximum number of iterations for linear solver.
277 nlsol_maxiter : int, optional
278 Maximum number of iterations for nonlinear solver.
280 Attributes
281 ----------
282 dx : float
283 Mesh grid width.
284 xs, xe : int
285 Define the ranges.
286 A : PETSc matrix object
287 Discretization matrix.
288 localX : PETSc vector object
289 Local vector for solution.
290 ksp : PETSc solver object
291 Linear solver object.
292 snes : PETSc solver object
293 Nonlinear solver object.
294 F : PETSc vector object
295 Global vector.
296 J : PETSc matrix object
297 Jacobi matrix.
299 References
300 ----------
301 .. [1] Z. Feng. Traveling wave behavior for a generalized fisher equation. Chaos Solitons Fractals 38(2),
302 481–488 (2008).
303 .. [2] PETSc Web page. Satish Balay et al. https://petsc.org/ (2023).
304 .. [3] Parallel distributed computing using Python. Lisandro D. Dalcin, Rodrigo R. Paz, Pablo A. Kler,
305 Alejandro Cosimo. Advances in Water Resources (2011).
306 """
308 dtype_u = petsc_vec
309 dtype_f = petsc_vec_comp2
311 def __init__(
312 self,
313 nvars,
314 lambda0,
315 nu,
316 interval,
317 comm=PETSc.COMM_WORLD,
318 lsol_tol=1e-10,
319 nlsol_tol=1e-10,
320 lsol_maxiter=None,
321 nlsol_maxiter=None,
322 ):
323 """Initialization routine"""
324 # create DMDA object which will be used for all grid operations
325 da = PETSc.DMDA().create([nvars], dof=1, stencil_width=1, comm=comm)
327 # invoke super init, passing number of dofs, dtype_u and dtype_f
328 super().__init__(init=da)
329 self._makeAttributeAndRegister(
330 'nvars',
331 'lambda0',
332 'nu',
333 'interval',
334 'comm',
335 'lsol_tol',
336 'nlsol_tol',
337 'lsol_maxiter',
338 'nlsol_maxiter',
339 localVars=locals(),
340 readOnly=True,
341 )
343 # compute dx and get local ranges
344 self.dx = (self.interval[1] - self.interval[0]) / (self.nvars - 1)
345 self.xs, self.xe = self.init.getRanges()[0]
347 # compute discretization matrix A and identity
348 self.A = self.__get_A()
349 self.localX = self.init.createLocalVec()
351 # setup linear solver
352 self.ksp = PETSc.KSP()
353 self.ksp.create(comm=self.comm)
354 self.ksp.setType('cg')
355 pc = self.ksp.getPC()
356 pc.setType('ilu')
357 self.ksp.setInitialGuessNonzero(True)
358 self.ksp.setFromOptions()
359 self.ksp.setTolerances(rtol=self.lsol_tol, atol=self.lsol_tol, max_it=self.lsol_maxiter)
360 self.ksp_itercount = 0
361 self.ksp_ncalls = 0
363 # setup nonlinear solver
364 self.snes = PETSc.SNES()
365 self.snes.create(comm=self.comm)
366 if self.nlsol_maxiter <= 1:
367 self.snes.setType('ksponly')
368 self.snes.getKSP().setType('cg')
369 pc = self.snes.getKSP().getPC()
370 pc.setType('ilu')
371 # self.snes.setType('ngmres')
372 self.snes.setFromOptions()
373 self.snes.setTolerances(
374 rtol=self.nlsol_tol,
375 atol=self.nlsol_tol,
376 stol=self.nlsol_tol,
377 max_it=self.nlsol_maxiter,
378 )
379 self.snes_itercount = 0
380 self.snes_ncalls = 0
381 self.F = self.init.createGlobalVec()
382 self.J = self.init.createMatrix()
384 def __get_A(self):
385 r"""
386 Helper function to assemble ``PETSc`` matrix A.
388 Returns
389 -------
390 A : PETSc matrix object
391 Discretization matrix.
392 """
393 # create matrix and set basic options
394 A = self.init.createMatrix()
395 A.setType('aij') # sparse
396 A.setFromOptions()
397 A.setPreallocationNNZ((3, 3))
398 A.setUp()
400 # fill matrix
401 A.zeroEntries()
402 row = PETSc.Mat.Stencil()
403 col = PETSc.Mat.Stencil()
404 mx = self.init.getSizes()[0]
405 xs, xe = self.init.getRanges()[0]
406 for i in range(xs, xe):
407 row.i = i
408 row.field = 0
409 if i == 0 or i == mx - 1:
410 A.setValueStencil(row, row, 1.0)
411 else:
412 diag = -2.0 / self.dx**2
413 for index, value in [
414 (i - 1, 1.0 / self.dx**2),
415 (i, diag),
416 (i + 1, 1.0 / self.dx**2),
417 ]:
418 col.i = index
419 col.field = 0
420 A.setValueStencil(row, col, value)
421 A.assemble()
422 return A
424 def get_sys_mat(self, factor):
425 """
426 Helper function to assemble the system matrix of the linear problem.
428 Parameters
429 ----------
430 factor : float
431 Factor to define the system matrix.
433 Returns
434 -------
435 A : PETSc matrix object
436 Matrix for the system to solve.
437 """
439 # create matrix and set basic options
440 A = self.init.createMatrix()
441 A.setType('aij') # sparse
442 A.setFromOptions()
443 A.setPreallocationNNZ((3, 3))
444 A.setUp()
446 # fill matrix
447 A.zeroEntries()
448 row = PETSc.Mat.Stencil()
449 col = PETSc.Mat.Stencil()
450 mx = self.init.getSizes()[0]
451 xs, xe = self.init.getRanges()[0]
452 for i in range(xs, xe):
453 row.i = i
454 row.field = 0
455 if i == 0 or i == mx - 1:
456 A.setValueStencil(row, row, 1.0)
457 else:
458 diag = 1.0 + factor * 2.0 / self.dx**2
459 for index, value in [
460 (i - 1, -factor / self.dx**2),
461 (i, diag),
462 (i + 1, -factor / self.dx**2),
463 ]:
464 col.i = index
465 col.field = 0
466 A.setValueStencil(row, col, value)
467 A.assemble()
468 return A
470 def eval_f(self, u, t):
471 """
472 Routine to evaluate the right-hand side of the problem.
474 Parameters
475 ----------
476 u : dtype_u
477 Current values of the numerical solution.
478 t : float
479 Current time the numerical solution is computed.
481 Returns
482 -------
483 f : dtype_f
484 Right-hand side of the problem.
485 """
487 f = self.dtype_f(self.init)
488 self.A.mult(u, f.comp1)
489 fa1 = self.init.getVecArray(f.comp1)
490 fa1[0] = 0
491 fa1[-1] = 0
493 fa2 = self.init.getVecArray(f.comp2)
494 xa = self.init.getVecArray(u)
495 for i in range(self.xs, self.xe):
496 fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
497 fa2[0] = 0
498 fa2[-1] = 0
500 return f
502 def solve_system_1(self, rhs, factor, u0, t):
503 r"""
504 Linear solver for :math:`(I - factor \cdot A)\vec{u} = \vec{rhs}`.
506 Parameters
507 ----------
508 rhs : dtype_f
509 Right-hand side for the linear system.
510 factor : float
511 Abbrev. for the local stepsize (or any other factor required).
512 u0 : dtype_u
513 Initial guess for the iterative solver.
514 t : float
515 Current time (e.g. for time-dependent BCs).
517 Returns
518 -------
519 me : dtype_u
520 Solution as mesh.
521 """
523 me = self.dtype_u(u0)
525 self.ksp.setOperators(self.get_sys_mat(factor))
526 self.ksp.solve(rhs, me)
528 self.ksp_itercount += self.ksp.getIterationNumber()
529 self.ksp_ncalls += 1
531 return me
533 def solve_system_2(self, rhs, factor, u0, t):
534 r"""
535 Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
537 Parameters
538 ----------
539 rhs : dtype_f
540 Right-hand side for the linear system.
541 factor : float
542 Abbrev. for the local stepsize (or any other factor required).
543 u0 : dtype_u
544 Initial guess for the iterative solver.
545 t : float
546 Current time (e.g. for time-dependent BCs).
548 Returns
549 -------
550 me : dtype_u
551 Solution as mesh.
552 """
554 me = self.dtype_u(u0)
555 target = Fisher_reaction(self.init, self, factor)
557 # assign residual function and Jacobian
558 F = self.init.createGlobalVec()
559 self.snes.setFunction(target.formFunction, F)
560 J = self.init.createMatrix()
561 self.snes.setJacobian(target.formJacobian, J)
563 self.snes.solve(rhs, me)
565 self.snes_itercount += self.snes.getIterationNumber()
566 self.snes_ncalls += 1
568 return me
570 def u_exact(self, t):
571 r"""
572 Routine to compute the exact solution at time :math:`t`.
574 Parameters
575 ----------
576 t : float
577 Time of the exact solution.
579 Returns
580 -------
581 me : dtype_u
582 Exact solution.
583 """
585 lam1 = self.lambda0 / 2.0 * ((self.nu / 2.0 + 1) ** 0.5 + (self.nu / 2.0 + 1) ** (-0.5))
586 sig1 = lam1 - np.sqrt(lam1**2 - self.lambda0**2)
587 me = self.dtype_u(self.init)
588 xa = self.init.getVecArray(me)
589 for i in range(self.xs, self.xe):
590 xa[i] = (
591 1
592 + (2 ** (self.nu / 2.0) - 1)
593 * np.exp(-self.nu / 2.0 * sig1 * (self.interval[0] + (i + 1) * self.dx + 2 * lam1 * t))
594 ) ** (-2.0 / self.nu)
596 return me
599class petsc_fisher_fullyimplicit(petsc_fisher_multiimplicit):
600 r"""
601 1D generalized Fisher equation with PETSc finite differences, fully implicit with PETSc's SNES.
603 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
604 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
605 problem [1]_ on a finite interval, where the first and last grid points keep their initial values (Dirichlet
606 boundary conditions)
608 .. math::
609 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
611 with exact solution
613 .. math::
614 u(x, 0) = \left[
615 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
616 \right]^{-2 / \nu}
618 for :math:`x \in \mathbb{R}`, and
620 .. math::
621 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
623 .. math::
624 \lambda_1 = \frac{\lambda_0}{2} \left[
625 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
626 \right].
628 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem is treated
629 *fully-implicitly*.
630 """
632 dtype_f = petsc_vec
634 def eval_f(self, u, t):
635 """
636 Routine to evaluate the right-hand side of the problem.
638 Parameters
639 ----------
640 u : dtype_u
641 Current values of the numerical solution.
642 t : float
643 Current time the numerical solution is computed.
645 Returns
646 -------
647 f : dtype_f
648 Right-hand side of the problem.
649 """
651 f = self.dtype_f(self.init)
652 self.A.mult(u, f)
654 fa2 = self.init.getVecArray(f)
655 xa = self.init.getVecArray(u)
656 for i in range(self.xs, self.xe):
657 fa2[i] += self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
658 fa2[0] = 0
659 fa2[-1] = 0
661 return f
663 def solve_system(self, rhs, factor, u0, t):
664 r"""
665 Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
667 Parameters
668 ----------
669 rhs : dtype_f
670 Right-hand side for the linear system.
671 factor : float
672 Abbrev. for the local stepsize (or any other factor required).
673 u0 : dtype_u
674 Initial guess for the iterative solver.
675 t : float
676 Current time (e.g. for time-dependent BCs).
678 Returns
679 -------
680 me : dtype_u
681 Solution as mesh.
682 """
684 me = self.dtype_u(u0)
685 target = Fisher_full(self.init, self, factor, self.dx)
687 # assign residual function and Jacobian
689 self.snes.setFunction(target.formFunction, self.F)
690 self.snes.setJacobian(target.formJacobian, self.J)
692 self.snes.solve(rhs, me)
694 self.snes_itercount += self.snes.getIterationNumber()
695 self.snes_ncalls += 1
697 return me
700class petsc_fisher_semiimplicit(petsc_fisher_multiimplicit):
701 r"""
702 1D generalized Fisher equation with PETSc finite differences, IMEX with diffusion implicit, reaction explicit.
704 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
705 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
706 problem [1]_ on a finite interval, where the first and last grid points keep their initial values (Dirichlet
707 boundary conditions)
709 .. math::
710 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
712 with exact solution
714 .. math::
715 u(x, 0) = \left[
716 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
717 \right]^{-2 / \nu}
719 for :math:`x \in \mathbb{R}`, and
721 .. math::
722 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
724 .. math::
725 \lambda_1 = \frac{\lambda_0}{2} \left[
726 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
727 \right].
729 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem here will be
730 solved in a *semi-implicit* way.
731 """
733 dtype_f = petsc_vec_imex
735 def eval_f(self, u, t):
736 """
737 Routine to evaluate the right-hand side of the problem.
739 Parameters
740 ----------
741 u : dtype_u
742 Current values of the numerical solution.
743 t : float
744 Current time the numerical solution is computed.
746 Returns
747 -------
748 f : dtype_f
749 Right-hand side of the problem.
750 """
752 f = self.dtype_f(self.init)
753 self.A.mult(u, f.impl)
754 fa1 = self.init.getVecArray(f.impl)
755 fa1[0] = 0
756 fa1[-1] = 0
758 fa2 = self.init.getVecArray(f.expl)
759 xa = self.init.getVecArray(u)
760 for i in range(self.xs, self.xe):
761 fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
762 fa2[0] = 0
763 fa2[-1] = 0
765 return f
767 def solve_system(self, rhs, factor, u0, t):
768 r"""
769 Simple linear solver for :math:`(I-factor \cdot A)\vec{u} = \vec{rhs}`.
771 Parameters
772 ----------
773 rhs : dtype_f
774 Right-hand side for the linear system.
775 factor : float
776 Abbrev. for the local stepsize (or any other factor required).
777 u0 : dtype_u
778 Initial guess for the iterative solver.
779 t : float
780 Current time (e.g. for time-dependent BCs).
782 Returns
783 -------
784 me : dtype_u
785 Solution as mesh.
786 """
788 me = self.dtype_u(u0)
790 self.ksp.setOperators(self.get_sys_mat(factor))
791 self.ksp.solve(rhs, me)
793 self.ksp_itercount += self.ksp.getIterationNumber()
794 self.ksp_ncalls += 1
796 return me