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

## 130 statements

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

1import numpy as np

2import cupy as cp

3import cupyx.scipy.sparse as csp

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

6from pySDC.core.errors import ProblemError

7from pySDC.core.problem import Problem

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

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

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

14 r"""

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

17 .. math::

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

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

22 .. math::

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

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.

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

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.

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

67 dtype_u = cupy_mesh

68 dtype_f = cupy_mesh

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,

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

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',

101 localVars=locals(),

103 )

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

110 self.newton_itercount = 0

111 self.lin_itercount = 0

112 self.newton_ncalls = 0

113 self.lin_ncalls = 0

115 @staticmethod

116 def __get_A(N, dx):

117 """

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

120 Parameters

121 ----------

122 N : list

123 Number of degrees of freedom.

124 dx : float

125 Distance between two spatial nodes.

127 Returns

128 -------

129 A : cupyx.scipy.sparse.csr_matrix

130 CuPy-matrix A in CSR format.

131 """

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

144 # noinspection PyTypeChecker

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

146 """

147 Simple Newton solver.

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

160 Returns

161 -------

162 u : dtype_u

163 The solution as mesh.

164 """

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

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

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

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

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

183 if res < self.newton_tol:

184 break

186 # assemble dg

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

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)

196 # if n == self.newton_maxiter:

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

199 me = self.dtype_u(self.init)

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

202 self.newton_ncalls += 1

203 self.newton_itercount += n

205 return me

207 def eval_f(self, u, t):

208 """

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

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

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)

227 return f

229 def u_exact(self, t):

230 r"""

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

233 Parameters

234 ----------

235 t : float

236 Time of the exact solution.

238 Returns

239 -------

240 me : dtype_u

241 The exact solution.

242 """

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

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

257 .. math::

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

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

262 .. math::

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

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.

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

270 """

272 dtype_f = imex_cupy_mesh

274 def eval_f(self, u, t):

275 """

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

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.

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)

295 return f

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

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

312 Returns

313 -------

314 me : dtype_u

315 The solution as mesh.

316 """

318 class context:

319 num_iter = 0

321 def callback(xk):

322 context.num_iter += 1

323 return context.num_iter

325 me = self.dtype_u(self.init)

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

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)

338 self.lin_ncalls += 1

339 self.lin_itercount += context.num_iter

341 return me

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

349 .. math::

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

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

354 .. math::

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

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.

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

363 """

365 dtype_f = imex_cupy_mesh

367 def eval_f(self, u, t):

368 """

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

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.

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)

388 return f

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

391 """

392 Simple Newton solver.

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

405 Returns

406 -------

407 me : dtype_u

408 The solution as mesh.

409 """

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

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

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

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

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

428 if res < self.newton_tol:

429 break

431 # assemble dg

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

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)

441 # if n == self.newton_maxiter:

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

444 me = self.dtype_u(self.init)

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

447 self.newton_ncalls += 1

448 self.newton_itercount += n

450 return me

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

458 .. math::

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

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

463 .. math::

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

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.

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

471 """

473 dtype_f = comp2_cupy_mesh

475 def eval_f(self, u, t):

476 """

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

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

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)

496 return f

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

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

513 Returns

514 -------

515 me : dtype_u

516 The solution as mesh.

517 """

519 class context:

520 num_iter = 0

522 def callback(xk):

523 context.num_iter += 1

524 return context.num_iter

526 me = self.dtype_u(self.init)

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

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)

539 self.lin_ncalls += 1

540 self.lin_itercount += context.num_iter

542 return me

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

545 """

546 Simple Newton solver.

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

559 Returns

560 -------

561 me : dtype_u

562 The solution as mesh.

563 """

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

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

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

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

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

582 if res < self.newton_tol:

583 break

585 # assemble dg

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

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)

595 # if n == self.newton_maxiter:

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

598 me = self.dtype_u(self.init)

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

601 self.newton_ncalls += 1

602 self.newton_itercount += n

604 return me

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

612 .. math::

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

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

617 .. math::

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

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.

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

626 """

628 dtype_f = comp2_cupy_mesh

630 def eval_f(self, u, t):

631 """

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

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

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)

651 return f

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

654 """

655 Simple Newton solver.

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

668 Returns

669 -------

670 me : dtype_u

671 The solution as mesh.

672 """

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

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

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

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

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

691 if res < self.newton_tol:

692 break

694 # assemble dg

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

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)

704 # if n == self.newton_maxiter:

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

707 me = self.dtype_u(self.init)

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

710 self.newton_ncalls += 1

711 self.newton_itercount += n

713 return me

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

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

730 Returns

731 -------

732 me : dtype_u

733 The solution as mesh.

734 """

736 me = self.dtype_u(self.init)

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

739 return me