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

## 166 statements

, created at 2024-09-19 09:13 +0000

1import numpy as np

3from pySDC.core.errors import ParameterError, ProblemError

4from pySDC.core.problem import Problem, WorkCounter

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

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)

17 .. math::

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

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

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.

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.

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

60 dtype_u = mesh

61 dtype_f = imex_mesh

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

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

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

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

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 )

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

108 self.A = self.switch_A[0]

110 self.t_switch = None

111 self.nswitches = 0

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:

118 :math:C_1 supplies energy if:

120 .. math::

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

123 :math:C_2 supplies energy if:

125 .. math::

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

128 :math:C_3 supplies energy if:

130 .. math::

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

133 :math:V_s supplies energy if:

135 .. math::

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

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.

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.

151 Returns

152 -------

153 f : dtype_f

154 The right-hand side of the problem.

155 """

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

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

160 if self.t_switch is not None:

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

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] else -1 for k in range(len(switch))])

168 if max_index == -1:

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

171 else:

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

174 return f

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

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

191 Returns

192 -------

193 me : dtype_u

194 The solution as mesh.

195 """

197 if self.t_switch is not None:

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

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] else -1 for k in range(len(switch))])

204 if max_index == -1:

205 self.A = self.switch_A[0]

207 else:

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

210 me = self.dtype_u(self.init)

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

212 return me

214 def u_exact(self, t):

215 r"""

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

218 Parameters

219 ----------

220 t : float

221 Time of the exact solution.

223 Returns

224 -------

225 me : dtype_u

226 The exact solution.

227 """

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

230 me = self.dtype_u(self.init)

232 me[0] = 0.0 # cL

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

234 return me

236 def get_switching_info(self, u, t):

237 """

238 Provides information about a discrete event for one subinterval.

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.

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

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

272 if break_flag:

273 break

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

277 return switch_detected, m_guess, state_function

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

285 self.nswitches += 1

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

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

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:

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

312 .. math::

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

315 .. math::

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

318 else:

320 .. math::

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

323 .. math::

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

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

328 Note

329 ----

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

331 """

333 dtype_f = imex_mesh

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)

339 def eval_f(self, u, t):

340 """

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

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.

350 Returns

351 -------

352 f : dtype_f

353 The right-hand side of the problem.

354 """

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

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

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

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

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

364 else:

365 f.expl[0] = 0

367 return f

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

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

384 Returns

385 -------

386 me : dtype_u

387 The solution as mesh.

388 """

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

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

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

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

396 else:

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

399 me = self.dtype_u(self.init)

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

401 return me

403 def u_exact(self, t):

404 r"""

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

407 Parameters

408 ----------

409 t : float

410 Time of the exact solution.

412 Returns

413 -------

414 me : dtype_u

415 The exact solution.

416 """

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

419 me = self.dtype_u(self.init)

421 me[0] = 0.0 # cL

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

424 return me

427class battery_implicit(battery):

428 r"""

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

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.

454 Attributes

455 ----------

456 work_counters : WorkCounter

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

458 """

460 dtype_f = mesh

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)

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

480 def eval_f(self, u, t):

481 """

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

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.

491 Returns

492 -------

493 f : dtype_f

494 The right-hand side of the problem.

495 """

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

498 non_f = np.zeros(2)

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

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

506 else:

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

508 non_f[0] = 0

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

511 return f

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

514 """

515 Simple Newton solver.

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

528 Returns

529 -------

530 me : dtype_u

531 The solution as mesh.

532 """

534 u = self.dtype_u(u0)

535 non_f = np.zeros(2)

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

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

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

544 else:

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

546 non_f[0] = 0

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)

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

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

558 if res < self.newton_tol:

559 break

561 # assemble dg

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

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

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

567 # increase iteration count

568 n += 1

569 self.work_counters['newton']()

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)

576 if n == self.newton_maxiter:

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

579 me = self.dtype_u(self.init)

580 me[:] = u[:]

582 return me