Coverage for pySDC/implementations/problem_classes/DiscontinuousTestODE.py: 100%
87 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +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 Parameters
27 ----------
28 newton_maxiter : int, optional
29 Maximum number of iterations for Newton's method.
30 newton_tol : float, optional
31 Absolute tolerance for Newton's method to terminate, applied to the maximum norm of the residual.
32 stop_at_nan : bool, optional
33 Indicates that the Newton solver should stop if ``nan`` values arise.
35 Attributes
36 ----------
37 t_switch_exact : float
38 Exact event time with :math:`t^* = \log(5)`.
39 t_switch: float
40 Time point of the discrete event found by switch estimation.
41 nswitches: int
42 Number of switches found by switch estimation.
44 Attributes
45 ----------
46 work_counters : WorkCounter
47 Counts different things, here: Number of Newton iterations is counted.
48 """
50 dtype_u = mesh
51 dtype_f = mesh
53 def __init__(self, newton_maxiter=100, newton_tol=1e-8, stop_at_nan=True):
54 """Initialization routine"""
55 nvars = 1
56 super().__init__(init=(nvars, None, np.dtype('float64')))
57 self._makeAttributeAndRegister('nvars', localVars=locals(), readOnly=True)
58 self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals())
60 if self.nvars != 1:
61 raise ParameterError('nvars has to be equal to 1!')
63 self.t_switch_exact = np.log(5)
64 self.t_switch = None
65 self.nswitches = 0
66 self.work_counters['newton'] = WorkCounter()
68 def eval_f(self, u, t):
69 """
70 Routine to evaluate the right-hand side of the problem.
72 Parameters
73 ----------
74 u : dtype_u
75 Current values of the numerical solution.
76 t : float
77 Current time of the numerical solution is computed.
79 Returns
80 -------
81 f : dtype_f
82 The right-hand side of the problem.
83 """
85 t_switch = np.inf if self.t_switch is None else self.t_switch
87 f = self.dtype_f(self.init, val=0.0)
88 h = u[0] - 5
89 if h >= 0 or t >= t_switch:
90 f[:] = 4 / self.t_switch_exact
91 else:
92 f[:] = u
93 return f
95 def solve_system(self, rhs, dt, u0, t):
96 r"""
97 Simple Newton solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
99 Parameters
100 ----------
101 rhs : dtype_f
102 Right-hand side for the linear system.
103 dt : float
104 Abbrev. for the local stepsize (or any other factor required).
105 u0 : dtype_u
106 Initial guess for the iterative solver.
107 t : float
108 Current time (e.g. for time-dependent BCs).
110 Returns
111 -------
112 me : dtype_u
113 The solution as mesh.
114 """
116 t_switch = np.inf if self.t_switch is None else self.t_switch
118 h = rhs[0] - 5
119 u = self.dtype_u(u0)
121 n = 0
122 res = 99
123 while n < self.newton_maxiter:
124 # form function g with g(u) = 0
125 if h >= 0 or t >= t_switch:
126 g = u - dt * (4 / self.t_switch_exact) - rhs
127 else:
128 g = u - dt * u - rhs
130 # if g is close to 0, then we are done
131 res = np.linalg.norm(g, np.inf)
133 if res < self.newton_tol:
134 break
136 if h >= 0 or t >= t_switch:
137 dg = 1
138 else:
139 dg = 1 - dt
141 # newton update
142 u -= 1.0 / dg * g
144 n += 1
145 self.work_counters['newton']()
147 if np.isnan(res) and self.stop_at_nan:
148 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
149 elif np.isnan(res):
150 self.logger.warning('Newton got nan after %i iterations...' % n)
152 if n == self.newton_maxiter:
153 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
155 me = self.dtype_u(self.init)
156 me[:] = u[:]
158 return me
160 def u_exact(self, t, u_init=None, t_init=None):
161 r"""
162 Routine to compute the exact solution at time :math:`t`.
164 Parameters
165 ----------
166 t : float
167 Time of the exact solution.
168 u_init : dtype_u
169 Initial conditions for getting the exact solution.
170 t_init : float
171 The starting time.
173 Returns
174 -------
175 me : dtype_u
176 The exact solution.
177 """
179 if t_init is not None and u_init is not None:
180 self.logger.warning(
181 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!'
182 )
184 me = self.dtype_u(self.init)
185 if t <= self.t_switch_exact:
186 me[:] = np.exp(t)
187 else:
188 me[:] = (4 * t) / self.t_switch_exact + 1
189 return me
191 def get_switching_info(self, u, t):
192 r"""
193 Provides information about the state function of the problem. When the state function changes its sign,
194 typically an event occurs. So the check for an event should be done in the way that the state function
195 is checked for a sign change. If this is the case, the intermediate value theorem states a root in this
196 step.
198 Parameters
199 ----------
200 u : dtype_u
201 Current values of the numerical solution at time :math:`t`.
202 t : float
203 Current time of the numerical solution.
205 Returns
206 -------
207 switch_detected : bool
208 Indicates whether a discrete event is found or not.
209 m_guess : int
210 The index before the sign changes.
211 state_function : list
212 Defines the values of the state function at collocation nodes where it changes the sign.
213 """
215 switch_detected = False
216 m_guess = -100
218 for m in range(1, len(u)):
219 h_prev_node = u[m - 1][0] - 5
220 h_curr_node = u[m][0] - 5
221 if h_prev_node < 0 and h_curr_node >= 0:
222 switch_detected = True
223 m_guess = m - 1
224 break
226 state_function = [u[m][0] - 5 for m in range(len(u))]
227 return switch_detected, m_guess, state_function
229 def count_switches(self):
230 """
231 Setter to update the number of switches if one is found.
232 """
233 self.nswitches += 1
236class ExactDiscontinuousTestODE(DiscontinuousTestODE):
237 r"""
238 Variant of ``DiscontinuousTestODE`` whose ``solve_system`` returns the exact solution, to test the
239 ``SwitchEstimator``.
240 The problem contains the exact dynamics
241 of the problem class ``DiscontinuousTestODE``.
243 Parameters
244 ----------
245 newton_maxiter : int, optional
246 Maximum number of iterations for Newton's method. Not used here, since ``solve_system`` returns the exact
247 solution.
248 newton_tol : float, optional
249 Tolerance for Newton's method to terminate. Not used here, since ``solve_system`` returns the exact solution.
250 """
252 def __init__(self, newton_maxiter=100, newton_tol=1e-8):
253 """Initialization routine"""
254 super().__init__(newton_maxiter, newton_tol)
256 def eval_f(self, u, t):
257 """
258 Derivative.
260 Parameters
261 ----------
262 u : dtype_u
263 Exact value of u.
264 t : float
265 Time :math:`t`.
267 Returns
268 -------
269 f : dtype_f
270 Derivative.
271 """
273 f = self.dtype_f(self.init)
275 t_switch = np.inf if self.t_switch is None else self.t_switch
276 h = u[0] - 5
277 if h >= 0 or t >= t_switch:
278 f[:] = 4 / self.t_switch_exact
279 else:
280 f[:] = np.exp(t)
281 return f
283 def solve_system(self, rhs, factor, u0, t):
284 """
285 Just return the exact solution...
287 Parameters
288 ----------
289 rhs : dtype_f
290 Right-hand side for the linear system.
291 factor : float
292 Abbrev. for the local stepsize (or any other factor required).
293 u0 : dtype_u
294 Initial guess for the iterative solver.
295 t : float
296 Current time (e.g. for time-dependent BCs).
298 Returns
299 -------
300 me : dtype_u
301 The solution as mesh.
302 """
304 return self.u_exact(t)