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

1#!/usr/bin/env python3 

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

3""" 

4Implementation of systems test problem ODEs. 

5 

6 

7Reference : 

8 

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

13 

14import numpy as np 

15 

16from pySDC.core.errors import ProblemError 

17from pySDC.core.problem import Problem, WorkCounter 

18from pySDC.implementations.datatype_classes.mesh import mesh 

19 

20 

21class ProtheroRobinsonAutonomous(Problem): 

22 r""" 

23 Implement the Prothero-Robinson problem into autonomous form: 

24 

25 .. math:: 

26 \begin{eqnarray*} 

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

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

29 \end{eqnarray*} 

30 

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

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

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

34 :math:`g(t)=\cos(t)`. 

35 

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

37 

38 .. math:: 

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

40 

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

42 and overload the `g`, `dg` and `dg2` methods. For instance, 

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

44 

45 >>> class MyProtheroRobinson(ProtheroRobinsonAutonomous): 

46 >>> 

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

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

49 >>> 

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

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

52 >>> 

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

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

55 

56 Parameters 

57 ---------- 

58 epsilon : float, optional 

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

60 nonLinear : bool, optional 

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

62 newton_maxiter : int, optional 

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

64 newton_tol : float, optional 

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

66 stop_at_nan : bool, optional 

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

68 

69 Reference 

70 --------- 

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

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

73 pp. 145–162. 

74 """ 

75 

76 dtype_u = mesh 

77 dtype_f = mesh 

78 

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

80 nvars = 2 

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

82 

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

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

85 self._makeAttributeAndRegister( 

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

87 ) 

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

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

90 

91 # ------------------------------------------------------------------------- 

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

93 # ------------------------------------------------------------------------- 

94 def g(self, t): 

95 return np.cos(t) 

96 

97 def dg(self, t): 

98 return -np.sin(t) 

99 

100 def dg2(self, t): 

101 return -np.cos(t) 

102 

103 # ------------------------------------------------------------------------- 

104 # f(u,t) and Jacobian functions 

105 # ------------------------------------------------------------------------- 

106 def f(self, u, t): 

107 raise NotImplementedError() 

108 

109 def f_LIN(self, u, t): 

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

111 

112 def f_NONLIN(self, u, t): 

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

114 

115 def dgInv(self, u, t): 

116 raise NotImplementedError() 

117 

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

119 e = self.epsilon 

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

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

122 

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

124 e = self.epsilon 

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

126 return np.array( 

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

128 ) 

129 

130 # ------------------------------------------------------------------------- 

131 # pySDC required methods 

132 # ------------------------------------------------------------------------- 

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

134 r""" 

135 Routine to return initial conditions or exact solutions. 

136 

137 Parameters 

138 ---------- 

139 t : float 

140 Time at which the exact solution is computed. 

141 u_init : dtype_u 

142 Initial conditions for getting the exact solution. 

143 t_init : float 

144 The starting time. 

145 

146 Returns 

147 ------- 

148 u : dtype_u 

149 The exact solution. 

150 """ 

151 u = self.dtype_u(self.init) 

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

153 u[1] = t 

154 return u 

155 

156 def eval_f(self, u, t): 

157 """ 

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

159 

160 Parameters 

161 ---------- 

162 u : dtype_u 

163 Current values of the numerical solution. 

164 t : float 

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

166 

167 Returns 

168 ------- 

169 f : dtype_f 

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

171 """ 

172 

173 f = self.dtype_f(self.init) 

174 u, t = u 

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

176 f[1] = 1 

177 self.work_counters['rhs']() 

178 return f 

179 

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

181 """ 

182 Simple Newton solver for the nonlinear equation 

183 

184 Parameters 

185 ---------- 

186 rhs : dtype_f 

187 Right-hand side for the nonlinear system. 

188 dt : float 

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

190 u0 : dtype_u 

191 Initial guess for the iterative solver. 

192 t : float 

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

194 

195 Returns 

196 ------- 

197 u : dtype_u 

198 The solution as mesh. 

199 """ 

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

201 u = self.dtype_u(u0) 

202 

203 # start newton iteration 

204 n, res = 0, np.inf 

205 while n < self.newton_maxiter: 

206 # evaluate RHS 

207 f = self.dtype_u(u) 

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

209 f[1] = 1 

210 

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

212 g = u - dt * f - rhs 

213 

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

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

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

217 break 

218 

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

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

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

222 u -= dgInv @ g 

223 

224 # increase iteration count and work counter 

225 n += 1 

226 self.work_counters['newton']() 

227 

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

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

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

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

232 

233 if n == self.newton_maxiter: 

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

235 

236 return u 

237 

238 

239class Kaps(Problem): 

240 r""" 

241 Implement the Kaps problem: 

242 

243 .. math:: 

244 \begin{eqnarray*} 

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

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

247 \end{eqnarray*} 

248 

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

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

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

252 

253 Parameters 

254 ---------- 

255 epsilon : float, optional 

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

257 newton_maxiter : int, optional 

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

259 newton_tol : float, optional 

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

261 stop_at_nan : bool, optional 

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

263 

264 Reference 

265 --------- 

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

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

268 12(5), 1000-1028. 

269 """ 

270 

271 dtype_u = mesh 

272 dtype_f = mesh 

273 

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

275 nvars = 2 

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

277 

278 self._makeAttributeAndRegister( 

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

280 ) 

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

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

283 

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

285 r""" 

286 Routine to return initial conditions or exact solutions. 

287 

288 Parameters 

289 ---------- 

290 t : float 

291 Time at which the exact solution is computed. 

292 u_init : dtype_u 

293 Initial conditions for getting the exact solution. 

294 t_init : float 

295 The starting time. 

296 

297 Returns 

298 ------- 

299 u : dtype_u 

300 The exact solution. 

301 """ 

302 u = self.dtype_u(self.init) 

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

304 return u 

305 

306 def eval_f(self, u, t): 

307 """ 

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

309 

310 Parameters 

311 ---------- 

312 u : dtype_u 

313 Current values of the numerical solution. 

314 t : float 

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

316 

317 Returns 

318 ------- 

319 f : dtype_f 

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

321 """ 

322 f = self.dtype_f(self.init) 

323 eps = self.epsilon 

324 x, y = u 

325 

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

327 self.work_counters['rhs']() 

328 return f 

329 

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

331 """ 

332 Simple Newton solver for the nonlinear equation 

333 

334 Parameters 

335 ---------- 

336 rhs : dtype_f 

337 Right-hand side for the nonlinear system. 

338 dt : float 

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

340 u0 : dtype_u 

341 Initial guess for the iterative solver. 

342 t : float 

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

344 

345 Returns 

346 ------- 

347 u : dtype_u 

348 The solution as mesh. 

349 """ 

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

351 u = self.dtype_u(u0) 

352 eps = self.epsilon 

353 

354 # start newton iteration 

355 n, res = 0, np.inf 

356 while n < self.newton_maxiter: 

357 x, y = u 

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

359 

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

361 g = u - dt * f - rhs 

362 

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

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

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

366 break 

367 

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

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

370 dgInv = ( 

371 1 

372 / prefactor 

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

374 ) 

375 

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

377 u -= dgInv @ g 

378 

379 # increase iteration count and work counter 

380 n += 1 

381 self.work_counters['newton']() 

382 

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

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

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

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

387 

388 if n == self.newton_maxiter: 

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

390 

391 return u 

392 

393 

394class ChemicalReaction3Var(Problem): 

395 r""" 

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

397 

398 .. math:: 

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

400 \begin{pmatrix} 

401 0.013+1000u_3 & 0 & 0 \\ 

402 0 & 2500u_3 0 \\ 

403 0.013 & 0 & 1000u_1 + 2500u_2 

404 \end{pmatrix} 

405 {\bf u}, 

406 

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

408 

409 Parameters 

410 ---------- 

411 newton_maxiter : int, optional 

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

413 newton_tol : float, optional 

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

415 stop_at_nan : bool, optional 

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

417 

418 Reference 

419 --------- 

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

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

422 12(5), 1000-1028. 

423 """ 

424 

425 dtype_u = mesh 

426 dtype_f = mesh 

427 

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

429 nvars = 3 

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

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

432 

433 self._makeAttributeAndRegister( 

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

435 ) 

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

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

438 

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

440 r""" 

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

442 

443 Parameters 

444 ---------- 

445 t : float 

446 Time at which the approximated exact solution is computed. 

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

448 Initial conditions for getting the exact solution. 

449 t_init : float 

450 The starting time. 

451 

452 Returns 

453 ------- 

454 me : dtype_u 

455 The approximated exact solution. 

456 """ 

457 

458 me = self.dtype_u(self.init) 

459 

460 if t > 0: 

461 

462 def eval_rhs(t, u): 

463 r""" 

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

465 

466 Args: 

467 t (float): Time 

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

469 

470 Returns: 

471 (numpy.ndarray): Right hand side evaluation 

472 """ 

473 return self.eval_f(u, t) 

474 

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

476 else: 

477 me[:] = self.u0 

478 return me 

479 

480 def eval_f(self, u, t): 

481 """ 

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

483 

484 Parameters 

485 ---------- 

486 u : dtype_u 

487 Current values of the numerical solution. 

488 t : float 

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

490 

491 Returns 

492 ------- 

493 f : dtype_f 

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

495 """ 

496 f = self.dtype_f(self.init) 

497 c1, c2, c3 = u 

498 

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

500 self.work_counters['rhs']() 

501 return f 

502 

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

504 """ 

505 Simple Newton solver for the nonlinear equation 

506 

507 Parameters 

508 ---------- 

509 rhs : dtype_f 

510 Right-hand side for the nonlinear system. 

511 dt : float 

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

513 u0 : dtype_u 

514 Initial guess for the iterative solver. 

515 t : float 

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

517 

518 Returns 

519 ------- 

520 u : dtype_u 

521 The solution as mesh. 

522 """ 

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

524 u = self.dtype_u(u0) 

525 

526 # start newton iteration 

527 n, res = 0, np.inf 

528 while n < self.newton_maxiter: 

529 c1, c2, c3 = u 

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

531 

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

533 g = u - dt * f - rhs 

534 

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

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

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

538 break 

539 

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

541 dgInv = np.array( 

542 [ 

543 [ 

544 ( 

545 2500000000.0 * c1 * c3**2 * dt**3 

546 + 32500.0 * c1 * c3 * dt**3 

547 + 3500000.0 * c1 * c3 * dt**2 

548 + 13.0 * c1 * dt**2 

549 + 1000.0 * c1 * dt 

550 + 2500000.0 * c2 * c3 * dt**2 

551 + 32.5 * c2 * dt**2 

552 + 2500.0 * c2 * dt 

553 + 2500000.0 * c3**2 * dt**2 

554 + 32.5 * c3 * dt**2 

555 + 3500.0 * c3 * dt 

556 + 0.013 * dt 

557 + 1.0 

558 ) 

559 / ( 

560 2500000000.0 * c1 * c3**2 * dt**3 

561 + 32500.0 * c1 * c3 * dt**3 

562 + 3500000.0 * c1 * c3 * dt**2 

563 + 13.0 * c1 * dt**2 

564 + 1000.0 * c1 * dt 

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

566 + 65000.0 * c2 * c3 * dt**3 

567 + 5000000.0 * c2 * c3 * dt**2 

568 + 0.4225 * c2 * dt**3 

569 + 65.0 * c2 * dt**2 

570 + 2500.0 * c2 * dt 

571 + 2500000000.0 * c3**3 * dt**3 

572 + 65000.0 * c3**2 * dt**3 

573 + 6000000.0 * c3**2 * dt**2 

574 + 0.4225 * c3 * dt**3 

575 + 91.0 * c3 * dt**2 

576 + 4500.0 * c3 * dt 

577 + 0.000169 * dt**2 

578 + 0.026 * dt 

579 + 1.0 

580 ), 

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

582 / ( 

583 2500000000.0 * c1 * c3**2 * dt**3 

584 + 32500.0 * c1 * c3 * dt**3 

585 + 3500000.0 * c1 * c3 * dt**2 

586 + 13.0 * c1 * dt**2 

587 + 1000.0 * c1 * dt 

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

589 + 65000.0 * c2 * c3 * dt**3 

590 + 5000000.0 * c2 * c3 * dt**2 

591 + 0.4225 * c2 * dt**3 

592 + 65.0 * c2 * dt**2 

593 + 2500.0 * c2 * dt 

594 + 2500000000.0 * c3**3 * dt**3 

595 + 65000.0 * c3**2 * dt**3 

596 + 6000000.0 * c3**2 * dt**2 

597 + 0.4225 * c3 * dt**3 

598 + 91.0 * c3 * dt**2 

599 + 4500.0 * c3 * dt 

600 + 0.000169 * dt**2 

601 + 0.026 * dt 

602 + 1.0 

603 ), 

604 ( 

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

606 - 32500.0 * c1 * c3 * dt**3 

607 - 3500000.0 * c1 * c3 * dt**2 

608 - 13.0 * c1 * dt**2 

609 - 1000.0 * c1 * dt 

610 ) 

611 / ( 

612 2500000000.0 * c1 * c3**2 * dt**3 

613 + 32500.0 * c1 * c3 * dt**3 

614 + 3500000.0 * c1 * c3 * dt**2 

615 + 13.0 * c1 * dt**2 

616 + 1000.0 * c1 * dt 

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

618 + 65000.0 * c2 * c3 * dt**3 

619 + 5000000.0 * c2 * c3 * dt**2 

620 + 0.4225 * c2 * dt**3 

621 + 65.0 * c2 * dt**2 

622 + 2500.0 * c2 * dt 

623 + 2500000000.0 * c3**3 * dt**3 

624 + 65000.0 * c3**2 * dt**3 

625 + 6000000.0 * c3**2 * dt**2 

626 + 0.4225 * c3 * dt**3 

627 + 91.0 * c3 * dt**2 

628 + 4500.0 * c3 * dt 

629 + 0.000169 * dt**2 

630 + 0.026 * dt 

631 + 1.0 

632 ), 

633 ], 

634 [ 

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

636 / ( 

637 6250000000.0 * c1 * c3 * dt**2 

638 + 2500000.0 * c1 * dt 

639 + 6250000000.0 * c2 * c3 * dt**2 

640 + 81250.0 * c2 * dt**2 

641 + 6250000.0 * c2 * dt 

642 + 6250000000.0 * c3**2 * dt**2 

643 + 81250.0 * c3 * dt**2 

644 + 8750000.0 * c3 * dt 

645 + 32.5 * dt 

646 + 2500.0 

647 ), 

648 ( 

649 2500000.0 * c1 * dt 

650 + 6250000000.0 * c2 * c3 * dt**2 

651 + 81250.0 * c2 * dt**2 

652 + 6250000.0 * c2 * dt 

653 + 2500000.0 * c3 * dt 

654 + 32.5 * dt 

655 + 2500.0 

656 ) 

657 / ( 

658 6250000000.0 * c1 * c3 * dt**2 

659 + 2500000.0 * c1 * dt 

660 + 6250000000.0 * c2 * c3 * dt**2 

661 + 81250.0 * c2 * dt**2 

662 + 6250000.0 * c2 * dt 

663 + 6250000000.0 * c3**2 * dt**2 

664 + 81250.0 * c3 * dt**2 

665 + 8750000.0 * c3 * dt 

666 + 32.5 * dt 

667 + 2500.0 

668 ), 

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

670 / ( 

671 6250000000.0 * c1 * c3 * dt**2 

672 + 2500000.0 * c1 * dt 

673 + 6250000000.0 * c2 * c3 * dt**2 

674 + 81250.0 * c2 * dt**2 

675 + 6250000.0 * c2 * dt 

676 + 6250000000.0 * c3**2 * dt**2 

677 + 81250.0 * c3 * dt**2 

678 + 8750000.0 * c3 * dt 

679 + 32.5 * dt 

680 + 2500.0 

681 ), 

682 ], 

683 [ 

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

685 / ( 

686 2500000.0 * c1 * c3 * dt**2 

687 + 1000.0 * c1 * dt 

688 + 2500000.0 * c2 * c3 * dt**2 

689 + 32.5 * c2 * dt**2 

690 + 2500.0 * c2 * dt 

691 + 2500000.0 * c3**2 * dt**2 

692 + 32.5 * c3 * dt**2 

693 + 3500.0 * c3 * dt 

694 + 0.013 * dt 

695 + 1.0 

696 ), 

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

698 / ( 

699 2500000.0 * c1 * c3 * dt**2 

700 + 1000.0 * c1 * dt 

701 + 2500000.0 * c2 * c3 * dt**2 

702 + 32.5 * c2 * dt**2 

703 + 2500.0 * c2 * dt 

704 + 2500000.0 * c3**2 * dt**2 

705 + 32.5 * c3 * dt**2 

706 + 3500.0 * c3 * dt 

707 + 0.013 * dt 

708 + 1.0 

709 ), 

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

711 / ( 

712 2500000.0 * c1 * c3 * dt**2 

713 + 1000.0 * c1 * dt 

714 + 2500000.0 * c2 * c3 * dt**2 

715 + 32.5 * c2 * dt**2 

716 + 2500.0 * c2 * dt 

717 + 2500000.0 * c3**2 * dt**2 

718 + 32.5 * c3 * dt**2 

719 + 3500.0 * c3 * dt 

720 + 0.013 * dt 

721 + 1.0 

722 ), 

723 ], 

724 ] 

725 ) 

726 

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

728 u -= dgInv @ g 

729 

730 # increase iteration count and work counter 

731 n += 1 

732 self.work_counters['newton']() 

733 

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

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

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

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

738 

739 if n == self.newton_maxiter: 

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

741 

742 return u 

743 

744 

745class JacobiElliptic(Problem): 

746 r""" 

747 Implement the Jacobi Elliptic non-linear problem: 

748 

749 .. math:: 

750 \begin{eqnarray*} 

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

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

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

754 \end{eqnarray*} 

755 

756 Parameters 

757 ---------- 

758 newton_maxiter : int, optional 

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

760 newton_tol : float, optional 

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

762 stop_at_nan : bool, optional 

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

764 

765 Reference 

766 --------- 

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

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

769 nonstiff initial value problems. Journal of computational and applied 

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

771 """ 

772 

773 dtype_u = mesh 

774 dtype_f = mesh 

775 

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

777 nvars = 3 

778 u0 = (0.0, 1.0, 1.0) 

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

780 

781 self._makeAttributeAndRegister( 

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

783 ) 

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

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

786 

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

788 r""" 

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

790 

791 Parameters 

792 ---------- 

793 t : float 

794 Time at which the approximated exact solution is computed. 

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

796 Initial conditions for getting the exact solution. 

797 t_init : float 

798 The starting time. 

799 

800 Returns 

801 ------- 

802 me : dtype_u 

803 The approximated exact solution. 

804 """ 

805 

806 me = self.dtype_u(self.init) 

807 

808 if t > 0: 

809 

810 def eval_rhs(t, u): 

811 r""" 

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

813 

814 Args: 

815 t (float): Time 

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

817 

818 Returns: 

819 (numpy.ndarray): Right hand side evaluation 

820 """ 

821 return self.eval_f(u, t) 

822 

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

824 else: 

825 me[:] = self.u0 

826 return me 

827 

828 def eval_f(self, u, t): 

829 """ 

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

831 

832 Parameters 

833 ---------- 

834 u : dtype_u 

835 Current values of the numerical solution. 

836 t : float 

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

838 

839 Returns 

840 ------- 

841 f : dtype_f 

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

843 """ 

844 f = self.dtype_f(self.init) 

845 u1, u2, u3 = u 

846 

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

848 self.work_counters['rhs']() 

849 return f 

850 

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

852 """ 

853 Simple Newton solver for the nonlinear equation 

854 

855 Parameters 

856 ---------- 

857 rhs : dtype_f 

858 Right-hand side for the nonlinear system. 

859 dt : float 

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

861 u0 : dtype_u 

862 Initial guess for the iterative solver. 

863 t : float