Coverage for pySDC/implementations/problem_classes/AllenCahn_1D_FD.py: 96%

1import numpy as np 

2import scipy.sparse as sp 

3from scipy.sparse.linalg import spsolve 

4 

5from pySDC.core.errors import ProblemError 

6from pySDC.core.problem import Problem, WorkCounter 

7from pySDC.helpers import problem_helper 

8from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh 

9 

10 

11class allencahn_front_fullyimplicit(Problem): 

12 r""" 

13 Example implementing the one-dimensional Allen-Cahn equation with driving force using inhomogeneous Dirichlet 

14 boundary conditions 

15 

16 .. math:: 

17 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

18 - 6 d_w u (1 - u) 

19 

20 for :math:`u \in [0, 1]`. The second order spatial derivative is approximated using centered finite differences. The 

21 exact solution is given by 

22 

23 .. math:: 

24 u(x, t)= 0.5 \left(1 + \tanh\left(\frac{x - vt}{\sqrt{2}\varepsilon}\right)\right) 

25 

26 with :math:`v = 3 \sqrt{2} \varepsilon d_w`. For time-stepping, this problem is implemented to be treated 

27 *fully-implicit* using Newton to solve the nonlinear system. 

28 

29 Parameters 

30 ---------- 

31 nvars : int 

32 Number of unknowns in the problem. 

33 dw : float 

34 Driving force. 

35 eps : float 

36 Scaling parameter :math:`\varepsilon`. 

37 newton_maxiter : int 

38 Maximum number of iterations for Newton's method. 

39 newton_tol : float 

40 Tolerance for Newton's method to terminate. 

41 interval : list 

42 Interval of spatial domain. 

43 stop_at_nan : bool, optional 

44 Indicates that the Newton solver should stop if ``nan`` values arise. 

45 

46 Attributes 

47 ---------- 

48 A : scipy.diags 

49 Second-order FD discretization of the 1D laplace operator. 

50 dx : float 

51 Distance between two spatial nodes. 

52 xvalues : np.1darray 

53 Spatial grid values. 

54 uext : dtype_u 

55 Contains additionally the external values of the boundary. 

56 work_counters : WorkCounter 

57 Counter for statistics. Here, number of Newton calls and number of evaluations 

58 of right-hand side are counted. 

59 """ 

60 

61 dtype_u = mesh 

62 dtype_f = mesh 

63 

64 def __init__( 

65 self, 

66 nvars=127, 

67 dw=-0.04, 

68 eps=0.04, 

69 newton_maxiter=100, 

70 newton_tol=1e-12, 

71 interval=(-0.5, 0.5), 

72 stop_at_nan=True, 

73 stop_at_maxiter=False, 

74 ): 

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

76 if (nvars + 1) % 2: 

77 raise ProblemError('setup requires nvars = 2^p - 1') 

78 

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

80 super().__init__((nvars, None, np.dtype('float64'))) 

81 self._makeAttributeAndRegister( 

82 'nvars', 

83 'dw', 

84 'eps', 

85 'newton_maxiter', 

86 'newton_tol', 

87 'interval', 

88 'stop_at_nan', 

89 'stop_at_maxiter', 

90 localVars=locals(), 

91 readOnly=True, 

92 ) 

93 

94 # compute dx and get discretization matrix A 

95 self.dx = (self.interval[1] - self.interval[0]) / (self.nvars + 1) 

96 self.xvalues = np.array([(i + 1 - (self.nvars + 1) / 2) * self.dx for i in range(self.nvars)]) 

97 

98 self.A, _ = problem_helper.get_finite_difference_matrix( 

99 derivative=2, 

100 order=2, 

101 stencil_type='center', 

102 dx=self.dx, 

103 size=self.nvars + 2, 

104 dim=1, 

105 bc='dirichlet-zero', 

106 ) 

107 self.uext = self.dtype_u((self.init[0] + 2, self.init[1], self.init[2]), val=0.0) 

108 

109 self.work_counters['newton'] = WorkCounter() 

110 self.work_counters['rhs'] = WorkCounter() 

111 

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

113 """ 

114 Simple Newton solver. 

115 

116 Parameters 

117 ---------- 

118 rhs : dtype_f 

119 Right-hand side for the nonlinear system. 

120 factor : float 

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

122 u0 : dtype_u 

123 Initial guess for the iterative solver. 

124 t : float 

125 Current time (required here for the BC). 

126 

127 Returns 

128 ------- 

129 me : dtype_u 

130 The solution as mesh. 

131 """ 

132 

133 u = self.dtype_u(u0) 

134 eps2 = self.eps**2 

135 dw = self.dw 

136 

137 Id = sp.eye(self.nvars) 

138 

139 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

140 self.uext[0] = 0.5 * (1 + np.tanh((self.interval[0] - v * t) / (np.sqrt(2) * self.eps))) 

141 self.uext[-1] = 0.5 * (1 + np.tanh((self.interval[1] - v * t) / (np.sqrt(2) * self.eps))) 

142 

143 A = self.A[1:-1, 1:-1] 

144 # start newton iteration 

145 n = 0 

146 res = 99 

147 while n < self.newton_maxiter: 

148 # print(n) 

149 # # form the function g(u), such that the solution to the nonlinear problem is a root of g 

150 self.uext[1:-1] = u[:] 

151 g = ( 

152 u 

153 - rhs 

154 - factor 

155 * ( 

156 self.A.dot(self.uext)[1:-1] 

157 - 2.0 / eps2 * u * (1.0 - u) * (1.0 - 2.0 * u) 

158 - 6.0 * dw * u * (1.0 - u) 

159 ) 

160 ) 

161 

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

163 res = np.linalg.norm(g, np.inf) 

164 

165 if res < self.newton_tol: 

166 break 

167 

168 # assemble dg 

169 dg = Id - factor * ( 

170 A 

171 - 2.0 

172 / eps2 

173 * sp.diags((1.0 - u) * (1.0 - 2.0 * u) - u * ((1.0 - 2.0 * u) + 2.0 * (1.0 - u)), offsets=0) 

174 - 6.0 * dw * sp.diags((1.0 - u) - u, offsets=0) 

175 ) 

176 

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

178 u -= spsolve(dg, g) 

179 # u -= gmres(dg, g, x0=z, rtol=self.lin_tol)[0] 

180 # increase iteration count 

181 n += 1 

182 self.work_counters['newton']() 

183 

184 if np.isnan(res) and self.stop_at_nan: 

185 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

186 elif np.isnan(res): 

187 self.logger.warning('Newton got nan after %i iterations...' % n) 

188 

189 if n == self.newton_maxiter: 

190 msg = 'Newton did not converge after %i iterations, error is %s' % (n, res) 

191 if self.stop_at_maxiter: 

192 raise ProblemError(msg) 

193 else: 

194 self.logger.warning(msg) 

195 

196 me = self.dtype_u(self.init) 

197 me[:] = u[:] 

198 

199 return me 

200 

201 def eval_f(self, u, t): 

202 """ 

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

204 

205 Parameters 

206 ---------- 

207 u : dtype_u 

208 Current values of the numerical solution. 

209 t : float 

210 Current time of the numerical solution is computed. 

211 

212 Returns 

213 ------- 

214 f : dtype_f 

215 The right-hand side of the problem. 

216 """ 

217 # set up boundary values to embed inner points 

218 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

219 self.uext[0] = 0.5 * (1 + np.tanh((self.interval[0] - v * t) / (np.sqrt(2) * self.eps))) 

220 self.uext[-1] = 0.5 * (1 + np.tanh((self.interval[1] - v * t) / (np.sqrt(2) * self.eps))) 

221 

222 self.uext[1:-1] = u[:] 

223 

224 f = self.dtype_f(self.init) 

225 f[:] = ( 

226 self.A.dot(self.uext)[1:-1] 

227 - 2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u) 

228 - 6.0 * self.dw * u * (1.0 - u) 

229 ) 

230 self.work_counters['rhs']() 

231 return f 

232 

233 def u_exact(self, t): 

234 r""" 

235 Routine to return initial condition or the exact solution 

236 

237 Parameters 

238 ---------- 

239 t : float 

240 Time at which the exact solution is computed. 

241 

242 Returns 

243 ------- 

244 me : dtype_u 

245 The exact solution (in space and time). 

246 """ 

247 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

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

249 me[:] = 0.5 * (1 + np.tanh((self.xvalues - v * t) / (np.sqrt(2) * self.eps))) 

250 return me 

251 

252 

253class allencahn_front_semiimplicit(allencahn_front_fullyimplicit): 

254 r""" 

255 This class implements the one-dimensional Allen-Cahn equation with driving force using inhomogeneous Dirichlet 

256 boundary conditions 

257 

258 .. math:: 

259 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

260 - 6 d_w u (1 - u) 

261 

262 for :math:`u \in [0, 1]`. Centered finite differences are used for discretization of the second order spatial derivative. 

263 The exact solution is given by 

264 

265 .. math:: 

266 u(x, t) = 0.5 \left(1 + \tanh\left(\frac{x - vt}{\sqrt{2}\varepsilon}\right)\right) 

267 

268 with :math:`v = 3 \sqrt{2} \varepsilon d_w`. For time-stepping, this problem will be treated in a 

269 *semi-implicit* way, i.e., the Laplacian is treated implicitly, and the rest of the right-hand side will be handled 

270 explicitly. 

271 """ 

272 

273 dtype_f = imex_mesh 

274 

275 def eval_f(self, u, t): 

276 """ 

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

278 

279 Parameters 

280 ---------- 

281 u : dtype_u 

282 Current values of the numerical solution. 

283 t : float 

284 Current time of the numerical solution is computed. 

285 

286 Returns 

287 ------- 

288 f : dtype_f 

289 The right-hand side of the problem. 

290 """ 

291 # set up boundary values to embed inner points 

292 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

293 self.uext[0] = 0.5 * (1 + np.tanh((self.interval[0] - v * t) / (np.sqrt(2) * self.eps))) 

294 self.uext[-1] = 0.5 * (1 + np.tanh((self.interval[1] - v * t) / (np.sqrt(2) * self.eps))) 

295 

296 self.uext[1:-1] = u[:] 

297 

298 f = self.dtype_f(self.init) 

299 f.impl[:] = self.A.dot(self.uext)[1:-1] 

300 f.expl[:] = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u) - 6.0 * self.dw * u * (1.0 - u) 

301 self.work_counters['rhs']() 

302 return f 

303 

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

305 r""" 

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

307 

308 Parameters 

309 ---------- 

310 rhs : dtype_f 

311 Right-hand side for the linear system. 

312 factor : float 

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

314 u0 : dtype_u 

315 Initial guess for the iterative solver. 

316 t : float 

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

318 

319 Returns 

320 ------- 

321 me : dtype_u 

322 The solution as mesh. 

323 """ 

324 

325 me = self.dtype_u(self.init) 

326 self.uext[0] = 0.0 

327 self.uext[-1] = 0.0 

328 self.uext[1:-1] = rhs[:] 

329 me[:] = spsolve(sp.eye(self.nvars + 2, format='csc') - factor * self.A, self.uext)[1:-1] 

330 return me 

331 

332 

333class allencahn_front_finel(allencahn_front_fullyimplicit): 

334 r""" 

335 This class implements the one-dimensional Allen-Cahn equation with driving force using inhomogeneous Dirichlet 

336 boundary conditions 

337 

338 .. math:: 

339 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

340 - 6 d_w u (1 - u) 

341 

342 for :math:`u \in [0, 1]`. Centered finite differences are used for discretization of the Laplacian. 

343 The exact solution is given by 

344 

345 .. math:: 

346 u(x, t) = 0.5 \left(1 + \tanh\left(\frac{x - vt}{\sqrt{2}\varepsilon}\right)\right) 

347 

348 with :math:`v = 3 \sqrt{2} \varepsilon d_w`. 

349 

350 Let :math:`A` denote the finite difference matrix to discretize :math:`\frac{\partial^2 u}{\partial x^2}`. Here, 

351 *Finel's trick* is used. Let 

352 

353 .. math:: 

354 a = \tanh\left(\frac{\Delta x}{\sqrt{2}\varepsilon}\right)^2, 

355 

356 then, the right-hand side of the problem can be written as 

357 

358 .. math:: 

359 \frac{\partial u}{\partial t} = A u - \frac{1}{\Delta x^2} \left[ 

360 \frac{1 - a}{1 - a (2u - 1)^2} - 1 

361 \right] (2u - 1). 

362 

363 For time-stepping, this problem will be treated in a *fully-implicit* way. The nonlinear system is solved using Newton. 

364 """ 

365 

366 # noinspection PyTypeChecker 

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

368 """ 

369 Simple Newton solver. 

370 

371 Parameters 

372 ---------- 

373 rhs : dtype_f 

374 Right-hand side for the nonlinear system. 

375 factor : float 

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

377 u0 : dtype_u 

378 Initial guess for the iterative solver. 

379 t : float 

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

381 

382 Returns 

383 ------- 

384 me : dtype_u 

385 The solution as mesh. 

386 """ 

387 

388 u = self.dtype_u(u0) 

389 dw = self.dw 

390 a2 = np.tanh(self.dx / (np.sqrt(2) * self.eps)) ** 2 

391 

392 Id = sp.eye(self.nvars) 

393 

394 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

395 self.uext[0] = 0.5 * (1 + np.tanh((self.interval[0] - v * t) / (np.sqrt(2) * self.eps))) 

396 self.uext[-1] = 0.5 * (1 + np.tanh((self.interval[1] - v * t) / (np.sqrt(2) * self.eps))) 

397 

398 A = self.A[1:-1, 1:-1] 

399 # start newton iteration 

400 n = 0 

401 res = 99 

402 while n < self.newton_maxiter: 

403 # form the function g(u), such that the solution to the nonlinear problem is a root of g 

404 self.uext[1:-1] = u[:] 

405 gprim = 1.0 / self.dx**2 * ((1.0 - a2) / (1.0 - a2 * (2.0 * u - 1.0) ** 2) - 1.0) * (2.0 * u - 1.0) 

406 g = u - rhs - factor * (self.A.dot(self.uext)[1:-1] - 1.0 * gprim - 6.0 * dw * u * (1.0 - u)) 

407 

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

409 res = np.linalg.norm(g, np.inf) 

410 

411 if res < self.newton_tol: 

412 break 

413 

414 # assemble dg 

415 dgprim = ( 

416 1.0 

417 / self.dx**2 

418 * ( 

419 2.0 * ((1.0 - a2) / (1.0 - a2 * (2.0 * u - 1.0) ** 2) - 1.0) 

420 + (2.0 * u - 1) ** 2 * (1.0 - a2) * 4 * a2 / (1.0 - a2 * (2.0 * u - 1.0) ** 2) ** 2 

421 ) 

422 ) 

423 

424 dg = Id - factor * (A - 1.0 * sp.diags(dgprim, offsets=0) - 6.0 * dw * sp.diags((1.0 - u) - u, offsets=0)) 

425 

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

427 u -= spsolve(dg, g) 

428 # For some reason, doing cg or gmres does not work so well here... 

429 # u -= cg(dg, g, x0=z, rtol=self.lin_tol)[0] 

430 # increase iteration count 

431 n += 1 

432 self.work_counters['newton']() 

433 

434 if np.isnan(res) and self.stop_at_nan: 

435 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

436 elif np.isnan(res): 

437 self.logger.warning('Newton got nan after %i iterations...' % n) 

438 

439 if n == self.newton_maxiter: 

440 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res)) 

441 

442 me = self.dtype_u(self.init) 

443 me[:] = u[:] 

444 

445 return me 

446 

447 def eval_f(self, u, t): 

448 """ 

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

450 

451 Parameters 

452 ---------- 

453 u : dtype_u 

454 Current values of the numerical solution. 

455 t : float 

456 Current time of the numerical solution is computed. 

457 

458 Returns 

459 ------- 

460 f : dtype_f 

461 The right-hand side of the problem. 

462 """ 

463 # set up boundary values to embed inner points 

464 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

465 self.uext[0] = 0.5 * (1 + np.tanh((self.interval[0] - v * t) / (np.sqrt(2) * self.eps))) 

466 self.uext[-1] = 0.5 * (1 + np.tanh((self.interval[1] - v * t) / (np.sqrt(2) * self.eps))) 

467 

468 self.uext[1:-1] = u[:] 

469 

470 a2 = np.tanh(self.dx / (np.sqrt(2) * self.eps)) ** 2 

471 gprim = 1.0 / self.dx**2 * ((1.0 - a2) / (1.0 - a2 * (2.0 * u - 1.0) ** 2) - 1) * (2.0 * u - 1.0) 

472 f = self.dtype_f(self.init) 

473 f[:] = self.A.dot(self.uext)[1:-1] - 1.0 * gprim - 6.0 * self.dw * u * (1.0 - u) 

474 self.work_counters['rhs']() 

475 return f 

476 

477 

478class allencahn_periodic_fullyimplicit(Problem): 

479 r""" 

480 Example implementing the one-dimensional Allen-Cahn equation with driving force and periodic boundary conditions 

481 

482 .. math:: 

483 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

484 - 6 d_w u (1 - u) 

485 

486 for :math:`u \in [0, 1]`. Centered finite differences are used for discretization of the Laplacian. 

487 The exact solution is 

488 

489 .. math:: 

490 u(x, t) = 0.5 \left(1 + \tanh\left(\frac{r - |x| - vt}{\sqrt{2}\varepsilon}\right)\right) 

491 

492 with :math:`v = 3 \sqrt{2} \varepsilon d_w` and radius :math:`r` of the circles. For time-stepping, the problem is treated 

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

494 

495 Parameters 

496 ---------- 

497 nvars : int 

498 Number of unknowns in the problem. 

499 dw : float 

500 Driving force. 

501 eps : float 

502 Scaling parameter :math:`\varepsilon`. 

503 newton_maxiter : int 

504 Maximum number of iterations for Newton's method. 

505 newton_tol : float 

506 Tolerance for Newton's method to terminate. 

507 interval : list 

508 Interval of spatial domain. 

509 radius : float 

510 Radius of the circles. 

511 stop_at_nan : bool, optional 

512 Indicates that the Newton solver should stop if nan values arise. 

513 

514 Attributes 

515 ---------- 

516 A : scipy.diags 

517 Second-order FD discretization of the 1D laplace operator. 

518 dx : float 

519 Distance between two spatial nodes. 

520 xvalues : np.1darray 

521 Spatial grid points. 

522 work_counters : WorkCounter 

523 Counter for statistics. Here, number of Newton calls and number of evaluations 

524 of right-hand side are counted. 

525 """ 

526 

527 dtype_u = mesh 

528 dtype_f = mesh 

529 

530 def __init__( 

531 self, 

532 nvars=128, 

533 dw=-0.04, 

534 eps=0.04, 

535 newton_maxiter=100, 

536 newton_tol=1e-12, 

537 interval=(-0.5, 0.5), 

538 radius=0.25, 

539 stop_at_nan=True, 

540 ): 

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

542 if (nvars) % 2: 

543 raise ProblemError('setup requires nvars = 2^p') 

544 

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

546 super().__init__((nvars, None, np.dtype('float64'))) 

547 self._makeAttributeAndRegister( 

548 'nvars', 

549 'dw', 

550 'eps', 

551 'newton_maxiter', 

552 'newton_tol', 

553 'interval', 

554 'radius', 

555 'stop_at_nan', 

556 localVars=locals(), 

557 readOnly=True, 

558 ) 

559 

560 # compute dx and get discretization matrix A 

561 self.dx = (self.interval[1] - self.interval[0]) / self.nvars 

562 self.xvalues = np.array([self.interval[0] + i * self.dx for i in range(self.nvars)]) 

563 

564 self.A, _ = problem_helper.get_finite_difference_matrix( 

565 derivative=2, 

566 order=2, 

567 stencil_type='center', 

568 dx=self.dx, 

569 size=self.nvars, 

570 dim=1, 

571 bc='periodic', 

572 ) 

573 

574 self.work_counters['newton'] = WorkCounter() 

575 self.work_counters['rhs'] = WorkCounter() 

576 

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

578 """ 

579 Simple Newton solver. 

580 

581 Parameters 

582 ---------- 

583 rhs : dtype_f 

584 Right-hand side for the nonlinear system. 

585 factor : float 

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

587 u0 : dtype_u 

588 Initial guess for the iterative solver. 

589 t : float 

590 Current time (required here for the BC). 

591 

592 Returns 

593 ------- 

594 u : dtype_u 

595 The solution as mesh. 

596 """ 

597 

598 u = self.dtype_u(u0) 

599 eps2 = self.eps**2 

600 dw = self.dw 

601 

602 Id = sp.eye(self.nvars) 

603 

604 # start newton iteration 

605 n = 0 

606 res = 99 

607 while n < self.newton_maxiter: 

608 # form the function g(u), such that the solution to the nonlinear problem is a root of g 

609 g = ( 

610 u 

611 - rhs 

612 - factor * (self.A.dot(u) - 2.0 / eps2 * u * (1.0 - u) * (1.0 - 2.0 * u) - 6.0 * dw * u * (1.0 - u)) 

613 ) 

614 

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

616 res = np.linalg.norm(g, np.inf) 

617 

618 if res < self.newton_tol: 

619 break 

620 

621 # assemble dg 

622 dg = Id - factor * ( 

623 self.A 

624 - 2.0 

625 / eps2 

626 * sp.diags((1.0 - u) * (1.0 - 2.0 * u) - u * ((1.0 - 2.0 * u) + 2.0 * (1.0 - u)), offsets=0) 

627 - 6.0 * dw * sp.diags((1.0 - u) - u, offsets=0) 

628 ) 

629 

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

631 u -= spsolve(dg, g) 

632 # u -= gmres(dg, g, x0=z, rtol=self.lin_tol)[0] 

633 # increase iteration count 

634 n += 1 

635 self.work_counters['newton']() 

636 

637 if np.isnan(res) and self.stop_at_nan: 

638 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

639 elif np.isnan(res): 

640 self.logger.warning('Newton got nan after %i iterations...' % n) 

641 

642 if n == self.newton_maxiter: 

643 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res)) 

644 

645 me = self.dtype_u(self.init) 

646 me[:] = u[:] 

647 

648 return me 

649 

650 def eval_f(self, u, t): 

651 """ 

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

653 

654 Parameters 

655 ---------- 

656 u : dtype_u 

657 Current values of the numerical solution. 

658 t : float 

659 Current time of the numerical solution is computed. 

660 

661 Returns 

662 ------- 

663 f : dtype_f 

664 The right-hand side of the problem. 

665 """ 

666 f = self.dtype_f(self.init) 

667 f[:] = self.A.dot(u) - 2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u) - 6.0 * self.dw * u * (1.0 - u) 

668 self.work_counters['rhs']() 

669 return f 

670 

671 def u_exact(self, t): 

672 r""" 

673 Routine to return initial condition or the exact solution. 

674 

675 Parameters 

676 ---------- 

677 t : float 

678 Time at which the approximated exact solution is computed. 

679 

680 Returns 

681 ------- 

682 me : dtype_u 

683 The approximated exact solution. 

684 """ 

685 v = 3.0 * np.sqrt(2) * self.eps * self.dw 

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

687 me[:] = 0.5 * (1 + np.tanh((self.radius - abs(self.xvalues) - v * t) / (np.sqrt(2) * self.eps))) 

688 return me 

689 

690 

691class allencahn_periodic_semiimplicit(allencahn_periodic_fullyimplicit): 

692 r""" 

693 This class implements the one-dimensional Allen-Cahn equation with driving force and periodic boundary conditions 

694 

695 .. math:: 

696 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

697 - 6 d_w u (1 - u) 

698 

699 for :math:`u \in [0, 1]`. For discretization of the Laplacian, centered finite differences are used. 

700 The exact solution is 

701 

702 .. math:: 

703 u(x, t) = 0.5 \left(1 + \tanh\left(\frac{r - |x| - vt}{\sqrt{2}\varepsilon}\right)\right) 

704 

705 with :math:`v = 3 \sqrt{2} \varepsilon d_w` and radius :math:`r` of the circles. For time-stepping, the problem is treated 

706 in *semi-implicit* way, i.e., the part containing the Laplacian is treated implicitly, and the rest of the right-hand 

707 side is only evaluated at each time. 

708 """ 

709 

710 dtype_f = imex_mesh 

711 

712 def __init__( 

713 self, 

714 nvars=128, 

715 dw=-0.04, 

716 eps=0.04, 

717 newton_maxiter=100, 

718 newton_tol=1e-12, 

719 interval=(-0.5, 0.5), 

720 radius=0.25, 

721 stop_at_nan=True, 

722 ): 

723 super().__init__(nvars, dw, eps, newton_maxiter, newton_tol, interval, radius, stop_at_nan) 

724 

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

726 r""" 

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

728 

729 Parameters 

730 ---------- 

731 rhs : dtype_f 

732 Right-hand side for the linear system. 

733 factor : float 

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

735 u0 : dtype_u 

736 Initial guess for the iterative solver. 

737 t : float 

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

739 

740 Returns 

741 ------- 

742 me : dtype_u 

743 The solution as mesh. 

744 """ 

745 

746 me = self.dtype_u(u0) 

747 me[:] = spsolve(sp.eye(self.nvars, format='csc') - factor * self.A, rhs) 

748 return me 

749 

750 def eval_f(self, u, t): 

751 """ 

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

753 

754 Parameters 

755 ---------- 

756 u : dtype_u 

757 Current values of the numerical solution. 

758 t : float 

759 Current time of the numerical solution is computed. 

760 

761 Returns 

762 ------- 

763 f : dtype_f 

764 The right-hand side of the problem. 

765 """ 

766 f = self.dtype_f(self.init) 

767 f.impl[:] = self.A.dot(u) 

768 f.expl[:] = ( 

769 -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u) - 6.0 * self.dw * u * (1.0 - u) + 0.0 / self.eps**2 * u 

770 ) 

771 self.work_counters['rhs']() 

772 return f 

773 

774 

775class allencahn_periodic_multiimplicit(allencahn_periodic_fullyimplicit): 

776 r""" 

777 This class implements the one-dimensional Allen-Cahn equation with driving force and periodic boundary conditions 

778 

779 .. math:: 

780 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u) 

781 - 6 d_w u (1 - u) 

782 

783 for :math:`u \in [0, 1]`. For discretization of the second order spatial derivative, centered finite differences are used. 

784 The exact solution is then given by 

785 

786 .. math:: 

787 u(x, t) = 0.5 \left(1 + \tanh\left(\frac{r - |x| - vt}{\sqrt{2}\varepsilon}\right)\right) 

788 

789 with :math:`v = 3 \sqrt{2} \varepsilon d_w` and radius :math:`r` of the circles. For time-stepping, the problem is treated 

790 in a *multi-implicit* fashion, i.e., the nonlinear system containing the part with the Laplacian is solved with a 

791 linear solver provided by a ``SciPy`` routine, and the nonlinear system including the rest of the right-hand side is solved by 

792 Newton's method. 

793 """ 

794 

795 dtype_f = comp2_mesh 

796 

797 def __init__( 

798 self, 

799 nvars=128, 

800 dw=-0.04, 

801 eps=0.04, 

802 newton_maxiter=100, 

803 newton_tol=1e-12, 

804 interval=(-0.5, 0.5), 

805 radius=0.25, 

806 stop_at_nan=True, 

807 ): 

808 super().__init__(nvars, dw, eps, newton_maxiter, newton_tol, interval, radius, stop_at_nan) 

809 

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

811 r""" 

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

813 

814 Parameters 

815 ---------- 

816 rhs : dtype_f 

817 Right-hand side for the linear system. 

818 factor : float 

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

820 u0 : dtype_u 

821 Initial guess for the iterative solver. 

822 t : float 

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

824 

825 Returns 

826 ------- 

827 me : dtype_u 

828 The solution as mesh. 

829 """ 

830 

831 me = self.dtype_u(u0) 

832 me[:] = spsolve(sp.eye(self.nvars, format='csc') - factor * self.A, rhs) 

833 return me 

834 

835 def eval_f(self, u, t): 

836 """ 

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

838 

839 Parameters 

840 ---------- 

841 u : dtype_u 

842 Current values of the numerical solution. 

843 t : float 

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

845 

846 Returns 

847 ------- 

848 f : dtype_f 

849 The right-hand side of the problem. 

850 """ 

851 f = self.dtype_f(self.init) 

852 f.comp1[:] = self.A.dot(u) 

853 f.comp2[:] = ( 

854 -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u) - 6.0 * self.dw * u * (1.0 - u) + 0.0 / self.eps**2 * u