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

1import numpy as np 

2from petsc4py import PETSc 

3 

4from pySDC.core.Problem import ptype 

5from pySDC.implementations.datatype_classes.petsc_vec import petsc_vec, petsc_vec_imex, petsc_vec_comp2 

6 

7 

8class GS_full(object): 

9 r""" 

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

11 

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. 

24 

25 Attributes 

26 ---------- 

27 localX : PETSc vector object 

28 Local vector for ``PETSc``. 

29 """ 

30 

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() 

40 

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. 

45 

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)`. 

54 

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 ) 

89 

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

91 """ 

92 Function to return the Jacobian matrix. 

93 

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. 

104 

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() 

116 

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) 

142 

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) 

172 

173 P.assemble() 

174 if J != P: 

175 J.assemble() # matrix-free operator 

176 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN 

177 

178 

179class GS_reaction(object): 

180 r""" 

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

182 

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`. 

191 

192 Attributes 

193 ---------- 

194 localX : PETSc vector object 

195 Local vector for ``PETSc``. 

196 """ 

197 

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() 

205 

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. 

210 

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)`. 

219 

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])) 

235 

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

237 """ 

238 Function to return the Jacobian matrix. 

239 

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. 

250 

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) 

278 

279 P.assemble() 

280 if J != P: 

281 J.assemble() # matrix-free operator 

282 return PETSc.Mat.Structure.SAME_NONZERO_PATTERN 

283 

284 

285class petsc_grayscott_multiimplicit(ptype): 

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 

291 

292 .. math:: 

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

294 

295 .. math:: 

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

297 

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. 

301 

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. 

324 

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. 

345 

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 """ 

354 

355 dtype_u = petsc_vec 

356 dtype_f = petsc_vec_comp2 

357 

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 ) 

380 

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(), 

395 readOnly=True, 

396 ) 

397 

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() 

402 

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() 

407 

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 

419 

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 

434 

435 def __get_A(self): 

436 r""" 

437 Helper function to assemble ``PETSc`` matrix A. 

438 

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() 

449 

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() 

495 

496 return A 

497 

498 def __get_Id(self): 

499 r""" 

500 Helper function to assemble ``PETSc`` identity matrix. 

501 

502 Returns 

503 ------- 

504 Id : PETSc matrix object 

505 Identity matrix. 

506 """ 

507 

508 Id = self.init.createMatrix() 

509 Id.setType('aij') # sparse 

510 Id.setFromOptions() 

511 Id.setPreallocationNNZ((1, 1)) 

512 Id.setUp() 

513 

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) 

524 

525 Id.assemble() 

526 

527 return Id 

528 

529 def eval_f(self, u, t): 

530 """ 

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

532 

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. 

539 

540 Returns 

541 ------- 

542 f : dtype_f 

543 Right-hand side of the problem. 

544 """ 

545 

546 f = self.dtype_f(self.init) 

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

548 

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] 

555 

556 return f 

557 

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}`. 

561 

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). 

572 

573 Returns 

574 ------- 

575 me : dtype_u 

576 Solution as mesh. 

577 """ 

578 

579 me = self.dtype_u(u0) 

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

581 self.ksp.solve(rhs, me) 

582 

583 self.ksp_ncalls += 1 

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

585 

586 return me 

587 

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}`. 

591 

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). 

602 

603 Returns 

604 ------- 

605 me : dtype_u 

606 Solution as mesh. 

607 """ 

608 

609 me = self.dtype_u(u0) 

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

611 

612 F = self.init.createGlobalVec() 

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

614 J = self.init.createMatrix() 

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

616 

617 self.snes.solve(rhs, me) 

618 

619 self.snes_ncalls += 1 

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

621 

622 return me 

623 

624 def u_exact(self, t): 

625 r""" 

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

627 

628 Parameters 

629 ---------- 

630 t : float 

631 Time of the exact solution. 

632 

633 Returns 

634 ------- 

635 me : dtype_u 

636 Exact solution. 

637 """ 

638 

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

640 

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 ) 

651 

652 return me 

653 

654 

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 

661 

662 .. math:: 

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

664 

665 .. math:: 

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

667 

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 """ 

672 

673 dtype_f = petsc_vec 

674 

675 def eval_f(self, u, t): 

676 """ 

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

678 

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. 

685 

686 Returns 

687 ------- 

688 f : dtype_f 

689 Right-hand side of the problem. 

690 """ 

691 

692 f = self.dtype_f(self.init) 

693 self.AMat.mult(u, f) 

694 

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] 

701 

702 return f 

703 

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

705 r""" 

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

707 

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). 

718 

719 Returns 

720 ------- 

721 me : dtype_u 

722 Solution as mesh. 

723 """ 

724 

725 me = self.dtype_u(u0) 

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

727 

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) 

733 

734 self.snes.solve(rhs, me) 

735 

736 self.snes_ncalls += 1 

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

738 

739 return me 

740 

741 

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 

748 

749 .. math:: 

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

751 

752 .. math:: 

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

754 

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 """ 

758 

759 dtype_f = petsc_vec_imex 

760 

761 def eval_f(self, u, t): 

762 """ 

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

764 

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. 

771 

772 Returns 

773 ------- 

774 f : dtype_f 

775 Right-hand side of the problem. 

776 """ 

777 

778 f = self.dtype_f(self.init) 

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

780 

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] 

787 

788 return f 

789 

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

791 r""" 

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

793 

794 Parameters 

795 ---------- 

796 rhs : dtype_f