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

## 55 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

3from pySDC.core.problem import WorkCounter

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

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:

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

14 .. math::

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

17 .. math::

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

20 else:

22 .. math::

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

25 .. math::

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

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

30 .. math::

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

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^*)).

37 Parameters

38 ----------

39 nvars : int

40 Number of unknowns of the system of DAEs.

41 newton_tol : float

42 Tolerance for Newton solver.

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.

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 """

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

60 """Initialization routine"""

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

63 self.t_switch_exact = np.arccosh(50)

64 self.t_switch = None

65 self.nswitches = 0

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

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

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.

81 Returns

82 -------

83 f : dtype_f

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

85 """

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

88 dy = du.diff[0]

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

92 h = 2 * y - 100

93 f = self.dtype_f(self.init)

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

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.

109 Parameters

110 ----------

111 t : float

112 Time of the exact solution.

114 Returns

115 -------

116 me : dtype_u

117 Exact solution.

118 """

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

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

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.

136 Parameters

137 ----------

138 t : float

139 Time of the exact solution.

141 Returns

142 -------

143 me : dtype_u

144 Derivative of exact solution.

145 """

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

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

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.

164 The state function for this problem is given by

166 .. math::

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

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.

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 """

186 switch_detected = False

187 m_guess = -100

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

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

198 return switch_detected, m_guess, state_function

200 def count_switches(self):

201 """

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

203 """

204 self.nswitches += 1