# Coverage for pySDC/implementations/problem_classes/GeneralizedFisher_1D_PETSc.py: 99%

## 243 statements

, created at 2024-09-20 17:10 +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(),

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