Coverage for pySDC/implementations/problem_classes/Lorenz.py: 94%
54 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-12 12:08 +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 Lorenz system of three chaotic ODEs, treated fully implicitly with a Newton solver.
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.
54 stop_at_nan : bool, optional
55 Indicates that the Newton solver should raise an error if ``nan`` values arise, instead of only warning.
57 Attributes
58 ----------
59 work_counter : dict
60 Counts the iterations/nfev (here for Newton's method and the nfev for the right-hand side).
61 """
63 dtype_u = mesh
64 dtype_f = mesh
66 def __init__(
67 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
68 ):
69 """Initialization routine"""
71 nvars = 3
72 # invoke super init, passing number of dofs, dtype_u and dtype_f
73 super().__init__(init=(nvars, None, np.dtype('float64')))
74 self._makeAttributeAndRegister('nvars', 'stop_at_nan', localVars=locals(), readOnly=True)
75 self._makeAttributeAndRegister(
76 'sigma', 'rho', 'beta', 'u0', 'newton_tol', 'newton_maxiter', localVars=locals(), readOnly=False
77 )
78 self.work_counters['newton'] = WorkCounter()
79 self.work_counters['rhs'] = WorkCounter()
81 def eval_f(self, u, t):
82 """
83 Routine to evaluate the right-hand side of the problem.
85 Parameters
86 ----------
87 u : dtype_u
88 Current values of the numerical solution.
89 t : float
90 Current time of the numerical solution is computed.
92 Returns
93 -------
94 f : dtype_f
95 The right-hand side of the problem.
96 """
97 # abbreviations
98 sigma = self.sigma
99 rho = self.rho
100 beta = self.beta
102 f = self.dtype_f(self.init)
104 f[0] = sigma * (u[1] - u[0])
105 f[1] = rho * u[0] - u[1] - u[0] * u[2]
106 f[2] = u[0] * u[1] - beta * u[2]
108 self.work_counters['rhs']()
109 return f
111 def solve_system(self, rhs, dt, u0, t):
112 """
113 Simple Newton solver for the nonlinear system.
115 Parameters
116 ----------
117 rhs : dtype_f
118 Right-hand side for the nonlinear system.
119 dt : float
120 Abbrev. for the local stepsize (or any other factor required).
121 u0 : dtype_u
122 Initial guess for the iterative solver
123 t : float
124 Current time (e.g. for time-dependent BCs).
126 Returns
127 -------
128 me : dtype_u
129 The solution as mesh.
130 """
132 # abbreviations
133 sigma = self.sigma
134 rho = self.rho
135 beta = self.beta
137 # start Newton iterations
138 u = self.dtype_u(u0)
139 res = np.inf
140 for _n in range(0, self.newton_maxiter):
141 # assemble G such that G(u) = 0 at the solution to the step
142 G = np.array(
143 [
144 u[0] - dt * sigma * (u[1] - u[0]) - rhs[0],
145 u[1] - dt * (rho * u[0] - u[1] - u[0] * u[2]) - rhs[1],
146 u[2] - dt * (u[0] * u[1] - beta * u[2]) - rhs[2],
147 ]
148 )
150 # compute the residual and determine if we are done
151 res = np.linalg.norm(G, np.inf)
152 if res <= self.newton_tol:
153 break
154 if np.isnan(res) and self.stop_at_nan:
155 self.logger.warning('Newton got nan after %i iterations...' % _n)
156 raise ConvergenceError('Newton got nan after %i iterations, aborting...' % _n)
157 elif np.isnan(res):
158 self.logger.warning('Newton got nan after %i iterations...' % _n)
159 break
161 # assemble inverse of Jacobian J of G
162 prefactor = 1.0 / (
163 dt**3 * sigma * (u[0] ** 2 + u[0] * u[1] + beta * (-rho + u[2] + 1))
164 + dt**2 * (beta * sigma + beta - rho * sigma + sigma + u[0] ** 2 + sigma * u[2])
165 + dt * (beta + sigma + 1)
166 + 1
167 )
168 J_inv = prefactor * np.array(
169 [
170 [
171 beta * dt**2 + dt**2 * u[0] ** 2 + beta * dt + dt + 1,
172 beta * dt**2 * sigma + dt * sigma,
173 -(dt**2) * sigma * u[0],
174 ],
175 [
176 beta * dt**2 * rho + dt**2 * (-u[0]) * u[1] - beta * dt**2 * u[2] + dt * rho - dt * u[2],
177 beta * dt**2 * sigma + beta * dt + dt * sigma + 1,
178 dt**2 * sigma * (-u[0]) - dt * u[0],
179 ],
180 [
181 dt**2 * rho * u[0] - dt**2 * u[0] * u[2] + dt**2 * u[1] + dt * u[1],
182 dt**2 * sigma * u[0] + dt**2 * sigma * u[1] + dt * u[0],
183 -(dt**2) * rho * sigma + dt**2 * sigma + dt**2 * sigma * u[2] + dt * sigma + dt + 1,
184 ],
185 ]
186 )
188 # solve the linear system for the Newton correction J delta = G
189 delta = J_inv @ G
191 # update solution
192 u = u - delta
193 self.work_counters['newton']()
195 return u
197 def u_exact(self, t, u_init=None, t_init=None):
198 r"""
199 Routine to return initial conditions or to approximate exact solution using ``SciPy``.
201 Parameters
202 ----------
203 t : float
204 Time at which the approximated exact solution is computed.
205 u_init : pySDC.implementations.problem_classes.Lorenz.dtype_u
206 Initial conditions for getting the exact solution.
207 t_init : float
208 The starting time.
210 Returns
211 -------
212 me : dtype_u
213 The approximated exact solution.
214 """
216 me = self.dtype_u(self.init)
218 if t > 0:
220 def eval_rhs(t, u):
221 r"""
222 Evaluate the right hand side, but switch the arguments for ``SciPy``.
224 Args:
225 t (float): Time
226 u (numpy.ndarray): Solution at time t
228 Returns:
229 (numpy.ndarray): Right hand side evaluation
230 """
231 return self.eval_f(u, t)
233 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
234 else:
235 me[:] = self.u0
236 return me