Coverage for pySDC/projects/DAE/problems/DiscontinuousTestDAE.py: 100%
46 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
3from pySDC.core.Problem import WorkCounter
4from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
7class DiscontinuousTestDAE(ptype_dae):
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 get_switching_info(self, u, t):
132 r"""
133 Provides information about the state function of the problem. A change in sign of the state function
134 indicates an event. If a sign change is detected, a root can be found within the step according to the
135 intermediate value theorem.
137 The state function for this problem is given by
139 .. math::
140 h(y(t)) = 2 y(t) - 100.
142 Parameters
143 ----------
144 u : dtype_u
145 Current values of the numerical solution at time :math:`t`.
146 t : float
147 Current time of the numerical solution.
149 Returns
150 -------
151 switch_detected : bool
152 Indicates whether a discrete event is found or not.
153 m_guess : int
154 The index before the sign changes.
155 state_function : list
156 Defines the values of the state function at collocation nodes where it changes the sign.
157 """
159 switch_detected = False
160 m_guess = -100
162 for m in range(1, len(u)):
163 h_prev_node = 2 * u[m - 1].diff[0] - 100
164 h_curr_node = 2 * u[m].diff[0] - 100
165 if h_prev_node < 0 and h_curr_node >= 0:
166 switch_detected = True
167 m_guess = m - 1
168 break
170 state_function = [2 * u[m].diff[0] - 100 for m in range(len(u))]
171 return switch_detected, m_guess, state_function
173 def count_switches(self):
174 """
175 Setter to update the number of switches if one is found.
176 """
177 self.nswitches += 1