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 Problem 

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

286 r""" 

287 2D periodic Gray-Scott system with PETSc finite differences, multi-implicit: diffusion (CG), reaction (SNES). 

288 

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

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

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

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

293 

294 .. math:: 

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

296 

297 .. math:: 

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

299 

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

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

302 solver. 

303 

304 Parameters 

305 ---------- 

306 nvars : tuple of int, optional 

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

308 Du : float, optional 

309 Diffusion rate for :math:`u`. 

310 Dv: float, optional 

311 Diffusion rate for :math:`v`. 

312 A : float, optional 

313 Feed rate for :math:`v`. 

314 B : float, optional 

315 Overall decay rate for :math:`u`. 

316 comm : PETSc.COMM_WORLD, optional 

317 Communicator for ``PETSc``. 

318 lsol_tol : float, optional 

319 Tolerance for linear solver to terminate. 

320 nlsol_tol : float, optional 

321 Tolerance for nonlinear solver to terminate. 

322 lsol_maxiter : int, optional 

323 Maximum number of iterations for linear solver. 

324 nlsol_maxiter : int, optional 

325 Maximum number of iterations for nonlinear solver. 

326 

327 Attributes 

328 ---------- 

329 dx : float 

330 Mesh grid width in x direction. 

331 dy : float 

332 Mesh grid width in y direction. 

333 AMat : PETSc matrix object 

334 Discretization matrix. 

335 Id : PETSc matrix object 

336 Identity matrix. 

337 localX : PETSc vector object 

338 Local vector for solution. 

339 ksp : PETSc solver object 

340 Linear solver object. 

341 snes : PETSc solver object 

342 Nonlinear solver object. 

343 snes_itercount : int 

344 Number of iterations done by nonlinear solver. 

345 snes_ncalls : int 

346 Number of calls of ``PETSc``'s nonlinear solver. 

347 

348 References 

349 ---------- 

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

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

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

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

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

355 """ 

356 

357 dtype_u = petsc_vec 

358 dtype_f = petsc_vec_comp2 

359 

360 def __init__( 

361 self, 

362 nvars, 

363 Du, 

364 Dv, 

365 A, 

366 B, 

367 comm=PETSc.COMM_WORLD, 

368 lsol_tol=1e-10, 

369 nlsol_tol=1e-10, 

370 lsol_maxiter=None, 

371 nlsol_maxiter=None, 

372 ): 

373 """Initialization routine""" 

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

375 da = PETSc.DMDA().create( 

376 [nvars[0], nvars[1]], 

377 dof=2, 

378 boundary_type=3, 

379 stencil_width=1, 

380 comm=comm, 

381 ) 

382 

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

384 super().__init__(init=da) 

385 self._makeAttributeAndRegister( 

386 'nvars', 

387 'Du', 

388 'Dv', 

389 'A', 

390 'B', 

391 'comm', 

392 'lsol_tol', 

393 'lsol_maxiter', 

394 'nlsol_tol', 

395 'nlsol_maxiter', 

396 localVars=locals(), 

397 readOnly=True, 

398 ) 

399 

400 # compute dx, dy and get local ranges 

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

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

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

404 

405 # compute discretization matrix A and identity 

406 self.AMat = self.__get_A() 

407 self.Id = self.__get_Id() 

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

409 

410 # setup linear solver 

411 self.ksp = PETSc.KSP() 

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

413 self.ksp.setType('cg') 

414 pc = self.ksp.getPC() 

415 pc.setType('none') 

416 self.ksp.setInitialGuessNonzero(True) 

417 self.ksp.setFromOptions() 

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

419 self.ksp_itercount = 0 

420 self.ksp_ncalls = 0 

421 

422 # setup nonlinear solver 

423 self.snes = PETSc.SNES() 

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

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

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

427 self.snes.setFromOptions() 

428 self.snes.setTolerances( 

429 rtol=self.nlsol_tol, 

430 atol=self.nlsol_tol, 

431 stol=self.nlsol_tol, 

432 max_it=self.nlsol_maxiter, 

433 ) 

434 self.snes_itercount = 0 

435 self.snes_ncalls = 0 

436 

437 def __get_A(self): 

438 r""" 

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

440 

441 Returns 

442 ------- 

443 A : PETSc matrix object 

444 Discretization matrix. 

445 """ 

446 A = self.init.createMatrix() 

447 A.setType('aij') # sparse 

448 A.setFromOptions() 

449 A.setPreallocationNNZ((5, 5)) 

450 A.setUp() 

451 

452 A.zeroEntries() 

453 row = PETSc.Mat.Stencil() 

454 col = PETSc.Mat.Stencil() 

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

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

457 for j in range(ys, ye): 

458 for i in range(xs, xe): 

459 row.index = (i, j) 

460 row.field = 0 

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

462 row.field = 1 

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

464 # if j > 0: 

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

466 col.field = 0 

467 row.field = 0 

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

469 col.field = 1 

470 row.field = 1 

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

472 # if j < my - 1: 

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

474 col.field = 0 

475 row.field = 0 

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

477 col.field = 1 

478 row.field = 1 

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

480 # if i > 0: 

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

482 col.field = 0 

483 row.field = 0 

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

485 col.field = 1 

486 row.field = 1 

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

488 # if i < mx - 1: 

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

490 col.field = 0 

491 row.field = 0 

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

493 col.field = 1 

494 row.field = 1 

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

496 A.assemble() 

497 

498 return A 

499 

500 def __get_Id(self): 

501 r""" 

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

503 

504 Returns 

505 ------- 

506 Id : PETSc matrix object 

507 Identity matrix. 

508 """ 

509 

510 Id = self.init.createMatrix() 

511 Id.setType('aij') # sparse 

512 Id.setFromOptions() 

513 Id.setPreallocationNNZ((1, 1)) 

514 Id.setUp() 

515 

516 Id.zeroEntries() 

517 row = PETSc.Mat.Stencil() 

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

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

520 for j in range(ys, ye): 

521 for i in range(xs, xe): 

522 for indx in [0, 1]: 

523 row.index = (i, j) 

524 row.field = indx 

525 Id.setValueStencil(row, row, 1.0) 

526 

527 Id.assemble() 

528 

529 return Id 

530 

531 def eval_f(self, u, t): 

532 """ 

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

534 

535 Parameters 

536 ---------- 

537 u : dtype_u 

538 Current values of the numerical solution. 

539 t : float 

540 Current time the numerical solution is computed. 

541 

542 Returns 

543 ------- 

544 f : dtype_f 

545 Right-hand side of the problem. 

546 """ 

547 

548 f = self.dtype_f(self.init) 

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

550 

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

552 xa = self.init.getVecArray(u) 

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

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

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

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

557 

558 return f 

559 

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

561 r""" 

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

563 

564 Parameters 

565 ---------- 

566 rhs : dtype_f 

567 Right-hand side for the linear system. 

568 factor : float 

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

570 u0 : dtype_u 

571 Initial guess for the iterative solver. 

572 t : float 

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

574 

575 Returns 

576 ------- 

577 me : dtype_u 

578 Solution as mesh. 

579 """ 

580 

581 me = self.dtype_u(u0) 

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

583 self.ksp.solve(rhs, me) 

584 

585 self.ksp_ncalls += 1 

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

587 

588 return me 

589 

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

591 r""" 

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

593 

594 Parameters 

595 ---------- 

596 rhs : dtype_f 

597 Right-hand side for the linear system. 

598 factor : float 

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

600 u0 : dtype_u 

601 Initial guess for the iterative solver. 

602 t : float 

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

604 

605 Returns 

606 ------- 

607 me : dtype_u 

608 Solution as mesh. 

609 """ 

610 

611 me = self.dtype_u(u0) 

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

613 

614 F = self.init.createGlobalVec() 

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

616 J = self.init.createMatrix() 

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

618 

619 self.snes.solve(rhs, me) 

620 

621 self.snes_ncalls += 1 

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

623 

624 return me 

625 

626 def u_exact(self, t): 

627 r""" 

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

629 

630 Parameters 

631 ---------- 

632 t : float 

633 Time of the exact solution. 

634 

635 Returns 

636 ------- 

637 me : dtype_u 

638 Exact solution. 

639 """ 

640 

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

642 

643 me = self.dtype_u(self.init) 

644 xa = self.init.getVecArray(me) 

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

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

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

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

649 ) 

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

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

652 ) 

653 

654 return me 

655 

656 

657class petsc_grayscott_fullyimplicit(petsc_grayscott_multiimplicit): 

658 r""" 

659 2D periodic Gray-Scott system with PETSc finite differences, fully implicit with PETSc's SNES. 

660 

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

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

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

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

665 

666 .. math:: 

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

668 

669 .. math:: 

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

671 

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

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

674 solved by PETSc's nonlinear solver. 

675 """ 

676 

677 dtype_f = petsc_vec 

678 

679 def eval_f(self, u, t): 

680 """ 

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

682 

683 Parameters 

684 ---------- 

685 u : dtype_u 

686 Current values of the numerical solution. 

687 t : float 

688 Current time the numerical solution is computed. 

689 

690 Returns 

691 ------- 

692 f : dtype_f 

693 Right-hand side of the problem. 

694 """ 

695 

696 f = self.dtype_f(self.init) 

697 self.AMat.mult(u, f) 

698 

699 fa = self.init.getVecArray(f) 

700 xa = self.init.getVecArray(u) 

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

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

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

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

705 

706 return f 

707 

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

709 r""" 

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

711 

712 Parameters 

713 ---------- 

714 rhs : dtype_f 

715 Right-hand side for the linear system. 

716 factor : float 

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

718 u0 : dtype_u 

719 Initial guess for the iterative solver. 

720 t : float 

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

722 

723 Returns 

724 ------- 

725 me : dtype_u 

726 Solution as mesh. 

727 """ 

728 

729 me = self.dtype_u(u0) 

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

731 

732 # assign residual function and Jacobian 

733 F = self.init.createGlobalVec() 

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

735 J = self.init.createMatrix() 

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

737 

738 self.snes.solve(rhs, me) 

739 

740 self.snes_ncalls += 1 

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

742 

743 return me 

744 

745 

746class petsc_grayscott_semiimplicit(petsc_grayscott_multiimplicit): 

747 r""" 

748 2D periodic Gray-Scott system with PETSc finite differences, IMEX with diffusion implicit, reaction explicit. 

749 

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

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

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

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

754 

755 .. math:: 

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

757 

758 .. math:: 

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

760 

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

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

763 """ 

764 

765 dtype_f = petsc_vec_imex 

766 

767 def eval_f(self, u, t): 

768 """ 

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

770 

771 Parameters 

772 ---------- 

773 u : dtype_u 

774 Current values of the numerical solution. 

775 t : float 

776 Current time the numerical solution is computed. 

777 

778 Returns 

779 ------- 

780 f : dtype_f 

781 Right-hand side of the problem. 

782 """ 

783 

784 f = self.dtype_f(self.init) 

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

786 

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

788 xa = self.init.getVecArray(u) 

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

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

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

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

793 

794 return f 

795 

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