Coverage for pySDC / implementations / problem_classes / odeScalar.py: 100%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Implementation of scalar test problem ODEs. 

5 

6 

7Reference : 

8 

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""" 

13 

14import numpy as np 

15 

16from pySDC.core.errors import ProblemError 

17from pySDC.core.problem import Problem, WorkCounter 

18from pySDC.implementations.datatype_classes.mesh import mesh 

19 

20 

21class ProtheroRobinson(Problem): 

22 r""" 

23 Implement the Prothero-Robinson problem: 

24 

25 .. math:: 

26 \frac{du}{dt} = -\frac{u-g(t)}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0)., 

27 

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)`. 

32 

33 Implement also the non-linear form of this problem: 

34 

35 .. math:: 

36 \frac{du}{dt} = -\frac{u^3-g(t)^3}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0). 

37 

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: 

41 

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) 

49 

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. 

62 

63 Reference 

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 """ 

69 

70 dtype_u = mesh 

71 dtype_f = mesh 

72 

73 def __init__(self, epsilon=1e-3, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True): 

74 nvars = 1 

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

76 

77 self.f = self.f_NONLIN if nonLinear else self.f_LIN 

78 self.jac = self.jac_NONLIN if nonLinear else self.jac_LIN 

79 self._makeAttributeAndRegister( 

80 'epsilon', 'nonLinear', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True 

81 ) 

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

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

84 

85 # ------------------------------------------------------------------------- 

86 # g function (analytical solution), and its first derivative 

87 # ------------------------------------------------------------------------- 

88 def g(self, t): 

89 return np.cos(t) 

90 

91 def dg(self, t): 

92 return -np.sin(t) 

93 

94 # ------------------------------------------------------------------------- 

95 # f(u,t) and Jacobian functions 

96 # ------------------------------------------------------------------------- 

97 def f(self, u, t): 

98 raise NotImplementedError() 

99 

100 def f_LIN(self, u, t): 

101 return -self.epsilon ** (-1) * (u - self.g(t)) + self.dg(t) 

102 

103 def f_NONLIN(self, u, t): 

104 return -self.epsilon ** (-1) * (u**3 - self.g(t) ** 3) + self.dg(t) 

105 

106 def jac(self, u, t): 

107 raise NotImplementedError() 

108 

109 def jac_LIN(self, u, t): 

110 return -self.epsilon ** (-1) 

111 

112 def jac_NONLIN(self, u, t): 

113 return -self.epsilon ** (-1) * 3 * u**2 

114 

115 # ------------------------------------------------------------------------- 

116 # pySDC required methods 

117 # ------------------------------------------------------------------------- 

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

119 r""" 

120 Routine to return initial conditions or exact solution. 

121 

122 Parameters 

123 ---------- 

124 t : float 

125 Time at which the exact solution is computed. 

126 u_init : dtype_u 

127 Initial conditions for getting the exact solution. 

128 t_init : float 

129 The starting time. 

130 

131 Returns 

132 ------- 

133 u : dtype_u 

134 The exact solution. 

135 """ 

136 u = self.dtype_u(self.init) 

137 u[:] = self.g(t) 

138 return u 

139 

140 def eval_f(self, u, t): 

141 """ 

142 Routine to evaluate the right-hand side of the problem. 

143 

144 Parameters 

145 ---------- 

146 u : dtype_u 

147 Current values of the numerical solution. 

148 t : float 

149 Current time of the numerical solution is computed (not used here). 

150 

151 Returns 

152 ------- 

153 f : dtype_f 

154 The right-hand side of the problem (one component). 

155 """ 

156 

157 f = self.dtype_f(self.init) 

158 f[:] = self.f(u, t) 

159 self.work_counters['rhs']() 

160 return f 

161 

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

163 """ 

164 Simple Newton solver for the nonlinear equation 

165 

166 Parameters 

167 ---------- 

168 rhs : dtype_f 

169 Right-hand side for the nonlinear system. 

170 dt : float 

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

172 u0 : dtype_u 

173 Initial guess for the iterative solver. 

174 t : float 

175 Time of the updated solution (e.g. for time-dependent BCs). 

176 

177 Returns 

178 ------- 

179 u : dtype_u 

180 The solution as mesh. 

181 """ 

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

183 u = self.dtype_u(u0) 

184 

185 # start newton iteration 

186 n, res = 0, np.inf 

187 while n < self.newton_maxiter: 

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

189 g = u - dt * self.f(u, t) - rhs 

190 

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

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

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

194 break 

195 

196 # assemble dg/du 

197 dg = 1 - dt * self.jac(u, t) 

198 

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

200 u -= dg ** (-1) * g 

201 

202 # increase iteration count and work counter 

203 n += 1 

204 self.work_counters['newton']() 

205 

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

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

208 elif np.isnan(res): # pragma: no cover 

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

210 

211 if n == self.newton_maxiter: 

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

213 

214 return u