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

1import numpy as np 

2 

3from pySDC.core.Errors import ParameterError, ProblemError 

4from pySDC.core.Problem import ptype, WorkCounter 

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

6 

7 

8class battery_n_capacitors(ptype): 

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 assertions_V_ref_2 = [ 

95 (alpha > V_ref[k] for k in range(n)), 

96 (V_ref[k] > 0 for k in range(n)), 

97 ] 

98 msg2 = "ERROR for V_ref: At least one of V_ref is less than zero and/or alpha!" 

99 assert all(assertions_V_ref_2), msg2 

100 

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

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

103 self._makeAttributeAndRegister( 

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

105 ) 

106 

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

108 self.A = self.switch_A[0] 

109 

110 self.t_switch = None 

111 self.nswitches = 0 

112 

113 def eval_f(self, u, t): 

114 r""" 

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

116 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: 

117 

118 :math:`C_1` supplies energy if: 

119 

120 .. math:: 

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

122 

123 :math:`C_2` supplies energy if: 

124 

125 .. math:: 

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

127 

128 :math:`C_3` supplies energy if: 

129 

130 .. math:: 

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

132 

133 :math:`V_s` supplies energy if: 

134 

135 .. math:: 

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

137 

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

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

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

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

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

143 

144 Parameters 

145 ---------- 

146 u : dtype_u 

147 Current values of the numerical solution. 

148 t : float 

149 Current time of the numerical solution is computed. 

150 

151 Returns 

152 ------- 

153 f : dtype_f 

154 The right-hand side of the problem. 

155 """ 

156 

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

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

159 

160 if self.t_switch is not None: 

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

162 

163 else: 

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

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

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

167 

168 if max_index == -1: 

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

170 

171 else: 

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

173 

174 return f 

175 

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

177 r""" 

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

179 

180 Parameters 

181 ---------- 

182 rhs : dtype_f 

183 Right-hand side for the linear system. 

184 factor : float 

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

186 u0 : dtype_u 

187 Initial guess for the iterative solver. 

188 t : float 

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

190 

191 Returns 

192 ------- 

193 me : dtype_u 

194 The solution as mesh. 

195 """ 

196 

197 if self.t_switch is not None: 

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

199 

200 else: 

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

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

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

204 if max_index == -1: 

205 self.A = self.switch_A[0] 

206 

207 else: 

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

209 

210 me = self.dtype_u(self.init) 

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

212 return me 

213 

214 def u_exact(self, t): 

215 r""" 

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

217 

218 Parameters 

219 ---------- 

220 t : float 

221 Time of the exact solution. 

222 

223 Returns 

224 ------- 

225 me : dtype_u 

226 The exact solution. 

227 """ 

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

229 

230 me = self.dtype_u(self.init) 

231 

232 me[0] = 0.0 # cL 

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

234 return me 

235 

236 def get_switching_info(self, u, t): 

237 """ 

238 Provides information about a discrete event for one subinterval. 

239 

240 Parameters 

241 ---------- 

242 u : dtype_u 

243 Current values of the numerical solution. 

244 t : float 

245 Current time of the numerical solution is computed. 

246 

247 Returns 

248 ------- 

249 switch_detected : bool 

250 Indicates if a switch is found or not. 

251 m_guess : int 

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

253 state_function : list 

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

255 """ 

256 

257 switch_detected = False 

258 m_guess = -100 

259 break_flag = False 

260 k_detected = 1 

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

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

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

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

265 if h_prev_node > 0 and h_curr_node <= 0: 

266 switch_detected = True 

267 m_guess = m - 1 

268 k_detected = k 

269 break_flag = True 

270 break 

271 

272 if break_flag: 

273 break 

274 

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

276 

277 return switch_detected, m_guess, state_function 

278 

279 def count_switches(self): 

280 """ 

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

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

283 """ 

284 

285 self.nswitches += 1 

286 

287 def get_problem_dict(self): 

288 """ 

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

290 """ 

291 

292 n = self.ncapacitors 

293 v = np.zeros(n + 1) 

294 v[0] = 1 

295 A, f = dict(), dict() 

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

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

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

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

300 return A, f 

301 

302 

303class battery(battery_n_capacitors): 

304 r""" 

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

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

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

308 the following equations: 

309 

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

311 

312 .. math:: 

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

314 

315 .. math:: 

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

317 

318 else: 

319 

320 .. math:: 

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

322 

323 .. math:: 

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

325 

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

327 

328 Note 

329 ---- 

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

331 """ 

332 

333 dtype_f = imex_mesh 

334 

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

336 """Initialization routine""" 

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

338 

339 def eval_f(self, u, t): 

340 """ 

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

342 

343 Parameters 

344 ---------- 

345 u : dtype_u 

346 Current values of the numerical solution. 

347 t : float 

348 Current time of the numerical solution is computed. 

349 

350 Returns 

351 ------- 

352 f : dtype_f 

353 The right-hand side of the problem. 

354 """ 

355 

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

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

358 

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

360 

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

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

363 

364 else: 

365 f.expl[0] = 0 

366 

367 return f 

368 

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

370 r""" 

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

372 

373 Parameters 

374 ---------- 

375 rhs : dtype_f 

376 Right-hand side for the linear system. 

377 factor : float 

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

379 u0 : dtype_u 

380 Initial guess for the iterative solver. 

381 t : float 

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

383 

384 Returns 

385 ------- 

386 me : dtype_u 

387 The solution as mesh. 

388 """ 

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

390 

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

392 

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

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

395 

396 else: 

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

398 

399 me = self.dtype_u(self.init) 

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

401 return me 

402 

403 def u_exact(self, t): 

404 r""" 

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

406 

407 Parameters 

408 ---------- 

409 t : float 

410 Time of the exact solution. 

411 

412 Returns 

413 ------- 

414 me : dtype_u 

415 The exact solution. 

416 """ 

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

418 

419 me = self.dtype_u(self.init) 

420 

421 me[0] = 0.0 # cL 

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

423 

424 return me 

425 

426 

427class battery_implicit(battery): 

428 r""" 

429 Example implementing the battery drain model as above. The method solve_system uses a fully-implicit computation. 

430 

431 Parameters 

432 ---------- 

433 ncapacitors : int, optional 

434 Number of capacitors in the circuit. 

435 Vs : float, optional 

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

437 Rs : float, optional 

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

439 C : np.1darray, optional 

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

441 R : float, optional 

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

443 L : float, optional 

444 Inductance of inductor :math:`L`. 

445 alpha : float, optional 

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

447 V_ref : float, optional 

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

449 newton_maxiter : int, optional 

450 Number of maximum iterations for the Newton solver. 

451 newton_tol : float, optional 

452 Tolerance for determination of the Newton solver. 

453 

454 Attributes 

455 ---------- 

456 work_counters : WorkCounter 

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

458 """ 

459 

460 dtype_f = mesh 

461 

462 def __init__( 

463 self, 

464 ncapacitors=1, 

465 Vs=5.0, 

466 Rs=0.5, 

467 C=None, 

468 R=1.0, 

469 L=1.0, 

470 alpha=1.2, 

471 V_ref=None, 

472 newton_maxiter=100, 

473 newton_tol=1e-11, 

474 ): 

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

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

477 

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

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. 

490 

491 Returns 

492 ------- 

493 f : dtype_f 

494 The right-hand side of the problem. 

495 """ 

496 

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

498 non_f = np.zeros(2) 

499 

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

501 

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

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

504 non_f[0] = self.Vs 

505 

506 else: 

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

508 non_f[0] = 0 

509 

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

511 return f 

512 

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

514 """ 

515 Simple Newton solver. 

516 

517 Parameters 

518 ---------- 

519 rhs : dtype_f 

520 Right-hand side for the nonlinear system. 

521 factor : float 

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

523 u0 : dtype_u 

524 Initial guess for the iterative solver 

525 t : float 

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

527 

528 Returns 

529 ------- 

530 me : dtype_u 

531 The solution as mesh. 

532 """ 

533 

534 u = self.dtype_u(u0) 

535 non_f = np.zeros(2) 

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

537 

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

539 

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

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

542 non_f[0] = self.Vs 

543 

544 else: 

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

546 non_f[0] = 0 

547 

548 # start newton iteration 

549 n = 0 

550 res = 99 

551 while n < self.newton_maxiter: 

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

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

554 

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

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

557 

558 if res < self.newton_tol: 

559 break 

560 

561 # assemble dg 

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

563 

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

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

566 

567 # increase iteration count 

568 n += 1 

569 self.work_counters['newton']() 

570 

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

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

573 elif np.isnan(res): 

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

575 

576 if n == self.newton_maxiter: 

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

578 

579 me = self.dtype_u(self.init) 

580 me[:] = u[:] 

581 

582 return me