# Coverage for pySDC/implementations/problem_classes/Van_der_Pol_implicit.py: 96%

## 57 statements

, created at 2024-09-20 17:10 +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')))

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