Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD_gpu.py: 0%

1import numpy as np 

2import cupy as cp 

3import cupyx.scipy.sparse as csp 

4from cupyx.scipy.sparse.linalg import cg # , spsolve, gmres, minres 

5 

6from pySDC.core.Errors import ProblemError 

7from pySDC.core.Problem import ptype 

8from pySDC.implementations.datatype_classes.cupy_mesh import cupy_mesh, imex_cupy_mesh, comp2_cupy_mesh 

9 

10# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf 

11 

12 

13class allencahn_fullyimplicit(ptype): # pragma: no cover 

14 r""" 

15 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

16 

17 .. math:: 

18 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) 

19 

20 for constant parameter :math:`\nu`. Initial condition are circles of the form 

21 

22 .. math:: 

23 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right) 

24 

25 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is 

26 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton. 

27 

28 This class is especially developed for solving it on GPUs using ``CuPy``. 

29 

30 Parameters 

31 ---------- 

32 nvars : tuple of int, optional 

33 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``. 

34 nu : float, optional 

35 Problem parameter :math:`\nu`. 

36 eps : float, optional 

37 Scaling parameter :math:`\varepsilon`. 

38 newton_maxiter : int, optional 

39 Maximum number of iterations for the Newton solver. 

40 newton_tol : float, optional 

41 Tolerance for Newton's method to terminate. 

42 lin_tol : float, optional 

43 Tolerance for linear solver to terminate. 

44 lin_maxiter : int, optional 

45 Maximum number of iterations for the linear solver. 

46 radius : float, optional 

47 Radius of the circles. 

48 

49 Attributes 

50 ---------- 

51 A : scipy.spdiags 

52 Second-order FD discretization of the 2D laplace operator. 

53 dx : float 

54 Distance between two spatial nodes (same for both directions). 

55 xvalues : np.1darray 

56 Spatial grid points, here both dimensions have the same grid points. 

57 newton_itercount : int 

58 Number of iterations of Newton solver. 

59 lin_itercount 

60 Number of iterations of linear solver. 

61 newton_ncalls : int 

62 Number of calls of Newton solver. 

63 lin_ncalls : int 

64 Number of calls of linear solver. 

65 """ 

66 

67 dtype_u = cupy_mesh 

68 dtype_f = cupy_mesh 

69 

70 def __init__( 

71 self, 

72 nvars=(128, 128), 

73 nu=2, 

74 eps=0.04, 

75 newton_maxiter=1e-12, 

76 newton_tol=100, 

77 lin_tol=1e-8, 

78 lin_maxiter=100, 

79 radius=0.25, 

80 ): 

81 """Initialization routine""" 

82 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on 

83 if len(nvars) != 2: 

84 raise ProblemError('this is a 2d example, got %s' % nvars) 

85 if nvars[0] != nvars[1]: 

86 raise ProblemError('need a square domain, got %s' % nvars) 

87 if nvars[0] % 2 != 0: 

88 raise ProblemError('the setup requires nvars = 2^p per dimension') 

89 

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

91 super().__init__((nvars, None, cp.dtype('float64'))) 

92 self._makeAttributeAndRegister( 

93 'nvars', 

94 'nu', 

95 'eps', 

96 'newton_maxiter', 

97 'newton_tol', 

98 'lin_tol', 

99 'lin_maxiter', 

100 'radius', 

101 localVars=locals(), 

102 readOnly=True, 

103 ) 

104 

105 # compute dx and get discretization matrix A 

106 self.dx = 1.0 / self.nvars[0] 

107 self.A = self.__get_A(self.nvars, self.dx) 

108 self.xvalues = cp.array([i * self.dx - 0.5 for i in range(self.nvars[0])]) 

109 

110 self.newton_itercount = 0 

111 self.lin_itercount = 0 

112 self.newton_ncalls = 0 

113 self.lin_ncalls = 0 

114 

115 @staticmethod 

116 def __get_A(N, dx): 

117 """ 

118 Helper function to assemble FD matrix A in sparse format. 

119 

120 Parameters 

121 ---------- 

122 N : list 

123 Number of degrees of freedom. 

124 dx : float 

125 Distance between two spatial nodes. 

126 

127 Returns 

128 ------- 

129 A : cupyx.scipy.sparse.csr_matrix 

130 CuPy-matrix A in CSR format. 

131 """ 

132 

133 stencil = cp.asarray([-2, 1]) 

134 A = stencil[0] * csp.eye(N[0], format='csr') 

135 for i in range(1, len(stencil)): 

136 A += stencil[i] * csp.eye(N[0], k=-i, format='csr') 

137 A += stencil[i] * csp.eye(N[0], k=+i, format='csr') 

138 A += stencil[i] * csp.eye(N[0], k=N[0] - i, format='csr') 

139 A += stencil[i] * csp.eye(N[0], k=-N[0] + i, format='csr') 

140 A = csp.kron(A, csp.eye(N[0])) + csp.kron(csp.eye(N[1]), A) 

141 A *= 1.0 / (dx**2) 

142 return A 

143 

144 # noinspection PyTypeChecker 

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

146 """ 

147 Simple Newton solver. 

148 

149 Parameters 

150 ---------- 

151 rhs : dtype_f 

152 Right-hand side for the nonlinear system. 

153 factor : float 

154 Abbrev. for the node-to-node stepsize (or any other factor required). 

155 u0 : dtype_u 

156 Initial guess for the iterative solver. 

157 t : float 

158 Current time (required here for the BC). 

159 

160 Returns 

161 ------- 

162 u : dtype_u 

163 The solution as mesh. 

164 """ 

165 

166 u = self.dtype_u(u0).flatten() 

167 z = self.dtype_u(self.init, val=0.0).flatten() 

168 nu = self.nu 

169 eps2 = self.eps**2 

170 

171 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

172 

173 # start newton iteration 

174 n = 0 

175 res = 99 

176 while n < self.newton_maxiter: 

177 # form the function g with g(u) = 0 

178 g = u - factor * (self.A.dot(u) + 1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten() 

179 

180 # if g is close to 0, then we are done 

181 res = cp.linalg.norm(g, np.inf) 

182 

183 if res < self.newton_tol: 

184 break 

185 

186 # assemble dg 

187 dg = Id - factor * (self.A + 1.0 / eps2 * csp.diags((1.0 - (nu + 1) * u**nu), offsets=0)) 

188 

189 # newton update: u1 = u0 - g/dg 

190 # u -= spsolve(dg, g) 

191 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0] 

192 # increase iteration count 

193 n += 1 

194 # print(n, res) 

195 

196 # if n == self.newton_maxiter: 

197 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

198 

199 me = self.dtype_u(self.init) 

200 me[:] = u.reshape(self.nvars) 

201 

202 self.newton_ncalls += 1 

203 self.newton_itercount += n 

204 

205 return me 

206 

207 def eval_f(self, u, t): 

208 """ 

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

210 

211 Parameters 

212 ---------- 

213 u : dtype_u 

214 Current values of the numerical solution. 

215 t : float 

216 Current time of the numerical solution is computed (not used here). 

217 

218 Returns 

219 ------- 

220 f : dtype_f 

221 The right-hand side of the problem. 

222 """ 

223 f = self.dtype_f(self.init) 

224 v = u.flatten() 

225 f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars) 

226 

227 return f 

228 

229 def u_exact(self, t): 

230 r""" 

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

232 

233 Parameters 

234 ---------- 

235 t : float 

236 Time of the exact solution. 

237 

238 Returns 

239 ------- 

240 me : dtype_u 

241 The exact solution. 

242 """ 

243 

244 assert t == 0, 'ERROR: u_exact only valid for t=0' 

245 me = self.dtype_u(self.init, val=0.0) 

246 mx, my = cp.meshgrid(self.xvalues, self.xvalues) 

247 me[:] = cp.tanh((self.radius - cp.sqrt(mx**2 + my**2)) / (cp.sqrt(2) * self.eps)) 

248 # print(type(me)) 

249 return me 

250 

251 

252# noinspection PyUnusedLocal 

253class allencahn_semiimplicit(allencahn_fullyimplicit): 

254 r""" 

255 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

256 

257 .. math:: 

258 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) 

259 

260 for constant parameter :math:`\nu`. Initial condition are circles of the form 

261 

262 .. math:: 

263 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right) 

264 

265 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is 

266 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the 

267 conjugate gradients method, and the system containing the rest of the right-hand side is only evaluated at each time. 

268 

269 This class is especially developed for solving it on GPUs using ``CuPy``. 

270 """ 

271 

272 dtype_f = imex_cupy_mesh 

273 

274 def eval_f(self, u, t): 

275 """ 

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

277 

278 Parameters 

279 ---------- 

280 u : dtype_u 

281 Current values of the numerical solution. 

282 t : float 

283 Current time of the numerical solution is computed. 

284 

285 Returns 

286 ------- 

287 f : dtype_f 

288 The right-hand side of the problem. 

289 """ 

290 f = self.dtype_f(self.init) 

291 v = u.flatten() 

292 f.impl[:] = self.A.dot(v).reshape(self.nvars) 

293 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars) 

294 

295 return f 

296 

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

298 r""" 

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

300 

301 Parameters 

302 ---------- 

303 rhs : dtype_f 

304 Right-hand side for the linear system. 

305 factor : float 

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

307 u0 : dtype_u 

308 Initial guess for the iterative solver. 

309 t : float 

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

311 

312 Returns 

313 ------- 

314 me : dtype_u 

315 The solution as mesh. 

316 """ 

317 

318 class context: 

319 num_iter = 0 

320 

321 def callback(xk): 

322 context.num_iter += 1 

323 return context.num_iter 

324 

325 me = self.dtype_u(self.init) 

326 

327 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

328 

329 me[:] = cg( 

330 Id - factor * self.A, 

331 rhs.flatten(), 

332 x0=u0.flatten(), 

333 tol=self.lin_tol, 

334 maxiter=self.lin_maxiter, 

335 callback=callback, 

336 )[0].reshape(self.nvars) 

337 

338 self.lin_ncalls += 1 

339 self.lin_itercount += context.num_iter 

340 

341 return me 

342 

343 

344# noinspection PyUnusedLocal 

345class allencahn_semiimplicit_v2(allencahn_fullyimplicit): 

346 r""" 

347 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

348 

349 .. math:: 

350 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) 

351 

352 for constant parameter :math:`\nu`. Initial condition are circles of the form 

353 

354 .. math:: 

355 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right) 

356 

357 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting 

358 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}` 

359 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u` 

360 is only evaluated at each time. 

361 

362 This class is especially developed for solving it on GPUs using ``CuPy``. 

363 """ 

364 

365 dtype_f = imex_cupy_mesh 

366 

367 def eval_f(self, u, t): 

368 """ 

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

370 

371 Parameters 

372 ---------- 

373 u : dtype_u 

374 Current values of the numerical solution. 

375 t : float 

376 Current time of the numerical solution is computed. 

377 

378 Returns 

379 ------- 

380 f : dtype_f 

381 The right-hand side of the problem. 

382 """ 

383 f = self.dtype_f(self.init) 

384 v = u.flatten() 

385 f.impl[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars) 

386 f.expl[:] = (1.0 / self.eps**2 * v).reshape(self.nvars) 

387 

388 return f 

389 

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

391 """ 

392 Simple Newton solver. 

393 

394 Parameters 

395 ---------- 

396 rhs : dtype_f 

397 Right-hand side for the nonlinear system. 

398 factor : float 

399 Abbrev. for the node-to-node stepsize (or any other factor required). 

400 u0 : dtype_u 

401 Initial guess for the iterative solver. 

402 t : float 

403 Current time (required here for the BC). 

404 

405 Returns 

406 ------- 

407 me : dtype_u 

408 The solution as mesh. 

409 """ 

410 

411 u = self.dtype_u(u0).flatten() 

412 z = self.dtype_u(self.init, val=0.0).flatten() 

413 nu = self.nu 

414 eps2 = self.eps**2 

415 

416 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

417 

418 # start newton iteration 

419 n = 0 

420 res = 99 

421 while n < self.newton_maxiter: 

422 # form the function g with g(u) = 0 

423 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten() 

424 

425 # if g is close to 0, then we are done 

426 res = cp.linalg.norm(g, np.inf) 

427 

428 if res < self.newton_tol: 

429 break 

430 

431 # assemble dg 

432 dg = Id - factor * (self.A - 1.0 / eps2 * csp.diags(((nu + 1) * u**nu), offsets=0)) 

433 

434 # newton update: u1 = u0 - g/dg 

435 # u -= spsolve(dg, g) 

436 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0] 

437 # increase iteration count 

438 n += 1 

439 # print(n, res) 

440 

441 # if n == self.newton_maxiter: 

442 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

443 

444 me = self.dtype_u(self.init) 

445 me[:] = u.reshape(self.nvars) 

446 

447 self.newton_ncalls += 1 

448 self.newton_itercount += n 

449 

450 return me 

451 

452 

453# noinspection PyUnusedLocal 

454class allencahn_multiimplicit(allencahn_fullyimplicit): 

455 r""" 

456 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

457 

458 .. math:: 

459 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) 

460 

461 for constant parameter :math:`\nu`. Initial condition are circles of the form 

462 

463 .. math:: 

464 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right) 

465 

466 for :math::`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is 

467 treated in *multi-implicit* fashion, i.e., the nonlinear system containing the second order term is solved by the 

468 conjugate gradients method, and the system containing the rest of the right-hand side will be solved by Newton's method. 

469 

470 This class is especially developed for solving it on GPUs using ``CuPy``. 

471 """ 

472 

473 dtype_f = comp2_cupy_mesh 

474 

475 def eval_f(self, u, t): 

476 """ 

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

478 

479 Parameters 

480 ---------- 

481 u : dtype_u 

482 Current values of the numerical solution. 

483 t : float 

484 Current time of the numerical solution is computed (not used here). 

485 

486 Returns 

487 ------- 

488 f : dtype_f 

489 The right-hand side of the problem. 

490 """ 

491 f = self.dtype_f(self.init) 

492 v = u.flatten() 

493 f.comp1[:] = self.A.dot(v).reshape(self.nvars) 

494 f.comp2[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars) 

495 

496 return f 

497 

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

499 r""" 

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

501 

502 Parameters 

503 ---------- 

504 rhs : dtype_f 

505 Right-hand side for the linear system. 

506 factor : float 

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

508 u0 : dtype_u 

509 Initial guess for the iterative solver. 

510 t : float 

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

512 

513 Returns 

514 ------- 

515 me : dtype_u 

516 The solution as mesh. 

517 """ 

518 

519 class context: 

520 num_iter = 0 

521 

522 def callback(xk): 

523 context.num_iter += 1 

524 return context.num_iter 

525 

526 me = self.dtype_u(self.init) 

527 

528 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

529 

530 me[:] = cg( 

531 Id - factor * self.A, 

532 rhs.flatten(), 

533 x0=u0.flatten(), 

534 tol=self.lin_tol, 

535 maxiter=self.lin_maxiter, 

536 callback=callback, 

537 )[0].reshape(self.nvars) 

538 

539 self.lin_ncalls += 1 

540 self.lin_itercount += context.num_iter 

541 

542 return me 

543 

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

545 """ 

546 Simple Newton solver. 

547 

548 Parameters 

549 ---------- 

550 rhs : dtype_f 

551 Right-hand side for the nonlinear system 

552 factor : float 

553 Abbrev. for the node-to-node stepsize (or any other factor required). 

554 u0 : dtype_u 

555 Initial guess for the iterative solver. 

556 t : float 

557 Current time (required here for the BC). 

558 

559 Returns 

560 ------- 

561 me : dtype_u 

562 The solution as mesh. 

563 """ 

564 

565 u = self.dtype_u(u0).flatten() 

566 z = self.dtype_u(self.init, val=0.0).flatten() 

567 nu = self.nu 

568 eps2 = self.eps**2 

569 

570 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

571 

572 # start newton iteration 

573 n = 0 

574 res = 99 

575 while n < self.newton_maxiter: 

576 # form the function g with g(u) = 0 

577 g = u - factor * (1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten() 

578 

579 # if g is close to 0, then we are done 

580 res = cp.linalg.norm(g, np.inf) 

581 

582 if res < self.newton_tol: 

583 break 

584 

585 # assemble dg 

586 dg = Id - factor * (1.0 / eps2 * csp.diags((1.0 - (nu + 1) * u**nu), offsets=0)) 

587 

588 # newton update: u1 = u0 - g/dg 

589 # u -= spsolve(dg, g) 

590 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0] 

591 # increase iteration count 

592 n += 1 

593 # print(n, res) 

594 

595 # if n == self.newton_maxiter: 

596 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

597 

598 me = self.dtype_u(self.init) 

599 me[:] = u.reshape(self.nvars) 

600 

601 self.newton_ncalls += 1 

602 self.newton_itercount += n 

603 

604 return me 

605 

606 

607# noinspection PyUnusedLocal 

608class allencahn_multiimplicit_v2(allencahn_fullyimplicit): 

609 r""" 

610 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2` 

611 

612 .. math:: 

613 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu) 

614 

615 for constant parameter :math:`\nu`. The initial condition has the form of circles 

616 

617 .. math:: 

618 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right) 

619 

620 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting 

621 is used here to get another kind of *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}` 

622 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u` 

623 is solved by a linear solver provided by a ``SciPy`` routine. 

624 

625 This class is especially developed for solving it on GPUs using ``CuPy``. 

626 """ 

627 

628 dtype_f = comp2_cupy_mesh 

629 

630 def eval_f(self, u, t): 

631 """ 

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

633 

634 Parameters 

635 ---------- 

636 u : dtype_u 

637 Current values of the numerical solution. 

638 t : float 

639 Current time of the numerical solution is computed (not used here). 

640 

641 Returns 

642 ------- 

643 f : dtype_f 

644 The right-hand side of the problem. 

645 """ 

646 f = self.dtype_f(self.init) 

647 v = u.flatten() 

648 f.comp1[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars) 

649 f.comp2[:] = (1.0 / self.eps**2 * v).reshape(self.nvars) 

650 

651 return f 

652 

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

654 """ 

655 Simple Newton solver. 

656 

657 Parameters 

658 ---------- 

659 rhs : dtype_f 

660 Right-hand side for the nonlinear system. 

661 factor : float 

662 Abbrev. for the node-to-node stepsize (or any other factor required). 

663 u0 : dtype_u 

664 Initial guess for the iterative solver. 

665 t : float 

666 Current time (required here for the BC). 

667 

668 Returns 

669 ------- 

670 me : dtype_u 

671 The solution as mesh. 

672 """ 

673 

674 u = self.dtype_u(u0).flatten() 

675 z = self.dtype_u(self.init, val=0.0).flatten() 

676 nu = self.nu 

677 eps2 = self.eps**2 

678 

679 Id = csp.eye(self.nvars[0] * self.nvars[1]) 

680 

681 # start newton iteration 

682 n = 0 

683 res = 99 

684 while n < self.newton_maxiter: 

685 # form the function g with g(u) = 0 

686 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten() 

687 

688 # if g is close to 0, then we are done 

689 res = cp.linalg.norm(g, np.inf) 

690 

691 if res < self.newton_tol: 

692 break 

693 

694 # assemble dg 

695 dg = Id - factor * (self.A - 1.0 / eps2 * csp.diags(((nu + 1) * u**nu), offsets=0)) 

696 

697 # newton update: u1 = u0 - g/dg 

698 # u -= spsolve(dg, g) 

699 u -= cg(dg, g, x0=z, tol=self.lin_tol)[0] 

700 # increase iteration count 

701 n += 1 

702 # print(n, res) 

703 

704 # if n == self.newton_maxiter: 

705 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

706 

707 me = self.dtype_u(self.init) 

708 me[:] = u.reshape(self.nvars) 

709 

710 self.newton_ncalls += 1 

711 self.newton_itercount += n 

712 

713 return me 

714 

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

716 r""" 

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

718 

719 Parameters 

720 ---------- 

721 rhs : dtype_f 

722 Right-hand side for the linear system. 

723 factor : float 

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

725 u0 : dtype_u 

726 Initial guess for the iterative solver. 

727 t : float 

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

729 

730 Returns 

731 ------- 

732 me : dtype_u 

733 The solution as mesh. 

734 """ 

735 

736 me = self.dtype_u(self.init) 

737 

738 me[:] = (1.0 / (1.0 - factor * 1.0 / self.eps**2) * rhs).reshape(self.nvars) 

739 return me