Coverage for physioblocks / library / blocks / valves.py: 100%

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

28Module describing Valve Blocks 

29""" 

30 

31from dataclasses import dataclass 

32from typing import Any 

33 

34import numpy as np 

35 

36from physioblocks.computing import ( 

37 Block, 

38 Expression, 

39 Quantity, 

40 diff, 

41 mid_alpha, 

42 mid_point, 

43) 

44from physioblocks.registers import register_type 

45from physioblocks.simulation import Time 

46 

47# Constant for the valve rl block type id 

48VALVE_RL_BLOCK_ID = "valve_rl_block" 

49 

50 

51# Flux through the block 

52VALVE_RL_BLOCK_FLUX_VAR_ID = "flux" 

53 

54# Pressure at local node 1 

55VALVE_RL_BLOCK_PRESSURE_1_DOF_ID = "pressure_1" 

56 

57# Pressure at local node 2 

58VALVE_RL_BLOCK_PRESSURE_2_DOF_ID = "pressure_2" 

59 

60 

61@dataclass 

62@register_type(VALVE_RL_BLOCK_ID) 

63class ValveRLBlock(Block): 

64 r""" 

65 Describes valve RL block quantities and flux definitions 

66 

67 .. tikz:: Valve RC Scheme 

68 

69 \draw (1,2) to[short, i=$Q_2$] (1,1) to[short, -*] (1,0) node[right]{$P_2$}; 

70 \draw (1,2) to[short, L=$L$] (1,4); 

71 \draw (0,4) -- (2,4); 

72 \draw (0,4) to[D] (0,6); 

73 \draw (0,6) to[R=$R_{\text{back}}$] (0,8); 

74 \draw (2,8) to[D] (2,6); 

75 \draw (2,6) to[R=$R$] (2,4); 

76 \draw (0,8) -- (2,8); 

77 \draw (1,9) to[short, i=$Q$] (1,8); 

78 \draw (1,9) to[short, i=$Q_1$] (1,10) to [short, -*] (1,10.5) 

79 node[right]{$P_1$}; 

80 

81 **Node 1:** 

82 

83 :math:`Q_1 = - Q` 

84 

85 **Node 2:** 

86 

87 :math:`Q_2 = Q` 

88 

89 **Internal Equations:** 

90 

91 .. math:: 

92 

93 L\ \dot{Q} + P_2 - P_1 + 

94 \begin{cases} 

95 RQ \text{ if Q $>$ 0 } \\ 

96 R_{\text{back}}Q \text{ else } 

97 \end{cases} = 0 

98 

99 **Discretisation:** 

100 

101 

102 .. math:: Q_1^{{n + \frac{1}{2}}} = - Q^{{n + \frac{1}{2}}} 

103 

104 .. math:: Q_2^{{n + \frac{1}{2}}} = Q^{{n + \frac{1}{2}}} 

105 

106 .. math:: 

107 

108 L\ \frac{Q^{n + 1} - Q^{n}}{\Delta t^n} 

109 + P_2^{n + \frac{1}{2}} - P_1^{n + \frac{1}{2}} 

110 + \begin{cases} 

111 RQ^{{n + \theta}} \text{ if } Q^{{n + \theta}} > 0 \\ 

112 R_{\text{back}}Q^{{n + \theta}} \text{ else } 

113 \end{cases} = 0 

114 

115 """ 

116 

117 flux: Quantity[np.float64] 

118 """Flux traversing the block""" 

119 

120 pressure_1: Quantity[np.float64] 

121 """Pressure quantity at local node 1""" 

122 

123 pressure_2: Quantity[np.float64] 

124 """Pressure quantity at local node 2""" 

125 

126 inductance: Quantity[np.float64] 

127 """Inductance quantity""" 

128 

129 conductance: Quantity[np.float64] 

130 """Conductance quantity for positive flux """ 

131 

132 backward_conductance: Quantity[np.float64] 

133 """Conductance quantity for negative flux""" 

134 

135 scheme_ts_flux: Quantity[np.float64] 

136 """Scheme time shift for flux""" 

137 

138 time: Time 

139 """Simulation time""" 

140 

141 def initialize(self) -> None: 

142 self.k_0 = ( 

143 2.0 

144 * (self.conductance.current * self.backward_conductance.current) 

145 / (self.conductance.current + self.backward_conductance.current) 

146 ) 

147 

148 def flux_residual(self) -> Any: 

149 """ 

150 Compute the residual giving the dynamics on the flux in the valve. 

151 

152 :return: the residual value 

153 :rtype: np.float64 

154 """ 

155 

156 q_mid_alpha = mid_alpha(self.flux, self.scheme_ts_flux.current) 

157 pressure_1_mid_point = mid_point(self.pressure_1) 

158 pressure_2_mid_point = mid_point(self.pressure_2) 

159 

160 conductance = self.conductance.current 

161 if q_mid_alpha < 0: 

162 conductance = self.backward_conductance.current 

163 

164 return ( 

165 self.inductance.current * self.time.inv_dt * diff(self.flux) 

166 - pressure_1_mid_point 

167 + pressure_2_mid_point 

168 + q_mid_alpha / conductance 

169 ) 

170 

171 def flux_residual_dflux(self) -> Any: 

172 """ 

173 Compute the residual partial derivative for ``flux`` 

174 

175 :return: the residual partial derivative for ``flux`` 

176 :rtype: np.float64 

177 """ 

178 q_mid_alpha = mid_alpha(self.flux, self.scheme_ts_flux.current) 

179 conductance = self.conductance.current 

180 if q_mid_alpha < 0: 

181 conductance = self.backward_conductance.current 

182 elif q_mid_alpha == 0: 

183 conductance = self.k_0 

184 

185 return ( 

186 self.inductance.current * self.time.inv_dt 

187 + (0.5 + self.scheme_ts_flux.current) / conductance 

188 ) 

189 

190 def flux_residual_dp1(self) -> Any: 

191 """ 

192 Compute the residual partial derivative for ``pressure_1`` 

193 

194 :return: the residual partial derivative for ``pressure_1`` 

195 :rtype: np.float64 

196 """ 

197 

198 return -0.5 

199 

200 def flux_residual_dp2(self) -> Any: 

201 """ 

202 Compute the residual partial derivative for ``pressure_2`` 

203 

204 :return: the residual partial derivative for ``pressure_2`` 

205 :rtype: np.float64 

206 """ 

207 

208 return 0.5 

209 

210 def flux_1(self) -> Any: 

211 """ 

212 Compute the block flux at node 1 

213 

214 

215 :return: the block flux at node 1 

216 :rtype: np.float64 

217 """ 

218 return -mid_point(self.flux) 

219 

220 def dflux_1_dflux(self) -> Any: 

221 """ 

222 Compute the block flux at node 1 partial derivative for ``flux`` 

223 

224 

225 :return: the block flux at node 1 partial derivative for ``flux`` 

226 :rtype: np.float64 

227 """ 

228 return -0.5 

229 

230 def flux_2(self) -> Any: 

231 """ 

232 Compute the block flux at node 2 

233 

234 :return: the block flux at node 2 

235 :rtype: np.float64 

236 """ 

237 return mid_point(self.flux) 

238 

239 def dflux_2_dflux(self) -> Any: 

240 """ 

241 Compute the block flux at node 2 partial derivative for ``flux`` 

242 

243 :return: the block flux at node 2 partial derivative for ``flux`` 

244 :rtype: np.float64 

245 """ 

246 return 0.5 

247 

248 

249# define the Valve RL internal variable residual expression 

250# giving the dynamics on the flux in the valve 

251_valve_rl_flux_residual_expr = Expression( 

252 1, 

253 ValveRLBlock.flux_residual, 

254 { 

255 VALVE_RL_BLOCK_FLUX_VAR_ID: ValveRLBlock.flux_residual_dflux, 

256 VALVE_RL_BLOCK_PRESSURE_1_DOF_ID: ValveRLBlock.flux_residual_dp1, 

257 VALVE_RL_BLOCK_PRESSURE_2_DOF_ID: ValveRLBlock.flux_residual_dp2, 

258 }, 

259) 

260 

261# define the ValveRL flux expression at node 1. 

262_valve_rl_flux_1_expr = Expression( 

263 1, ValveRLBlock.flux_1, {VALVE_RL_BLOCK_FLUX_VAR_ID: ValveRLBlock.dflux_1_dflux} 

264) 

265 

266# define the ValveRL flux expression at node 2. 

267_valve_rl_flux_2_expr = Expression( 

268 1, 

269 ValveRLBlock.flux_2, 

270 {VALVE_RL_BLOCK_FLUX_VAR_ID: ValveRLBlock.dflux_2_dflux}, 

271) 

272 

273ValveRLBlock.declares_internal_expression( 

274 VALVE_RL_BLOCK_FLUX_VAR_ID, _valve_rl_flux_residual_expr 

275) 

276ValveRLBlock.declares_flux_expression( 

277 1, VALVE_RL_BLOCK_PRESSURE_1_DOF_ID, _valve_rl_flux_1_expr 

278) 

279ValveRLBlock.declares_flux_expression( 

280 2, VALVE_RL_BLOCK_PRESSURE_2_DOF_ID, _valve_rl_flux_2_expr 

281)