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

## 305 statements

, created at 2024-09-09 14:59 +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 GS_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.

22 dy : float

23 Grid spacing in y direction.

25 Attributes

26 ----------

27 localX : PETSc vector object

28 Local vector for PETSc.

29 """

31 def __init__(self, da, prob, factor, dx, dy):

32 """Initialization routine"""

33 assert da.getDim() == 2

34 self.da = da

35 self.prob = prob

36 self.factor = factor

37 self.dx = dx

38 self.dy = dy

39 self.localX = da.createLocalVec()

41 def formFunction(self, snes, X, F):

42 r"""

43 Function to evaluate the residual for the Newton solver. This function should be equal to the RHS

44 in the solution.

46 Parameters

47 ----------

48 snes : PETSc solver object

49 Nonlinear solver object.

50 X : PETSc vector object

51 Input vector.

52 F : PETSc vector object

53 Output vector :math:F(X).

55 Returns

56 -------

57 None

58 Overwrites F.

59 """

60 self.da.globalToLocal(X, self.localX)

61 x = self.da.getVecArray(self.localX)

62 f = self.da.getVecArray(F)

63 (xs, xe), (ys, ye) = self.da.getRanges()

64 for j in range(ys, ye):

65 for i in range(xs, xe):

66 u = x[i, j] # center

67 u_e = x[i + 1, j] # east

68 u_w = x[i - 1, j] # west

69 u_s = x[i, j + 1] # south

70 u_n = x[i, j - 1] # north

71 u_xx = u_e - 2 * u + u_w

72 u_yy = u_n - 2 * u + u_s

73 f[i, j, 0] = x[i, j, 0] - (

74 self.factor

75 * (

76 self.prob.Du * (u_xx[0] / self.dx**2 + u_yy[0] / self.dy**2)

77 - x[i, j, 0] * x[i, j, 1] ** 2

78 + self.prob.A * (1 - x[i, j, 0])

79 )

80 )

81 f[i, j, 1] = x[i, j, 1] - (

82 self.factor

83 * (

84 self.prob.Dv * (u_xx[1] / self.dx**2 + u_yy[1] / self.dy**2)

85 + x[i, j, 0] * x[i, j, 1] ** 2

86 - self.prob.B * x[i, j, 1]

87 )

88 )

90 def formJacobian(self, snes, X, J, P):

91 """

92 Function to return the Jacobian matrix.

94 Parameters

95 ----------

96 snes : PETSc solver object

97 Nonlinear solver object.

98 X : PETSc vector object

99 Input vector.

100 J : PETSc matrix object

101 Current Jacobian matrix.

102 P : PETSc matrix object

103 New Jacobian matrix.

105 Returns

106 -------

107 PETSc.Mat.Structure.SAME_NONZERO_PATTERN

108 Matrix status.

109 """

110 self.da.globalToLocal(X, self.localX)

111 x = self.da.getVecArray(self.localX)

112 P.zeroEntries()

113 row = PETSc.Mat.Stencil()

114 col = PETSc.Mat.Stencil()

115 (xs, xe), (ys, ye) = self.da.getRanges()

117 for j in range(ys, ye):

118 for i in range(xs, xe):

119 # diagnoal 2-by-2 block (for u and v)

120 row.index = (i, j)

121 col.index = (i, j)

122 row.field = 0

123 col.field = 0

124 val = 1.0 - self.factor * (

125 self.prob.Du * (-2.0 / self.dx**2 - 2.0 / self.dy**2) - x[i, j, 1] ** 2 - self.prob.A

126 )

127 P.setValueStencil(row, col, val)

128 row.field = 0

129 col.field = 1

130 val = self.factor * 2.0 * x[i, j, 0] * x[i, j, 1]

131 P.setValueStencil(row, col, val)

132 row.field = 1

133 col.field = 1

134 val = 1.0 - self.factor * (

135 self.prob.Dv * (-2.0 / self.dx**2 - 2.0 / self.dy**2) + 2.0 * x[i, j, 0] * x[i, j, 1] - self.prob.B

136 )

137 P.setValueStencil(row, col, val)

138 row.field = 1

139 col.field = 0

140 val = -self.factor * x[i, j, 1] ** 2

141 P.setValueStencil(row, col, val)

143 # coupling through finite difference part

144 col.index = (i, j - 1)

145 col.field = 0

146 row.field = 0

147 P.setValueStencil(row, col, -self.factor * self.prob.Du / self.dx**2)

148 col.field = 1

149 row.field = 1

150 P.setValueStencil(row, col, -self.factor * self.prob.Dv / self.dy**2)

151 col.index = (i, j + 1)

152 col.field = 0

153 row.field = 0

154 P.setValueStencil(row, col, -self.factor * self.prob.Du / self.dx**2)

155 col.field = 1

156 row.field = 1

157 P.setValueStencil(row, col, -self.factor * self.prob.Dv / self.dy**2)

158 col.index = (i - 1, j)

159 col.field = 0

160 row.field = 0

161 P.setValueStencil(row, col, -self.factor * self.prob.Du / self.dx**2)

162 col.field = 1

163 row.field = 1

164 P.setValueStencil(row, col, -self.factor * self.prob.Dv / self.dy**2)

165 col.index = (i + 1, j)

166 col.field = 0

167 row.field = 0

168 P.setValueStencil(row, col, -self.factor * self.prob.Du / self.dx**2)

169 col.field = 1

170 row.field = 1

171 P.setValueStencil(row, col, -self.factor * self.prob.Dv / self.dy**2)

173 P.assemble()

174 if J != P:

175 J.assemble() # matrix-free operator

176 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

179class GS_reaction(object):

180 r"""

181 Helper class to generate residual and Jacobian matrix for PETSc's nonlinear solver SNES.

183 Parameters

184 ----------

185 da : DMDA object

186 Object of PETSc.

187 prob : problem instance

188 Contains problem information for PETSc.

189 factor : float

190 Temporal factor :math:\Delta t Q_\Delta.

192 Attributes

193 ----------

194 localX : PETSc vector object

195 Local vector for PETSc.

196 """

198 def __init__(self, da, prob, factor):

199 """Initialization routine"""

200 assert da.getDim() == 2

201 self.da = da

202 self.prob = prob

203 self.factor = factor

204 self.localX = da.createLocalVec()

206 def formFunction(self, snes, X, F):

207 r"""

208 Function to evaluate the residual for the Newton solver. This function should be equal to the RHS

209 in the solution.

211 Parameters

212 ----------

213 snes : PETSc solver object

214 Nonlinear solver object.

215 X : PETSc vector object

216 Input vector.

217 F : PETSc vector object

218 Output vector :math:F(X).

220 Returns

221 -------

222 None

223 Overwrites F.

224 """

225 self.da.globalToLocal(X, self.localX)

226 x = self.da.getVecArray(self.localX)

227 f = self.da.getVecArray(F)

228 (xs, xe), (ys, ye) = self.da.getRanges()

229 for j in range(ys, ye):

230 for i in range(xs, xe):

231 f[i, j, 0] = x[i, j, 0] - (

232 self.factor * (-x[i, j, 0] * x[i, j, 1] ** 2 + self.prob.A * (1 - x[i, j, 0]))

233 )

234 f[i, j, 1] = x[i, j, 1] - (self.factor * (x[i, j, 0] * x[i, j, 1] ** 2 - self.prob.B * x[i, j, 1]))

236 def formJacobian(self, snes, X, J, P):

237 """

238 Function to return the Jacobian matrix.

240 Parameters

241 ----------

242 snes : PETSc solver object

243 Nonlinear solver object.

244 X : PETSc vector object

245 Input vector.

246 J : PETSc matrix object

247 Current Jacobian matrix.

248 P : PETSc matrix object

249 New Jacobian matrix.

251 Returns

252 -------

253 PETSc.Mat.Structure.SAME_NONZERO_PATTERN

254 Matrix status.

255 """

256 self.da.globalToLocal(X, self.localX)

257 x = self.da.getVecArray(self.localX)

258 P.zeroEntries()

259 row = PETSc.Mat.Stencil()

260 col = PETSc.Mat.Stencil()

261 (xs, xe), (ys, ye) = self.da.getRanges()

262 for j in range(ys, ye):

263 for i in range(xs, xe):

264 row.index = (i, j)

265 col.index = (i, j)

266 row.field = 0

267 col.field = 0

268 P.setValueStencil(row, col, 1.0 - self.factor * (-x[i, j, 1] ** 2 - self.prob.A))

269 row.field = 0

270 col.field = 1

271 P.setValueStencil(row, col, self.factor * 2.0 * x[i, j, 0] * x[i, j, 1])

272 row.field = 1

273 col.field = 1

274 P.setValueStencil(row, col, 1.0 - self.factor * (2.0 * x[i, j, 0] * x[i, j, 1] - self.prob.B))

275 row.field = 1

276 col.field = 0

277 P.setValueStencil(row, col, -self.factor * x[i, j, 1] ** 2)

279 P.assemble()

280 if J != P:

281 J.assemble() # matrix-free operator

282 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

285class petsc_grayscott_multiimplicit(Problem):

286 r"""

287 The Gray-Scott system [1]_ describes a reaction-diffusion process of two substances :math:u and :math:v,

288 where they diffuse over time. During the reaction :math:u is used up with overall decay rate :math:B,

289 whereas :math:v is produced with feed rate :math:A. :math:D_u,\, D_v are the diffusion rates for

290 :math:u,\, v. This process is described by the two-dimensional model using periodic boundary conditions

292 .. math::

293 \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),

295 .. math::

296 \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u

298 for :math:x \in \Omega:=[0, 100]. The spatial solve of the problem is realized by PETSc [2]_, [3]_. For time-stepping,

299 the diffusion part is solved by one of PETSc's linear solver, whereas the reaction part will be solved by a nonlinear

300 solver.

302 Parameters

303 ----------

304 nvars : tuple of int, optional

305 Spatial resolution, i.e., number of degrees of freedom in space, e.g. nvars=(256, 256).

306 Du : float, optional

307 Diffusion rate for :math:u.

308 Dv: float, optional

309 Diffusion rate for :math:v.

310 A : float, optional

311 Feed rate for :math:v.

312 B : float, optional

313 Overall decay rate for :math:u.

314 comm : PETSc.COMM_WORLD, optional

315 Communicator for PETSc.

316 lsol_tol : float, optional

317 Tolerance for linear solver to terminate.

318 nlsol_tol : float, optional

319 Tolerance for nonlinear solver to terminate.

320 lsol_maxiter : int, optional

321 Maximum number of iterations for linear solver.

322 nlsol_maxiter : int, optional

323 Maximum number of iterations for nonlinear solver.

325 Attributes

326 ----------

327 dx : float

328 Mesh grid width in x direction.

329 dy : float

330 Mesh grid width in y direction.

331 AMat : PETSc matrix object

332 Discretization matrix.

333 Id : PETSc matrix object

334 Identity matrix.

335 localX : PETSc vector object

336 Local vector for solution.

337 ksp : PETSc solver object

338 Linear solver object.

339 snes : PETSc solver object

340 Nonlinear solver object.

341 snes_itercount : int

342 Number of iterations done by nonlinear solver.

343 snes_ncalls : int

344 Number of calls of PETSc's nonlinear solver.

346 References

347 ----------

348 .. [1] Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms

349 of multistability. P. Gray, S. K. Scott. Chem. Eng. Sci. 38, 1 (1983).

350 .. [2] PETSc Web page. Satish Balay et al. https://petsc.org/ (2023).

351 .. [3] Parallel distributed computing using Python. Lisandro D. Dalcin, Rodrigo R. Paz, Pablo A. Kler,

352 Alejandro Cosimo. Advances in Water Resources (2011).

353 """

355 dtype_u = petsc_vec

356 dtype_f = petsc_vec_comp2

358 def __init__(

359 self,

360 nvars,

361 Du,

362 Dv,

363 A,

364 B,

365 comm=PETSc.COMM_WORLD,

366 lsol_tol=1e-10,

367 nlsol_tol=1e-10,

368 lsol_maxiter=None,

369 nlsol_maxiter=None,

370 ):

371 """Initialization routine"""

372 # create DMDA object which will be used for all grid operations (boundary_type=3 -> periodic BC)

373 da = PETSc.DMDA().create(

374 [nvars[0], nvars[1]],

375 dof=2,

376 boundary_type=3,

377 stencil_width=1,

378 comm=comm,

379 )

381 # invoke super init, passing number of dofs, dtype_u and dtype_f

382 super().__init__(init=da)

383 self._makeAttributeAndRegister(

384 'nvars',

385 'Du',

386 'Dv',

387 'A',

388 'B',

389 'comm',

390 'lsol_tol',

391 'lsol_maxiter',

392 'nlsol_tol',

393 'nlsol_maxiter',

394 localVars=locals(),

396 )

398 # compute dx, dy and get local ranges

399 self.dx = 100.0 / (self.nvars[0])

400 self.dy = 100.0 / (self.nvars[1])

401 (self.xs, self.xe), (self.ys, self.ye) = self.init.getRanges()

403 # compute discretization matrix A and identity

404 self.AMat = self.__get_A()

405 self.Id = self.__get_Id()

406 self.localX = self.init.createLocalVec()

408 # setup linear solver

409 self.ksp = PETSc.KSP()

410 self.ksp.create(comm=self.comm)

411 self.ksp.setType('cg')

412 pc = self.ksp.getPC()

413 pc.setType('none')

414 self.ksp.setInitialGuessNonzero(True)

415 self.ksp.setFromOptions()

416 self.ksp.setTolerances(rtol=self.lsol_tol, atol=self.lsol_tol, max_it=self.lsol_maxiter)

417 self.ksp_itercount = 0

418 self.ksp_ncalls = 0

420 # setup nonlinear solver

421 self.snes = PETSc.SNES()

422 self.snes.create(comm=self.comm)

423 # self.snes.getKSP().setType('cg')

424 # self.snes.setType('ngmres')

425 self.snes.setFromOptions()

426 self.snes.setTolerances(

427 rtol=self.nlsol_tol,

428 atol=self.nlsol_tol,

429 stol=self.nlsol_tol,

430 max_it=self.nlsol_maxiter,

431 )

432 self.snes_itercount = 0

433 self.snes_ncalls = 0

435 def __get_A(self):

436 r"""

437 Helper function to assemble PETSc matrix A.

439 Returns

440 -------

441 A : PETSc matrix object

442 Discretization matrix.

443 """

444 A = self.init.createMatrix()

445 A.setType('aij') # sparse

446 A.setFromOptions()

447 A.setPreallocationNNZ((5, 5))

448 A.setUp()

450 A.zeroEntries()

451 row = PETSc.Mat.Stencil()

452 col = PETSc.Mat.Stencil()

453 mx, my = self.init.getSizes()

454 (xs, xe), (ys, ye) = self.init.getRanges()

455 for j in range(ys, ye):

456 for i in range(xs, xe):

457 row.index = (i, j)

458 row.field = 0

459 A.setValueStencil(row, row, self.Du * (-2.0 / self.dx**2 - 2.0 / self.dy**2))

460 row.field = 1

461 A.setValueStencil(row, row, self.Dv * (-2.0 / self.dx**2 - 2.0 / self.dy**2))

462 # if j > 0:

463 col.index = (i, j - 1)

464 col.field = 0

465 row.field = 0

466 A.setValueStencil(row, col, self.Du / self.dy**2)

467 col.field = 1

468 row.field = 1

469 A.setValueStencil(row, col, self.Dv / self.dy**2)

470 # if j < my - 1:

471 col.index = (i, j + 1)

472 col.field = 0

473 row.field = 0

474 A.setValueStencil(row, col, self.Du / self.dy**2)

475 col.field = 1

476 row.field = 1

477 A.setValueStencil(row, col, self.Dv / self.dy**2)

478 # if i > 0:

479 col.index = (i - 1, j)

480 col.field = 0

481 row.field = 0

482 A.setValueStencil(row, col, self.Du / self.dx**2)

483 col.field = 1

484 row.field = 1

485 A.setValueStencil(row, col, self.Dv / self.dx**2)

486 # if i < mx - 1:

487 col.index = (i + 1, j)

488 col.field = 0

489 row.field = 0

490 A.setValueStencil(row, col, self.Du / self.dx**2)

491 col.field = 1

492 row.field = 1

493 A.setValueStencil(row, col, self.Dv / self.dx**2)

494 A.assemble()

496 return A

498 def __get_Id(self):

499 r"""

500 Helper function to assemble PETSc identity matrix.

502 Returns

503 -------

504 Id : PETSc matrix object

505 Identity matrix.

506 """

508 Id = self.init.createMatrix()

509 Id.setType('aij') # sparse

510 Id.setFromOptions()

511 Id.setPreallocationNNZ((1, 1))

512 Id.setUp()

514 Id.zeroEntries()

515 row = PETSc.Mat.Stencil()

516 mx, my = self.init.getSizes()

517 (xs, xe), (ys, ye) = self.init.getRanges()

518 for j in range(ys, ye):

519 for i in range(xs, xe):

520 for indx in [0, 1]:

521 row.index = (i, j)

522 row.field = indx

523 Id.setValueStencil(row, row, 1.0)

525 Id.assemble()

527 return Id

529 def eval_f(self, u, t):

530 """

531 Routine to evaluate the right-hand side of the problem.

533 Parameters

534 ----------

535 u : dtype_u

536 Current values of the numerical solution.

537 t : float

538 Current time the numerical solution is computed.

540 Returns

541 -------

542 f : dtype_f

543 Right-hand side of the problem.

544 """

546 f = self.dtype_f(self.init)

547 self.AMat.mult(u, f.comp1)

549 fa = self.init.getVecArray(f.comp2)

550 xa = self.init.getVecArray(u)

551 for i in range(self.xs, self.xe):

552 for j in range(self.ys, self.ye):

553 fa[i, j, 0] = -xa[i, j, 0] * xa[i, j, 1] ** 2 + self.A * (1 - xa[i, j, 0])

554 fa[i, j, 1] = xa[i, j, 0] * xa[i, j, 1] ** 2 - self.B * xa[i, j, 1]

556 return f

558 def solve_system_1(self, rhs, factor, u0, t):

559 r"""

560 Linear solver for :math:(I - factor \cdot A)\vec{u} = \vec{rhs}.

562 Parameters

563 ----------

564 rhs : dtype_f

565 Right-hand side for the linear system.

566 factor : float

567 Abbrev. for the local stepsize (or any other factor required).

568 u0 : dtype_u

569 Initial guess for the iterative solver.

570 t : float

571 Current time (e.g. for time-dependent BCs).

573 Returns

574 -------

575 me : dtype_u

576 Solution as mesh.

577 """

579 me = self.dtype_u(u0)

580 self.ksp.setOperators(self.Id - factor * self.AMat)

581 self.ksp.solve(rhs, me)

583 self.ksp_ncalls += 1

584 self.ksp_itercount += self.ksp.getIterationNumber()

586 return me

588 def solve_system_2(self, rhs, factor, u0, t):

589 r"""

590 Nonlinear solver for :math:(I - factor \cdot F)(\vec{u}) = \vec{rhs}.

592 Parameters

593 ----------

594 rhs : dtype_f

595 Right-hand side for the linear system.

596 factor : float

597 Abbrev. for the local stepsize (or any other factor required).

598 u0 : dtype_u

599 Initial guess for the iterative solver.

600 t : float

601 Current time (e.g. for time-dependent BCs).

603 Returns

604 -------

605 me : dtype_u

606 Solution as mesh.

607 """

609 me = self.dtype_u(u0)

610 target = GS_reaction(self.init, self, factor)

612 F = self.init.createGlobalVec()

613 self.snes.setFunction(target.formFunction, F)

614 J = self.init.createMatrix()

615 self.snes.setJacobian(target.formJacobian, J)

617 self.snes.solve(rhs, me)

619 self.snes_ncalls += 1

620 self.snes_itercount += self.snes.getIterationNumber()

622 return me

624 def u_exact(self, t):

625 r"""

626 Routine to compute the exact solution at time :math:t.

628 Parameters

629 ----------

630 t : float

631 Time of the exact solution.

633 Returns

634 -------

635 me : dtype_u

636 Exact solution.

637 """

639 assert t == 0, 'ERROR: u_exact is only valid for the initial solution'

641 me = self.dtype_u(self.init)

642 xa = self.init.getVecArray(me)

643 for i in range(self.xs, self.xe):

644 for j in range(self.ys, self.ye):

645 xa[i, j, 0] = 1.0 - 0.5 * np.power(

646 np.sin(np.pi * i * self.dx / 100) * np.sin(np.pi * j * self.dy / 100), 100

647 )

648 xa[i, j, 1] = 0.25 * np.power(

649 np.sin(np.pi * i * self.dx / 100) * np.sin(np.pi * j * self.dy / 100), 100

650 )

652 return me

655class petsc_grayscott_fullyimplicit(petsc_grayscott_multiimplicit):

656 r"""

657 The Gray-Scott system [1]_ describes a reaction-diffusion process of two substances :math:u and :math:v,

658 where they diffuse over time. During the reaction :math:u is used up with overall decay rate :math:B,

659 whereas :math:v is produced with feed rate :math:A. :math:D_u,\, D_v are the diffusion rates for

660 :math:u,\, v. This process is described by the two-dimensional model using periodic boundary conditions

662 .. math::

663 \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),

665 .. math::

666 \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u

668 for :math:x \in \Omega:=[0, 100]. The spatial solve of the problem is realized by PETSc [2]_, [3]_. For time-stepping, the

669 problem is handled in a *fully-implicit* way, i.e., the nonlinear system containing the full right-hand side will be

670 solved by PETSc's nonlinear solver.

671 """

673 dtype_f = petsc_vec

675 def eval_f(self, u, t):

676 """

677 Routine to evaluate the right-hand side of the problem.

679 Parameters

680 ----------

681 u : dtype_u

682 Current values of the numerical solution.

683 t : float

684 Current time the numerical solution is computed.

686 Returns

687 -------

688 f : dtype_f

689 Right-hand side of the problem.

690 """

692 f = self.dtype_f(self.init)

693 self.AMat.mult(u, f)

695 fa = self.init.getVecArray(f)

696 xa = self.init.getVecArray(u)

697 for i in range(self.xs, self.xe):

698 for j in range(self.ys, self.ye):

699 fa[i, j, 0] += -xa[i, j, 0] * xa[i, j, 1] ** 2 + self.A * (1 - xa[i, j, 0])

700 fa[i, j, 1] += xa[i, j, 0] * xa[i, j, 1] ** 2 - self.B * xa[i, j, 1]

702 return f

704 def solve_system(self, rhs, factor, u0, t):

705 r"""

706 Nonlinear solver for :math:(I - factor \cdot F)(\vec{u}) = \vec{rhs}.

708 Parameters

709 ----------

710 rhs : dtype_f

711 Right-hand side for the linear system.

712 factor : float

713 Abbrev. for the local stepsize (or any other factor required).

714 u0 : dtype_u

715 Initial guess for the iterative solver.

716 t : float

717 Current time (e.g. for time-dependent BCs).

719 Returns

720 -------

721 me : dtype_u

722 Solution as mesh.

723 """

725 me = self.dtype_u(u0)

726 target = GS_full(self.init, self, factor, self.dx, self.dy)

728 # assign residual function and Jacobian

729 F = self.init.createGlobalVec()

730 self.snes.setFunction(target.formFunction, F)

731 J = self.init.createMatrix()

732 self.snes.setJacobian(target.formJacobian, J)

734 self.snes.solve(rhs, me)

736 self.snes_ncalls += 1

737 self.snes_itercount += self.snes.getIterationNumber()

739 return me

742class petsc_grayscott_semiimplicit(petsc_grayscott_multiimplicit):

743 r"""

744 The Gray-Scott system [1]_ describes a reaction-diffusion process of two substances :math:u and :math:v,

745 where they diffuse over time. During the reaction :math:u is used up with overall decay rate :math:B,

746 whereas :math:v is produced with feed rate :math:A. :math:D_u,\, D_v are the diffusion rates for

747 :math:u,\, v. This process is described by the two-dimensional model using periodic boundary conditions

749 .. math::

750 \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),

752 .. math::

753 \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u

755 for :math:x \in \Omega:=[0, 100]. The spatial solve of the problem is realized by PETSc [2]_, [3]_. For time-stepping,

756 the problem is treated *semi-implicitly*, i.e., the system with diffusion part is solved by PETSc's linear solver.

757 """

759 dtype_f = petsc_vec_imex

761 def eval_f(self, u, t):

762 """

763 Routine to evaluate the right-hand side of the problem.

765 Parameters

766 ----------

767 u : dtype_u

768 Current values of the numerical solution.

769 t : float

770 Current time the numerical solution is computed.

772 Returns

773 -------

774 f : dtype_f

775 Right-hand side of the problem.

776 """

778 f = self.dtype_f(self.init)

779 self.AMat.mult(u, f.impl)

781 fa = self.init.getVecArray(f.expl)

782 xa = self.init.getVecArray(u)

783 for i in range(self.xs, self.xe):

784 for j in range(self.ys, self.ye):

785 fa[i, j, 0] = -xa[i, j, 0] * xa[i, j, 1] ** 2 + self.A * (1 - xa[i, j, 0])

786 fa[i, j, 1] = xa[i, j, 0] * xa[i, j, 1] ** 2 - self.B * xa[i, j, 1]

788 return f

790 def solve_system(self, rhs, factor, u0, t):

791 r"""

792 Linear solver for :math:(I - factor \cdot A)\vec{u} = \vec{rhs}.

794 Parameters

795 ----------

796 rhs : dtype_f

797 Right-hand side for the linear system.

798 factor : float

799 Abbrev. for the local stepsize (or any other factor required).

800 u0 : dtype_u

801 Initial guess for the iterative solver.

802 t : float

803 Current time (e.g. for time-dependent BCs).

805 Returns

806 -------

807 me : dtype_u

808 Solution as mesh.

809 """

811 me = self.dtype_u(u0)

812 self.ksp.setOperators(self.Id - factor * self.AMat)

813 self.ksp.solve(rhs, me)

815 self.ksp_ncalls += 1

816 self.ksp_itercount += self.ksp.getIterationNumber()

818 return me