Coverage for pySDC/implementations/problem_classes/DiscontinuousTestODE.py: 100%
87 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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')))
48 self._makeAttributeAndRegister('nvars', localVars=locals(), readOnly=True)
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)