Coverage for pySDC/implementations/problem_classes/Van_der_Pol_implicit.py: 96%
57 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 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['rhs'] = WorkCounter()
74 def u_exact(self, t, u_init=None, t_init=None):
75 r"""
76 Routine to approximate the exact solution at time t by ``SciPy`` or give initial conditions when called at :math:`t=0`.
78 Parameters
79 ----------
80 t : float
81 Current time.
82 u_init : pySDC.problem.vanderpol.dtype_u
83 Initial conditions for getting the exact solution.
84 t_init : float
85 The starting time.
87 Returns
88 -------
89 me : dtype_u
90 Approximate exact solution.
91 """
93 me = self.dtype_u(self.init)
95 if t > 0.0:
97 def eval_rhs(t, u):
98 return self.eval_f(u, t)
100 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
101 else:
102 me[:] = self.u0
103 return me
105 def eval_f(self, u, t):
106 """
107 Routine to compute the right-hand side for both components simultaneously.
109 Parameters
110 ----------
111 u : dtype_u
112 Current values of the numerical solution.
113 t : float
114 Current time at which the numerical solution is computed (not used here).
116 Returns
117 -------
118 f : dtype_f
119 The right-hand side (contains 2 components).
120 """
122 x1 = u[0]
123 x2 = u[1]
124 f = self.f_init
125 f[0] = x2
126 f[1] = self.mu * (1 - x1**2) * x2 - x1
127 self.work_counters['rhs']()
128 return f
130 def solve_system(self, rhs, dt, u0, t):
131 """
132 Simple Newton solver for the nonlinear system.
134 Parameters
135 ----------
136 rhs : dtype_f
137 Right-hand side for the nonlinear system.
138 dt : float
139 Abbrev. for the node-to-node stepsize (or any other factor required).
140 u0 : dtype_u
141 Initial guess for the iterative solver.
142 t : float
143 Current time (e.g. for time-dependent BCs).
145 Returns
146 -------
147 u : dtype_u
148 The solution u.
149 """
151 mu = self.mu
153 # create new mesh object from u0 and set initial values for iteration
154 u = self.dtype_u(u0)
155 x1 = u[0]
156 x2 = u[1]
158 # start newton iteration
159 n = 0
160 res = 99
161 while n < self.newton_maxiter:
162 # form the function g with g(u) = 0
163 g = np.array([x1 - dt * x2 - rhs[0], x2 - dt * (mu * (1 - x1**2) * x2 - x1) - rhs[1]])
165 # if g is close to 0, then we are done
166 res = np.linalg.norm(g, np.inf) / (abs(u) if self.relative_tolerance else 1.0)
167 if res < self.newton_tol or np.isnan(res):
168 break
170 # prefactor for dg/du
171 c = 1.0 / (-2 * dt**2 * mu * x1 * x2 - dt**2 - 1 + dt * mu * (1 - x1**2))
172 # assemble dg/du
173 dg = c * np.array([[dt * mu * (1 - x1**2) - 1, -dt], [2 * dt * mu * x1 * x2 + dt, -1]])
175 # newton update: u1 = u0 - g/dg
176 u -= np.dot(dg, g)
178 # set new values and increase iteration count
179 x1 = u[0]
180 x2 = u[1]
181 n += 1
182 self.work_counters['newton']()
184 if np.isnan(res) and self.stop_at_nan:
185 self.logger.warning('Newton got nan after %i iterations...' % n)
186 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
187 elif np.isnan(res):
188 self.logger.warning('Newton got nan after %i iterations...' % n)
190 if n == self.newton_maxiter and self.crash_at_maxiter:
191 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
193 return u