Coverage for pySDC/implementations/problem_classes/odeSystem.py: 100%

187 statements

, created at 2024-09-20 17:10 +0000

1#!/usr/bin/env python3

2# -*- coding: utf-8 -*-

3"""

4Implementation of systems test problem ODEs.

7Reference :

9Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods

10on parallel computers. SIAM journal on scientific and statistical computing,

1112(5), 1000-1028.

12"""

13import numpy as np

15from pySDC.core.errors import ProblemError

16from pySDC.core.problem import Problem, WorkCounter

17from pySDC.implementations.datatype_classes.mesh import mesh

20class ProtheroRobinsonAutonomous(Problem):

21 r"""

22 Implement the Prothero-Robinson problem into autonomous form:

24 .. math::

25 \begin{eqnarray*}

26 \frac{du}{dt} &=& -\frac{u^3-g(v)^3}{\epsilon} + \frac{dg}{dv}, &\quad u(0) = g(0),\\

27 \frac{dv}{dt} &=& 1, &\quad v(0) = 0,

28 \end{eqnarray*}

30 with :math:\epsilon a stiffness parameter, that makes the problem more stiff

31 the smaller it is (usual taken value is :math:\epsilon=1e^{-3}).

32 Exact solution is given by :math:u(t)=g(t),\;v(t)=t, and this implementation uses

33 :math:g(t)=\cos(t).

35 Implement also the non-linear form of this problem:

37 .. math::

38 \frac{du}{dt} = -\frac{u^3-g(v)^3}{\epsilon} + \frac{dg}{dv}, \quad u(0) = g(0).

40 To use an other exact solution, one just have to derivate this class

41 and overload the g, dg and dg2 methods. For instance,

42 to use :math:g(t)=e^{-0.2t}, define and use the following class:

44 >>> class MyProtheroRobinson(ProtheroRobinsonAutonomous):

45 >>>

46 >>> def g(self, t):

47 >>> return np.exp(-0.2 * t)

48 >>>

49 >>> def dg(self, t):

50 >>> return (-0.2) * np.exp(-0.2 * t)

51 >>>

52 >>> def dg2(self, t):

53 >>> return (-0.2) ** 2 * np.exp(-0.2 * t)

55 Parameters

56 ----------

57 epsilon : float, optional

58 Stiffness parameter. The default is 1e-3.

59 nonLinear : bool, optional

60 Wether or not to use the non-linear form of the problem. The default is False.

61 newton_maxiter : int, optional

62 Maximum number of Newton iteration in solve_system. The default is 200.

63 newton_tol : float, optional

64 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

65 stop_at_nan : bool, optional

66 Wheter to stop or not solve_system when getting NAN. The default is True.

68 Reference

69 ---------

70 A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving

71 stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974),

72 pp. 145–162.

73 """

75 dtype_u = mesh

76 dtype_f = mesh

78 def __init__(self, epsilon=1e-3, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

79 nvars = 2

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

82 self.f = self.f_NONLIN if nonLinear else self.f_LIN

83 self.dgInv = self.dgInv_NONLIN if nonLinear else self.dgInv_LIN

84 self._makeAttributeAndRegister(

85 'epsilon', 'nonLinear', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

86 )

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

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

90 # -------------------------------------------------------------------------

91 # g function (analytical solution), and its first and second derivative

92 # -------------------------------------------------------------------------

93 def g(self, t):

94 return np.cos(t)

96 def dg(self, t):

97 return -np.sin(t)

99 def dg2(self, t):

100 return -np.cos(t)

102 # -------------------------------------------------------------------------

103 # f(u,t) and Jacobian functions

104 # -------------------------------------------------------------------------

105 def f(self, u, t):

106 raise NotImplementedError()

108 def f_LIN(self, u, t):

109 return -self.epsilon ** (-1) * (u - self.g(t)) + self.dg(t)

111 def f_NONLIN(self, u, t):

112 return -self.epsilon ** (-1) * (u**3 - self.g(t) ** 3) + self.dg(t)

114 def dgInv(self, u, t):

115 raise NotImplementedError()

117 def dgInv_LIN(self, u, t, dt):

118 e = self.epsilon

119 g1, g2 = self.dg(t), self.dg2(t)

120 return np.array([[1 / (dt / e + 1), (dt * g2 + dt * g1 / e) / (dt / e + 1)], [0, 1]])

122 def dgInv_NONLIN(self, u, t, dt):

123 e = self.epsilon

124 g, g1, g2 = self.g(t), self.dg(t), self.dg2(t)

125 return np.array(

126 [[1 / (3 * dt * u**2 / e + 1), (dt * g2 + 3 * dt * g**2 * g1 / e) / (3 * dt * u**2 / e + 1)], [0, 1]]

127 )

129 # -------------------------------------------------------------------------

130 # pySDC required methods

131 # -------------------------------------------------------------------------

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

133 r"""

134 Routine to return initial conditions or exact solutions.

136 Parameters

137 ----------

138 t : float

139 Time at which the exact solution is computed.

140 u_init : dtype_u

141 Initial conditions for getting the exact solution.

142 t_init : float

143 The starting time.

145 Returns

146 -------

147 u : dtype_u

148 The exact solution.

149 """

150 u = self.dtype_u(self.init)

151 u[0] = self.g(t)

152 u[1] = t

153 return u

155 def eval_f(self, u, t):

156 """

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

159 Parameters

160 ----------

161 u : dtype_u

162 Current values of the numerical solution.

163 t : float

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

166 Returns

167 -------

168 f : dtype_f

169 The right-hand side of the problem (one component).

170 """

172 f = self.dtype_f(self.init)

173 u, t = u

174 f[0] = self.f(u, t)

175 f[1] = 1

176 self.work_counters['rhs']()

177 return f

179 def solve_system(self, rhs, dt, u0, t):

180 """

181 Simple Newton solver for the nonlinear equation

183 Parameters

184 ----------

185 rhs : dtype_f

186 Right-hand side for the nonlinear system.

187 dt : float

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

189 u0 : dtype_u

190 Initial guess for the iterative solver.

191 t : float

192 Time of the updated solution (e.g. for time-dependent BCs).

194 Returns

195 -------

196 u : dtype_u

197 The solution as mesh.

198 """

199 # create new mesh object from u0 and set initial values for iteration

200 u = self.dtype_u(u0)

202 # start newton iteration

203 n, res = 0, np.inf

204 while n < self.newton_maxiter:

205 # evaluate RHS

206 f = self.dtype_u(u)

207 f[0] = self.f(*u)

208 f[1] = 1

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

211 g = u - dt * f - rhs

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

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

215 if res < self.newton_tol or np.isnan(res):

216 break

218 # assemble (dg/du)^{-1}

219 dgInv = self.dgInv(u[0], u[1], dt)

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

221 u -= dgInv @ g

223 # increase iteration count and work counter

224 n += 1

225 self.work_counters['newton']()

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

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

229 elif np.isnan(res): # pragma: no cover

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

232 if n == self.newton_maxiter:

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

235 return u

238class Kaps(Problem):

239 r"""

240 Implement the Kaps problem:

242 .. math::

243 \begin{eqnarray*}

244 \frac{du}{dt} &=& -(2+\epsilon^{-1})u + \frac{v^2}{\epsilon}, &\quad u(0) = 1,\\

245 \frac{dv}{dt} &=& u - v(1+v), &\quad v(0) = 1,

246 \end{eqnarray*}

248 with :math:\epsilon a stiffness parameter, that makes the problem more stiff

249 the smaller it is (usual taken value is :math:\epsilon=1e^{-3}).

250 Exact solution is given by :math:u(t)=e^{-2t},\;v(t)=e^{-t}.

252 Parameters

253 ----------

254 epsilon : float, optional

255 Stiffness parameter. The default is 1e-3.

256 newton_maxiter : int, optional

257 Maximum number of Newton iteration in solve_system. The default is 200.

258 newton_tol : float, optional

259 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

260 stop_at_nan : bool, optional

261 Wheter to stop or not solve_system when getting NAN. The default is True.

263 Reference

264 ---------

265 Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods

266 on parallel computers. SIAM journal on scientific and statistical computing,

267 12(5), 1000-1028.

268 """

270 dtype_u = mesh

271 dtype_f = mesh

273 def __init__(self, epsilon=1e-3, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

274 nvars = 2

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

277 self._makeAttributeAndRegister(

278 'epsilon', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

279 )

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

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

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

284 r"""

285 Routine to return initial conditions or exact solutions.

287 Parameters

288 ----------

289 t : float

290 Time at which the exact solution is computed.

291 u_init : dtype_u

292 Initial conditions for getting the exact solution.

293 t_init : float

294 The starting time.

296 Returns

297 -------

298 u : dtype_u

299 The exact solution.

300 """

301 u = self.dtype_u(self.init)

302 u[:] = [np.exp(-2 * t), np.exp(-t)]

303 return u

305 def eval_f(self, u, t):

306 """

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

309 Parameters

310 ----------

311 u : dtype_u

312 Current values of the numerical solution.

313 t : float

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

316 Returns

317 -------

318 f : dtype_f

319 The right-hand side of the problem (one component).

320 """

321 f = self.dtype_f(self.init)

322 eps = self.epsilon

323 x, y = u

325 f[:] = [-(2 + 1 / eps) * x + y**2 / eps, x - y * (1 + y)]

326 self.work_counters['rhs']()

327 return f

329 def solve_system(self, rhs, dt, u0, t):

330 """

331 Simple Newton solver for the nonlinear equation

333 Parameters

334 ----------

335 rhs : dtype_f

336 Right-hand side for the nonlinear system.

337 dt : float

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

339 u0 : dtype_u

340 Initial guess for the iterative solver.

341 t : float

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

344 Returns

345 -------

346 u : dtype_u

347 The solution as mesh.

348 """

349 # create new mesh object from u0 and set initial values for iteration

350 u = self.dtype_u(u0)

351 eps = self.epsilon

353 # start newton iteration

354 n, res = 0, np.inf

355 while n < self.newton_maxiter:

356 x, y = u

357 f = np.array([-(2 + 1 / eps) * x + y**2 / eps, x - y * (1 + y)])

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

360 g = u - dt * f - rhs

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

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

364 if res < self.newton_tol or np.isnan(res):

365 break

367 # assemble (dg/du)^(-1)

368 prefactor = 4 * dt**2 * eps * y + 2 * dt**2 * eps + dt**2 + 2 * dt * eps * y + 3 * dt * eps + dt + eps

369 dgInv = (

370 1

371 / prefactor

372 * np.array([[2 * dt * eps * y + dt * eps + eps, 2 * dt * y], [dt * eps, 2 * dt * eps + dt + eps]])

373 )

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

376 u -= dgInv @ g

378 # increase iteration count and work counter

379 n += 1

380 self.work_counters['newton']()

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

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

384 elif np.isnan(res): # pragma: no cover

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

387 if n == self.newton_maxiter:

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

390 return u

393class ChemicalReaction3Var(Problem):

394 r"""

395 Chemical reaction with three components, modeled by the non-linear system:

397 .. math::

398 \frac{d{\bf u}}{dt} =

399 \begin{pmatrix}

400 0.013+1000u_3 & 0 & 0 \\

401 0 & 2500u_3 0 \\

402 0.013 & 0 & 1000u_1 + 2500u_2

403 \end{pmatrix}

404 {\bf u},

406 with initial solution :math:u(0)=(0.990731920827, 1.009264413846, -0.366532612659e-5).

408 Parameters

409 ----------

410 newton_maxiter : int, optional

411 Maximum number of Newton iteration in solve_system. The default is 200.

412 newton_tol : float, optional

413 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

414 stop_at_nan : bool, optional

415 Wheter to stop or not solve_system when getting NAN. The default is True.

417 Reference

418 ---------

419 Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods

420 on parallel computers. SIAM journal on scientific and statistical computing,

421 12(5), 1000-1028.

422 """

424 dtype_u = mesh

425 dtype_f = mesh

427 def __init__(self, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

428 nvars = 3

429 u0 = (0.990731920827, 1.009264413846, -0.366532612659e-5)

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

432 self._makeAttributeAndRegister(

433 'u0', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

434 )

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

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

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

439 r"""

440 Routine to return initial conditions or to approximate exact solution using SciPy.

442 Parameters

443 ----------

444 t : float

445 Time at which the approximated exact solution is computed.

446 u_init : pySDC.implementations.problem_classes.Lorenz.dtype_u

447 Initial conditions for getting the exact solution.

448 t_init : float

449 The starting time.

451 Returns

452 -------

453 me : dtype_u

454 The approximated exact solution.

455 """

457 me = self.dtype_u(self.init)

459 if t > 0:

461 def eval_rhs(t, u):

462 r"""

463 Evaluate the right hand side, but switch the arguments for SciPy.

465 Args:

466 t (float): Time

467 u (numpy.ndarray): Solution at time t

469 Returns:

470 (numpy.ndarray): Right hand side evaluation

471 """

472 return self.eval_f(u, t)

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

475 else:

476 me[:] = self.u0

477 return me

479 def eval_f(self, u, t):

480 """

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

483 Parameters

484 ----------

485 u : dtype_u

486 Current values of the numerical solution.

487 t : float

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

490 Returns

491 -------

492 f : dtype_f

493 The right-hand side of the problem (one component).

494 """

495 f = self.dtype_f(self.init)

496 c1, c2, c3 = u

498 f[:] = -np.array([0.013 * c1 + 1000 * c3 * c1, 2500 * c3 * c2, 0.013 * c1 + 1000 * c1 * c3 + 2500 * c2 * c3])

499 self.work_counters['rhs']()

500 return f

502 def solve_system(self, rhs, dt, u0, t):

503 """

504 Simple Newton solver for the nonlinear equation

506 Parameters

507 ----------

508 rhs : dtype_f

509 Right-hand side for the nonlinear system.

510 dt : float

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

512 u0 : dtype_u

513 Initial guess for the iterative solver.

514 t : float

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

517 Returns

518 -------

519 u : dtype_u

520 The solution as mesh.

521 """

522 # create new mesh object from u0 and set initial values for iteration

523 u = self.dtype_u(u0)

525 # start newton iteration

526 n, res = 0, np.inf

527 while n < self.newton_maxiter:

528 c1, c2, c3 = u

529 f = -np.array([0.013 * c1 + 1000 * c3 * c1, 2500 * c3 * c2, 0.013 * c1 + 1000 * c1 * c3 + 2500 * c2 * c3])

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

532 g = u - dt * f - rhs

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

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

536 if res < self.newton_tol or np.isnan(res):

537 break

539 # assemble (dg/du)^(-1)

540 dgInv = np.array(

541 [

542 [

543 (

544 2500000000.0 * c1 * c3**2 * dt**3

545 + 32500.0 * c1 * c3 * dt**3

546 + 3500000.0 * c1 * c3 * dt**2

547 + 13.0 * c1 * dt**2

548 + 1000.0 * c1 * dt

549 + 2500000.0 * c2 * c3 * dt**2

550 + 32.5 * c2 * dt**2

551 + 2500.0 * c2 * dt

552 + 2500000.0 * c3**2 * dt**2

553 + 32.5 * c3 * dt**2

554 + 3500.0 * c3 * dt

555 + 0.013 * dt

556 + 1.0

557 )

558 / (

559 2500000000.0 * c1 * c3**2 * dt**3

560 + 32500.0 * c1 * c3 * dt**3

561 + 3500000.0 * c1 * c3 * dt**2

562 + 13.0 * c1 * dt**2

563 + 1000.0 * c1 * dt

564 + 2500000000.0 * c2 * c3**2 * dt**3

565 + 65000.0 * c2 * c3 * dt**3

566 + 5000000.0 * c2 * c3 * dt**2

567 + 0.4225 * c2 * dt**3

568 + 65.0 * c2 * dt**2

569 + 2500.0 * c2 * dt

570 + 2500000000.0 * c3**3 * dt**3

571 + 65000.0 * c3**2 * dt**3

572 + 6000000.0 * c3**2 * dt**2

573 + 0.4225 * c3 * dt**3

574 + 91.0 * c3 * dt**2

575 + 4500.0 * c3 * dt

576 + 0.000169 * dt**2

577 + 0.026 * dt

578 + 1.0

579 ),

580 (2500000000.0 * c1 * c3**2 * dt**3 + 32500.0 * c1 * c3 * dt**3 + 2500000.0 * c1 * c3 * dt**2)

581 / (

582 2500000000.0 * c1 * c3**2 * dt**3

583 + 32500.0 * c1 * c3 * dt**3

584 + 3500000.0 * c1 * c3 * dt**2

585 + 13.0 * c1 * dt**2

586 + 1000.0 * c1 * dt

587 + 2500000000.0 * c2 * c3**2 * dt**3

588 + 65000.0 * c2 * c3 * dt**3

589 + 5000000.0 * c2 * c3 * dt**2

590 + 0.4225 * c2 * dt**3

591 + 65.0 * c2 * dt**2

592 + 2500.0 * c2 * dt

593 + 2500000000.0 * c3**3 * dt**3

594 + 65000.0 * c3**2 * dt**3

595 + 6000000.0 * c3**2 * dt**2

596 + 0.4225 * c3 * dt**3

597 + 91.0 * c3 * dt**2

598 + 4500.0 * c3 * dt

599 + 0.000169 * dt**2

600 + 0.026 * dt

601 + 1.0

602 ),

603 (

604 -2500000000.0 * c1 * c3**2 * dt**3

605 - 32500.0 * c1 * c3 * dt**3

606 - 3500000.0 * c1 * c3 * dt**2

607 - 13.0 * c1 * dt**2

608 - 1000.0 * c1 * dt

609 )

610 / (

611 2500000000.0 * c1 * c3**2 * dt**3

612 + 32500.0 * c1 * c3 * dt**3

613 + 3500000.0 * c1 * c3 * dt**2

614 + 13.0 * c1 * dt**2

615 + 1000.0 * c1 * dt

616 + 2500000000.0 * c2 * c3**2 * dt**3

617 + 65000.0 * c2 * c3 * dt**3

618 + 5000000.0 * c2 * c3 * dt**2

619 + 0.4225 * c2 * dt**3

620 + 65.0 * c2 * dt**2

621 + 2500.0 * c2 * dt

622 + 2500000000.0 * c3**3 * dt**3

623 + 65000.0 * c3**2 * dt**3

624 + 6000000.0 * c3**2 * dt**2

625 + 0.4225 * c3 * dt**3

626 + 91.0 * c3 * dt**2

627 + 4500.0 * c3 * dt

628 + 0.000169 * dt**2

629 + 0.026 * dt

630 + 1.0

631 ),

632 ],

633 [

634 (6250000000.0 * c2 * c3 * dt**2 + 81250.0 * c2 * dt**2)

635 / (

636 6250000000.0 * c1 * c3 * dt**2

637 + 2500000.0 * c1 * dt

638 + 6250000000.0 * c2 * c3 * dt**2

639 + 81250.0 * c2 * dt**2

640 + 6250000.0 * c2 * dt

641 + 6250000000.0 * c3**2 * dt**2

642 + 81250.0 * c3 * dt**2

643 + 8750000.0 * c3 * dt

644 + 32.5 * dt

645 + 2500.0

646 ),

647 (

648 2500000.0 * c1 * dt

649 + 6250000000.0 * c2 * c3 * dt**2

650 + 81250.0 * c2 * dt**2

651 + 6250000.0 * c2 * dt

652 + 2500000.0 * c3 * dt

653 + 32.5 * dt

654 + 2500.0

655 )

656 / (

657 6250000000.0 * c1 * c3 * dt**2

658 + 2500000.0 * c1 * dt

659 + 6250000000.0 * c2 * c3 * dt**2

660 + 81250.0 * c2 * dt**2

661 + 6250000.0 * c2 * dt

662 + 6250000000.0 * c3**2 * dt**2

663 + 81250.0 * c3 * dt**2

664 + 8750000.0 * c3 * dt

665 + 32.5 * dt

666 + 2500.0

667 ),

668 (-6250000000.0 * c2 * c3 * dt**2 - 81250.0 * c2 * dt**2 - 6250000.0 * c2 * dt)

669 / (

670 6250000000.0 * c1 * c3 * dt**2

671 + 2500000.0 * c1 * dt

672 + 6250000000.0 * c2 * c3 * dt**2

673 + 81250.0 * c2 * dt**2

674 + 6250000.0 * c2 * dt

675 + 6250000000.0 * c3**2 * dt**2

676 + 81250.0 * c3 * dt**2

677 + 8750000.0 * c3 * dt

678 + 32.5 * dt

679 + 2500.0

680 ),

681 ],

682 [

683 (-2500000.0 * c3**2 * dt**2 - 32.5 * c3 * dt**2 - 1000.0 * c3 * dt - 0.013 * dt)

684 / (

685 2500000.0 * c1 * c3 * dt**2

686 + 1000.0 * c1 * dt

687 + 2500000.0 * c2 * c3 * dt**2

688 + 32.5 * c2 * dt**2

689 + 2500.0 * c2 * dt

690 + 2500000.0 * c3**2 * dt**2

691 + 32.5 * c3 * dt**2

692 + 3500.0 * c3 * dt

693 + 0.013 * dt

694 + 1.0

695 ),

696 (-2500000.0 * c3**2 * dt**2 - 32.5 * c3 * dt**2 - 2500.0 * c3 * dt)

697 / (

698 2500000.0 * c1 * c3 * dt**2

699 + 1000.0 * c1 * dt

700 + 2500000.0 * c2 * c3 * dt**2

701 + 32.5 * c2 * dt**2

702 + 2500.0 * c2 * dt

703 + 2500000.0 * c3**2 * dt**2

704 + 32.5 * c3 * dt**2

705 + 3500.0 * c3 * dt

706 + 0.013 * dt

707 + 1.0

708 ),

709 (2500000.0 * c3**2 * dt**2 + 32.5 * c3 * dt**2 + 3500.0 * c3 * dt + 0.013 * dt + 1.0)

710 / (

711 2500000.0 * c1 * c3 * dt**2

712 + 1000.0 * c1 * dt

713 + 2500000.0 * c2 * c3 * dt**2

714 + 32.5 * c2 * dt**2

715 + 2500.0 * c2 * dt

716 + 2500000.0 * c3**2 * dt**2

717 + 32.5 * c3 * dt**2

718 + 3500.0 * c3 * dt

719 + 0.013 * dt

720 + 1.0

721 ),

722 ],

723 ]

724 )

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

727 u -= dgInv @ g

729 # increase iteration count and work counter

730 n += 1

731 self.work_counters['newton']()

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

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

735 elif np.isnan(res): # pragma: no cover

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

738 if n == self.newton_maxiter:

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

741 return u

744class JacobiElliptic(Problem):

745 r"""

746 Implement the Jacobi Elliptic non-linear problem:

748 .. math::

749 \begin{eqnarray*}

750 \frac{du}{dt} &=& vw, &\quad u(0) = 0, \\

751 \frac{dv}{dt} &=& -uw, &\quad v(0) = 1, \\

752 \frac{dw}{dt} &=& -0.51uv, &\quad w(0) = 1.

753 \end{eqnarray*}

755 Parameters

756 ----------

757 newton_maxiter : int, optional

758 Maximum number of Newton iteration in solve_system. The default is 200.

759 newton_tol : float, optional

760 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

761 stop_at_nan : bool, optional

762 Wheter to stop or not solve_system when getting NAN. The default is True.

764 Reference

765 ---------

766 Van Der Houwen, P. J., Sommeijer, B. P., & Van Der Veen, W. A. (1995).

767 Parallel iteration across the steps of high-order Runge-Kutta methods for

768 nonstiff initial value problems. Journal of computational and applied

769 mathematics, 60(3), 309-329.

770 """

772 dtype_u = mesh

773 dtype_f = mesh

775 def __init__(self, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

776 nvars = 3

777 u0 = (0.0, 1.0, 1.0)

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

780 self._makeAttributeAndRegister(

781 'u0', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

782 )

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

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

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

787 r"""

788 Routine to return initial conditions or to approximate exact solution using SciPy.

790 Parameters

791 ----------

792 t : float

793 Time at which the approximated exact solution is computed.

794 u_init : pySDC.implementations.problem_classes.Lorenz.dtype_u

795 Initial conditions for getting the exact solution.

796 t_init : float

797 The starting time.

799 Returns

800 -------

801 me : dtype_u

802 The approximated exact solution.

803 """

805 me = self.dtype_u(self.init)

807 if t > 0:

809 def eval_rhs(t, u):

810 r"""

811 Evaluate the right hand side, but switch the arguments for SciPy.

813 Args:

814 t (float): Time

815 u (numpy.ndarray): Solution at time t

817 Returns:

818 (numpy.ndarray): Right hand side evaluation

819 """

820 return self.eval_f(u, t)

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

823 else:

824 me[:] = self.u0

825 return me

827 def eval_f(self, u, t):

828 """

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

831 Parameters

832 ----------

833 u : dtype_u

834 Current values of the numerical solution.

835 t : float

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

838 Returns

839 -------

840 f : dtype_f

841 The right-hand side of the problem (one component).

842 """

843 f = self.dtype_f(self.init)

844 u1, u2, u3 = u

846 f[:] = np.array([u2 * u3, -u1 * u3, -0.51 * u1 * u2])

847 self.work_counters['rhs']()

848 return f

850 def solve_system(self, rhs, dt, u0, t):

851 """

852 Simple Newton solver for the nonlinear equation

854 Parameters

855 ----------

856 rhs : dtype_f

857 Right-hand side for the nonlinear system.

858 dt : float

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

860 u0 : dtype_u

861 Initial guess for the iterative solver.

862 t : float

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

865 Returns

866 -------

867 u : dtype_u

868 The solution as mesh.

869 """

870 # create new mesh object from u0 and set initial values for iteration

871 u = self.dtype_u(u0)

873 # start newton iteration

874 n, res = 0, np.inf

875 while n < self.newton_maxiter:

876 u1, u2, u3 = u

877 f = np.array([u2 * u3, -u1 * u3, -0.51 * u1 * u2])

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

880 g = u - dt * f - rhs

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

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

884 if res < self.newton_tol or np.isnan(res):

885 break

887 # assemble (dg/du)^(-1)

888 dgInv = np.array(

889 [

890 [

891 0.51 * dt**2 * u1**2 - 1.0,

892 0.51 * dt**2 * u1 * u2 - 1.0 * dt * u3,

893 1.0 * dt**2 * u1 * u3 - 1.0 * dt * u2,

894 ],

895 [

896 -0.51 * dt**2 * u1 * u2 + 1.0 * dt * u3,

897 -0.51 * dt**2 * u2**2 - 1.0,

898 1.0 * dt**2 * u2 * u3 + 1.0 * dt * u1,

899 ],

900 [

901 -0.51 * dt**2 * u1 * u3 + 0.51 * dt * u2,

902 0.51 * dt**2 * u2 * u3 + 0.51 * dt * u1,

903 -1.0 * dt**2 * u3**2 - 1.0,

904 ],

905 ]

906 )

907 dgInv /= (

908 1.02 * dt**3 * u1 * u2 * u3 + 0.51 * dt**2 * u1**2 - 0.51 * dt**2 * u2**2 - 1.0 * dt**2 * u3**2 - 1.0

909 )

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

912 u -= dgInv @ g

914 # increase iteration count and work counter

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): # pragma: no cover

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

923 if n == self.newton_maxiter:

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

926 return u