# Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD.py: 96%

## 201 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

2import scipy.sparse as sp

3from scipy.sparse.linalg import cg

5from pySDC.core.errors import ParameterError, 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

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

14# noinspection PyUnusedLocal

15class allencahn_fullyimplicit(Problem):

16 r"""

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

19 .. math::

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

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

24 .. math::

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

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

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

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.

48 order : int, optional

49 Order of the finite difference matrix.

51 Attributes

52 ----------

53 A : scipy.spdiags

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

55 dx : float

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

57 xvalues : np.1darray

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

59 newton_itercount : int

60 Number of iterations of Newton solver.

61 lin_itercount

62 Number of iterations of linear solver.

63 newton_ncalls : int

64 Number of calls of Newton solver.

65 lin_ncalls : int

66 Number of calls of linear solver.

67 """

69 dtype_u = mesh

70 dtype_f = mesh

72 def __init__(

73 self,

74 nvars=(128, 128),

75 nu=2,

76 eps=0.04,

77 newton_maxiter=200,

78 newton_tol=1e-12,

79 lin_tol=1e-8,

80 lin_maxiter=100,

81 inexact_linear_ratio=None,

83 order=2,

84 ):

85 """Initialization routine"""

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

87 if len(nvars) != 2:

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

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

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

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

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

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

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

96 self._makeAttributeAndRegister(

97 'nvars',

98 'nu',

99 'eps',

101 'order',

102 localVars=locals(),

104 )

105 self._makeAttributeAndRegister(

106 'newton_maxiter',

107 'newton_tol',

108 'lin_tol',

109 'lin_maxiter',

110 'inexact_linear_ratio',

111 localVars=locals(),

113 )

115 # compute dx and get discretization matrix A

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

117 self.A, _ = problem_helper.get_finite_difference_matrix(

118 derivative=2,

119 order=self.order,

120 stencil_type='center',

121 dx=self.dx,

122 size=self.nvars[0],

123 dim=2,

124 bc='periodic',

125 )

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

128 self.newton_itercount = 0

129 self.lin_itercount = 0

130 self.newton_ncalls = 0

131 self.lin_ncalls = 0

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

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

135 self.work_counters['linear'] = WorkCounter()

137 # noinspection PyTypeChecker

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

139 """

140 Simple Newton solver.

142 Parameters

143 ----------

144 rhs : dtype_f

145 Right-hand side for the nonlinear system

146 factor : float

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

148 u0 : dtype_u

149 Initial guess for the iterative solver.

150 t : float

151 Current time (required here for the BC).

153 Returns

154 -------

155 me : dtype_u

156 The solution as mesh.

157 """

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

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

161 nu = self.nu

162 eps2 = self.eps**2

164 Id = sp.eye(self.nvars[0] * self.nvars[1])

166 # start newton iteration

167 n = 0

168 res = 99

169 while n < self.newton_maxiter:

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

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

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

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

176 # do inexactness in the linear solver

177 if self.inexact_linear_ratio:

178 self.lin_tol = res * self.inexact_linear_ratio

180 if res < self.newton_tol:

181 break

183 # assemble dg

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

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

187 # u -= spsolve(dg, g)

188 u -= cg(

189 dg, g, x0=z, rtol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']

190 )[0]

191 # increase iteration count

192 n += 1

193 # print(n, res)

195 self.work_counters['newton']()

197 # if n == self.newton_maxiter:

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

200 me = self.dtype_u(self.init)

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

203 self.newton_ncalls += 1

204 self.newton_itercount += n

206 return me

208 def eval_f(self, u, t):

209 """

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

212 Parameters

213 ----------

214 u : dtype_u

215 Current values of the numerical solution.

216 t : float

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

219 Returns

220 -------

221 f : dtype_f

222 The right-hand side of the problem.

223 """

224 f = self.dtype_f(self.init)

225 v = u.flatten()

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

228 self.work_counters['rhs']()

229 return f

231 def u_exact(self, t, u_init=None, t_init=None):

232 r"""

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

235 Parameters

236 ----------

237 t : float

238 Time of the exact solution.

240 Returns

241 -------

242 me : dtype_u

243 The exact solution.

244 """

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

246 if t > 0:

248 def eval_rhs(t, u):

249 return self.eval_f(u.reshape(self.init[0]), t).flatten()

251 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)

253 else:

254 X, Y = np.meshgrid(self.xvalues, self.xvalues)

255 r2 = X**2 + Y**2

256 me[:] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))

258 return me

261# noinspection PyUnusedLocal

262class allencahn_semiimplicit(allencahn_fullyimplicit):

263 r"""

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

266 .. math::

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

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

271 .. math::

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

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

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

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

277 """

279 dtype_f = imex_mesh

281 def eval_f(self, u, t):

282 """

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

285 Parameters

286 ----------

287 u : dtype_u

288 Current values of the numerical solution.

289 t : float

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

292 Returns

293 -------

294 f : dtype_f

295 The right-hand side of the problem.

296 """

297 f = self.dtype_f(self.init)

298 v = u.flatten()

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

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

302 self.work_counters['rhs']()

303 return f

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

306 r"""

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

309 Parameters

310 ----------

311 rhs : dtype_f

312 Right-hand side for the linear system.

313 factor : float

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

315 u0 : dtype_u

316 Initial guess for the iterative solver.

317 t : float

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

320 Returns

321 -------

322 me : dtype_u

323 The solution as mesh.

324 """

326 class context:

327 num_iter = 0

329 def callback(xk):

330 context.num_iter += 1

331 self.work_counters['linear']()

332 return context.num_iter

334 me = self.dtype_u(self.init)

336 Id = sp.eye(self.nvars[0] * self.nvars[1])

338 me[:] = cg(

339 Id - factor * self.A,

340 rhs.flatten(),

341 x0=u0.flatten(),

342 rtol=self.lin_tol,

343 maxiter=self.lin_maxiter,

344 atol=0,

345 callback=callback,

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

348 self.lin_ncalls += 1

349 self.lin_itercount += context.num_iter

351 return me

353 def u_exact(self, t, u_init=None, t_init=None):

354 """

355 Routine to compute the exact solution at time t.

357 Parameters

358 ----------

359 t : float

360 Time of the exact solution.

362 Returns

363 -------

364 me : dtype_u

365 The exact solution.

366 """

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

368 if t > 0:

370 def eval_rhs(t, u):

371 f = self.eval_f(u.reshape(self.init[0]), t)

372 return (f.impl + f.expl).flatten()

374 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)

375 else:

376 me[:] = super().u_exact(t, u_init, t_init)

377 return me

380# noinspection PyUnusedLocal

381class allencahn_semiimplicit_v2(allencahn_fullyimplicit):

382 r"""

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

385 .. math::

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

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

390 .. math::

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

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

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

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

396 is only evaluated at each time.

397 """

399 dtype_f = imex_mesh

401 def eval_f(self, u, t):

402 """

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

405 Parameters

406 ----------

407 u : dtype_u

408 Current values of the numerical solution.

409 t : float

410 Current time of the numerical solution is computed.

412 Returns

413 -------

414 f : dtype_f

415 The right-hand side of the problem.

416 """

417 f = self.dtype_f(self.init)

418 v = u.flatten()

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

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

422 return f

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

425 """

426 Simple Newton solver.

428 Parameters

429 ----------

430 rhs : dtype_f

431 Right-hand side for the nonlinear system.

432 factor : float

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

434 u0 : dtype_u

435 Initial guess for the iterative solver.

436 t : float

437 Current time (required here for the BC).

439 Returns

440 -------

441 me : dtype_u

442 The solution as mesh.

443 """

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

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

447 nu = self.nu

448 eps2 = self.eps**2

450 Id = sp.eye(self.nvars[0] * self.nvars[1])

452 # start newton iteration

453 n = 0

454 res = 99

455 while n < self.newton_maxiter:

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

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

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

460 # res = np.linalg.norm(g, np.inf)

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

463 if res < self.newton_tol:

464 break

466 # assemble dg

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

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

470 # u -= spsolve(dg, g)

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

472 # increase iteration count

473 n += 1

474 # print(n, res)

476 # if n == self.newton_maxiter:

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

479 me = self.dtype_u(self.init)

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

482 self.newton_ncalls += 1

483 self.newton_itercount += n

485 return me

488# noinspection PyUnusedLocal

489class allencahn_multiimplicit(allencahn_fullyimplicit):

490 r"""

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

493 .. math::

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

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

498 .. math::

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

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

502 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients

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

504 """

506 dtype_f = comp2_mesh

508 def eval_f(self, u, t):

509 """

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

512 Parameters

513 ----------

514 u : dtype_u

515 Current values of the numerical solution.

516 t : float

517 Current time of the numerical solution is computed.

519 Returns

520 -------

521 f : dtype_f

522 The right-hand side of the problem.

523 """

524 f = self.dtype_f(self.init)

525 v = u.flatten()

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

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

529 return f

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

532 r"""

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

535 Parameters

536 ----------

537 rhs : dtype_f

538 Right-hand side for the linear system.

539 factor : float

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

541 u0 : dtype_u

542 Initial guess for the iterative solver.

543 t : float

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

546 Returns

547 -------

548 me : dtype_u

549 The solution as mesh.

550 """

552 class context:

553 num_iter = 0

555 def callback(xk):

556 context.num_iter += 1

557 return context.num_iter

559 me = self.dtype_u(self.init)

561 Id = sp.eye(self.nvars[0] * self.nvars[1])

563 me[:] = cg(

564 Id - factor * self.A,

565 rhs.flatten(),

566 x0=u0.flatten(),

567 rtol=self.lin_tol,

568 maxiter=self.lin_maxiter,

569 atol=0,

570 callback=callback,

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

573 self.lin_ncalls += 1

574 self.lin_itercount += context.num_iter

576 return me

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

579 """

580 Simple Newton solver.

582 Parameters

583 ----------

584 rhs : dtype_f

585 Right-hand side for the nonlinear system.

586 factor : float

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

588 u0 : dtype_u

589 Initial guess for the iterative solver.

590 t : float

591 Current time (required here for the BC).

593 Returns

594 -------

595 me : dtype_u

596 The solution as mesh.

597 """

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

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

601 nu = self.nu

602 eps2 = self.eps**2

604 Id = sp.eye(self.nvars[0] * self.nvars[1])

606 # start newton iteration

607 n = 0

608 res = 99

609 while n < self.newton_maxiter:

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

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

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

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

616 if res < self.newton_tol:

617 break

619 # assemble dg

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

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

623 # u -= spsolve(dg, g)

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

625 # increase iteration count

626 n += 1

627 # print(n, res)

629 # if n == self.newton_maxiter:

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

632 me = self.dtype_u(self.init)

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

635 self.newton_ncalls += 1

636 self.newton_itercount += n

638 return me

641# noinspection PyUnusedLocal

642class allencahn_multiimplicit_v2(allencahn_fullyimplicit):

643 r"""

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

646 .. math::

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

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

651 .. math::

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

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

655 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}

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

657 is solved by a linear solver provided by a SciPy routine.

658 """

660 dtype_f = comp2_mesh

662 def eval_f(self, u, t):

663 """

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

666 Parameters

667 ----------

668 u : dtype_u

669 Current values of the numerical solution.

670 t : float

671 Current time of the numerical solution is computed.

673 Returns

674 -------

675 f : dtype_f

676 The right-hand side of the problem.

677 """

678 f = self.dtype_f(self.init)

679 v = u.flatten()

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

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

683 return f

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

686 """

687 Simple Newton solver.

689 Parameters

690 ----------

691 rhs : dtype_f

692 Right-hand side for the nonlinear system.

693 factor : float

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

695 u0 : dtype_u

696 Initial guess for the iterative solver.

697 t : float

698 Current time (required here for the BC).

700 Returns

701 ------

702 me : dtype_u

703 The solution as mesh.

704 """

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

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

708 nu = self.nu

709 eps2 = self.eps**2

711 Id = sp.eye(self.nvars[0] * self.nvars[1])

713 # start newton iteration

714 n = 0

715 res = 99

716 while n < self.newton_maxiter:

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

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

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

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

723 if res < self.newton_tol:

724 break

726 # assemble dg

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

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

730 # u -= spsolve(dg, g)

731 u -= cg(

732 dg,

733 g,

734 x0=z,

735 rtol=self.lin_tol,

736 atol=0,

737 )[0]

738 # increase iteration count

739 n += 1

740 # print(n, res)

742 # if n == self.newton_maxiter:

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

745 me = self.dtype_u(self.init)

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

748 self.newton_ncalls += 1

749 self.newton_itercount += n

751 return me

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

754 r"""

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

757 Parameters

758 ----------

759 rhs : dtype_f

760 Right-hand side for the linear system.

761 factor : float

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

763 u0 : dtype_u

764 Initial guess for the iterative solver.

765 t : float

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

768 Returns

769 -------

770 me : dtype_u

771 The solution as mesh.

772 """

774 me = self.dtype_u(self.init)

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

777 return me