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

13import numpy as np 

14 

15from pySDC.core.errors import ProblemError 

16from pySDC.core.problem import Problem, WorkCounter 

17from pySDC.implementations.datatype_classes.mesh import mesh 

18 

19 

20class ProtheroRobinson(Problem): 

21 r""" 

22 Implement the Prothero-Robinson problem: 

23 

24 .. math:: 

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

26 

27 with :math:`\epsilon` a stiffness parameter, that makes the problem more stiff 

28 the smaller it is (usual taken value is :math:`\epsilon=1e^{-3}`). 

29 Exact solution is given by :math:`u(t)=g(t)`, and this implementation uses 

30 :math:`g(t)=\cos(t)`. 

31 

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

33 

34 .. math:: 

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

36 

37 To use an other exact solution, one just have to derivate this class 

38 and overload the `g` and `dg` methods. For instance, 

39 to use :math:`g(t)=e^{-0.2*t}`, define and use the following class: 

40 

41 >>> class MyProtheroRobinson(ProtheroRobinson): 

42 >>> 

43 >>> def g(self, t): 

44 >>> return np.exp(-0.2 * t) 

45 >>> 

46 >>> def dg(self, t): 

47 >>> return (-0.2) * np.exp(-0.2 * t) 

48 

49 Parameters 

50 ---------- 

51 epsilon : float, optional 

52 Stiffness parameter. The default is 1e-3. 

53 nonLinear : bool, optional 

54 Wether or not to use the non-linear form of the problem. The default is False. 

55 newton_maxiter : int, optional 

56 Maximum number of Newton iteration in solve_system. The default is 200. 

57 newton_tol : float, optional 

58 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11. 

59 stop_at_nan : bool, optional 

60 Wheter to stop or not solve_system when getting NAN. The default is True. 

61 

62 Reference 

63 --------- 

64 A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving 

65 stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974), 

66 pp. 145–162. 

67 """ 

68 

69 dtype_u = mesh 

70 dtype_f = mesh 

71 

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

73 nvars = 1 

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

75 

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

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

78 self._makeAttributeAndRegister( 

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

80 ) 

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

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

83 

84 # ------------------------------------------------------------------------- 

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

86 # ------------------------------------------------------------------------- 

87 def g(self, t): 

88 return np.cos(t) 

89 

90 def dg(self, t): 

91 return -np.sin(t) 

92 

93 # ------------------------------------------------------------------------- 

94 # f(u,t) and Jacobian functions 

95 # ------------------------------------------------------------------------- 

96 def f(self, u, t): 

97 raise NotImplementedError() 

98 

99 def f_LIN(self, u, t): 

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

101 

102 def f_NONLIN(self, u, t): 

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

104 

105 def jac(self, u, t): 

106 raise NotImplementedError() 

107 

108 def jac_LIN(self, u, t): 

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

110 

111 def jac_NONLIN(self, u, t): 

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

113 

114 # ------------------------------------------------------------------------- 

115 # pySDC required methods 

116 # ------------------------------------------------------------------------- 

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

118 r""" 

119 Routine to return initial conditions or exact solution. 

120 

121 Parameters 

122 ---------- 

123 t : float 

124 Time at which the exact solution is computed. 

125 u_init : dtype_u 

126 Initial conditions for getting the exact solution. 

127 t_init : float 

128 The starting time. 

129 

130 Returns 

131 ------- 

132 u : dtype_u 

133 The exact solution. 

134 """ 

135 u = self.dtype_u(self.init) 

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

137 return u 

138 

139 def eval_f(self, u, t): 

140 """ 

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

142 

143 Parameters 

144 ---------- 

145 u : dtype_u 

146 Current values of the numerical solution. 

147 t : float 

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

149 

150 Returns 

151 ------- 

152 f : dtype_f 

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

154 """ 

155 

156 f = self.dtype_f(self.init) 

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

158 self.work_counters['rhs']() 

159 return f 

160 

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

162 """ 

163 Simple Newton solver for the nonlinear equation 

164 

165 Parameters 

166 ---------- 

167 rhs : dtype_f 

168 Right-hand side for the nonlinear system. 

169 dt : float 

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

171 u0 : dtype_u 

172 Initial guess for the iterative solver. 

173 t : float 

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

175 

176 Returns 

177 ------- 

178 u : dtype_u 

179 The solution as mesh. 

180 """ 

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

182 u = self.dtype_u(u0) 

183 

184 # start newton iteration 

185 n, res = 0, np.inf 

186 while n < self.newton_maxiter: 

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

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

189 

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

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

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

193 break 

194 

195 # assemble dg/du 

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

197 

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

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

200 

201 # increase iteration count and work counter 

202 n += 1 

203 self.work_counters['newton']() 

204 

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

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

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

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

209 

210 if n == self.newton_maxiter: 

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

212 

213 return u