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

13import numpy as np 

14 

15from pySDC.core.Errors import ProblemError 

16from pySDC.core.Problem import ptype, WorkCounter 

17from pySDC.implementations.datatype_classes.mesh import mesh 

18 

19 

20class ProtheroRobinsonAutonomous(ptype): 

21 r""" 

22 Implement the Prothero-Robinson problem into autonomous form: 

23 

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

29 

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

34 

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

36 

37 .. math:: 

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

39 

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: 

43 

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) 

54 

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. 

67 

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

74 

75 dtype_u = mesh 

76 dtype_f = mesh 

77 

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

81 

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

89 

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

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

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

93 def g(self, t): 

94 return np.cos(t) 

95 

96 def dg(self, t): 

97 return -np.sin(t) 

98 

99 def dg2(self, t): 

100 return -np.cos(t) 

101 

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

103 # f(u,t) and Jacobian functions 

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

105 def f(self, u, t): 

106 raise NotImplementedError() 

107 

108 def f_LIN(self, u, t): 

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

110 

111 def f_NONLIN(self, u, t): 

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

113 

114 def dgInv(self, u, t): 

115 raise NotImplementedError() 

116 

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

121 

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 ) 

128 

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. 

135 

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. 

144 

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 

154 

155 def eval_f(self, u, t): 

156 """ 

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

158 

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

165 

166 Returns 

167 ------- 

168 f : dtype_f 

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

170 """ 

171 

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 

178 

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

180 """ 

181 Simple Newton solver for the nonlinear equation 

182 

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

193 

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) 

201 

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 

209 

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

211 g = u - dt * f - rhs 

212 

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 

217 

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 

222 

223 # increase iteration count and work counter 

224 n += 1 

225 self.work_counters['newton']() 

226 

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) 

231 

232 if n == self.newton_maxiter: 

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

234 

235 return u 

236 

237 

238class Kaps(ptype): 

239 r""" 

240 Implement the Kaps problem: 

241 

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

247 

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

251 

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. 

262 

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

269 

270 dtype_u = mesh 

271 dtype_f = mesh 

272 

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

276 

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

282 

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

284 r""" 

285 Routine to return initial conditions or exact solutions. 

286 

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. 

295 

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 

304 

305 def eval_f(self, u, t): 

306 """ 

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

308 

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

315 

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 

324 

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

326 self.work_counters['rhs']() 

327 return f 

328 

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

330 """ 

331 Simple Newton solver for the nonlinear equation 

332 

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

343 

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 

352 

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

358 

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

360 g = u - dt * f - rhs 

361 

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 

366 

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 ) 

374 

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

376 u -= dgInv @ g 

377 

378 # increase iteration count and work counter 

379 n += 1 

380 self.work_counters['newton']() 

381 

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) 

386 

387 if n == self.newton_maxiter: 

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

389 

390 return u 

391 

392 

393class ChemicalReaction3Var(ptype): 

394 r""" 

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

396 

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

405 

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

407 

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. 

416 

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

423 

424 dtype_u = mesh 

425 dtype_f = mesh 

426 

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

431 

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

437 

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

441 

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. 

450 

451 Returns 

452 ------- 

453 me : dtype_u 

454 The approximated exact solution. 

455 """ 

456 

457 me = self.dtype_u(self.init) 

458 

459 if t > 0: 

460 

461 def eval_rhs(t, u): 

462 r""" 

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

464 

465 Args: 

466 t (float): Time 

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

468 

469 Returns: 

470 (numpy.ndarray): Right hand side evaluation 

471 """ 

472 return self.eval_f(u, t) 

473 

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

475 else: 

476 me[:] = self.u0 

477 return me 

478 

479 def eval_f(self, u, t): 

480 """ 

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

482 

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

489 

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 

497 

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 

501 

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

503 """ 

504 Simple Newton solver for the nonlinear equation 

505 

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

516 

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) 

524 

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

530 

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

532 g = u - dt * f - rhs 

533 

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 

538 

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 ) 

725 

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

727 u -= dgInv @ g 

728 

729 # increase iteration count and work counter 

730 n += 1 

731 self.work_counters['newton']() 

732 

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) 

737 

738 if n == self.newton_maxiter: 

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

740 

741 return u 

742 

743 

744class JacobiElliptic(ptype): 

745 r""" 

746 Implement the Jacobi Elliptic non-linear problem: 

747 

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

754 

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. 

763 

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

771 

772 dtype_u = mesh 

773 dtype_f = mesh 

774 

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

779 

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

785 

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

789 

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. 

798 

799 Returns 

800 ------- 

801 me : dtype_u 

802 The approximated exact solution. 

803 """ 

804 

805 me = self.dtype_u(self.init) 

806 

807 if t > 0: 

808 

809 def eval_rhs(t, u): 

810 r""" 

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

812 

813 Args: 

814 t (float): Time 

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

816 

817 Returns: 

818 (numpy.ndarray): Right hand side evaluation 

819 """ 

820 return self.eval_f(u, t) 

821 

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

823 else: 

824 me[:] = self.u0 

825 return me 

826 

827 def eval_f(self, u, t): 

828 """ 

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

830 

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

837 

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 

845 

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

847 self.work_counters['rhs']() 

848 return f 

849 

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

851 """ 

852 Simple Newton solver for the nonlinear equation 

853 

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