Coverage for pySDC/projects/DAE/problems/pendulum2D.py: 88%

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

2 

3 

4class Pendulum2D(ProblemDAE): 

5 r""" 

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

7 The DAE system is given by the equations 

8 

9 .. math:: 

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

11 

12 .. math:: 

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

14 

15 .. math:: 

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

17 

18 .. math:: 

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

20 

21 .. math:: 

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

23 

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

25 of [1]_. 

26 

27 Parameters 

28 ---------- 

29 nvars : int 

30 Number of unknowns of the system of DAEs. 

31 newton_tol : float 

32 Tolerance for Newton solver. 

33 

34 Attributes 

35 ---------- 

36 t_end: float 

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

38 

39 References 

40 ---------- 

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

42 Lect. Notes Math. (1989). 

43 """ 

44 

45 def __init__(self, newton_tol): 

46 """Initialization routine""" 

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

48 # load reference solution 

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

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

51 # t = data[:, 0] 

52 # solution = data[:, 1:] 

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

54 self.t_end = 0.0 

55 

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

57 r""" 

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

59 

60 Parameters 

61 ---------- 

62 u : dtype_u 

63 Current values of the numerical solution at time t. 

64 du : dtype_u 

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

66 t : float 

67 Current time of the numerical solution. 

68 

69 Returns 

70 ------- 

71 f : dtype_f 

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

73 """ 

74 g = 9.8 

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

76 # weight somehow 

77 f = self.dtype_f(self.init) 

78 f.diff[:4] = ( 

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

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

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

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

83 ) 

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

85 self.work_counters['rhs']() 

86 return f 

87 

88 def u_exact(self, t): 

89 """ 

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

91 

92 Parameters 

93 ---------- 

94 t : float 

95 The time of the reference solution. 

96 

97 Returns 

98 ------- 

99 me : dtype_u 

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

101 """ 

102 me = self.dtype_u(self.init) 

103 if t == 0: 

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

105 me.alg[0] = 0 

106 elif t < self.t_end: 

107 u_ref = self.u_ref(t) 

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

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

110 else: 

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

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

113 me.alg[0] = 0 

114 return me