Coverage for pySDC/projects/DAE/problems/simple_DAE.py: 95%

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 

9class pendulum_2d(ptype_dae): 

10 r""" 

11 Example implementing the well known 2D pendulum as a first order differential-algebraic equation (DAE) of index 3. 

12 The DAE system is given by the equations 

13 

14 .. math:: 

15 \frac{dp}{dt} = u, 

16 

17 .. math:: 

18 \frac{dq}{dt} = v, 

19 

20 .. math:: 

21 m\frac{du}{dt} = -p \lambda, 

22 

23 .. math:: 

24 m\frac{dv}{dt} = -q \lambda - g, 

25 

26 .. math:: 

27 0 = p^2 + q^2 - l^2 

28 

29 for :math:`l=1` and :math:`m=1`. The pendulum is used in most introductory literature on DAEs, for example on page 8 

30 of [1]_. 

31 

32 Parameters 

33 ---------- 

34 nvars : int 

35 Number of unknowns of the system of DAEs. 

36 newton_tol : float 

37 Tolerance for Newton solver. 

38 

39 Attributes 

40 ---------- 

41 t_end: float 

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

43 

44 References 

45 ---------- 

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

47 Lect. Notes Math. (1989). 

48 """ 

49 

50 def __init__(self, newton_tol): 

51 """Initialization routine""" 

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

53 # load reference solution 

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

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

56 # t = data[:, 0] 

57 # solution = data[:, 1:] 

58 # self.u_ref = interp1d(t, solution, kind='cubic', axis=0, fill_value='extrapolate') 

59 self.t_end = 0.0 

60 

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

62 r""" 

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

64 

65 Parameters 

66 ---------- 

67 u : dtype_u 

68 Current values of the numerical solution at time t. 

69 du : dtype_u 

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

71 t : float 

72 Current time of the numerical solution. 

73 

74 Returns 

75 ------- 

76 f : dtype_f 

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

78 """ 

79 g = 9.8 

80 # The last element of u is a Lagrange multiplier. Not sure if this needs to be time dependent, but must model the 

81 # weight somehow 

82 f = self.dtype_f(self.init) 

83 f.diff[:4] = ( 

84 du.diff[0] - u.diff[2], 

85 du.diff[1] - u.diff[3], 

86 du.diff[2] + u.alg[0] * u.diff[0], 

87 du.diff[3] + u.alg[0] * u.diff[1] + g, 

88 ) 

89 f.alg[0] = u.diff[0] ** 2 + u.diff[1] ** 2 - 1 

90 self.work_counters['rhs']() 

91 return f 

92 

93 def u_exact(self, t): 

94 """ 

95 Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution. 

96 

97 Parameters 

98 ---------- 

99 t : float 

100 The time of the reference solution. 

101 

102 Returns 

103 ------- 

104 me : dtype_u 

105 The reference solution as mesh object. It contains fixed initial conditions at initial time. 

106 """ 

107 me = self.dtype_u(self.init) 

108 if t == 0: 

109 me.diff[:4] = (-1, 0, 0, 0) 

110 me.alg[0] = 0 

111 elif t < self.t_end: 

112 u_ref = self.u_ref(t) 

113 me.diff[:4] = u_ref[:4] 

114 me.alg[0] = u_ref[5] 

115 else: 

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

117 me.diff[:4] = (0, 0, 0, 0) 

118 me.alg[0] = 0 

119 return me 

120 

121 

122class simple_dae_1(ptype_dae): 

123 r""" 

124 Example implementing a smooth linear index-2 differential-algebraic equation (DAE) with known analytical solution. 

125 The DAE system is given by 

126 

127 .. math:: 

128 \frac{d u_1 (t)}{dt} = (\alpha - \frac{1}{2 - t}) u_1 (t) + (2-t) \alpha z (t) + \frac{3 - t}{2 - t}, 

129 

130 .. math:: 

131 \frac{d u_2 (t)}{dt} = \frac{1 - \alpha}{t - 2} u_1 (t) - u_2 (t) + (\alpha - 1) z (t) + 2 e^{t}, 

132 

133 .. math:: 

134 0 = (t + 2) u_1 (t) + (t^{2} - 4) u_2 (t) - (t^{2} + t - 2) e^{t}. 

135 

136 The exact solution of this system is 

137 

138 .. math:: 

139 u_1 (t) = u_2 (t) = e^{t}, 

140 

141 .. math:: 

142 z (t) = -\frac{e^{t}}{2 - t}. 

143 

144 This example is commonly used to test that numerical implementations are functioning correctly. See, for example, 

145 page 267 of [1]_. 

146 

147 Parameters 

148 ---------- 

149 nvars : int 

150 Number of unknowns of the system of DAEs. 

151 newton_tol : float 

152 Tolerance for Newton solver. 

153 

154 References 

155 ---------- 

156 .. [1] U. Ascher, L. R. Petzold. Computer method for ordinary differential equations and differential-algebraic 

157 equations. Society for Industrial and Applied Mathematics (1998). 

158 """ 

159 

160 def __init__(self, newton_tol=1e-10): 

161 """Initialization routine""" 

162 super().__init__(nvars=3, newton_tol=newton_tol) 

163 

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

165 r""" 

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

167 

168 Parameters 

169 ---------- 

170 u : dtype_u 

171 Current values of the numerical solution at time t. 

172 du : dtype_u 

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

174 t : float 

175 Current time of the numerical solution. 

176 

177 Returns 

178 ------- 

179 f : dtype_f 

180 Current value of the right-hand side of f (which includes three components). 

181 """ 

182 # Smooth index-2 DAE pg. 267 Ascher and Petzold (also the first example in KDC Minion paper) 

183 a = 10.0 

184 f = self.dtype_f(self.init) 

185 

186 f.diff[:2] = ( 

187 -du.diff[0] + (a - 1 / (2 - t)) * u.diff[0] + (2 - t) * a * u.alg[0] + (3 - t) / (2 - t) * np.exp(t), 

188 -du.diff[1] + (1 - a) / (t - 2) * u.diff[0] - u.diff[1] + (a - 1) * u.alg[0] + 2 * np.exp(t), 

189 ) 

190 f.alg[0] = (t + 2) * u.diff[0] + (t**2 - 4) * u.diff[1] - (t**2 + t - 2) * np.exp(t) 

191 self.work_counters['rhs']() 

192 return f 

193 

194 def u_exact(self, t): 

195 """ 

196 Routine for the exact solution. 

197 

198 Parameters 

199 ---------- 

200 t : float 

201 The time of the reference solution. 

202 

203 Returns 

204 ------- 

205 me : dtype_u 

206 The reference solution as mesh object containing three components. 

207 """ 

208 me = self.dtype_u(self.init) 

209 me.diff[:2] = (np.exp(t), np.exp(t)) 

210 me.alg[0] = -np.exp(t) / (2 - t) 

211 return me 

212 

213 

214class problematic_f(ptype_dae): 

215 r""" 

216 Standard example of a very simple fully implicit index-2 differential algebraic equation (DAE) that is not 

217 numerically solvable for certain choices of the parameter :math:`\eta`. The DAE system is given by 

218 

219 .. math:: 

220 \frac{d y(t)}{dt} + \eta t \frac{d z(t)}{dt} + (1 + \eta) z (t) = g (t). 

221 

222 .. math:: 

223 y (t) + \eta t z (t) = f(t), 

224 

225 See, for example, page 264 of [1]_. 

226 

227 Parameters 

228 ---------- 

229 nvars : int 

230 Number of unknowns of the system of DAEs. 

231 newton_tol : float 

232 Tolerance for Newton solver. 

233 

234 Attributes 

235 ---------- 

236 eta : float 

237 Specific parameter of the problem. 

238 

239 References 

240 ---------- 

241 .. [1] U. Ascher, L. R. Petzold. Computer method for ordinary differential equations and differential-algebraic 

242 equations. Society for Industrial and Applied Mathematics (1998). 

243 """ 

244 

245 dtype_u = mesh 

246 dtype_f = mesh 

247 

248 def __init__(self, newton_tol, eta=1): 

249 """Initialization routine""" 

250 super().__init__(nvars=2, newton_tol=newton_tol) 

251 self._makeAttributeAndRegister('eta', localVars=locals()) 

252 

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

254 r""" 

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

256 

257 Parameters 

258 ---------- 

259 u : dtype_u 

260 Current values of the numerical solution at time t. 

261 du : dtype_u 

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

263 t : float 

264 Current time of the numerical solution. 

265 

266 Returns 

267 ------- 

268 f : dtype_f 

269 Current value of the right-hand side of f (which includes two components). 

270 """ 

271 f = self.dtype_f(self.init) 

272 f[:] = ( 

273 u[0] + self.eta * t * u[1] - np.sin(t), 

274 du[0] + self.eta * t * du[1] + (1 + self.eta) * u[1] - np.cos(t), 

275 ) 

276 self.work_counters['rhs']() 

277 return f 

278 

279 def u_exact(self, t): 

280 """ 

281 Routine for the exact solution. 

282 

283 Parameters 

284 ---------- 

285 t : float 

286 The time of the reference solution. 

287 

288 Returns 

289 ------- 

290 me : dtype_u 

291 The reference solution as mesh object containing two components. 

292 """ 

293 me = self.dtype_u(self.init) 

294 me[:] = (np.sin(t), 0) 

295 return me