Coverage for pySDC/implementations/problem_classes/Van_der_Pol_implicit.py: 98%
64 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-21 06:49 +0000
1import numpy as np
3from pySDC.core.errors import ProblemError
4from pySDC.core.problem import Problem, WorkCounter
5from pySDC.implementations.datatype_classes.mesh import mesh
8# noinspection PyUnusedLocal
9class vanderpol(Problem):
10 r"""
11 This class implements the stiff Van der Pol oscillator given by the equation
13 .. math::
14 \frac{d^2 u(t)}{d t^2} - \mu (1 - u(t)^2) \frac{d u(t)}{dt} + u(t) = 0.
16 Parameters
17 ----------
18 u0 : sequence of array_like, optional
19 Initial condition.
20 mu : float, optional
21 Stiff parameter :math:`\mu`.
22 newton_maxiter : int, optional
23 Maximum number of iterations for Newton's method to terminate.
24 newton_tol : float, optional
25 Tolerance for Newton to terminate.
26 stop_at_nan : bool, optional
27 Indicate whether Newton's method should stop if ``nan`` values arise.
28 crash_at_maxiter : bool, optional
29 Indicates whether Newton's method should stop if maximum number of iterations
30 ``newton_maxiter`` is reached.
31 relative_tolerance : bool, optional
32 Use a relative or absolute tolerance for the Newton solver
34 Attributes
35 ----------
36 work_counters : WorkCounter
37 Counts different things, here: Number of evaluations of the right-hand side in ``eval_f``
38 and number of Newton calls in each Newton iterations are counted.
39 """
41 dtype_u = mesh
42 dtype_f = mesh
44 def __init__(
45 self,
46 u0=None,
47 mu=5.0,
48 newton_maxiter=100,
49 newton_tol=1e-9,
50 stop_at_nan=True,
51 crash_at_maxiter=True,
52 relative_tolerance=False,
53 ):
54 """Initialization routine"""
55 nvars = 2
57 if u0 is None:
58 u0 = [2.0, 0.0]
60 super().__init__((nvars, None, np.dtype('float64')))
61 self._makeAttributeAndRegister('nvars', 'u0', localVars=locals(), readOnly=True)
62 self._makeAttributeAndRegister(
63 'mu',
64 'newton_maxiter',
65 'newton_tol',
66 'stop_at_nan',
67 'crash_at_maxiter',
68 'relative_tolerance',
69 localVars=locals(),
70 )
71 self.work_counters['newton'] = WorkCounter()
72 self.work_counters['jacobian_solves'] = WorkCounter()
73 self.work_counters['rhs'] = WorkCounter()
75 def u_exact(self, t, u_init=None, t_init=None):
76 r"""
77 Routine to approximate the exact solution at time t by ``SciPy`` or give initial conditions when called at :math:`t=0`.
79 Parameters
80 ----------
81 t : float
82 Current time.
83 u_init : pySDC.problem.vanderpol.dtype_u
84 Initial conditions for getting the exact solution.
85 t_init : float
86 The starting time.
88 Returns
89 -------
90 me : dtype_u
91 Approximate exact solution.
92 """
94 me = self.dtype_u(self.init)
96 if t > 0.0:
98 def eval_rhs(t, u):
99 return self.eval_f(u, t)
101 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
102 else:
103 me[:] = self.u0
104 return me
106 def eval_f(self, u, t):
107 """
108 Routine to compute the right-hand side for both components simultaneously.
110 Parameters
111 ----------
112 u : dtype_u
113 Current values of the numerical solution.
114 t : float
115 Current time at which the numerical solution is computed (not used here).
117 Returns
118 -------
119 f : dtype_f
120 The right-hand side (contains 2 components).
121 """
123 x1 = u[0]
124 x2 = u[1]
125 f = self.f_init
126 f[0] = x2
127 f[1] = self.mu * (1 - x1**2) * x2 - x1
128 self.work_counters['rhs']()
129 return f
131 def solve_system(self, rhs, dt, u0, t):
132 """
133 Simple Newton solver for the nonlinear system.
135 Parameters
136 ----------
137 rhs : dtype_f
138 Right-hand side for the nonlinear system.
139 dt : float
140 Abbrev. for the node-to-node stepsize (or any other factor required).
141 u0 : dtype_u
142 Initial guess for the iterative solver.
143 t : float
144 Current time (e.g. for time-dependent BCs).
146 Returns
147 -------
148 u : dtype_u
149 The solution u.
150 """
152 mu = self.mu
154 # create new mesh object from u0 and set initial values for iteration
155 u = self.dtype_u(u0)
156 x1 = u[0]
157 x2 = u[1]
159 # start newton iteration
160 n = 0
161 res = 99
162 while n < self.newton_maxiter:
163 # form the function g with g(u) = 0
164 g = np.array([x1 - dt * x2 - rhs[0], x2 - dt * (mu * (1 - x1**2) * x2 - x1) - rhs[1]])
166 # if g is close to 0, then we are done
167 res = np.linalg.norm(g, np.inf) / (abs(u) if self.relative_tolerance else 1.0)
168 if res < self.newton_tol or np.isnan(res):
169 break
171 u -= self.solve_jacobian(g, dt, u)
173 # set new values and increase iteration count
174 x1 = u[0]
175 x2 = u[1]
176 n += 1
177 self.work_counters['newton']()
179 if np.isnan(res) and self.stop_at_nan:
180 self.logger.warning('Newton got nan after %i iterations...' % n)
181 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
182 elif np.isnan(res):
183 self.logger.warning('Newton got nan after %i iterations...' % n)
185 if n == self.newton_maxiter and self.crash_at_maxiter:
186 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
188 return u
190 def solve_jacobian(self, rhs, dt, u, **kwargs):
191 mu = self.mu
192 u1 = u[0]
193 u2 = u[1]
195 # assemble prefactor
196 c = 1.0 / (-2 * dt**2 * mu * u1 * u2 - dt**2 - 1 + dt * mu * (1 - u1**2))
197 # assemble (dg/du)^-1
198 dg = c * np.array([[dt * mu * (1 - u1**2) - 1, -dt], [2 * dt * mu * u1 * u2 + dt, -1]])
200 self.work_counters['jacobian_solves']()
201 return np.dot(dg, rhs)