Coverage for physioblocks / library / model_components / velocity_law.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 a discretized Velocity Law Block. 

29 

30It defines two residual giving a dynamics on the velocity and acceleration. 

31""" 

32 

33from dataclasses import dataclass 

34 

35import numpy as np 

36from numpy.typing import NDArray 

37 

38from physioblocks.computing import ( 

39 Expression, 

40 ModelComponent, 

41 Quantity, 

42 diff, 

43 mid_alpha, 

44 mid_point, 

45) 

46from physioblocks.registers import register_type 

47from physioblocks.simulation import Time 

48 

49# Constant for the velocity law hht type name 

50VELOCITY_LAW_HHT_TYPE_ID = "velocity_law_hht" 

51 

52 

53@dataclass 

54@register_type(VELOCITY_LAW_HHT_TYPE_ID) 

55class VelocityLawHHTModelComponent(ModelComponent): 

56 r""" 

57 Velocity law HHT quantities and equation definitions 

58 

59 Implement the extended version of the **HHT time integration scheme**. 

60 Velocity and acceleration unknowns are introduced and solved for by 

61 

62 **Discretised Internal Equations:** 

63 

64 .. math:: 

65 

66 \frac{\dot{y}^{n+1}-\dot{y}^n}{\Delta t^n} 

67 - (\tfrac{1}{2}+\alpha)\, \ddot{y}^{n+1} 

68 - (\tfrac{1}{2}-\alpha)\, \ddot{y}^n = 0 

69 

70 .. math:: 

71 

72 \frac{y^{n+1}-y^n}{\Delta t^n} 

73 - \dot{y}^{n + \frac{1}{2}} 

74 - \frac{\alpha^2}{4}\Delta t^n\, (\ddot{y}^{n+1}-\ddot{y}^n) = 0 

75 """ 

76 

77 scheme_ts_hht: Quantity[np.float64] 

78 """:math:`\\alpha` the time shift scheme""" 

79 

80 accel: Quantity[np.float64] 

81 """:math:`\\ddot{y}` the acceleration""" 

82 

83 vel: Quantity[np.float64] 

84 """:math:`\\dot{y}` the velocity""" 

85 

86 disp: Quantity[np.float64] 

87 """:math:`y` the displacement""" 

88 

89 time: Time 

90 """:math:`t` the time""" 

91 

92 def velocity_law_residual(self) -> NDArray[np.float64]: 

93 """ 

94 Compute the velocity law residual 

95 

96 :return: the residual value 

97 :rtype: NDArray[np.float64] 

98 """ 

99 accel_mid_alpha = mid_alpha(self.accel, self.scheme_ts_hht.current) 

100 vel_mid_point = mid_point(self.vel) 

101 

102 return np.array( 

103 [ 

104 (self.time.inv_dt * diff(self.vel) - accel_mid_alpha), 

105 ( 

106 self.time.inv_dt * diff(self.disp) 

107 - vel_mid_point 

108 - np.power(self.scheme_ts_hht.current, 2) 

109 * self.time.dt 

110 / 4.0 

111 * diff(self.accel) 

112 ), 

113 ], 

114 ) 

115 

116 def velocity_law_residual_daccel(self) -> NDArray[np.float64]: 

117 """ 

118 Compute the velocity law residual partial derivative for ``accel`` 

119 

120 :return: the partial derivative for accel 

121 :rtype: NDArray[np.float64] 

122 """ 

123 return np.array( 

124 [ 

125 -(0.5 + self.scheme_ts_hht.current), 

126 -np.power(self.scheme_ts_hht.current, 2) * self.time.dt * 0.25, 

127 ], 

128 ) 

129 

130 def velocity_law_residual_dvel(self) -> NDArray[np.float64]: 

131 """ 

132 Compute the velocity law residual partial derivative for ``vel`` 

133 

134 :return: the partial derivative for vel 

135 :rtype: NDArray[np.float64] 

136 """ 

137 return np.array( 

138 [self.time.inv_dt, -0.5], 

139 ) 

140 

141 def velocity_law_residual_ddisp(self) -> NDArray[np.float64]: 

142 """ 

143 Compute the velocity law residual partial derivative for ``disp`` 

144 

145 :return: the partial derivative for disp 

146 :rtype: NDArray[np.float64] 

147 """ 

148 return np.array( 

149 [0.0, self.time.inv_dt], 

150 ) 

151 

152 

153# Define the residual expression of the velocity law and its partial derivatives 

154_velocity_law_hht_expression = Expression( 

155 2, 

156 VelocityLawHHTModelComponent.velocity_law_residual, 

157 { 

158 "accel": VelocityLawHHTModelComponent.velocity_law_residual_daccel, 

159 "vel": VelocityLawHHTModelComponent.velocity_law_residual_dvel, 

160 "disp": VelocityLawHHTModelComponent.velocity_law_residual_ddisp, 

161 }, 

162) 

163 

164VelocityLawHHTModelComponent.declares_internal_expression( 

165 "vel", _velocity_law_hht_expression, 1, 0 

166) 

167VelocityLawHHTModelComponent.declares_internal_expression( 

168 "accel", _velocity_law_hht_expression, 1, 1 

169)