# Coverage for pySDC/implementations/problem_classes/DiscontinuousTestODE.py: 100%

## 87 statements

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

1import numpy as np

3from pySDC.core.errors import ParameterError, ProblemError

4from pySDC.core.problem import Problem, WorkCounter

5from pySDC.implementations.datatype_classes.mesh import mesh

8class DiscontinuousTestODE(Problem):

9 r"""

10 This class implements a very simple test case of a ordinary differential equation consisting of one discrete event. The dynamics of

11 the solution changes when the state function :math:h(u) := u - 5 changes the sign. The problem is defined by:

13 if :math:u - 5 < 0:

15 .. math::

16 \frac{d u}{dt} = u

18 else:

20 .. math::

21 \frac{d u}{dt} = \frac{4}{t^*},

23 where :math:t^* = \log(5) \approx 1.6094379. For :math:h(u) < 0, i.e. :math:t \leq t^*, the exact solution is

24 :math:u(t) = \exp(t); for :math:h(u) \geq 0, i.e. :math:t \geq t^*, the exact solution is :math:u(t) = \frac{4 t}{t^*} + 1.

26 Attributes

27 ----------

28 t_switch_exact : float

29 Exact event time with :math:t^* = \log(5).

30 t_switch: float

31 Time point of the discrete event found by switch estimation.

32 nswitches: int

33 Number of switches found by switch estimation.

35 Attributes

36 ----------

37 work_counters : WorkCounter

38 Counts different things, here: Number of Newton iterations is counted.

39 """

41 dtype_u = mesh

42 dtype_f = mesh

44 def __init__(self, newton_maxiter=100, newton_tol=1e-8, stop_at_nan=True):

45 """Initialization routine"""

46 nvars = 1

47 super().__init__(init=(nvars, None, np.dtype('float64')))

49 self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals())

51 if self.nvars != 1:

52 raise ParameterError('nvars has to be equal to 1!')

54 self.t_switch_exact = np.log(5)

55 self.t_switch = None

56 self.nswitches = 0

57 self.work_counters['newton'] = WorkCounter()

59 def eval_f(self, u, t):

60 """

61 Routine to evaluate the right-hand side of the problem.

63 Parameters

64 ----------

65 u : dtype_u

66 Current values of the numerical solution.

67 t : float

68 Current time of the numerical solution is computed.

70 Returns

71 -------

72 f : dtype_f

73 The right-hand side of the problem.

74 """

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

78 f = self.dtype_f(self.init, val=0.0)

79 h = u[0] - 5

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

81 f[:] = 4 / self.t_switch_exact

82 else:

83 f[:] = u

84 return f

86 def solve_system(self, rhs, dt, u0, t):

87 r"""

88 Simple Newton solver for :math:(I-factor\cdot A)\vec{u}=\vec{rhs}.

90 Parameters

91 ----------

92 rhs : dtype_f

93 Right-hand side for the linear system.

94 dt : float

95 Abbrev. for the local stepsize (or any other factor required).

96 u0 : dtype_u

97 Initial guess for the iterative solver.

98 t : float

99 Current time (e.g. for time-dependent BCs).

101 Returns

102 -------

103 me : dtype_u

104 The solution as mesh.

105 """

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

109 h = rhs[0] - 5

110 u = self.dtype_u(u0)

112 n = 0

113 res = 99

114 while n < self.newton_maxiter:

115 # form function g with g(u) = 0

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

117 g = u - dt * (4 / self.t_switch_exact) - rhs

118 else:

119 g = u - dt * u - rhs

121 # if g is close to 0, then we are done

122 res = np.linalg.norm(g, np.inf)

124 if res < self.newton_tol:

125 break

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

128 dg = 1

129 else:

130 dg = 1 - dt

132 # newton update

133 u -= 1.0 / dg * g

135 n += 1

136 self.work_counters['newton']()

138 if np.isnan(res) and self.stop_at_nan:

139 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)

140 elif np.isnan(res):

141 self.logger.warning('Newton got nan after %i iterations...' % n)

143 if n == self.newton_maxiter:

144 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))

146 me = self.dtype_u(self.init)

147 me[:] = u[:]

149 return me

151 def u_exact(self, t, u_init=None, t_init=None):

152 r"""

153 Routine to compute the exact solution at time :math:t.

155 Parameters

156 ----------

157 t : float

158 Time of the exact solution.

159 u_init : dtype_u

160 Initial conditions for getting the exact solution.

161 t_init : float

162 The starting time.

164 Returns

165 -------

166 me : dtype_u

167 The exact solution.

168 """

170 if t_init is not None and u_init is not None:

171 self.logger.warning(

172 f'{type(self).__name__} uses an analytic exact solution from t=0. If you try to compute the local error, you will get the global error instead!'

173 )

175 me = self.dtype_u(self.init)

176 if t <= self.t_switch_exact:

177 me[:] = np.exp(t)

178 else:

179 me[:] = (4 * t) / self.t_switch_exact + 1

180 return me

182 def get_switching_info(self, u, t):

183 r"""

184 Provides information about the state function of the problem. When the state function changes its sign,

185 typically an event occurs. So the check for an event should be done in the way that the state function

186 is checked for a sign change. If this is the case, the intermediate value theorem states a root in this

187 step.

189 Parameters

190 ----------

191 u : dtype_u

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

193 t : float

194 Current time of the numerical solution.

196 Returns

197 -------

198 switch_detected : bool

199 Indicates whether a discrete event is found or not.

200 m_guess : int

201 The index before the sign changes.

202 state_function : list

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

204 """

206 switch_detected = False

207 m_guess = -100

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

210 h_prev_node = u[m - 1][0] - 5

211 h_curr_node = u[m][0] - 5

212 if h_prev_node < 0 and h_curr_node >= 0:

213 switch_detected = True

214 m_guess = m - 1

215 break

217 state_function = [u[m][0] - 5 for m in range(len(u))]

218 return switch_detected, m_guess, state_function

220 def count_switches(self):

221 """

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

223 """

224 self.nswitches += 1

227class ExactDiscontinuousTestODE(DiscontinuousTestODE):

228 r"""

229 Dummy ODE problem for testing the SwitchEstimator class. The problem contains the exact dynamics

230 of the problem class DiscontinuousTestODE.

231 """

233 def __init__(self, newton_maxiter=100, newton_tol=1e-8):

234 """Initialization routine"""

235 super().__init__(newton_maxiter, newton_tol)

237 def eval_f(self, u, t):

238 """

239 Derivative.

241 Parameters

242 ----------

243 u : dtype_u

244 Exact value of u.

245 t : float

246 Time :math:t.

248 Returns

249 -------

250 f : dtype_f

251 Derivative.

252 """

254 f = self.dtype_f(self.init)

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

257 h = u[0] - 5

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

259 f[:] = 1

260 else:

261 f[:] = np.exp(t)

262 return f

264 def solve_system(self, rhs, factor, u0, t):

265 """

266 Just return the exact solution...

268 Parameters

269 ----------

270 rhs : dtype_f

271 Right-hand side for the linear system.

272 factor : float

273 Abbrev. for the local stepsize (or any other factor required).

274 u0 : dtype_u

275 Initial guess for the iterative solver.

276 t : float

277 Current time (e.g. for time-dependent BCs).

279 Returns

280 -------

281 me : dtype_u

282 The solution as mesh.

283 """

285 return self.u_exact(t)