Coverage for pySDC/projects/DAE/problems/discontinuousTestDAE.py: 96%
56 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +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^* = 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 After the event, :math:`h` stays at zero instead of changing its sign. A switch estimator, which looks for a sign
38 change, then only sees the overshoot of the iterates. So with ``switch_on_state=False``, the dynamics switch only
39 at the event time :math:`t_{switch}` the estimator has located, and :math:`h` changes its sign across the event.
41 Parameters
42 ----------
43 nvars : int
44 Number of unknowns of the system of DAEs.
45 newton_tol : float
46 Tolerance for Newton solver.
47 switch_on_state : bool, optional
48 Whether the dynamics switch where :math:`h(y) \geq 0`, or only after the located event time
49 :math:`t_{switch}`. Set it to False when the event is located by a switch estimator.
51 Attributes
52 ----------
53 t_switch_exact: float
54 Exact time of the event :math:`t^*`.
55 t_switch: float
56 Time point of the event found by switch estimation.
57 nswitches: int
58 Number of switches found by switch estimation.
60 References
61 ----------
62 .. [1] L. Lopez, S. Maset. Numerical event location techniques in discontinuous differential algebraic equations.
63 Appl. Numer. Math. 178, 98-122 (2022).
64 """
66 def __init__(self, newton_tol=1e-12, switch_on_state=True):
67 """Initialization routine"""
68 super().__init__(nvars=2, newton_tol=newton_tol)
69 self._makeAttributeAndRegister('switch_on_state', localVars=locals(), readOnly=True)
71 self.t_switch_exact = np.arccosh(50)
72 self.t_switch = None
73 self.nswitches = 0
74 self.work_counters['rhs'] = WorkCounter()
76 def eval_f(self, u, du, t):
77 r"""
78 Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`.
80 Parameters
81 ----------
82 u : dtype_u
83 Current values of the numerical solution at time t.
84 du : dtype_u
85 Current values of the derivative of the numerical solution at time t.
86 t : float
87 Current time of the numerical solution.
89 Returns
90 -------
91 f : dtype_f
92 The right-hand side of f (contains two components).
93 """
95 y, z = u.diff[0], u.alg[0]
96 dy = du.diff[0]
98 t_switch = np.inf if self.t_switch is None else self.t_switch
100 h = 2 * y - 100
101 f = self.dtype_f(self.init)
103 if (self.switch_on_state and h >= 0) or t >= t_switch:
104 f.diff[0] = dy
105 f.alg[0] = y**2 - z**2 - 1
106 else:
107 f.diff[0] = dy - z
108 f.alg[0] = y**2 - z**2 - 1
109 self.work_counters['rhs']()
110 return f
112 def u_exact(self, t, **kwargs):
113 r"""
114 Routine for the exact solution at time :math:`t \leq 1`. For this problem, the exact
115 solution is piecewise.
117 Parameters
118 ----------
119 t : float
120 Time of the exact solution.
122 Returns
123 -------
124 me : dtype_u
125 Exact solution.
126 """
128 assert t >= 1, 'ERROR: u_exact only available for t>=1'
130 me = self.dtype_u(self.init)
131 if t <= self.t_switch_exact:
132 me.diff[0] = np.cosh(t)
133 me.alg[0] = np.sinh(t)
134 else:
135 me.diff[0] = np.cosh(self.t_switch_exact)
136 me.alg[0] = np.sinh(self.t_switch_exact)
137 return me
139 def du_exact(self, t, **kwargs):
140 r"""
141 Routine for the derivative of the exact solution at time :math:`t \leq 1`.
142 For this problem, the exact solution is piecewise.
144 Parameters
145 ----------
146 t : float
147 Time of the exact solution.
149 Returns
150 -------
151 me : dtype_u
152 Derivative of exact solution.
153 """
155 assert t >= 1, 'ERROR: u_exact only available for t>=1'
157 me = self.dtype_u(self.init)
158 if t <= self.t_switch_exact:
159 me.diff[0] = np.sinh(t)
160 me.alg[0] = np.cosh(t)
161 else:
162 me.diff[0] = np.sinh(self.t_switch_exact)
163 me.alg[0] = np.cosh(self.t_switch_exact)
164 return me
166 def get_switching_info(self, u, t):
167 r"""
168 Provides information about the state function of the problem. A change in sign of the state function
169 indicates an event. If a sign change is detected, a root can be found within the step according to the
170 intermediate value theorem.
172 The state function for this problem is given by
174 .. math::
175 h(y(t)) = 2 y(t) - 100.
177 Parameters
178 ----------
179 u : dtype_u
180 Current values of the numerical solution at time :math:`t`.
181 t : float
182 Current time of the numerical solution.
184 Returns
185 -------
186 switch_detected : bool
187 Indicates whether a discrete event is found or not.
188 m_guess : int
189 The index before the sign changes.
190 state_function : list
191 Defines the values of the state function at collocation nodes where it changes the sign.
192 """
194 switch_detected = False
195 m_guess = -100
197 for m in range(1, len(u)):
198 h_prev_node = 2 * u[m - 1].diff[0] - 100
199 h_curr_node = 2 * u[m].diff[0] - 100
200 if h_prev_node < 0 and h_curr_node >= 0:
201 switch_detected = True
202 m_guess = m - 1
203 break
205 state_function = [2 * u[m].diff[0] - 100 for m in range(len(u))]
206 return switch_detected, m_guess, state_function
208 def count_switches(self):
209 """
210 Setter to update the number of switches if one is found.
211 """
212 self.nswitches += 1