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

## 54 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

2from pySDC.core.problem import Problem, WorkCounter

3from pySDC.implementations.datatype_classes.mesh import mesh

4from pySDC.core.errors import ConvergenceError

7class LorenzAttractor(Problem):

8 r"""

9 Simple script to run a Lorenz attractor problem.

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.

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.

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

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

24 The system of ODEs is given by

26 .. math::

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

29 .. math::

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

32 .. math::

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

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.

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.

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

61 dtype_u = mesh

62 dtype_f = mesh

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

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')))

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()

79 def eval_f(self, u, t):

80 """

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

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.

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

100 f = self.dtype_f(self.init)

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]

106 self.work_counters['rhs']()

107 return f

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

110 """

111 Simple Newton solver for the nonlinear system.

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

124 Returns

125 -------

126 me : dtype_u

127 The solution as mesh.

128 """

130 # abbreviations

131 sigma = self.sigma

132 rho = self.rho

133 beta = self.beta

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 )

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

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 )

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

187 delta = J_inv @ G

189 # update solution

190 u = u - delta

191 self.work_counters['newton']()

193 return u

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.

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.

208 Returns

209 -------

210 me : dtype_u

211 The approximated exact solution.

212 """

214 me = self.dtype_u(self.init)

216 if t > 0:

218 def eval_rhs(t, u):

219 r"""

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

222 Args:

223 t (float): Time

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

226 Returns:

227 (numpy.ndarray): Right hand side evaluation

228 """

229 return self.eval_f(u, t)

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

232 else:

233 me[:] = self.u0

234 return me