Coverage for pySDC/implementations/problem_classes/odeScalar.py: 100%
53 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Implementation of scalar test problem ODEs.
7Reference :
9Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods
10on parallel computers. SIAM journal on scientific and statistical computing,
1112(5), 1000-1028.
12"""
14import numpy as np
16from pySDC.core.errors import ProblemError
17from pySDC.core.problem import Problem, WorkCounter
18from pySDC.implementations.datatype_classes.mesh import mesh
21class ProtheroRobinson(Problem):
22 r"""
23 Implement the Prothero-Robinson problem:
25 .. math::
26 \frac{du}{dt} = -\frac{u-g(t)}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).,
28 with :math:`\epsilon` a stiffness parameter, that makes the problem more stiff
29 the smaller it is (usual taken value is :math:`\epsilon=1e^{-3}`).
30 Exact solution is given by :math:`u(t)=g(t)`, and this implementation uses
31 :math:`g(t)=\cos(t)`.
33 Implement also the non-linear form of this problem:
35 .. math::
36 \frac{du}{dt} = -\frac{u^3-g(t)^3}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).
38 To use an other exact solution, one just have to derivate this class
39 and overload the `g` and `dg` methods. For instance,
40 to use :math:`g(t)=e^{-0.2*t}`, define and use the following class:
42 >>> class MyProtheroRobinson(ProtheroRobinson):
43 >>>
44 >>> def g(self, t):
45 >>> return np.exp(-0.2 * t)
46 >>>
47 >>> def dg(self, t):
48 >>> return (-0.2) * np.exp(-0.2 * t)
50 Parameters
51 ----------
52 epsilon : float, optional
53 Stiffness parameter. The default is 1e-3.
54 nonLinear : bool, optional
55 Wether or not to use the non-linear form of the problem. The default is False.
56 newton_maxiter : int, optional
57 Maximum number of Newton iteration in solve_system. The default is 200.
58 newton_tol : float, optional
59 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.
60 stop_at_nan : bool, optional
61 Wheter to stop or not solve_system when getting NAN. The default is True.
63 References
64 ----------
65 A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving
66 stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974),
67 pp. 145–162.
68 """
70 dtype_u = mesh
71 dtype_f = mesh
73 def __init__(self, epsilon=1e-3, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):
74 """Initialization routine"""
75 nvars = 1
76 super().__init__((nvars, None, np.dtype('float64')))
78 self.f = self.f_NONLIN if nonLinear else self.f_LIN
79 self.jac = self.jac_NONLIN if nonLinear else self.jac_LIN
80 self._makeAttributeAndRegister(
81 'epsilon', 'nonLinear', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True
82 )
83 self.work_counters['newton'] = WorkCounter()
84 self.work_counters['rhs'] = WorkCounter()
86 # -------------------------------------------------------------------------
87 # g function (analytical solution), and its first derivative
88 # -------------------------------------------------------------------------
89 def g(self, t):
90 r"""
91 Function :math:`g(t) = \cos(t)`, which is the exact solution.
93 Parameters
94 ----------
95 t : float
96 Time.
98 Returns
99 -------
100 float
101 The value of :math:`g(t)`.
102 """
103 return np.cos(t)
105 def dg(self, t):
106 r"""
107 First derivative :math:`g'(t) = -\sin(t)` of :math:`g`.
109 Parameters
110 ----------
111 t : float
112 Time.
114 Returns
115 -------
116 float
117 The value of :math:`g'(t)`.
118 """
119 return -np.sin(t)
121 # -------------------------------------------------------------------------
122 # f(u,t) and Jacobian functions
123 # -------------------------------------------------------------------------
124 def f(self, u, t):
125 """
126 Right-hand side, replaced in ``__init__`` by ``f_LIN`` or ``f_NONLIN``.
128 Parameters
129 ----------
130 u : dtype_u
131 Current values of the numerical solution.
132 t : float
133 Current time.
135 Raises
136 ------
137 NotImplementedError
138 Always.
139 """
140 raise NotImplementedError()
142 def f_LIN(self, u, t):
143 r"""
144 Right-hand side :math:`-(u - g(t)) / \epsilon + g'(t)` of the linear form.
146 Parameters
147 ----------
148 u : dtype_u
149 Current values of the numerical solution.
150 t : float
151 Current time.
153 Returns
154 -------
155 dtype_u
156 The right-hand side.
157 """
158 return -self.epsilon ** (-1) * (u - self.g(t)) + self.dg(t)
160 def f_NONLIN(self, u, t):
161 r"""
162 Right-hand side :math:`-(u^3 - g(t)^3) / \epsilon + g'(t)` of the nonlinear form.
164 Parameters
165 ----------
166 u : dtype_u
167 Current values of the numerical solution.
168 t : float
169 Current time.
171 Returns
172 -------
173 dtype_u
174 The right-hand side.
175 """
176 return -self.epsilon ** (-1) * (u**3 - self.g(t) ** 3) + self.dg(t)
178 def jac(self, u, t):
179 """
180 Jacobian of the right-hand side, replaced in ``__init__`` by ``jac_LIN`` or ``jac_NONLIN``.
182 Parameters
183 ----------
184 u : dtype_u
185 Current values of the numerical solution.
186 t : float
187 Current time.
189 Raises
190 ------
191 NotImplementedError
192 Always.
193 """
194 raise NotImplementedError()
196 def jac_LIN(self, u, t):
197 r"""
198 Jacobian :math:`-1/\epsilon` of the right-hand side of the linear form.
200 Parameters
201 ----------
202 u : dtype_u
203 Current values of the numerical solution.
204 t : float
205 Current time.
207 Returns
208 -------
209 float
210 The Jacobian.
211 """
212 return -self.epsilon ** (-1)
214 def jac_NONLIN(self, u, t):
215 r"""
216 Jacobian :math:`-3 u^2 / \epsilon` of the right-hand side of the nonlinear form.
218 Parameters
219 ----------
220 u : dtype_u
221 Current values of the numerical solution.
222 t : float
223 Current time.
225 Returns
226 -------
227 dtype_u
228 The Jacobian.
229 """
230 return -self.epsilon ** (-1) * 3 * u**2
232 # -------------------------------------------------------------------------
233 # pySDC required methods
234 # -------------------------------------------------------------------------
235 def u_exact(self, t, u_init=None, t_init=None):
236 r"""
237 Routine to return initial conditions or exact solution.
239 Parameters
240 ----------
241 t : float
242 Time at which the exact solution is computed.
243 u_init : dtype_u
244 Initial conditions for getting the exact solution.
245 t_init : float
246 The starting time.
248 Returns
249 -------
250 u : dtype_u
251 The exact solution.
252 """
253 u = self.dtype_u(self.init)
254 u[:] = self.g(t)
255 return u
257 def eval_f(self, u, t):
258 """
259 Routine to evaluate the right-hand side of the problem.
261 Parameters
262 ----------
263 u : dtype_u
264 Current values of the numerical solution.
265 t : float
266 Current time of the numerical solution is computed (not used here).
268 Returns
269 -------
270 f : dtype_f
271 The right-hand side of the problem (one component).
272 """
274 f = self.dtype_f(self.init)
275 f[:] = self.f(u, t)
276 self.work_counters['rhs']()
277 return f
279 def solve_system(self, rhs, dt, u0, t):
280 """
281 Simple Newton solver for the nonlinear equation
283 Parameters
284 ----------
285 rhs : dtype_f
286 Right-hand side for the nonlinear system.
287 dt : float
288 Abbrev. for the node-to-node stepsize (or any other factor required).
289 u0 : dtype_u
290 Initial guess for the iterative solver.
291 t : float
292 Time of the updated solution (e.g. for time-dependent BCs).
294 Returns
295 -------
296 u : dtype_u
297 The solution as mesh.
298 """
299 # create new mesh object from u0 and set initial values for iteration
300 u = self.dtype_u(u0)
302 # start newton iteration
303 n, res = 0, np.inf
304 while n < self.newton_maxiter:
305 # form the function g with g(u) = 0
306 g = u - dt * self.f(u, t) - rhs
308 # if g is close to 0, then we are done
309 res = np.linalg.norm(g, np.inf)
310 if res < self.newton_tol or np.isnan(res):
311 break
313 # assemble dg/du
314 dg = 1 - dt * self.jac(u, t)
316 # newton update: u1 = u0 - g/dg
317 u -= dg ** (-1) * g
319 # increase iteration count and work counter
320 n += 1
321 self.work_counters['newton']()
323 if np.isnan(res) and self.stop_at_nan:
324 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
325 elif np.isnan(res): # pragma: no cover
326 self.logger.warning('Newton got nan after %i iterations...' % n)
328 if n == self.newton_maxiter:
329 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
331 return u