Coverage for physioblocks / library / blocks / capacitances.py: 99%

1# SPDX-FileCopyrightText: Copyright INRIA 

2# 

3# SPDX-License-Identifier: LGPL-3.0-only 

4# 

5# Copyright INRIA 

6# 

7# This file is part of PhysioBlocks, a library mostly developed by the 

8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA. 

9# 

10# Authors: 

11# - Colin Drieu 

12# - Dominique Chapelle 

13# - François Kimmig 

14# - Philippe Moireau 

15# 

16# PhysioBlocks is free software: you can redistribute it and/or modify it under the 

17# terms of the GNU Lesser General Public License as published by the Free Software 

18# Foundation, version 3 of the License. 

19# 

20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY 

21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 

22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along with 

25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>. 

26 

27"""Module describing capacitance blocks""" 

28 

29from dataclasses import dataclass 

30from typing import Any 

31 

32import numpy as np 

33 

34from physioblocks.computing import Block, Expression, Quantity, diff, mid_point 

35from physioblocks.registers import register_type 

36from physioblocks.simulation import Time 

37 

38# C BLOCK Definition 

39 

40 

41# Constant for the c block type id 

42C_BLOCK_TYPE_ID = "c_block" 

43 

44# Constant for the c block dof id 

45C_BLOCK_PRESSURE_DOF_ID = "pressure" 

46 

47 

48@register_type(C_BLOCK_TYPE_ID) 

49@dataclass 

50class CBlock(Block): 

51 r""" 

52 C Block quantities and flux definitions 

53 

54 .. tikz:: C Block scheme 

55 

56 \draw (0,1) to[short, -*, i=$Q$] (0,3) node[right]{$P$} 

57 (0,1) to[C=$C$] (0,0) node[ground]{}; 

58 

59 **Node 1:** 

60 

61 :math:`Q = - C\dot{P}` 

62 

63 **Discretisation:** 

64 

65 .. math:: 

66 

67 Q^{n + \frac{1}{2}} = - C\ \frac{P^{n + 1} - P^{n}}{\Delta t^n} 

68 

69 """ 

70 

71 pressure: Quantity[np.float64] 

72 """Pressure quantity""" 

73 

74 capacitance: Quantity[np.float64] 

75 """Block capacitance quantity""" 

76 

77 time: Time 

78 """Simulation time""" 

79 

80 def flux(self) -> Any: 

81 """ 

82 Compute the flux at local node 1 

83 

84 :return: the flux 

85 :rtype: np.float64 

86 """ 

87 return -self.capacitance.current * diff(self.pressure) * self.time.inv_dt 

88 

89 def dflux_dpressure(self) -> Any: 

90 """ 

91 Compute the flux at local node 1 partial derivative for pressure 

92 

93 :return: the flux derivative 

94 :rtype: np.float64 

95 """ 

96 return -self.capacitance.current * self.time.inv_dt 

97 

98 

99_c_block_flux_expression = Expression( 

100 1, CBlock.flux, {C_BLOCK_PRESSURE_DOF_ID: CBlock.dflux_dpressure} 

101) 

102CBlock.declares_flux_expression(1, C_BLOCK_PRESSURE_DOF_ID, _c_block_flux_expression) 

103 

104 

105# RC BLOCK Definition 

106 

107# Constant for the rc block type id 

108RC_BLOCK_TYPE_ID = "rc_block" 

109 

110# Constant for the pressure in local id 

111RC_BLOCK_PRESSURE_1_DOF_ID = "pressure_1" 

112 

113# Constant for the pressure out local id 

114RC_BLOCK_PRESSURE_2_DOF_ID = "pressure_2" 

115 

116 

117@register_type(RC_BLOCK_TYPE_ID) 

118@dataclass 

119class RCBlock(Block): 

120 r""" 

121 RC Block quantities and fluxes definitions 

122 

123 .. tikz:: RC Block scheme 

124 

125 \draw (0,3) node[below]{$P_{1}$} to[short, *-] (1,3); 

126 \draw (2, 3) to[short, i=$Q_{1}$] (1, 3); 

127 \draw (2, 3) to[R=$R$] (4,3); 

128 \draw (4,3) to[C=$C$] (4,0) node[ground]{}; 

129 \draw (4,3) to[short, -*, i=$Q_{2}$] (6,3) node[below]{$P_{2}$}; 

130 

131 **Node 1:** 

132 

133 :math:`Q_1 = \frac{P_2 - P_1}{R}` 

134 

135 **Node 2:** 

136 

137 :math:`Q_2 = \frac{P_1 - P_2}{R} - C\dot{P_2}` 

138 

139 **Discretisation:** 

140 

141 .. math:: 

142 

143 Q_1^{n + \frac{1}{2}} = \frac{P_2^{n + \frac{1}{2}} - P_1^{n + \frac{1}{2}}}{R} 

144 

145 .. math:: 

146 

147 Q_2^{n + \frac{1}{2}} = \frac{P_1^{n + \frac{1}{2}} - P_2^{n + \frac{1}{2}}}{R} 

148 - C\ \frac{P_2^{n + 1} - P_2^{n}}{\Delta t^n} 

149 

150 """ 

151 

152 pressure_1: Quantity[np.float64] 

153 """Pressure at local node 1 of the block""" 

154 

155 pressure_2: Quantity[np.float64] 

156 """Pressure at local node 2 of the block""" 

157 

158 resistance: Quantity[np.float64] 

159 """Resistance value of the block""" 

160 

161 capacitance: Quantity[np.float64] 

162 """Capacitor value of the block""" 

163 

164 time: Time 

165 """The simulation time""" 

166 

167 def flux_1(self) -> Any: 

168 """ 

169 Computes the outlet flux at local node 1. 

170 

171 :return: the flux value for current block values 

172 :rtype: np.float64 

173 """ 

174 pressure_1 = mid_point(self.pressure_1) 

175 pressure_2 = mid_point(self.pressure_2) 

176 return (pressure_2 - pressure_1) / self.resistance.current 

177 

178 def dflux_1_dpressure_1(self) -> Any: 

179 """ 

180 Computes the outlet flux at node 1 derivative for pressure_1. 

181 

182 :return: the flux derivative for pressure_1 

183 :rtype: np.float64 

184 """ 

185 return -0.5 / self.resistance.current 

186 

187 def dflux_1_dpressure_2(self) -> Any: 

188 """ 

189 Computes the outlet flux at node 1 derivative for pressure_2. 

190 

191 :return: the flux derivative for pressure_2 

192 :rtype: np.float64 

193 """ 

194 return 0.5 / self.resistance.current 

195 

196 def flux_2(self) -> Any: 

197 """ 

198 Computes the outlet flux at node 2. 

199 

200 :return: the flux value for current block values 

201 :rtype: np.float64 

202 """ 

203 pressure_1 = mid_point(self.pressure_1) 

204 pressure_2 = mid_point(self.pressure_2) 

205 dpressure_2 = diff(self.pressure_2) 

206 return ( 

207 (pressure_1 - pressure_2) / self.resistance.current 

208 - self.capacitance.current * self.time.inv_dt * dpressure_2 

209 ) 

210 

211 def dflux_2_dpressure_1(self) -> Any: 

212 """ 

213 Computes the outlet flux at node 2 derivative for pressure_1. 

214 

215 :return: the flux derivative for pressure_1 

216 :rtype: np.float64 

217 """ 

218 return 0.5 / self.resistance.current 

219 

220 def dflux_2_dpressure_2(self) -> Any: 

221 """ 

222 Computes the outlet flux at node 2 derivative for pressure_2. 

223 

224 :return: the flux derivative for pressure_2 

225 :rtype: np.float64 

226 """ 

227 return ( 

228 -0.5 / self.resistance.current - self.capacitance.current * self.time.inv_dt 

229 ) 

230 

231 

232# Define the flux expression going in the input node for rc_block 

233_rc_block_flux_1_expr = Expression( 

234 1, 

235 RCBlock.flux_1, 

236 { 

237 RC_BLOCK_PRESSURE_1_DOF_ID: RCBlock.dflux_1_dpressure_1, 

238 RC_BLOCK_PRESSURE_2_DOF_ID: RCBlock.dflux_1_dpressure_2, 

239 }, 

240) 

241 

242# Define the flux expression going in the output node for rc_block 

243_rc_block_flux_2_expr = Expression( 

244 1, 

245 RCBlock.flux_2, 

246 { 

247 RC_BLOCK_PRESSURE_1_DOF_ID: RCBlock.dflux_2_dpressure_1, 

248 RC_BLOCK_PRESSURE_2_DOF_ID: RCBlock.dflux_2_dpressure_2, 

249 }, 

250) 

251 

252RCBlock.declares_flux_expression(1, RC_BLOCK_PRESSURE_1_DOF_ID, _rc_block_flux_1_expr) 

253RCBlock.declares_flux_expression(2, RC_BLOCK_PRESSURE_2_DOF_ID, _rc_block_flux_2_expr) 

254 

255# RCR BLOCK Definition 

256 

257# Constant for the rcr block type id 

258RCR_BLOCK_TYPE_ID = "rcr_block" 

259 

260# Constant for the rcr block volume saved quantity 

261RCR_BLOCK_VOLUME_OUTPUT_ID = "volume_stored" 

262 

263# Constant for the rcr block pressure at the mid point 

264RCR_BLOCK_PRESSURE_MID_ID = "pressure_mid" 

265 

266# Constant for the rcr block pressure at local node 1 

267RCR_BLOCK_PRESSURE_1_ID = "pressure_1" 

268 

269# Constant for the rcr block pressure at local node 2 

270RCR_BLOCK_PRESSURE_2_ID = "pressure_2" 

271 

272 

273@register_type(RCR_BLOCK_TYPE_ID) 

274@dataclass 

275class RCRBlock(Block): 

276 r""" 

277 RCR Block quantities and flux definitions 

278 

279 .. tikz:: 

280 

281 \draw (-5,3) node[above]{$P_1$} to[short, *-] (-4,3); 

282 \draw (-3, 3) to[short, i=$Q_1$] (-4,3); 

283 \draw (-3, 3) to [R=$R_1$] (0, 3); 

284 \draw (0, 3) node[above]{$P_{mid}$} to [short, *-, R=$R_2$] (3, 3); 

285 \draw (3, 3) to[short, i=$Q_2$] (4, 3); 

286 \draw (4,3) to[short, -*] (5,3) node[above]{$P_2$}; 

287 \draw (0,3) to[C=$C$] (0,0) node[ground]{}; 

288 

289 **Node 1:** 

290 

291 :math:`Q_1 = \frac{P_{mid} - P_1}{R_1}` 

292 

293 **Node 2:** 

294 

295 :math:`Q_2 = \frac{P_{mid} - P_2}{R_2}` 

296 

297 **Internal equation:** 

298 

299 .. math:: 

300 

301 \frac{P_1 - P_{mid}}{R_1} + \frac{P_2 - P_{mid}}{R_2} - C\dot{P}_{mid} = 0 

302 

303 **Discretisation:** 

304 

305 .. math:: 

306 

307 Q_1^{n + \frac{1}{2}} = \frac{P_{mid}^{n + \frac{1}{2}} 

308 - P_1^{n + \frac{1}{2}}}{R_1} 

309 

310 .. math:: 

311 

312 Q_2^{n + \frac{1}{2}} = \frac{P_{mid}^{n + \frac{1}{2}} 

313 - P_2^{n + \frac{1}{2}}}{R_2} 

314 

315 .. math:: 

316 

317 \frac{P_1^{n + \frac{1}{2}} - P_{mid}^{n + \frac{1}{2}}}{R_1} 

318 + \frac{P_2^{n + \frac{1}{2}} - P_{mid}^{n + \frac{1}{2}}}{R_2} 

319 - C\ \frac{P_{mid}^{n + 1} - P_{mid}^{n}}{\Delta t^n} = 0 

320 

321 """ 

322 

323 pressure_1: Quantity[np.float64] 

324 """Pressure at the input of the block""" 

325 

326 pressure_2: Quantity[np.float64] 

327 """Pressure at the output of the block""" 

328 

329 pressure_mid: Quantity[np.float64] 

330 """Pressure at the output of the block""" 

331 

332 resistance_1: Quantity[np.float64] 

333 """Resistance in value of the block""" 

334 

335 capacitance: Quantity[np.float64] 

336 """Capacitor value of the block""" 

337 

338 resistance_2: Quantity[np.float64] 

339 """Resistance out value of the block""" 

340 

341 time: Time 

342 """The simulation time""" 

343 

344 def flux_1(self) -> Any: 

345 """ 

346 Computes the outlet flux at node 1. 

347 

348 :return: the flux value 

349 :rtype: np.float64 

350 """ 

351 

352 pressure_1_discr = mid_point(self.pressure_1) 

353 pressure_mid_discr = mid_point(self.pressure_mid) 

354 

355 return (pressure_mid_discr - pressure_1_discr) / self.resistance_1.current 

356 

357 def dflux_1_dp_1(self) -> Any: 

358 """ 

359 Computes the outlet flux at node 1 derivative for pressure_1. 

360 

361 :return: flux derivative for pressure_1 

362 :rtype: np.float64 

363 """ 

364 

365 return -0.5 / self.resistance_1.current 

366 

367 def dflux_1_dp_mid(self) -> Any: 

368 """ 

369 Computes the outlet flux at node 1 derivative for pressure_mid. 

370 

371 :return: flux derivative for pressure_2 

372 :rtype: np.float64 

373 """ 

374 

375 return 0.5 / self.resistance_1.current 

376 

377 def flux_2(self) -> Any: 

378 """ 

379 Computes the flux at node 2. 

380 

381 :return: the flux value 

382 :rtype: np.float64 

383 """ 

384 pressure_2_discr = mid_point(self.pressure_2) 

385 pressure_mid_discr = mid_point(self.pressure_mid) 

386 

387 return (pressure_mid_discr - pressure_2_discr) / self.resistance_2.current 

388 

389 def dflux_2_dp_mid(self) -> Any: 

390 """ 

391 Computes the outlet flux at node 2 derivative for pressure_mid. 

392 

393 :return: flux derivative for pressure_mid 

394 :rtype: np.float64 

395 """ 

396 

397 return 0.5 / self.resistance_2.current 

398 

399 def dflux_2_dp_2(self) -> Any: 

400 """ 

401 Computes the outlet flux at node 2 derivative for pressure_2. 

402 

403 :return: flux derivative for pressure_2 

404 :rtype: np.float64 

405 """ 

406 

407 return -0.5 / self.resistance_2.current 

408 

409 def pressure_mid_residual(self) -> Any: 

410 """ 

411 Compute the residual representing dynamics of the mid node pressure. 

412 

413 :return: the residual value 

414 :rtype: np.float64 

415 """ 

416 

417 p_1_discr = mid_point(self.pressure_1) 

418 p_mid_discr = mid_point(self.pressure_mid) 

419 p_2_discr = mid_point(self.pressure_2) 

420 

421 return ( 

422 +(p_1_discr - p_mid_discr) / self.resistance_1.current 

423 + (p_2_discr - p_mid_discr) / self.resistance_2.current 

424 - self.capacitance.current * self.time.inv_dt * diff(self.pressure_mid) 

425 ) 

426 

427 def pressure_mid_residual_dp_1(self) -> Any: 

428 """ 

429 Compute the residual derivative for pressure_1 

430 

431 :return: the residual derivative for pressure_1 

432 :rtype: np.float64 

433 """ 

434 

435 return 0.5 / self.resistance_1.current 

436 

437 def pressure_mid_residual_dp_2(self) -> Any: 

438 """ 

439 Compute the residual derivative for pressure_2 

440 

441 :return: the residual derivative for pressure_2 

442 :rtype: np.float64 

443 """ 

444 

445 return 0.5 / self.resistance_2.current 

446 

447 def pressure_mid_residual_dp_mid(self) -> Any: 

448 """ 

449 Compute the residual derivative for pressure_mid 

450 

451 :return: the residual derivative for pressure_mid 

452 :rtype: np.float64 

453 """ 

454 

455 return ( 

456 -self.capacitance.current * self.time.inv_dt 

457 - 0.5 / self.resistance_1.current 

458 - 0.5 / self.resistance_2.current 

459 ) 

460 

461 def compute_volume_stored(self) -> Any: 

462 """ 

463 Computes volume stored in the capacitance. 

464 

465 :return: volume stored in the capacitance 

466 :rtype: np.float64 

467 """ 

468 

469 return self.capacitance.current * self.pressure_mid.current 

470 

471 

472# Define the flux expression going in node 1 for rcr block 

473_rcr_block_flux_1_expr = Expression( 

474 1, 

475 RCRBlock.flux_1, 

476 { 

477 RCR_BLOCK_PRESSURE_1_ID: RCRBlock.dflux_1_dp_1, 

478 RCR_BLOCK_PRESSURE_MID_ID: RCRBlock.dflux_1_dp_mid, 

479 }, 

480) 

481 

482 

483# Define the flux expression going in node 2 for rcr block 

484_rcr_block_flux_2_expr = Expression( 

485 1, 

486 RCRBlock.flux_2, 

487 { 

488 RCR_BLOCK_PRESSURE_MID_ID: RCRBlock.dflux_2_dp_mid, 

489 RCR_BLOCK_PRESSURE_2_ID: RCRBlock.dflux_2_dp_2, 

490 }, 

491) 

492 

493# Define the residual expression giving the pressure at the mid node 

494_rcr_block_pressure_mid_residual_expr = Expression( 

495 1, 

496 RCRBlock.pressure_mid_residual, 

497 { 

498 RCR_BLOCK_PRESSURE_1_ID: RCRBlock.pressure_mid_residual_dp_1, 

499 RCR_BLOCK_PRESSURE_MID_ID: RCRBlock.pressure_mid_residual_dp_mid, 

500 RCR_BLOCK_PRESSURE_2_ID: RCRBlock.pressure_mid_residual_dp_2, 

501 }, 

502) 

503 

504# Derfine the stored volume saved quantity expression. 

505_rcr_volume_stored_expr = Expression(1, RCRBlock.compute_volume_stored) 

506 

507 

508RCRBlock.declares_internal_expression( 

509 RCR_BLOCK_PRESSURE_MID_ID, _rcr_block_pressure_mid_residual_expr 

510) 

511 

512RCRBlock.declares_flux_expression(1, RCR_BLOCK_PRESSURE_1_ID, _rcr_block_flux_1_expr) 

513RCRBlock.declares_flux_expression(2, RCR_BLOCK_PRESSURE_2_ID, _rcr_block_flux_2_expr) 

514RCRBlock.declares_saved_quantity_expression( 

515 RCR_BLOCK_VOLUME_OUTPUT_ID, _rcr_volume_stored_expr 

516)