Coverage for pySDC/projects/DAE/problems/discontinuousTestDAE.py: 96%

1import numpy as np 

2 

3from pySDC.core.problem import WorkCounter 

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

5 

6 

7class DiscontinuousTestDAE(ProblemDAE): 

8 r""" 

9 This class implements a scalar test discontinuous differential-algebraic equation (DAE) similar to [1]_. The event function 

10 is defined as :math:`h(y):= 2y - 100`. Then, the discontinuous DAE model reads: 

11 

12 - if :math:`h(y) \leq 0`: 

13 

14 .. math:: 

15 \dfrac{d y(t)}{dt} = z(t), 

16 

17 .. math:: 

18 y(t)^2 - z(t)^2 - 1 = 0, 

19 

20 else: 

21 

22 .. math:: 

23 \dfrac{d y(t)}{dt} = 0, 

24 

25 .. math:: 

26 y(t)^2 - z(t)^2 - 1 = 0, 

27 

28 for :math:`t \geq 1`. If :math:`h(y) < 0`, the solution is given by 

29 

30 .. math:: 

31 (y(t), z(t)) = (cosh(t), sinh(t)), 

32 

33 and the event takes place at :math:`t^* = 0.5 * arccosh(50) = 4.60507` and event point :math:`(cosh(t^*), sinh(t^*))`. 

34 As soon as :math:`t \geq t^*`, i.e., for :math:`h(y) \geq 0`, the solution is given by the constant value 

35 :math:`(cosh(t^*), sinh(t^*))`. 

36 

37 Parameters 

38 ---------- 

39 nvars : int 

40 Number of unknowns of the system of DAEs. 

41 newton_tol : float 

42 Tolerance for Newton solver. 

43 

44 Attributes 

45 ---------- 

46 t_switch_exact: float 

47 Exact time of the event :math:`t^*`. 

48 t_switch: float 

49 Time point of the event found by switch estimation. 

50 nswitches: int 

51 Number of switches found by switch estimation. 

52 

53 References 

54 ---------- 

55 .. [1] L. Lopez, S. Maset. Numerical event location techniques in discontinuous differential algebraic equations. 

56 Appl. Numer. Math. 178, 98-122 (2022). 

57 """ 

58 

59 def __init__(self, newton_tol=1e-12): 

60 """Initialization routine""" 

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

62 

63 self.t_switch_exact = np.arccosh(50) 

64 self.t_switch = None 

65 self.nswitches = 0 

66 self.work_counters['rhs'] = WorkCounter() 

67 

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

69 r""" 

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

71 

72 Parameters 

73 ---------- 

74 u : dtype_u 

75 Current values of the numerical solution at time t. 

76 du : dtype_u 

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

78 t : float 

79 Current time of the numerical solution. 

80 

81 Returns 

82 ------- 

83 f : dtype_f 

84 The right-hand side of f (contains two components). 

85 """ 

86 

87 y, z = u.diff[0], u.alg[0] 

88 dy = du.diff[0] 

89 

90 t_switch = np.inf if self.t_switch is None else self.t_switch 

91 

92 h = 2 * y - 100 

93 f = self.dtype_f(self.init) 

94 

95 if h >= 0 or t >= t_switch: 

96 f.diff[0] = dy 

97 f.alg[0] = y**2 - z**2 - 1 

98 else: 

99 f.diff[0] = dy - z 

100 f.alg[0] = y**2 - z**2 - 1 

101 self.work_counters['rhs']() 

102 return f 

103 

104 def u_exact(self, t, **kwargs): 

105 r""" 

106 Routine for the exact solution at time :math:`t \leq 1`. For this problem, the exact 

107 solution is piecewise. 

108 

109 Parameters 

110 ---------- 

111 t : float 

112 Time of the exact solution. 

113 

114 Returns 

115 ------- 

116 me : dtype_u 

117 Exact solution. 

118 """ 

119 

120 assert t >= 1, 'ERROR: u_exact only available for t>=1' 

121 

122 me = self.dtype_u(self.init) 

123 if t <= self.t_switch_exact: 

124 me.diff[0] = np.cosh(t) 

125 me.alg[0] = np.sinh(t) 

126 else: 

127 me.diff[0] = np.cosh(self.t_switch_exact) 

128 me.alg[0] = np.sinh(self.t_switch_exact) 

129 return me 

130 

131 def du_exact(self, t, **kwargs): 

132 r""" 

133 Routine for the derivative of the exact solution at time :math:`t \leq 1`. 

134 For this problem, the exact solution is piecewise. 

135 

136 Parameters 

137 ---------- 

138 t : float 

139 Time of the exact solution. 

140 

141 Returns 

142 ------- 

143 me : dtype_u 

144 Derivative of exact solution. 

145 """ 

146 

147 assert t >= 1, 'ERROR: u_exact only available for t>=1' 

148 

149 me = self.dtype_u(self.init) 

150 if t <= self.t_switch_exact: 

151 me.diff[0] = np.sinh(t) 

152 me.alg[0] = np.cosh(t) 

153 else: 

154 me.diff[0] = np.sinh(self.t_switch_exact) 

155 me.alg[0] = np.cosh(self.t_switch_exact) 

156 return me 

157 

158 def get_switching_info(self, u, t): 

159 r""" 

160 Provides information about the state function of the problem. A change in sign of the state function 

161 indicates an event. If a sign change is detected, a root can be found within the step according to the 

162 intermediate value theorem. 

163 

164 The state function for this problem is given by 

165 

166 .. math:: 

167 h(y(t)) = 2 y(t) - 100. 

168 

169 Parameters 

170 ---------- 

171 u : dtype_u 

172 Current values of the numerical solution at time :math:`t`. 

173 t : float 

174 Current time of the numerical solution. 

175 

176 Returns 

177 ------- 

178 switch_detected : bool 

179 Indicates whether a discrete event is found or not. 

180 m_guess : int 

181 The index before the sign changes. 

182 state_function : list 

183 Defines the values of the state function at collocation nodes where it changes the sign. 

184 """ 

185 

186 switch_detected = False 

187 m_guess = -100 

188 

189 for m in range(1, len(u)): 

190 h_prev_node = 2 * u[m - 1].diff[0] - 100 

191 h_curr_node = 2 * u[m].diff[0] - 100 

192 if h_prev_node < 0 and h_curr_node >= 0: 

193 switch_detected = True 

194 m_guess = m - 1 

195 break 

196 

197 state_function = [2 * u[m].diff[0] - 100 for m in range(len(u))] 

198 return switch_detected, m_guess, state_function 

199 

200 def count_switches(self): 

201 """ 

202 Setter to update the number of switches if one is found. 

203 """ 

204 self.nswitches += 1