Coverage for pySDC/projects/DAE/problems/transistor_amplifier.py: 97%

1import warnings 

2import numpy as np 

3from scipy.interpolate import interp1d 

4 

5from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae 

6from pySDC.implementations.datatype_classes.mesh import mesh 

7 

8 

9# Helper function 

10def _transistor(u_in): 

11 return 1e-6 * (np.exp(u_in / 0.026) - 1) 

12 

13 

14class one_transistor_amplifier(ptype_dae): 

15 r""" 

16 The one transistor amplifier example from pg. 377 in [1]_. The problem is an index-1 differential-algebraic equation 

17 (DAE) having the equations 

18 

19 .. math:: 

20 \frac{U_e (t)}{R_0} - \frac{U_1 (t)}{R_0} + C_1 (\frac{d U_2 (t)}{dt} - \frac{d U_1 (t)}{dt}) = 0, 

21 

22 .. math:: 

23 \frac{U_b}{R_2} - U_2 (t) (\frac{1}{R_1} + \frac{1}{R_2}) + C_1 (\frac{d U_1 (t)}{dt} - \frac{d U_2 (t)}{dt}) - 0.01 f(U_2 (t) - U_3 (t)) = 0, 

24 

25 .. math:: 

26 f(U_2 (t) - U_3 (t)) - \frac{U_3 (t)}{R_3} - C_2 \frac{d U_3 (t)}{dt} = 0, 

27 

28 .. math:: 

29 \frac{U_b}{R_4} - \frac{U_4 (t)}{R_4} + C_3 (\frac{d U_5 (t)}{dt} - \frac{d U_4 (t)}{dt}) - 0.99 f(U_2 (t) - U_3 (t)) = 0, 

30 

31 .. math:: 

32 -\frac{U_5 (t)}{R_5} + C_3 (\frac{d U_4 (t)}{dt} - \frac{d U_5 (t)}{dt}) = 0, 

33 

34 with 

35 

36 .. math:: 

37 f(U(t)) = 10^{-6} (exp(\frac{U (t)}{0.026}) - 1). 

38 

39 The initial signal :math:`U_e (t)` is defined as 

40 

41 .. math:: 

42 U_e (t) = 0.4 \sin(200 \pi t). 

43 

44 Constants are fixed as :math:`U_b = 6`, :math:`R_0 = 1000`, :math:`R_k = 9000` for :math:`k=1,..,5`, 

45 `C_j = j \cdot 10^{-6}` for :math:`j=1,2,3`.They are also defined in the method `eval_f`. 

46 

47 Parameters 

48 ---------- 

49 nvars : int 

50 Number of unknowns of the system of DAEs. 

51 newton_tol : float 

52 Tolerance for Newton solver. 

53 

54 Attributes 

55 ---------- 

56 t_end: float 

57 The end time at which the reference solution is determined. 

58 

59 References 

60 ---------- 

61 .. [1] E. Hairer, G. Wanner. Solving ordinary differential equations II: Stiff and differential-algebraic problems. 

62 Springer (2009). 

63 """ 

64 

65 dtype_u = mesh 

66 dtype_f = mesh 

67 

68 def __init__(self, newton_tol): 

69 super().__init__(nvars=5, newton_tol=newton_tol) 

70 # load reference solution 

71 # data file must be generated and stored under misc/data and self.t_end = t[-1] 

72 # data = np.load(r'pySDC/projects/DAE/misc/data/one_trans_amp.npy') 

73 # x = data[:, 0] 

74 # # The last column contains the input signal 

75 # y = data[:, 1:-1] 

76 # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') 

77 self.t_end = 0.0 

78 

79 def eval_f(self, u, du, t): 

80 r""" 

81 Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. 

82 

83 Parameters 

84 ---------- 

85 u : dtype_u 

86 Current values of the numerical solution at time t. 

87 du : dtype_u 

88 Current values of the derivative of the numerical solution at time t. 

89 t : float 

90 Current time of the numerical solution. 

91 

92 Returns 

93 ------- 

94 f : dtype_f 

95 Current value of the right-hand side of f (which includes five components). 

96 """ 

97 u_b = 6.0 

98 u_e = 0.4 * np.sin(200 * np.pi * t) 

99 alpha = 0.99 

100 r_0 = 1000 

101 r_k = 9000 

102 c_1, c_2, c_3 = 1e-6, 2e-6, 3e-6 

103 f = self.dtype_f(self.init) 

104 f[:] = ( 

105 (u_e - u[0]) / r_0 + c_1 * (du[1] - du[0]), 

106 (u_b - u[1]) / r_k - u[1] / r_k + c_1 * (du[0] - du[1]) - (1 - alpha) * _transistor(u[1] - u[2]), 

107 _transistor(u[1] - u[2]) - u[2] / r_k - c_2 * du[2], 

108 (u_b - u[3]) / r_k + c_3 * (du[4] - du[3]) - alpha * _transistor(u[1] - u[2]), 

109 -u[4] / r_k + c_3 * (du[3] - du[4]), 

110 ) 

111 self.work_counters['rhs']() 

112 return f 

113 

114 def u_exact(self, t): 

115 """ 

116 Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical 

117 reference solution. 

118 

119 Parameters 

120 ---------- 

121 t : float 

122 The time of the reference solution. 

123 

124 Returns 

125 ------- 

126 me : dtype_u 

127 The reference solution as mesh object containing five components and fixed initial conditions. 

128 """ 

129 me = self.dtype_u(self.init) 

130 

131 if t == 0: 

132 me[:] = (0, 3, 3, 6, 0) 

133 elif t < self.t_end: 

134 me[:] = self.u_ref(t) 

135 else: 

136 self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.") 

137 me[:] = (0, 0, 0, 0, 0) 

138 return me 

139 

140 

141class two_transistor_amplifier(ptype_dae): 

142 r""" 

143 The two transistor amplifier example from page 108 in [1]_. The problem is an index-1 differential-algebraic equation 

144 (DAE) having the equations 

145 

146 .. math:: 

147 \frac{U_e (t)}{R_0} - \frac{U_1 (t)}{R_0} + C_1 (\frac{d U_2 (t)}{dt} - \frac{d U_1 (t)}{dt}) = 0, 

148 

149 .. math:: 

150 \frac{U_b}{R_2} - U_2 (t) (\frac{1}{R_1} + \frac{1}{R_2}) + C_1 (\frac{d U_1 (t)}{dt} - \frac{d U_2 (t)}{dt}) - (\alpha - 1) f(U_2 (t) - U_3 (t)) = 0, 

151 

152 .. math:: 

153 f(U_2 (t) - U_3 (t)) - \frac{U_3 (t)}{R_3} - C_2 \frac{d U_3 (t)}{dt} = 0, 

154 

155 .. math:: 

156 \frac{U_b}{R_4} - \frac{U_4 (t)}{R_4} + C_3 (\frac{d U_5 (t)}{dt} - \frac{d U_4 (t)}{dt}) - \alpha f(U_2 (t) - U_3 (t)) = 0, 

157 

158 .. math:: 

159 \frac{U_b}{R_6} - U_5 (t) (\frac{1}{R_5} + \frac{1}{R_6}) + C_3 (\frac{d U_4 (t)}{dt} - \frac{d U_5 (t)}{dt}) + (\alpha - 1) f(U_5 (t) - U_6 (t)) = 0, 

160 

161 .. math:: 

162 f(U_5 (t) - U_6 (t)) - \frac{U_6 (t)}{R_7} - C_4 \frac{d U_6 (t)}{dt} = 0, 

163 

164 .. math:: 

165 \frac{U_b}{R_8} - \frac{U_7 (t)}{R_8} - C_5 (\frac{d U_7 (t)}{dt} - \frac{d U_8 (t)}{dt}) - \alpha f(U_5 (t) - U_6 (t)) = 0, 

166 

167 .. math:: 

168 \frac{U_8 (t)}{R_9} - C_5 (\frac{d U_7 (t)}{dt} - \frac{d U_7 (t)}{dt}) = 0, 

169 

170 with 

171 

172 .. math:: 

173 f(U_i (t) - U_j (t)) = \beta (\exp(\frac{U_i (t) - U_j (t)}{U_F}) - 1). 

174 

175 The initial signal :math:`U_e (t)` is defined as 

176 

177 .. math:: 

178 U_e (t) = 0.1 \sin(200 \pi t). 

179 

180 Constants are fixed as :math:`U_b = 6`, :math:`U_F = 0.026`, :math:`\alpha = 0.99`, :math:`\beta = 10^{-6}`, :math:`R_0 = 1000`, 

181 :math:`R_k = 9000` for :math:`k=1,..,9`, `C_j = j \cdot 10^{-6}` for :math:`j=1,..,5`. They are also defined in the 

182 method `eval_f`. 

183 

184 Parameters 

185 ---------- 

186 nvars : int 

187 Number of unknowns of the system of DAEs. 

188 newton_tol : float 

189 Tolerance for Newton solver. 

190 

191 Attributes 

192 ---------- 

193 t_end: float 

194 The end time at which the reference solution is determined. 

195 

196 References 

197 ---------- 

198 .. [1] E. Hairer, C. Lubich, M. Roche. The numerical solution of differential-algebraic systems by Runge-Kutta methods. 

199 Lect. Notes Math. (1989). 

200 """ 

201 

202 dtype_u = mesh 

203 dtype_f = mesh 

204 

205 def __init__(self, newton_tol): 

206 super().__init__(nvars=8, newton_tol=newton_tol) 

207 # load reference solution 

208 # data file must be generated and stored under misc/data and self.t_end = t[-1] 

209 # data = np.load(r'pySDC/projects/DAE/misc/data/two_trans_amp.npy') 

210 # x = data[:, 0] 

211 # The last column contains the input signal 

212 # y = data[:, 1:-1] 

213 # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') 

214 self.t_end = 0.0 

215 

216 def eval_f(self, u, du, t): 

217 r""" 

218 Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. 

219 

220 Parameters 

221 ---------- 

222 u : dtype_u 

223 Current values of the numerical solution at time t. 

224 du : dtype_u 

225 Current values of the derivative of the numerical solution at time t. 

226 t : float 

227 Current time of the numerical solution. 

228 

229 Returns 

230 ------- 

231 f : dtype_f 

232 Current value of the right-hand side of f (which includes eight components). 

233 """ 

234 u_b = 6.0 

235 u_e = 0.1 * np.sin(200 * np.pi * t) 

236 alpha = 0.99 

237 r_0 = 1000.0 

238 r_k = 9000.0 

239 c_1, c_2, c_3, c_4, c_5 = 1e-6, 2e-6, 3e-6, 4e-6, 5e-6 

240 f = self.dtype_f(self.init) 

241 f[:] = ( 

242 (u_e - u[0]) / r_0 - c_1 * (du[0] - du[1]), 

243 (u_b - u[1]) / r_k - u[1] / r_k + c_1 * (du[0] - du[1]) + (alpha - 1) * _transistor(u[1] - u[2]), 

244 _transistor(u[1] - u[2]) - u[2] / r_k - c_2 * du[2], 

245 (u_b - u[3]) / r_k - c_3 * (du[3] - du[4]) - alpha * _transistor(u[1] - u[2]), 

246 (u_b - u[4]) / r_k - u[4] / r_k + c_3 * (du[3] - du[4]) + (alpha - 1) * _transistor(u[4] - u[5]), 

247 _transistor(u[4] - u[5]) - u[5] / r_k - c_4 * du[5], 

248 (u_b - u[6]) / r_k - c_5 * (du[6] - du[7]) - alpha * _transistor(u[4] - u[5]), 

249 -u[7] / r_k + c_5 * (du[6] - du[7]), 

250 ) 

251 self.work_counters['rhs']() 

252 return f 

253 

254 def u_exact(self, t): 

255 """ 

256 Dummy exact solution that should only be used to get initial conditions for the problem. This makes 

257 initialisation compatible with problems that have a known analytical solution. Could be used to output a 

258 reference solution if generated/available. 

259 

260 Parameters 

261 ---------- 

262 t : float 

263 The time of the reference solution. 

264 

265 Returns 

266 ------- 

267 me : dtype_u 

268 The reference solution as mesh object containing eight components and fixed initial conditions. 

269 """ 

270 me = self.dtype_u(self.init) 

271 

272 if t == 0: 

273 me[:] = (0, 3, 3, 6, 3, 3, 6, 0) 

274 elif t < self.t_end: 

275 me[:] = self.u_ref(t) 

276 else: 

277 self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.") 

278 me[:] = (0, 0, 0, 0, 0, 0, 0, 0) 

279 return me