Coverage for pySDC/implementations/problem_classes/GeneralizedFisher_1D_PETSc.py: 99%
243 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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`.
142 dx : float
143 Grid spacing in x direction.
145 Attributes
146 ----------
147 localX : PETSc vector object
148 Local vector for ``PETSc``.
149 """
151 def __init__(self, da, prob, factor):
152 """Initialization routine"""
153 assert da.getDim() == 1
154 self.da = da
155 self.prob = prob
156 self.factor = factor
157 self.localX = da.createLocalVec()
159 def formFunction(self, snes, X, F):
160 r"""
161 Function to evaluate the residual for the Newton solver. This function should be equal to the RHS
162 in the solution.
164 Parameters
165 ----------
166 snes : PETSc solver object
167 Nonlinear solver object.
168 X : PETSc vector object
169 Input vector.
170 F : PETSc vector object
171 Output vector :math:`F(X)`.
173 Returns
174 -------
175 None
176 Overwrites F.
177 """
178 self.da.globalToLocal(X, self.localX)
179 x = self.da.getVecArray(self.localX)
180 f = self.da.getVecArray(F)
181 mx = self.da.getSizes()[0]
182 (xs, xe) = self.da.getRanges()[0]
183 for i in range(xs, xe):
184 if i == 0 or i == mx - 1:
185 f[i] = x[i]
186 else:
187 f[i] = x[i] - self.factor * self.prob.lambda0**2 * x[i] * (1 - x[i] ** self.prob.nu)
189 def formJacobian(self, snes, X, J, P):
190 """
191 Function to return the Jacobian matrix.
193 Parameters
194 ----------
195 snes : PETSc solver object
196 Nonlinear solver object.
197 X : PETSc vector object
198 Input vector.
199 J : PETSc matrix object
200 Current Jacobian matrix.
201 P : PETSc matrix object
202 New Jacobian matrix.
204 Returns
205 -------
206 PETSc.Mat.Structure.SAME_NONZERO_PATTERN
207 Matrix status.
208 """
209 self.da.globalToLocal(X, self.localX)
210 x = self.da.getVecArray(self.localX)
211 P.zeroEntries()
212 row = PETSc.Mat.Stencil()
213 mx = self.da.getSizes()[0]
214 (xs, xe) = self.da.getRanges()[0]
215 for i in range(xs, xe):
216 row.i = i
217 row.field = 0
218 if i == 0 or i == mx - 1:
219 P.setValueStencil(row, row, 1.0)
220 else:
221 diag = 1.0 - self.factor * self.prob.lambda0**2 * (1.0 - (self.prob.nu + 1) * x[i] ** self.prob.nu)
222 P.setValueStencil(row, row, diag)
223 P.assemble()
224 if J != P:
225 J.assemble() # matrix-free operator
226 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN
229class petsc_fisher_multiimplicit(Problem):
230 r"""
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]_ using periodic boundary conditions
235 .. math::
236 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
238 with exact solution
240 .. math::
241 u(x, 0) = \left[
242 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
243 \right]^{-2 / \nu}
245 for :math:`x \in \mathbb{R}`, and
247 .. math::
248 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
250 .. math::
251 \lambda_1 = \frac{\lambda_0}{2} \left[
252 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
253 \right].
255 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem will be solved
256 *multi-implicitly*.
258 Parameters
259 ----------
260 nvars : int, optional
261 Spatial resolution.
262 lambda0 : float, optional
263 Problem parameter : math:`\lambda_0`.
264 nu : float, optional
265 Problem parameter :math:`\nu`.
266 interval : tuple, optional
267 Defines the spatial domain.
268 comm : PETSc.COMM_WORLD, optional
269 Communicator for PETSc.
270 lsol_tol : float, optional
271 Tolerance for linear solver to terminate.
272 nlsol_tol : float, optional
273 Tolerance for nonlinear solver to terminate.
274 lsol_maxiter : int, optional
275 Maximum number of iterations for linear solver.
276 nlsol_maxiter : int, optional
277 Maximum number of iterations for nonlinear solver.
279 Attributes
280 ----------
281 dx : float
282 Mesh grid width.
283 xs, xe : int
284 Define the ranges.
285 A : PETSc matrix object
286 Discretization matrix.
287 localX : PETSc vector object
288 Local vector for solution.
289 ksp : PETSc solver object
290 Linear solver object.
291 snes : PETSc solver object
292 Nonlinear solver object.
293 F : PETSc vector object
294 Global vector.
295 J : PETSc matrix object
296 Jacobi matrix.
298 References
299 ----------
300 .. [1] Z. Feng. Traveling wave behavior for a generalized fisher equation. Chaos Solitons Fractals 38(2),
301 481–488 (2008).
302 .. [2] PETSc Web page. Satish Balay et al. https://petsc.org/ (2023).
303 .. [3] Parallel distributed computing using Python. Lisandro D. Dalcin, Rodrigo R. Paz, Pablo A. Kler,
304 Alejandro Cosimo. Advances in Water Resources (2011).
305 """
307 dtype_u = petsc_vec
308 dtype_f = petsc_vec_comp2
310 def __init__(
311 self,
312 nvars,
313 lambda0,
314 nu,
315 interval,
316 comm=PETSc.COMM_WORLD,
317 lsol_tol=1e-10,
318 nlsol_tol=1e-10,
319 lsol_maxiter=None,
320 nlsol_maxiter=None,
321 ):
322 """Initialization routine"""
323 # create DMDA object which will be used for all grid operations
324 da = PETSc.DMDA().create([nvars], dof=1, stencil_width=1, comm=comm)
326 # invoke super init, passing number of dofs, dtype_u and dtype_f
327 super().__init__(init=da)
328 self._makeAttributeAndRegister(
329 'nvars',
330 'lambda0',
331 'nu',
332 'interval',
333 'comm',
334 'lsol_tol',
335 'nlsol_tol',
336 'lsol_maxiter',
337 'nlsol_maxiter',
338 localVars=locals(),
339 readOnly=True,
340 )
342 # compute dx and get local ranges
343 self.dx = (self.interval[1] - self.interval[0]) / (self.nvars - 1)
344 (self.xs, self.xe) = self.init.getRanges()[0]
346 # compute discretization matrix A and identity
347 self.A = self.__get_A()
348 self.localX = self.init.createLocalVec()
350 # setup linear solver
351 self.ksp = PETSc.KSP()
352 self.ksp.create(comm=self.comm)
353 self.ksp.setType('cg')
354 pc = self.ksp.getPC()
355 pc.setType('ilu')
356 self.ksp.setInitialGuessNonzero(True)
357 self.ksp.setFromOptions()
358 self.ksp.setTolerances(rtol=self.lsol_tol, atol=self.lsol_tol, max_it=self.lsol_maxiter)
359 self.ksp_itercount = 0
360 self.ksp_ncalls = 0
362 # setup nonlinear solver
363 self.snes = PETSc.SNES()
364 self.snes.create(comm=self.comm)
365 if self.nlsol_maxiter <= 1:
366 self.snes.setType('ksponly')
367 self.snes.getKSP().setType('cg')
368 pc = self.snes.getKSP().getPC()
369 pc.setType('ilu')
370 # self.snes.setType('ngmres')
371 self.snes.setFromOptions()
372 self.snes.setTolerances(
373 rtol=self.nlsol_tol,
374 atol=self.nlsol_tol,
375 stol=self.nlsol_tol,
376 max_it=self.nlsol_maxiter,
377 )
378 self.snes_itercount = 0
379 self.snes_ncalls = 0
380 self.F = self.init.createGlobalVec()
381 self.J = self.init.createMatrix()
383 def __get_A(self):
384 r"""
385 Helper function to assemble ``PETSc`` matrix A.
387 Returns
388 -------
389 A : PETSc matrix object
390 Discretization matrix.
391 """
392 # create matrix and set basic options
393 A = self.init.createMatrix()
394 A.setType('aij') # sparse
395 A.setFromOptions()
396 A.setPreallocationNNZ((3, 3))
397 A.setUp()
399 # fill matrix
400 A.zeroEntries()
401 row = PETSc.Mat.Stencil()
402 col = PETSc.Mat.Stencil()
403 mx = self.init.getSizes()[0]
404 (xs, xe) = self.init.getRanges()[0]
405 for i in range(xs, xe):
406 row.i = i
407 row.field = 0
408 if i == 0 or i == mx - 1:
409 A.setValueStencil(row, row, 1.0)
410 else:
411 diag = -2.0 / self.dx**2
412 for index, value in [
413 (i - 1, 1.0 / self.dx**2),
414 (i, diag),
415 (i + 1, 1.0 / self.dx**2),
416 ]:
417 col.i = index
418 col.field = 0
419 A.setValueStencil(row, col, value)
420 A.assemble()
421 return A
423 def get_sys_mat(self, factor):
424 """
425 Helper function to assemble the system matrix of the linear problem.
427 Parameters
428 ----------
429 factor : float
430 Factor to define the system matrix.
432 Returns
433 -------
434 A : PETSc matrix object
435 Matrix for the system to solve.
436 """
438 # create matrix and set basic options
439 A = self.init.createMatrix()
440 A.setType('aij') # sparse
441 A.setFromOptions()
442 A.setPreallocationNNZ((3, 3))
443 A.setUp()
445 # fill matrix
446 A.zeroEntries()
447 row = PETSc.Mat.Stencil()
448 col = PETSc.Mat.Stencil()
449 mx = self.init.getSizes()[0]
450 (xs, xe) = self.init.getRanges()[0]
451 for i in range(xs, xe):
452 row.i = i
453 row.field = 0
454 if i == 0 or i == mx - 1:
455 A.setValueStencil(row, row, 1.0)
456 else:
457 diag = 1.0 + factor * 2.0 / self.dx**2
458 for index, value in [
459 (i - 1, -factor / self.dx**2),
460 (i, diag),
461 (i + 1, -factor / self.dx**2),
462 ]:
463 col.i = index
464 col.field = 0
465 A.setValueStencil(row, col, value)
466 A.assemble()
467 return A
469 def eval_f(self, u, t):
470 """
471 Routine to evaluate the right-hand side of the problem.
473 Parameters
474 ----------
475 u : dtype_u
476 Current values of the numerical solution.
477 t : float
478 Current time the numerical solution is computed.
480 Returns
481 -------
482 f : dtype_f
483 Right-hand side of the problem.
484 """
486 f = self.dtype_f(self.init)
487 self.A.mult(u, f.comp1)
488 fa1 = self.init.getVecArray(f.comp1)
489 fa1[0] = 0
490 fa1[-1] = 0
492 fa2 = self.init.getVecArray(f.comp2)
493 xa = self.init.getVecArray(u)
494 for i in range(self.xs, self.xe):
495 fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
496 fa2[0] = 0
497 fa2[-1] = 0
499 return f
501 def solve_system_1(self, rhs, factor, u0, t):
502 r"""
503 Linear solver for :math:`(I - factor \cdot A)\vec{u} = \vec{rhs}`.
505 Parameters
506 ----------
507 rhs : dtype_f
508 Right-hand side for the linear system.
509 factor : float
510 Abbrev. for the local stepsize (or any other factor required).
511 u0 : dtype_u
512 Initial guess for the iterative solver.
513 t : float
514 Current time (e.g. for time-dependent BCs).
516 Returns
517 -------
518 me : dtype_u
519 Solution as mesh.
520 """
522 me = self.dtype_u(u0)
524 self.ksp.setOperators(self.get_sys_mat(factor))
525 self.ksp.solve(rhs, me)
527 self.ksp_itercount += self.ksp.getIterationNumber()
528 self.ksp_ncalls += 1
530 return me
532 def solve_system_2(self, rhs, factor, u0, t):
533 r"""
534 Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
536 Parameters
537 ----------
538 rhs : dtype_f
539 Right-hand side for the linear system.
540 factor : float
541 Abbrev. for the local stepsize (or any other factor required).
542 u0 : dtype_u
543 Initial guess for the iterative solver.
544 t : float
545 Current time (e.g. for time-dependent BCs).
547 Returns
548 -------
549 me : dtype_u
550 Solution as mesh.
551 """
553 me = self.dtype_u(u0)
554 target = Fisher_reaction(self.init, self, factor)
556 # assign residual function and Jacobian
557 F = self.init.createGlobalVec()
558 self.snes.setFunction(target.formFunction, F)
559 J = self.init.createMatrix()
560 self.snes.setJacobian(target.formJacobian, J)
562 self.snes.solve(rhs, me)
564 self.snes_itercount += self.snes.getIterationNumber()
565 self.snes_ncalls += 1
567 return me
569 def u_exact(self, t):
570 r"""
571 Routine to compute the exact solution at time :math:`t`.
573 Parameters
574 ----------
575 t : float
576 Time of the exact solution.
578 Returns
579 -------
580 me : dtype_u
581 Exact solution.
582 """
584 lam1 = self.lambda0 / 2.0 * ((self.nu / 2.0 + 1) ** 0.5 + (self.nu / 2.0 + 1) ** (-0.5))
585 sig1 = lam1 - np.sqrt(lam1**2 - self.lambda0**2)
586 me = self.dtype_u(self.init)
587 xa = self.init.getVecArray(me)
588 for i in range(self.xs, self.xe):
589 xa[i] = (
590 1
591 + (2 ** (self.nu / 2.0) - 1)
592 * np.exp(-self.nu / 2.0 * sig1 * (self.interval[0] + (i + 1) * self.dx + 2 * lam1 * t))
593 ) ** (-2.0 / self.nu)
595 return me
598class petsc_fisher_fullyimplicit(petsc_fisher_multiimplicit):
599 r"""
600 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
601 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
602 problem [1]_ using periodic boundary conditions
604 .. math::
605 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
607 with exact solution
609 .. math::
610 u(x, 0) = \left[
611 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
612 \right]^{-2 / \nu}
614 for :math:`x \in \mathbb{R}`, and
616 .. math::
617 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
619 .. math::
620 \lambda_1 = \frac{\lambda_0}{2} \left[
621 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
622 \right].
624 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem is treated
625 *fully-implicitly*.
626 """
628 dtype_f = petsc_vec
630 def eval_f(self, u, t):
631 """
632 Routine to evaluate the right-hand side of the problem.
634 Parameters
635 ----------
636 u : dtype_u
637 Current values of the numerical solution.
638 t : float
639 Current time the numerical solution is computed.
641 Returns
642 -------
643 f : dtype_f
644 Right-hand side of the problem.
645 """
647 f = self.dtype_f(self.init)
648 self.A.mult(u, f)
650 fa2 = self.init.getVecArray(f)
651 xa = self.init.getVecArray(u)
652 for i in range(self.xs, self.xe):
653 fa2[i] += self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
654 fa2[0] = 0
655 fa2[-1] = 0
657 return f
659 def solve_system(self, rhs, factor, u0, t):
660 r"""
661 Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
663 Parameters
664 ----------
665 rhs : dtype_f
666 Right-hand side for the linear system.
667 factor : float
668 Abbrev. for the local stepsize (or any other factor required).
669 u0 : dtype_u
670 Initial guess for the iterative solver.
671 t : float
672 Current time (e.g. for time-dependent BCs).
674 Returns
675 -------
676 me : dtype_u
677 Solution as mesh.
678 """
680 me = self.dtype_u(u0)
681 target = Fisher_full(self.init, self, factor, self.dx)
683 # assign residual function and Jacobian
685 self.snes.setFunction(target.formFunction, self.F)
686 self.snes.setJacobian(target.formJacobian, self.J)
688 self.snes.solve(rhs, me)
690 self.snes_itercount += self.snes.getIterationNumber()
691 self.snes_ncalls += 1
693 return me
696class petsc_fisher_semiimplicit(petsc_fisher_multiimplicit):
697 r"""
698 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
699 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
700 problem [1]_ using periodic boundary conditions
702 .. math::
703 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
705 with exact solution
707 .. math::
708 u(x, 0) = \left[
709 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
710 \right]^{-2 / \nu}
712 for :math:`x \in \mathbb{R}`, and
714 .. math::
715 \sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
717 .. math::
718 \lambda_1 = \frac{\lambda_0}{2} \left[
719 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
720 \right].
722 This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem here will be
723 solved in a *semi-implicit* way.
724 """
726 dtype_f = petsc_vec_imex
728 def eval_f(self, u, t):
729 """
730 Routine to evaluate the right-hand side of the problem.
732 Parameters
733 ----------
734 u : dtype_u
735 Current values of the numerical solution.
736 t : float
737 Current time the numerical solution is computed.
739 Returns
740 -------
741 f : dtype_f
742 Right-hand side of the problem.
743 """
745 f = self.dtype_f(self.init)
746 self.A.mult(u, f.impl)
747 fa1 = self.init.getVecArray(f.impl)
748 fa1[0] = 0
749 fa1[-1] = 0
751 fa2 = self.init.getVecArray(f.expl)
752 xa = self.init.getVecArray(u)
753 for i in range(self.xs, self.xe):
754 fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
755 fa2[0] = 0
756 fa2[-1] = 0
758 return f
760 def solve_system(self, rhs, factor, u0, t):
761 r"""
762 Simple linear solver for :math:`(I-factor \cdot A)\vec{u} = \vec{rhs}`.
764 Parameters
765 ----------
766 rhs : dtype_f
767 Right-hand side for the linear system.
768 factor : float
769 Abbrev. for the local stepsize (or any other factor required).
770 u0 : dtype_u
771 Initial guess for the iterative solver.
772 t : float
773 Current time (e.g. for time-dependent BCs).
775 Returns
776 -------
777 me : dtype_u
778 Solution as mesh.
779 """
781 me = self.dtype_u(u0)
783 self.ksp.setOperators(self.get_sys_mat(factor))
784 self.ksp.solve(rhs, me)
786 self.ksp_itercount += self.ksp.getIterationNumber()
787 self.ksp_ncalls += 1
789 return me