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

1import warnings

2import numpy as np

3from scipy.interpolate import interp1d

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

6from pySDC.implementations.datatype_classes.mesh import mesh

9# Helper function

10def _transistor(u_in):

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

14class OneTransistorAmplifier(ProblemDAE):

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

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,

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,

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,

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,

31 .. math::

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

34 with

36 .. math::

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

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

41 .. math::

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

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

47 Parameters

48 ----------

49 nvars : int

50 Number of unknowns of the system of DAEs.

51 newton_tol : float

52 Tolerance for Newton solver.

54 Attributes

55 ----------

56 t_end: float

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

59 References

60 ----------

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

62 Springer (2009).

63 """

65 dtype_u = mesh

66 dtype_f = mesh

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

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

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.

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 TwoTransistorAmplifier(ProblemDAE):

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