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

## 225 statements

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

1import numpy as np

2import scipy.sparse as sp

3from scipy.sparse.linalg import spsolve

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

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

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)

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

23 .. math::

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

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.

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.

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

61 dtype_u = mesh

62 dtype_f = mesh

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

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

92 )

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

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)

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

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

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

113 """

114 Simple Newton solver.

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

127 Returns

128 -------

129 me : dtype_u

130 The solution as mesh.

131 """

133 u = self.dtype_u(u0)

134 eps2 = self.eps**2

135 dw = self.dw

137 Id = sp.eye(self.nvars)

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

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 )

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

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

165 if res < self.newton_tol:

166 break

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 )

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

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)

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)

196 me = self.dtype_u(self.init)

197 me[:] = u[:]

199 return me

201 def eval_f(self, u, t):

202 """

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

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.

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

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

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

233 def u_exact(self, t):

234 r"""

235 Routine to return initial condition or the exact solution

237 Parameters

238 ----------

239 t : float

240 Time at which the exact solution is computed.

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

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

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)

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

265 .. math::

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

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

273 dtype_f = imex_mesh

275 def eval_f(self, u, t):

276 """

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

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.

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

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

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

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

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

319 Returns

320 -------

321 me : dtype_u

322 The solution as mesh.

323 """

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

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

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)

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

343 The exact solution is given by

345 .. math::

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

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

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

353 .. math::

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

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

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

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

364 """

366 # noinspection PyTypeChecker

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

368 """

369 Simple Newton solver.

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

382 Returns

383 -------

384 me : dtype_u

385 The solution as mesh.

386 """

388 u = self.dtype_u(u0)

389 dw = self.dw

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

392 Id = sp.eye(self.nvars)

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

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

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

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

411 if res < self.newton_tol:

412 break

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 )

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

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

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)

439 if n == self.newton_maxiter:

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

442 me = self.dtype_u(self.init)

443 me[:] = u[:]

445 return me

447 def eval_f(self, u, t):

448 """

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

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.

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

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

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

478class allencahn_periodic_fullyimplicit(Problem):

479 r"""

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

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)

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

487 The exact solution is

489 .. math::

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

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.

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.

511 stop_at_nan : bool, optional

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

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

527 dtype_u = mesh

528 dtype_f = mesh

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

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

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

555 'stop_at_nan',

556 localVars=locals(),

558 )

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

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 )

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

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

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

578 """

579 Simple Newton solver.

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

592 Returns

593 -------

594 u : dtype_u

595 The solution as mesh.

596 """

598 u = self.dtype_u(u0)

599 eps2 = self.eps**2

600 dw = self.dw

602 Id = sp.eye(self.nvars)

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 )

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

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

618 if res < self.newton_tol:

619 break

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 )

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

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)

642 if n == self.newton_maxiter:

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

645 me = self.dtype_u(self.init)

646 me[:] = u[:]

648 return me

650 def eval_f(self, u, t):

651 """

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

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.

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

671 def u_exact(self, t):

672 r"""

673 Routine to return initial condition or the exact solution.

675 Parameters

676 ----------

677 t : float

678 Time at which the approximated exact solution is computed.

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

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

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)

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

700 The exact solution is

702 .. math::

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

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

710 dtype_f = imex_mesh

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

721 stop_at_nan=True,

722 ):

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

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

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

740 Returns

741 -------

742 me : dtype_u

743 The solution as mesh.

744 """

746 me = self.dtype_u(u0)

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

748 return me

750 def eval_f(self, u, t):

751 """

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

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.

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

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

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)

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

786 .. math::

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

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

795 dtype_f = comp2_mesh

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

806 stop_at_nan=True,

807 ):

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

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

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

825 Returns

826 -------

827 me : dtype_u

828 The solution as mesh.

829 """

831 me = self.dtype_u(u0)

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

833 return me

835 def eval_f(self, u, t):

836 """

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

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

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

855 )

856 self.work_counters['rhs']()

857 return f

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

860 r"""

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

863 Parameters

864 ----------

865 rhs : dtype_f

866 Right-hand side for the linear system.

867 factor : float

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

869 u0 : dtype_u

870 Initial guess for the iterative solver.

871 t : float

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

874 Returns

875 -------

876 me : dtype_u

877 The solution as mesh.

878 """

880 u = self.dtype_u(u0)

881 eps2 = self.eps**2

882 dw = self.dw

884 Id = sp.eye(self.nvars)

886 # start newton iteration

887 n = 0

888 res = 99

889 while n < self.newton_maxiter:

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

891 g = (

892 u

893 - rhs

894 - factor

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

896 )

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

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

901 if res < self.newton_tol:

902 break

904 # assemble dg

905 dg = Id - factor * (

906 -2.0 / eps2 * sp.diags((1.0 - u) * (1.0 - 2.0 * u) - u * ((1.0 - 2.0 * u) + 2.0 * (1.0 - u)), offsets=0)

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

908 + 0.0 / self.eps**2 * Id

909 )

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

912 u -= spsolve(dg, g)

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

914 # increase iteration count

915 n += 1

916 self.work_counters['newton']()

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

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

920 elif np.isnan(res):

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

923 if n == self.newton_maxiter:

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

926 me = self.dtype_u(self.init)

927 me[:] = u[:]

929 return me