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

1import numpy as np 

2 

3from pySDC.core.Errors import ProblemError 

4from pySDC.core.Problem import ptype, WorkCounter 

5from pySDC.implementations.datatype_classes.mesh import mesh 

6 

7 

8# noinspection PyUnusedLocal 

9class vanderpol(ptype): 

10 r""" 

11 This class implements the stiff Van der Pol oscillator given by the equation 

12 

13 .. math:: 

14 \frac{d^2 u(t)}{d t^2} - \mu (1 - u(t)^2) \frac{d u(t)}{dt} + u(t) = 0. 

15 

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 

32 Attributes 

33 ---------- 

34 work_counters : WorkCounter 

35 Counts different things, here: Number of evaluations of the right-hand side in ``eval_f`` 

36 and number of Newton calls in each Newton iterations are counted. 

37 """ 

38 

39 dtype_u = mesh 

40 dtype_f = mesh 

41 

42 def __init__(self, u0=None, mu=5.0, newton_maxiter=100, newton_tol=1e-9, stop_at_nan=True, crash_at_maxiter=True): 

43 """Initialization routine""" 

44 nvars = 2 

45 

46 if u0 is None: 

47 u0 = [2.0, 0.0] 

48 

49 super().__init__((nvars, None, np.dtype('float64'))) 

50 self._makeAttributeAndRegister('nvars', 'u0', localVars=locals(), readOnly=True) 

51 self._makeAttributeAndRegister( 

52 'mu', 'newton_maxiter', 'newton_tol', 'stop_at_nan', 'crash_at_maxiter', localVars=locals() 

53 ) 

54 self.work_counters['newton'] = WorkCounter() 

55 self.work_counters['rhs'] = WorkCounter() 

56 

57 def u_exact(self, t, u_init=None, t_init=None): 

58 r""" 

59 Routine to approximate the exact solution at time t by ``SciPy`` or give initial conditions when called at :math:`t=0`. 

60 

61 Parameters 

62 ---------- 

63 t : float 

64 Current time. 

65 u_init : pySDC.problem.vanderpol.dtype_u 

66 Initial conditions for getting the exact solution. 

67 t_init : float 

68 The starting time. 

69 

70 Returns 

71 ------- 

72 me : dtype_u 

73 Approximate exact solution. 

74 """ 

75 

76 me = self.dtype_u(self.init) 

77 

78 if t > 0.0: 

79 

80 def eval_rhs(t, u): 

81 return self.eval_f(u, t) 

82 

83 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init) 

84 else: 

85 me[:] = self.u0 

86 return me 

87 

88 def eval_f(self, u, t): 

89 """ 

90 Routine to compute the right-hand side for both components simultaneously. 

91 

92 Parameters 

93 ---------- 

94 u : dtype_u 

95 Current values of the numerical solution. 

96 t : float 

97 Current time at which the numerical solution is computed (not used here). 

98 

99 Returns 

100 ------- 

101 f : dtype_f 

102 The right-hand side (contains 2 components). 

103 """ 

104 

105 x1 = u[0] 

106 x2 = u[1] 

107 f = self.f_init 

108 f[0] = x2 

109 f[1] = self.mu * (1 - x1**2) * x2 - x1 

110 self.work_counters['rhs']() 

111 return f 

112 

113 def solve_system(self, rhs, dt, u0, t): 

114 """ 

115 Simple Newton solver for the nonlinear system. 

116 

117 Parameters 

118 ---------- 

119 rhs : dtype_f 

120 Right-hand side for the nonlinear system. 

121 dt : float 

122 Abbrev. for the node-to-node stepsize (or any other factor required). 

123 u0 : dtype_u 

124 Initial guess for the iterative solver. 

125 t : float 

126 Current time (e.g. for time-dependent BCs). 

127 

128 Returns 

129 ------- 

130 u : dtype_u 

131 The solution u. 

132 """ 

133 

134 mu = self.mu 

135 

136 # create new mesh object from u0 and set initial values for iteration 

137 u = self.dtype_u(u0) 

138 x1 = u[0] 

139 x2 = u[1] 

140 

141 # start newton iteration 

142 n = 0 

143 res = 99 

144 while n < self.newton_maxiter: 

145 # form the function g with g(u) = 0 

146 g = np.array([x1 - dt * x2 - rhs[0], x2 - dt * (mu * (1 - x1**2) * x2 - x1) - rhs[1]]) 

147 

148 # if g is close to 0, then we are done 

149 res = np.linalg.norm(g, np.inf) 

150 if res < self.newton_tol or np.isnan(res): 

151 break 

152 

153 # prefactor for dg/du 

154 c = 1.0 / (-2 * dt**2 * mu * x1 * x2 - dt**2 - 1 + dt * mu * (1 - x1**2)) 

155 # assemble dg/du 

156 dg = c * np.array([[dt * mu * (1 - x1**2) - 1, -dt], [2 * dt * mu * x1 * x2 + dt, -1]]) 

157 

158 # newton update: u1 = u0 - g/dg 

159 u -= np.dot(dg, g) 

160 

161 # set new values and increase iteration count 

162 x1 = u[0] 

163 x2 = u[1] 

164 n += 1 

165 self.work_counters['newton']() 

166 

167 if np.isnan(res) and self.stop_at_nan: 

168 self.logger.warning('Newton got nan after %i iterations...' % n) 

169 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

170 elif np.isnan(res): 

171 self.logger.warning('Newton got nan after %i iterations...' % n) 

172 

173 if n == self.newton_maxiter and self.crash_at_maxiter: 

174 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

175 

176 return u