Coverage for pySDC/implementations/problem_classes/Lorenz.py: 94%

1import numpy as np 

2from pySDC.core.Problem import ptype, WorkCounter 

3from pySDC.implementations.datatype_classes.mesh import mesh 

4from pySDC.core.Errors import ConvergenceError 

5 

6 

7class LorenzAttractor(ptype): 

8 r""" 

9 Simple script to run a Lorenz attractor problem. 

10 

11 The Lorenz attractor is a system of three ordinary differential equations (ODEs) that exhibits some chaotic behaviour. 

12 It is well known for the "Butterfly Effect", because the solution looks like a butterfly (solve to :math:`T_{end} = 100` 

13 or so to see this with these initial conditions) and because of the chaotic nature. 

14 

15 Lorenz developed this system from equations modelling convection in a layer of fluid with the top and bottom surfaces 

16 kept at different temperatures. In the notation used here, the first component of u is proportional to the convective 

17 motion, the second component is proportional to the temperature difference between the surfaces and the third component 

18 is proportional to the distortion of the vertical temperature profile from linearity. 

19 

20 See doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2 for the original publication. 

21 

22 Since the problem is non-linear, we need to use a Newton solver. 

23 

24 The system of ODEs is given by 

25 

26 .. math:: 

27 \frac{d y_1(t)}{dt} = \sigma (y_2 (t) - y_1 (t)), 

28 

29 .. math:: 

30 \frac{d y_2(t)}{dt} = \rho y_1 (t) - y_2 (t) - y_1 (t) y_3 (t), 

31 

32 .. math:: 

33 \frac{d y_3(t)}{dt} = y_1 (t) y_2 (t) - \beta y_3 (t) 

34 

35 with initial condition :math:`(y_1(0), y_2(0), y_3(0))^T = (1, 1, 1)^T` (default) for :math:`t \in [0, 1]`. 

36 The problem parameters for this problem are :math:`\sigma = 10`, :math:`\rho = 28` and :math:`\beta = 8/3`. 

37 Lorenz chose these parameters such that the Reynolds number :math:`\rho` is slightly supercritical 

38 as to provoke instability of steady convection. 

39 

40 Parameters 

41 ---------- 

42 sigma : float, optional 

43 Parameter :math:`\sigma` of the problem. 

44 rho : float, optional 

45 Parameter :math:`\rho` of the problem. 

46 beta : float, optional 

47 Parameter :math:`\beta` of the problem. 

48 u0 : tuple, optional 

49 Initial solution :math:`u_0` of the problem. 

50 newton_tol : float, optional 

51 Tolerance for Newton for termination. 

52 newton_maxiter : int, optional 

53 Maximum number of iterations for Newton's method. 

54 

55 Attributes 

56 ---------- 

57 work_counter : dict 

58 Counts the iterations/nfev (here for Newton's method and the nfev for the right-hand side). 

59 """ 

60 

61 dtype_u = mesh 

62 dtype_f = mesh 

63 

64 def __init__( 

65 self, sigma=10.0, rho=28.0, beta=8.0 / 3.0, u0=(1, 1, 1), newton_tol=1e-9, newton_maxiter=99, stop_at_nan=True 

66 ): 

67 """Initialization routine""" 

68 

69 nvars = 3 

70 # invoke super init, passing number of dofs, dtype_u and dtype_f 

71 super().__init__(init=(nvars, None, np.dtype('float64'))) 

72 self._makeAttributeAndRegister('nvars', 'stop_at_nan', localVars=locals(), readOnly=True) 

73 self._makeAttributeAndRegister( 

74 'sigma', 'rho', 'beta', 'u0', 'newton_tol', 'newton_maxiter', localVars=locals(), readOnly=False 

75 ) 

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

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

78 

79 def eval_f(self, u, t): 

80 """ 

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

82 

83 Parameters 

84 ---------- 

85 u : dtype_u 

86 Current values of the numerical solution. 

87 t : float 

88 Current time of the numerical solution is computed. 

89 

90 Returns 

91 ------- 

92 f : dtype_f 

93 The right-hand side of the problem. 

94 """ 

95 # abbreviations 

96 sigma = self.sigma 

97 rho = self.rho 

98 beta = self.beta 

99 

100 f = self.dtype_f(self.init) 

101 

102 f[0] = sigma * (u[1] - u[0]) 

103 f[1] = rho * u[0] - u[1] - u[0] * u[2] 

104 f[2] = u[0] * u[1] - beta * u[2] 

105 

106 self.work_counters['rhs']() 

107 return f 

108 

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

110 """ 

111 Simple Newton solver for the nonlinear system. 

112 

113 Parameters 

114 ---------- 

115 rhs : dtype_f 

116 Right-hand side for the nonlinear system. 

117 factor : float 

118 Abbrev. for the local stepsize (or any other factor required). 

119 u0 : dtype_u 

120 Initial guess for the iterative solver 

121 t : float 

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

123 

124 Returns 

125 ------- 

126 me : dtype_u 

127 The solution as mesh. 

128 """ 

129 

130 # abbreviations 

131 sigma = self.sigma 

132 rho = self.rho 

133 beta = self.beta 

134 

135 # start Newton iterations 

136 u = self.dtype_u(u0) 

137 res = np.inf 

138 for _n in range(0, self.newton_maxiter): 

139 # assemble G such that G(u) = 0 at the solution to the step 

140 G = np.array( 

141 [ 

142 u[0] - dt * sigma * (u[1] - u[0]) - rhs[0], 

143 u[1] - dt * (rho * u[0] - u[1] - u[0] * u[2]) - rhs[1], 

144 u[2] - dt * (u[0] * u[1] - beta * u[2]) - rhs[2], 

145 ] 

146 ) 

147 

148 # compute the residual and determine if we are done 

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

150 if res <= self.newton_tol: 

151 break 

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

153 self.logger.warning('Newton got nan after %i iterations...' % _n) 

154 raise ConvergenceError('Newton got nan after %i iterations, aborting...' % _n) 

155 elif np.isnan(res): 

156 self.logger.warning('Newton got nan after %i iterations...' % _n) 

157 break 

158 

159 # assemble inverse of Jacobian J of G 

160 prefactor = 1.0 / ( 

161 dt**3 * sigma * (u[0] ** 2 + u[0] * u[1] + beta * (-rho + u[2] + 1)) 

162 + dt**2 * (beta * sigma + beta - rho * sigma + sigma + u[0] ** 2 + sigma * u[2]) 

163 + dt * (beta + sigma + 1) 

164 + 1 

165 ) 

166 J_inv = prefactor * np.array( 

167 [ 

168 [ 

169 beta * dt**2 + dt**2 * u[0] ** 2 + beta * dt + dt + 1, 

170 beta * dt**2 * sigma + dt * sigma, 

171 -(dt**2) * sigma * u[0], 

172 ], 

173 [ 

174 beta * dt**2 * rho + dt**2 * (-u[0]) * u[1] - beta * dt**2 * u[2] + dt * rho - dt * u[2], 

175 beta * dt**2 * sigma + beta * dt + dt * sigma + 1, 

176 dt**2 * sigma * (-u[0]) - dt * u[0], 

177 ], 

178 [ 

179 dt**2 * rho * u[0] - dt**2 * u[0] * u[2] + dt**2 * u[1] + dt * u[1], 

180 dt**2 * sigma * u[0] + dt**2 * sigma * u[1] + dt * u[0], 

181 -(dt**2) * rho * sigma + dt**2 * sigma + dt**2 * sigma * u[2] + dt * sigma + dt + 1, 

182 ], 

183 ] 

184 ) 

185 

186 # solve the linear system for the Newton correction J delta = G 

187 delta = J_inv @ G 

188 

189 # update solution 

190 u = u - delta 

191 self.work_counters['newton']() 

192 

193 return u 

194 

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

196 r""" 

197 Routine to return initial conditions or to approximate exact solution using ``SciPy``. 

198 

199 Parameters 

200 ---------- 

201 t : float 

202 Time at which the approximated exact solution is computed. 

203 u_init : pySDC.implementations.problem_classes.Lorenz.dtype_u 

204 Initial conditions for getting the exact solution. 

205 t_init : float 

206 The starting time. 

207 

208 Returns 

209 ------- 

210 me : dtype_u 

211 The approximated exact solution. 

212 """ 

213 

214 me = self.dtype_u(self.init) 

215 

216 if t > 0: 

217 

218 def eval_rhs(t, u): 

219 r""" 

220 Evaluate the right hand side, but switch the arguments for ``SciPy``. 

221 

222 Args: 

223 t (float): Time 

224 u (numpy.ndarray): Solution at time t 

225 

226 Returns: 

227 (numpy.ndarray): Right hand side evaluation 

228 """ 

229 return self.eval_f(u, t) 

230 

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

232 else: 

233 me[:] = self.u0 

234 return me