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 1D Allen-Cahn front with driving force, Dirichlet BCs and finite differences, fully implicit with Newton. 

14 

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

16 boundary conditions 

17 

18 .. math:: 

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

20 - 6 d_w u (1 - u) 

21 

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

23 exact solution is given by 

24 

25 .. math:: 

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

27 

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

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

30 

31 Parameters 

32 ---------- 

33 nvars : int 

34 Number of unknowns in the problem. 

35 dw : float 

36 Driving force. 

37 eps : float 

38 Scaling parameter :math:`\varepsilon`. 

39 newton_maxiter : int 

40 Maximum number of iterations for Newton's method. 

41 newton_tol : float 

42 Tolerance for Newton's method to terminate. 

43 interval : list 

44 Interval of spatial domain. 

45 stop_at_nan : bool, optional 

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

47 stop_at_maxiter : bool, optional 

48 Indicates that the Newton solver should raise an error instead of only warning if it has not converged after 

49 ``newton_maxiter`` iterations. 

50 

51 Attributes 

52 ---------- 

53 A : scipy.diags 

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

55 dx : float 

56 Distance between two spatial nodes. 

57 xvalues : np.1darray 

58 Spatial grid values. 

59 uext : dtype_u 

60 Contains additionally the external values of the boundary. 

61 work_counters : WorkCounter 

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

63 of right-hand side are counted. 

64 """ 

65 

66 dtype_u = mesh 

67 dtype_f = mesh 

68 

69 def __init__( 

70 self, 

71 nvars=127, 

72 dw=-0.04, 

73 eps=0.04, 

74 newton_maxiter=100, 

75 newton_tol=1e-12, 

76 interval=(-0.5, 0.5), 

77 stop_at_nan=True, 

78 stop_at_maxiter=False, 

79 ): 

80 """Initialization routine""" 

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

82 if (nvars + 1) % 2: 

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

84 

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

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

87 self._makeAttributeAndRegister( 

88 'nvars', 

89 'dw', 

90 'eps', 

91 'newton_maxiter', 

92 'newton_tol', 

93 'interval', 

94 'stop_at_nan', 

95 'stop_at_maxiter', 

96 localVars=locals(), 

97 readOnly=True, 

98 ) 

99 

100 # compute dx and get discretization matrix A 

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

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

103 

104 self.A, _ = problem_helper.get_finite_difference_matrix( 

105 derivative=2, 

106 order=2, 

107 stencil_type='center', 

108 dx=self.dx, 

109 size=self.nvars + 2, 

110 dim=1, 

111 bc='dirichlet-zero', 

112 ) 

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

114 

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

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

117 

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

119 """ 

120 Simple Newton solver. 

121 

122 Parameters 

123 ---------- 

124 rhs : dtype_f 

125 Right-hand side for the nonlinear system. 

126 factor : float 

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

128 u0 : dtype_u 

129 Initial guess for the iterative solver. 

130 t : float 

131 Current time (required here for the BC). 

132 

133 Returns 

134 ------- 

135 me : dtype_u 

136 The solution as mesh. 

137 """ 

138 

139 u = self.dtype_u(u0) 

140 eps2 = self.eps**2 

141 dw = self.dw 

142 

143 Id = sp.eye(self.nvars) 

144 

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

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

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

148 

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

150 # start newton iteration 

151 n = 0 

152 res = 99 

153 while n < self.newton_maxiter: 

154 # print(n) 

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

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

157 g = ( 

158 u 

159 - rhs 

160 - factor 

161 * ( 

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

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

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

165 ) 

166 ) 

167 

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

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

170 

171 if res < self.newton_tol: 

172 break 

173 

174 # assemble dg 

175 dg = Id - factor * ( 

176 A 

177 - 2.0 

178 / eps2 

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

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

181 ) 

182 

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

184 u -= spsolve(dg, g) 

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

186 # increase iteration count 

187 n += 1 

188 self.work_counters['newton']() 

189 

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

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

192 elif np.isnan(res): 

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

194 

195 if n == self.newton_maxiter: 

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

197 if self.stop_at_maxiter: 

198 raise ProblemError(msg) 

199 else: 

200 self.logger.warning(msg) 

201 

202 me = self.dtype_u(self.init) 

203 me[:] = u[:] 

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. 

217 

218 Returns 

219 ------- 

220 f : dtype_f 

221 The right-hand side of the problem. 

222 """ 

223 # set up boundary values to embed inner points 

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

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

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

227 

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

229 

230 f = self.dtype_f(self.init) 

231 f[:] = ( 

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

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

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

235 ) 

236 self.work_counters['rhs']() 

237 return f 

238 

239 def u_exact(self, t): 

240 r""" 

241 Routine to return initial condition or the exact solution 

242 

243 Parameters 

244 ---------- 

245 t : float 

246 Time at which the exact solution is computed. 

247 

248 Returns 

249 ------- 

250 me : dtype_u 

251 The exact solution (in space and time). 

252 """ 

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

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

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

256 return me 

257 

258 

259class allencahn_front_semiimplicit(allencahn_front_fullyimplicit): 

260 r""" 

261 1D Allen-Cahn front with driving force, Dirichlet BCs and finite differences, IMEX with Laplacian implicit. 

262 

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

264 boundary conditions 

265 

266 .. math:: 

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

268 - 6 d_w u (1 - u) 

269 

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

271 The exact solution is given by 

272 

273 .. math:: 

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

275 

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

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

278 explicitly. 

279 """ 

280 

281 dtype_f = imex_mesh 

282 

283 def eval_f(self, u, t): 

284 """ 

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

286 

287 Parameters 

288 ---------- 

289 u : dtype_u 

290 Current values of the numerical solution. 

291 t : float 

292 Current time of the numerical solution is computed. 

293 

294 Returns 

295 ------- 

296 f : dtype_f 

297 The right-hand side of the problem. 

298 """ 

299 # set up boundary values to embed inner points 

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

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

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

303 

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

305 

306 f = self.dtype_f(self.init) 

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

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

309 self.work_counters['rhs']() 

310 return f 

311 

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

313 r""" 

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

315 

316 Parameters 

317 ---------- 

318 rhs : dtype_f 

319 Right-hand side for the linear system. 

320 factor : float 

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

322 u0 : dtype_u 

323 Initial guess for the iterative solver. 

324 t : float 

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

326 

327 Returns 

328 ------- 

329 me : dtype_u 

330 The solution as mesh. 

331 """ 

332 

333 me = self.dtype_u(self.init) 

334 self.uext[0] = 0.0 

335 self.uext[-1] = 0.0 

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

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

338 return me 

339 

340 

341class allencahn_front_finel(allencahn_front_fullyimplicit): 

342 r""" 

343 1D Allen-Cahn front with driving force and Dirichlet BCs, fully implicit with Newton, using Finel's trick. 

344 

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

346 boundary conditions 

347 

348 .. math:: 

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

350 - 6 d_w u (1 - u) 

351 

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

353 The exact solution is given by 

354 

355 .. math:: 

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

357 

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

359 

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

361 *Finel's trick* is used. Let 

362 

363 .. math:: 

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

365 

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

367 

368 .. math:: 

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

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

371 \right] (2u - 1) - 6 d_w u (1 - u). 

372 

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

374 """ 

375 

376 # noinspection PyTypeChecker 

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

378 """ 

379 Simple Newton solver. 

380 

381 Parameters 

382 ---------- 

383 rhs : dtype_f 

384 Right-hand side for the nonlinear system. 

385 factor : float 

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

387 u0 : dtype_u 

388 Initial guess for the iterative solver. 

389 t : float 

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

391 

392 Returns 

393 ------- 

394 me : dtype_u 

395 The solution as mesh. 

396 """ 

397 

398 u = self.dtype_u(u0) 

399 dw = self.dw 

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

401 

402 Id = sp.eye(self.nvars) 

403 

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

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

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

407 

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

409 # start newton iteration 

410 n = 0 

411 res = 99 

412 while n < self.newton_maxiter: 

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

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

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

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

417 

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

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

420 

421 if res < self.newton_tol: 

422 break 

423 

424 # assemble dg 

425 dgprim = ( 

426 1.0 

427 / self.dx**2 

428 * ( 

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

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

431 ) 

432 ) 

433 

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

435 

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

437 u -= spsolve(dg, g) 

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

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

440 # increase iteration count 

441 n += 1 

442 self.work_counters['newton']() 

443 

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

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

446 elif np.isnan(res): 

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

448 

449 if n == self.newton_maxiter: 

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

451 

452 me = self.dtype_u(self.init) 

453 me[:] = u[:] 

454 

455 return me 

456 

457 def eval_f(self, u, t): 

458 """ 

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

460 

461 Parameters 

462 ---------- 

463 u : dtype_u 

464 Current values of the numerical solution. 

465 t : float 

466 Current time of the numerical solution is computed. 

467 

468 Returns 

469 ------- 

470 f : dtype_f 

471 The right-hand side of the problem. 

472 """ 

473 # set up boundary values to embed inner points 

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

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

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

477 

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

479 

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

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

482 f = self.dtype_f(self.init) 

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

484 self.work_counters['rhs']() 

485 return f 

486 

487 

488class allencahn_periodic_fullyimplicit(Problem): 

489 r""" 

490 1D Allen-Cahn equation with driving force, periodic BCs and finite differences, fully implicit with Newton. 

491 

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

493 

494 .. math:: 

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

496 - 6 d_w u (1 - u) 

497 

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

499 The exact solution is 

500 

501 .. math:: 

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

503 

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

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

506 

507 Parameters 

508 ---------- 

509 nvars : int 

510 Number of unknowns in the problem. 

511 dw : float 

512 Driving force. 

513 eps : float 

514 Scaling parameter :math:`\varepsilon`. 

515 newton_maxiter : int 

516 Maximum number of iterations for Newton's method. 

517 newton_tol : float 

518 Tolerance for Newton's method to terminate. 

519 interval : list 

520 Interval of spatial domain. 

521 radius : float 

522 Radius of the circles. 

523 stop_at_nan : bool, optional 

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

525 

526 Attributes 

527 ---------- 

528 A : scipy.diags 

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

530 dx : float 

531 Distance between two spatial nodes. 

532 xvalues : np.1darray 

533 Spatial grid points. 

534 work_counters : WorkCounter 

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

536 of right-hand side are counted. 

537 """ 

538 

539 dtype_u = mesh 

540 dtype_f = mesh 

541 

542 def __init__( 

543 self, 

544 nvars=128, 

545 dw=-0.04, 

546 eps=0.04, 

547 newton_maxiter=100, 

548 newton_tol=1e-12, 

549 interval=(-0.5, 0.5), 

550 radius=0.25, 

551 stop_at_nan=True, 

552 ): 

553 """Initialization routine""" 

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

555 if (nvars) % 2: 

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

557 

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

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

560 self._makeAttributeAndRegister( 

561 'nvars', 

562 'dw', 

563 'eps', 

564 'newton_maxiter', 

565 'newton_tol', 

566 'interval', 

567 'radius', 

568 'stop_at_nan', 

569 localVars=locals(), 

570 readOnly=True, 

571 ) 

572 

573 # compute dx and get discretization matrix A 

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

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

576 

577 self.A, _ = problem_helper.get_finite_difference_matrix( 

578 derivative=2, 

579 order=2, 

580 stencil_type='center', 

581 dx=self.dx, 

582 size=self.nvars, 

583 dim=1, 

584 bc='periodic', 

585 ) 

586 

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

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

589 

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

591 """ 

592 Simple Newton solver. 

593 

594 Parameters 

595 ---------- 

596 rhs : dtype_f 

597 Right-hand side for the nonlinear system. 

598 factor : float 

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

600 u0 : dtype_u 

601 Initial guess for the iterative solver. 

602 t : float 

603 Current time (required here for the BC). 

604 

605 Returns 

606 ------- 

607 u : dtype_u 

608 The solution as mesh. 

609 """ 

610 

611 u = self.dtype_u(u0) 

612 eps2 = self.eps**2 

613 dw = self.dw 

614 

615 Id = sp.eye(self.nvars) 

616 

617 # start newton iteration 

618 n = 0 

619 res = 99 

620 while n < self.newton_maxiter: 

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

622 g = ( 

623 u 

624 - rhs 

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

626 ) 

627 

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

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

630 

631 if res < self.newton_tol: 

632 break 

633 

634 # assemble dg 

635 dg = Id - factor * ( 

636 self.A 

637 - 2.0 

638 / eps2 

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

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

641 ) 

642 

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

644 u -= spsolve(dg, g) 

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

646 # increase iteration count 

647 n += 1 

648 self.work_counters['newton']() 

649 

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

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

652 elif np.isnan(res): 

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

654 

655 if n == self.newton_maxiter: 

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

657 

658 me = self.dtype_u(self.init) 

659 me[:] = u[:] 

660 

661 return me 

662 

663 def eval_f(self, u, t): 

664 """ 

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

666 

667 Parameters 

668 ---------- 

669 u : dtype_u 

670 Current values of the numerical solution. 

671 t : float 

672 Current time of the numerical solution is computed. 

673 

674 Returns 

675 ------- 

676 f : dtype_f 

677 The right-hand side of the problem. 

678 """ 

679 f = self.dtype_f(self.init) 

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

681 self.work_counters['rhs']() 

682 return f 

683 

684 def u_exact(self, t): 

685 r""" 

686 Routine to return initial condition or the exact solution. 

687 

688 Parameters 

689 ---------- 

690 t : float 

691 Time at which the approximated exact solution is computed. 

692 

693 Returns 

694 ------- 

695 me : dtype_u 

696 The approximated exact solution. 

697 """ 

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

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

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

701 return me 

702 

703 

704class allencahn_periodic_semiimplicit(allencahn_periodic_fullyimplicit): 

705 r""" 

706 1D Allen-Cahn equation with driving force, periodic BCs and finite differences, IMEX with Laplacian implicit. 

707 

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

709 

710 .. math:: 

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

712 - 6 d_w u (1 - u) 

713 

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

715 The exact solution is 

716 

717 .. math:: 

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

719 

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

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

722 side is only evaluated at each time. 

723 

724 Parameters 

725 ---------- 

726 nvars : int, optional 

727 Number of unknowns in the problem, i.e., points of the periodic grid. Has to be even. 

728 dw : float, optional 

729 Driving force :math:`d_w`. 

730 eps : float, optional 

731 Scaling parameter :math:`\varepsilon`. 

732 newton_maxiter : int, optional 

733 Maximum number of iterations for Newton's method. Not used here, since the implicit part is linear. 

734 newton_tol : float, optional 

735 Tolerance for Newton's method to terminate. Not used here, since the implicit part is linear. 

736 interval : list, optional 

737 Interval of spatial domain. 

738 radius : float, optional 

739 Radius of the circles. 

740 stop_at_nan : bool, optional 

741 Indicates that the Newton solver should stop if ``nan`` values arise. Not used here. 

742 """ 

743 

744 dtype_f = imex_mesh 

745 

746 def __init__( 

747 self, 

748 nvars=128, 

749 dw=-0.04, 

750 eps=0.04, 

751 newton_maxiter=100, 

752 newton_tol=1e-12, 

753 interval=(-0.5, 0.5), 

754 radius=0.25, 

755 stop_at_nan=True, 

756 ): 

757 """Initialization routine""" 

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

759 

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

761 r""" 

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

763 

764 Parameters 

765 ---------- 

766 rhs : dtype_f 

767 Right-hand side for the linear system. 

768 factor : float 

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

770 u0 : dtype_u 

771 Initial guess for the iterative solver. 

772 t : float 

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

774 

775 Returns 

776 ------- 

777 me : dtype_u 

778 The solution as mesh. 

779 """ 

780 

781 me = self.dtype_u(u0) 

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

783 return me 

784 

785 def eval_f(self, u, t): 

786 """ 

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

788 

789 Parameters 

790 ---------- 

791 u : dtype_u 

792 Current values of the numerical solution. 

793 t : float 

794 Current time of the numerical solution is computed. 

795 

796 Returns 

797 ------- 

798 f : dtype_f 

799 The right-hand side of the problem. 

800 """ 

801 f = self.dtype_f(self.init) 

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

803 f.expl[:] = ( 

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

805 ) 

806 self.work_counters['rhs']() 

807 return f 

808 

809 

810class allencahn_periodic_multiimplicit(allencahn_periodic_fullyimplicit): 

811 r""" 

812 1D periodic Allen-Cahn equation with driving force, multi-implicit: linear solve for Laplacian, Newton for rest. 

813 

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

815 

816 .. math:: 

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

818 - 6 d_w u (1 - u) 

819 

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

821 The exact solution is then given by 

822 

823 .. math:: 

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

825 

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

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

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

829 Newton's method. 

830 

831 Parameters 

832 ---------- 

833 nvars : int, optional 

834 Number of unknowns in the problem, i.e., points of the periodic grid. Has to be even. 

835 dw : float, optional 

836 Driving force :math:`d_w`. 

837 eps : float, optional 

838 Scaling parameter :math:`\varepsilon`. 

839 newton_maxiter : int, optional 

840 Maximum number of iterations for Newton's method in ``solve_system_2``. 

841 newton_tol : float, optional 

842 Absolute tolerance for Newton's method to terminate, applied to the maximum norm of the residual. 

843 interval : list, optional 

844 Interval of spatial domain. 

845 radius : float, optional 

846 Radius of the circles. 

847 stop_at_nan : bool, optional 

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

849 """ 

850 

851 dtype_f = comp2_mesh 

852 

853 def __init__( 

854 self, 

855 nvars=128, 

856 dw=-0.04, 

857 eps=0.04, 

858 newton_maxiter=100, 

859 newton_tol=1e-12, 

860 interval=(-0.5, 0.5), 

861 radius=0.25, 

862 stop_at_nan=True, 

863 ): 

864 """Initialization routine""" 

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

866 

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

868 r""" 

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