Coverage for pySDC/implementations/problem_classes/Battery.py: 99%

1import numpy as np 

2 

3from pySDC.core.errors import ParameterError, ProblemError 

4from pySDC.core.problem import Problem, WorkCounter 

5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh 

6 

7 

8class battery_n_capacitors(Problem): 

9 r""" 

10 Example implementing the battery drain model with :math:`N` capacitors, where :math:`N` is an arbitrary integer greater than zero. 

11 First, the capacitor :math:`C` serves as a battery and provides energy. When the voltage of the capacitor :math:`v_{C_n}` for 

12 :math:`n=1,..,N` drops below their reference value :math:`V_{ref,n-1}`, the circuit switches to the next capacitor. If all capacitors 

13 has dropped below their reference value, the voltage source :math:`V_s` provides further energy. The problem of simulating the 

14 battery draining has :math:`N + 1` different states. Each of this state can be expressed as a nonhomogeneous linear system of 

15 ordinary differential equations (ODEs) 

16 

17 .. math:: 

18 \frac{d u(t)}{dt} = A_k u(t) + f_k (t) 

19 

20 for :math:`k=1,..,N+1` using an initial condition. 

21 

22 Parameters 

23 ---------- 

24 ncapacitors : int, optional 

25 Number of capacitors :math:`n_{capacitors}` in the circuit. 

26 Vs : float, optional 

27 Voltage at the voltage source :math:`V_s`. 

28 Rs : float, optional 

29 Resistance of the resistor :math:`R_s` at the voltage source. 

30 C : np.1darray, optional 

31 Capacitances of the capacitors :math:`C_n`. 

32 R : float, optional 

33 Resistance for the load :math:`R_\ell`. 

34 L : float, optional 

35 Inductance of inductor :math:`L`. 

36 alpha : float, optional 

37 Factor greater than zero to describe the storage of the capacitor(s). 

38 V_ref : np.1darray, optional 

39 Array contains the reference values greater than zero for each capacitor to switch to the next energy source. 

40 

41 Attributes 

42 ---------- 

43 A : matrix 

44 Coefficients matrix of the linear system of ordinary differential equations (ODEs). 

45 switch_A : dict 

46 Dictionary that contains the coefficients for the coefficient matrix A. 

47 switch_f : dict 

48 Dictionary that contains the coefficients of the right-hand side f of the ODE system. 

49 t_switch : float 

50 Time point of the discrete event found by switch estimation. 

51 nswitches : int 

52 Number of switches found by switch estimation. 

53 

54 Note 

55 ---- 

56 The array containing the capacitances :math:`C_n` and the array containing the reference values :math:`V_{ref, n-1}` 

57 for each capacitor must be equal to the number of capacitors :math:`n_{capacitors}`. 

58 """ 

59 

60 dtype_u = mesh 

61 dtype_f = imex_mesh 

62 

63 def __init__(self, ncapacitors=2, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.2, V_ref=None): 

64 """Initialization routine""" 

65 n = ncapacitors 

66 nvars = n + 1 

67 

68 if C is None: 

69 if ncapacitors == 1: 

70 C = np.array([1.0]) 

71 elif ncapacitors == 2: 

72 C = np.array([1.0, 1.0]) 

73 else: 

74 raise ParameterError(f"No default value for C is set up for ncapacitors={ncapacitors}") 

75 else: 

76 msgC = "ERROR for capacitance C: C has to be an np.ndarray and/or length of array needs to be equal to number of capacitances" 

77 assert all([type(C) == np.ndarray, np.shape(C)[0] == n]), msgC 

78 

79 if V_ref is None: 

80 if ncapacitors == 1: 

81 V_ref = np.array([1.0]) 

82 elif ncapacitors == 2: 

83 V_ref = np.array([1.0, 1.0]) 

84 else: 

85 raise ParameterError(f"No default value for V_ref is set up for ncapacitors={ncapacitors}") 

86 else: 

87 assertions_V_ref_1 = [ 

88 type(V_ref) == np.ndarray, 

89 np.shape(V_ref)[0] == n, 

90 ] 

91 msg1 = "ERROR for reference value V_ref: V_ref has to be an np.ndarray and/or length of array needs to be equal to number of capacitances" 

92 assert all(assertions_V_ref_1), msg1 

93 

94 msg2 = "ERROR for V_ref: every entry of V_ref has to be positive and smaller than alpha!" 

95 assert np.all(V_ref > 0) and np.all(V_ref < alpha), msg2 

96 

97 # invoke super init, passing number of dofs, dtype_u and dtype_f 

98 super().__init__(init=(nvars, None, np.dtype('float64'))) 

99 self._makeAttributeAndRegister( 

100 'nvars', 'ncapacitors', 'Vs', 'Rs', 'C', 'R', 'L', 'alpha', 'V_ref', localVars=locals(), readOnly=True 

101 ) 

102 

103 self.switch_A, self.switch_f = self.get_problem_dict() 

104 self.A = self.switch_A[0] 

105 

106 self.t_switch = None 

107 self.nswitches = 0 

108 

109 def eval_f(self, u, t): 

110 r""" 

111 Routine to evaluate the right-hand side of the problem. Let :math:`v_k:=v_{C_k}` be the voltage of capacitor :math:`C_k` for :math:`k=1,..,N` 

112 with reference value :math:`V_{ref, k-1}`. No switch estimator is used: For :math:`N = 3` there are :math:`N + 1 = 4` different states of the battery: 

113 

114 :math:`C_1` supplies energy if: 

115 

116 .. math:: 

117 v_1 > V_{ref,0}, v_2 > V_{ref,1}, v_3 > V_{ref,2}, 

118 

119 :math:`C_2` supplies energy if: 

120 

121 .. math:: 

122 v_1 \leq V_{ref,0}, v_2 > V_{ref,1}, v_3 > V_{ref,2}, 

123 

124 :math:`C_3` supplies energy if: 

125 

126 .. math:: 

127 v_1 \leq V_{ref,0}, v_2 \leq V_{ref,1}, v_3 > V_{ref,2}, 

128 

129 :math:`V_s` supplies energy if: 

130 

131 .. math:: 

132 v_1 \leq V_{ref,0}, v_2 \leq V_{ref,1}, v_3 \leq V_{ref,2}. 

133 

134 :math:`max_{index}` is initialized to :math:`-1`. List "switch" contains a True if :math:`v_k \leq V_{ref,k-1}` is satisfied. 

135 

136 - Is no True there (i.e., :math:`max_{index}=-1`), we are in the first case. 

137 - :math:`max_{index}=k\geq 0` means we are in the :math:`(k+1)`-th case. 

138 So, the actual RHS has key :math:`max_{index}`-1 in the dictionary self.switch_f. 

139 

140 In case of using the switch estimator, we count the number of switches which illustrates in which case of voltage source we are. 

141 

142 Parameters 

143 ---------- 

144 u : dtype_u 

145 Current values of the numerical solution. 

146 t : float 

147 Current time of the numerical solution is computed. 

148 

149 Returns 

150 ------- 

151 f : dtype_f 

152 The right-hand side of the problem. 

153 """ 

154 

155 f = self.dtype_f(self.init, val=0.0) 

156 f.impl[:] = self.A.dot(u) 

157 

158 if self.t_switch is not None: 

159 f.expl[:] = self.switch_f[self.nswitches] 

160 

161 else: 

162 # proof all switching conditions and find largest index where it drops below V_ref 

163 switch = [True if u[k] <= self.V_ref[k - 1] else False for k in range(1, len(u))] 

164 max_index = max([k if switch[k] else -1 for k in range(len(switch))]) 

165 

166 if max_index == -1: 

167 f.expl[:] = self.switch_f[0] 

168 

169 else: 

170 f.expl[:] = self.switch_f[max_index + 1] 

171 

172 return f 

173 

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

175 r""" 

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

177 

178 Parameters 

179 ---------- 

180 rhs : dtype_f 

181 Right-hand side for the linear system. 

182 factor : float 

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

184 u0 : dtype_u 

185 Initial guess for the iterative solver. 

186 t : float 

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

188 

189 Returns 

190 ------- 

191 me : dtype_u 

192 The solution as mesh. 

193 """ 

194 

195 if self.t_switch is not None: 

196 self.A = self.switch_A[self.nswitches] 

197 

198 else: 

199 # proof all switching conditions and find largest index where it drops below V_ref 

200 switch = [True if rhs[k] <= self.V_ref[k - 1] else False for k in range(1, len(rhs))] 

201 max_index = max([k if switch[k] else -1 for k in range(len(switch))]) 

202 if max_index == -1: 

203 self.A = self.switch_A[0] 

204 

205 else: 

206 self.A = self.switch_A[max_index + 1] 

207 

208 me = self.dtype_u(self.init) 

209 me[:] = np.linalg.solve(np.eye(self.nvars) - factor * self.A, rhs) 

210 return me 

211 

212 def u_exact(self, t): 

213 r""" 

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

215 

216 Parameters 

217 ---------- 

218 t : float 

219 Time of the exact solution. 

220 

221 Returns 

222 ------- 

223 me : dtype_u 

224 The exact solution. 

225 """ 

226 assert t == 0, 'ERROR: u_exact only valid for t=0' 

227 

228 me = self.dtype_u(self.init) 

229 

230 me[0] = 0.0 # cL 

231 me[1:] = self.alpha * self.V_ref # vC's 

232 return me 

233 

234 def get_switching_info(self, u, t): 

235 """ 

236 Provides information about a discrete event for one subinterval. 

237 

238 Parameters 

239 ---------- 

240 u : dtype_u 

241 Current values of the numerical solution. 

242 t : float 

243 Current time of the numerical solution is computed. 

244 

245 Returns 

246 ------- 

247 switch_detected : bool 

248 Indicates if a switch is found or not. 

249 m_guess : int 

250 Index of collocation node inside one subinterval of where the discrete event was found. 

251 state_function : list 

252 Contains values of the state function (for interpolation). 

253 """ 

254 

255 switch_detected = False 

256 m_guess = -100 

257 break_flag = False 

258 k_detected = 1 

259 for m in range(1, len(u)): 

260 for k in range(1, self.nvars): 

261 h_prev_node = u[m - 1][k] - self.V_ref[k - 1] 

262 h_curr_node = u[m][k] - self.V_ref[k - 1] 

263 if h_prev_node > 0 and h_curr_node <= 0: 

264 switch_detected = True 

265 m_guess = m - 1 

266 k_detected = k 

267 break_flag = True 

268 break 

269 

270 if break_flag: 

271 break 

272 

273 state_function = [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))] 

274 

275 return switch_detected, m_guess, state_function 

276 

277 def count_switches(self): 

278 """ 

279 Counts the number of switches. This function is called when a switch is found inside the range of tolerance 

280 (in pySDC/projects/PinTSimE/switch_estimator.py) 

281 """ 

282 

283 self.nswitches += 1 

284 

285 def get_problem_dict(self): 

286 """ 

287 Helper to create dictionaries for both the coefficent matrix of the ODE system and the nonhomogeneous part. 

288 """ 

289 

290 n = self.ncapacitors 

291 v = np.zeros(n + 1) 

292 v[0] = 1 

293 A, f = dict(), dict() 

294 A = {k: np.diag(-1 / (self.C[k] * self.R) * np.roll(v, k + 1)) for k in range(n)} 

295 A.update({n: np.diag(-(self.Rs + self.R) / self.L * v)}) 

296 f = {k: np.zeros(n + 1) for k in range(n)} 

297 f.update({n: self.Vs / self.L * v}) 

298 return A, f 

299 

300 

301class battery(battery_n_capacitors): 

302 r""" 

303 Example implementing the battery drain model with :math:`N=1` capacitor, inherits from ``battery_n_capacitors``. This model is an example 

304 of a discontinuous problem. The state function :math:`decides` which differential equation is solved. When the state function has 

305 a sign change the dynamics of the solution changes by changing the differential equation. The ODE system of this model is given by 

306 the following equations: 

307 

308 If :math:`h(v_1) := v_1 - V_{ref, 0} > 0:` 

309 

310 .. math:: 

311 \frac{d i_L (t)}{dt} = 0, 

312 

313 .. math:: 

314 \frac{d v_1 (t)}{dt} = -\frac{1}{CR}v_1 (t), 

315 

316 else: 

317 

318 .. math:: 

319 \frac{d i_L(t)}{dt} = -\frac{R_s + R}{L}i_L (t) + \frac{1}{L} V_s, 

320 

321 .. math:: 

322 \frac{d v_1(t)}{dt} = 0, 

323 

324 where :math:`i_L` denotes the function of the current over time :math:`t`. 

325 

326 Parameters 

327 ---------- 

328 ncapacitors : int, optional 

329 Number of capacitors in the circuit. Has to be one here, since the right-hand side and the solver are 

330 implemented for a single capacitor only. 

331 Vs : float, optional 

332 Voltage at the voltage source :math:`V_s`. 

333 Rs : float, optional 

334 Resistance of the resistor :math:`R_s` at the voltage source. 

335 C : np.1darray, optional 

336 Capacitance of the capacitor as an array of length one. ``None`` means ``np.array([1.0])``. 

337 R : float, optional 

338 Resistance for the load :math:`R_\ell`. 

339 L : float, optional 

340 Inductance of inductor :math:`L`. 

341 alpha : float, optional 

342 Factor greater than zero to describe the storage of the capacitor. The initial voltage of the capacitor is 

343 :math:`\alpha V_{ref, 0}`. 

344 V_ref : np.1darray, optional 

345 Reference value :math:`V_{ref, 0}` greater than zero for the battery to switch to the voltage source, as an 

346 array of length one. ``None`` means ``np.array([1.0])``. 

347 

348 Note 

349 ---- 

350 This class has the same attributes as the class it inherits from. 

351 """ 

352 

353 dtype_f = imex_mesh 

354 

355 def __init__(self, ncapacitors=1, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.2, V_ref=None): 

356 """Initialization routine""" 

357 super().__init__(ncapacitors, Vs, Rs, C, R, L, alpha, V_ref) 

358 

359 def eval_f(self, u, t): 

360 """ 

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

362 

363 Parameters 

364 ---------- 

365 u : dtype_u 

366 Current values of the numerical solution. 

367 t : float 

368 Current time of the numerical solution is computed. 

369 

370 Returns 

371 ------- 

372 f : dtype_f 

373 The right-hand side of the problem. 

374 """ 

375 

376 f = self.dtype_f(self.init, val=0.0) 

377 f.impl[:] = self.A.dot(u) 

378 

379 t_switch = np.inf if self.t_switch is None else self.t_switch 

380 

381 if u[1] - self.V_ref[0] <= 0 or t >= t_switch: 

382 f.expl[0] = self.Vs / self.L 

383 

384 else: 

385 f.expl[0] = 0 

386 

387 return f 

388 

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

390 r""" 

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

392 

393 Parameters 

394 ---------- 

395 rhs : dtype_f 

396 Right-hand side for the linear system. 

397 factor : float 

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

399 u0 : dtype_u 

400 Initial guess for the iterative solver. 

401 t : float 

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

403 

404 Returns 

405 ------- 

406 me : dtype_u 

407 The solution as mesh. 

408 """ 

409 self.A = np.zeros((2, 2)) 

410 

411 t_switch = np.inf if self.t_switch is None else self.t_switch 

412 

413 if rhs[1] - self.V_ref[0] <= 0 or t >= t_switch: 

414 self.A[0, 0] = -(self.Rs + self.R) / self.L 

415 

416 else: 

417 self.A[1, 1] = -1 / (self.C[0] * self.R) 

418 

419 me = self.dtype_u(self.init) 

420 me[:] = np.linalg.solve(np.eye(self.nvars) - factor * self.A, rhs) 

421 return me 

422 

423 def u_exact(self, t): 

424 r""" 

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

426 

427 Parameters 

428 ---------- 

429 t : float 

430 Time of the exact solution. 

431 

432 Returns 

433 ------- 

434 me : dtype_u 

435 The exact solution. 

436 """ 

437 assert t == 0, 'ERROR: u_exact only valid for t=0' 

438 

439 me = self.dtype_u(self.init) 

440 

441 me[0] = 0.0 # cL 

442 me[1] = self.alpha * self.V_ref[0] # vC 

443 

444 return me 

445 

446 

447class battery_implicit(battery): 

448 r""" 

449 Battery drain model with one capacitor, treated fully implicitly with Newton instead of the IMEX splitting. 

450 The method solve_system uses a fully-implicit computation. 

451 

452 Parameters 

453 ---------- 

454 ncapacitors : int, optional 

455 Number of capacitors in the circuit. 

456 Vs : float, optional 

457 Voltage at the voltage source :math:`V_s`. 

458 Rs : float, optional 

459 Resistance of the resistor :math:`R_s` at the voltage source. 

460 C : np.1darray, optional 

461 Capacitances of the capacitors. Length of array must equal to number of capacitors. 

462 R : float, optional 

463 Resistance for the load :math:`R_\ell`. 

464 L : float, optional 

465 Inductance of inductor :math:`L`. 

466 alpha : float, optional 

467 Factor greater than zero to describe the storage of the capacitor(s). 

468 V_ref : float, optional 

469 Reference value greater than zero for the battery to switch to the voltage source. 

470 newton_maxiter : int, optional 

471 Number of maximum iterations for the Newton solver. 

472 newton_tol : float, optional 

473 Tolerance for determination of the Newton solver. 

474 stop_at_nan : bool, optional 

475 Raise a ``ProblemError`` if Newton's method produces ``nan``, instead of only logging a warning. 

476 

477 Attributes 

478 ---------- 

479 work_counters : WorkCounter 

480 Counts different things, here: Number of Newton iterations is counted. 

481 """ 

482 

483 dtype_f = mesh 

484 

485 def __init__( 

486 self, 

487 ncapacitors=1, 

488 Vs=5.0, 

489 Rs=0.5, 

490 C=None, 

491 R=1.0, 

492 L=1.0, 

493 alpha=1.2, 

494 V_ref=None, 

495 newton_maxiter=100, 

496 newton_tol=1e-11, 

497 stop_at_nan=True, 

498 ): 

499 """Initialization routine""" 

500 super().__init__(ncapacitors, Vs, Rs, C, R, L, alpha, V_ref) 

501 self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True) 

502 

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

504 

505 def eval_f(self, u, t): 

506 """ 

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

508 

509 Parameters 

510 ---------- 

511 u : dtype_u 

512 Current values of the numerical solution. 

513 t : float 

514 Current time of the numerical solution is computed. 

515 

516 Returns 

517 ------- 

518 f : dtype_f 

519 The right-hand side of the problem. 

520 """ 

521 

522 f = self.dtype_f(self.init, val=0.0) 

523 non_f = np.zeros(2) 

524 

525 t_switch = np.inf if self.t_switch is None else self.t_switch 

526 

527 if u[1] - self.V_ref[0] <= 0 or t >= t_switch: 

528 self.A[0, 0] = -(self.Rs + self.R) / self.L 

529 non_f[0] = self.Vs / self.L 

530 

531 else: 

532 self.A[1, 1] = -1 / (self.C[0] * self.R) 

533 non_f[0] = 0 

534 

535 f[:] = self.A.dot(u) + non_f 

536 return f 

537 

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

539 """ 

540 Simple Newton solver. 

541 

542 Parameters 

543 ---------- 

544 rhs : dtype_f 

545 Right-hand side for the nonlinear system. 

546 factor : float 

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

548 u0 : dtype_u 

549 Initial guess for the iterative solver 

550 t : float 

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

552 

553 Returns 

554 ------- 

555 me : dtype_u 

556 The solution as mesh. 

557 """ 

558 

559 u = self.dtype_u(u0) 

560 non_f = np.zeros(2) 

561 self.A = np.zeros((2, 2)) 

562 

563 t_switch = np.inf if self.t_switch is None else self.t_switch 

564 

565 if rhs[1] - self.V_ref[0] <= 0 or t >= t_switch: 

566 self.A[0, 0] = -(self.Rs + self.R) / self.L 

567 non_f[0] = self.Vs / self.L 

568 

569 else: 

570 self.A[1, 1] = -1 / (self.C[0] * self.R) 

571 non_f[0] = 0 

572 

573 # start newton iteration 

574 n = 0 

575 res = 99 

576 while n < self.newton_maxiter: 

577 # form function g with g(u) = 0 

578 g = u - rhs - factor * (self.A.dot(u) + non_f) 

579 

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

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

582 

583 if res < self.newton_tol: 

584 break 

585 

586 # assemble dg 

587 dg = np.eye(self.nvars) - factor * self.A 

588 

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

590 u -= np.linalg.solve(dg, g) 

591 

592 # increase iteration count 

593 n += 1 

594 self.work_counters['newton']() 

595 

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

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

598 elif np.isnan(res): 

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

600 

601 if n == self.newton_maxiter: 

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

603 

604 me = self.dtype_u(self.init) 

605 me[:] = u[:] 

606 

607 return me