Coverage for pySDC/implementations/problem_classes/Lorenz.py: 94%
54 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +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')))
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()
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