Coverage for pySDC/implementations/problem_classes/Burgers.py: 100%
96 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 10:50 +0000
1import numpy as np
3from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
4from pySDC.implementations.problem_classes.generic_spectral import GenericSpectralLinear
7class Burgers1D(GenericSpectralLinear):
8 """
9 1D viscous Burgers equation with Dirichlet BCs and a Chebychev method, IMEX with diffusion implicit, advection
10 explicit.
12 See https://en.wikipedia.org/wiki/Burgers'_equation for the equation that is solved.
13 Discretization is done with a Chebychev method, which requires a first order derivative formulation.
14 Feel free to do a more efficient implementation using an ultraspherical method to avoid the first order business.
16 Parameters:
17 N (int): Spatial resolution
18 epsilon (float): viscosity
19 BCl (float): Value at left boundary
20 BCr (float): Value at right boundary
21 f (int): Frequency of the initial conditions
22 mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
23 """
25 dtype_u = mesh
26 dtype_f = imex_mesh
28 def __init__(self, N=64, epsilon=0.1, BCl=1, BCr=-1, f=0, mode='T2U', **kwargs):
29 """
30 Constructor. `kwargs` are forwarded to parent class constructor.
32 Args:
33 N (int): Spatial resolution
34 epsilon (float): viscosity
35 BCl (float): Value at left boundary
36 BCr (float): Value at right boundary
37 f (int): Frequency of the initial conditions
38 mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
39 """
40 self._makeAttributeAndRegister('N', 'epsilon', 'BCl', 'BCr', 'f', 'mode', localVars=locals(), readOnly=True)
42 bases = [{'base': 'cheby', 'N': N}]
43 components = ['u', 'ux']
45 super().__init__(bases=bases, components=components, spectral_space=False, **kwargs)
47 self.x = self.get_grid()[0]
49 # prepare matrices
50 Dx = self.get_differentiation_matrix(axes=(0,))
51 I = self.get_Id()
53 T2U = self.get_basis_change_matrix(conv=mode)
55 self.Dx = Dx
57 # construct linear operator
58 L_lhs = {'u': {'ux': -epsilon * (T2U @ Dx)}, 'ux': {'u': -T2U @ Dx, 'ux': T2U @ I}}
59 self.setup_L(L_lhs)
61 # construct mass matrix
62 M_lhs = {'u': {'u': T2U @ I}}
63 self.setup_M(M_lhs)
65 # boundary conditions
66 self.add_BC(component='u', equation='u', axis=0, x=1, v=BCr, kind='Dirichlet')
67 self.add_BC(component='u', equation='ux', axis=0, x=-1, v=BCl, kind='Dirichlet')
68 self.setup_BCs()
70 def u_exact(self, t=0, *args, **kwargs):
71 r"""
72 Routine to return the initial conditions at :math:`t=0` or the stationary solution for :math:`t=\infty`.
74 At :math:`t=0`, ``u`` is the linear function through the boundary values multiplied by :math:`\cos(\pi f x)`,
75 and ``ux`` is its derivative. For ``t=np.inf``, which requires ``f=0`` and ``BCl=-BCr``, ``u`` is the stationary
76 viscous shock connecting the boundary values and ``ux`` is left zero. Other times are not implemented.
78 Parameters
79 ----------
80 t : float, optional
81 Time, either ``0`` or ``np.inf``.
82 *args, **kwargs
83 Not used.
85 Returns
86 -------
87 me : dtype_u
88 The solution in physical space.
89 """
90 me = self.u_init
92 # x = (self.x + 1) / 2
93 # g = 4 * (1 + np.exp(-(4 * x + t)/self.epsilon/32))
94 # g_x = 4 * np.exp(-(4 * x + t)/self.epsilon/32) * (-4/self.epsilon/32)
96 # me[0] = 3./4. - 1./g
97 # me[1] = 1/g**2 * g_x
99 # return me
101 if t == 0:
102 me[self.index('u')][:] = ((self.BCr + self.BCl) / 2 + (self.BCr - self.BCl) / 2 * self.x) * np.cos(
103 self.x * np.pi * self.f
104 )
105 me[self.index('ux')][:] = (self.BCr - self.BCl) / 2 * np.cos(self.x * np.pi * self.f) + (
106 (self.BCr + self.BCl) / 2 + (self.BCr - self.BCl) / 2 * self.x
107 ) * self.f * np.pi * -np.sin(self.x * np.pi * self.f)
108 elif t == np.inf and self.f == 0 and self.BCl == -self.BCr:
109 me[0] = (self.BCl * np.exp((self.BCr - self.BCl) / (2 * self.epsilon) * self.x) + self.BCr) / (
110 np.exp((self.BCr - self.BCl) / (2 * self.epsilon) * self.x) + 1
111 )
112 else:
113 raise NotImplementedError
115 return me
117 def eval_f(self, u, *args, **kwargs):
118 r"""
119 Routine to evaluate the right-hand side, split into diffusion and advection.
121 The implicit part of the ``u`` equation is the diffusion :math:`\varepsilon \partial_x u_x`, computed by
122 applying the Chebychev differentiation matrix to the auxiliary component ``ux``. The explicit part is the
123 advection :math:`-u u_x`, computed pointwise. The algebraic ``ux`` equation gets zero in both parts.
125 Parameters
126 ----------
127 u : dtype_u
128 Current values of the numerical solution, in physical space.
129 *args, **kwargs
130 Not used, the right-hand side does not depend on time.
132 Returns
133 -------
134 f : dtype_f
135 The right-hand side, with ``impl`` and ``expl`` parts.
136 """
137 f = self.f_init
138 iu, iux = self.index('u'), self.index('ux')
140 u_hat = self.transform(u)
142 Dx_u_hat = self.u_init_forward
143 Dx_u_hat[iux] = (self.Dx @ u_hat[iux].flatten()).reshape(u_hat[iu].shape)
145 f.impl[iu] = self.epsilon * self.itransform(Dx_u_hat)[iux].real
146 f.expl[iu] = -u[iu] * u[iux]
147 return f
149 def get_fig(self): # pragma: no cover
150 """
151 Get a figure suitable to plot the solution of this problem.
153 Returns
154 -------
155 self.fig : matplotlib.pyplot.figure.Figure
156 """
157 import matplotlib.pyplot as plt
159 self.fig, axs = plt.subplots(constrained_layout=True)
160 return self.fig
162 def plot(self, u, t=None, fig=None, comp='u'): # pragma: no cover
163 r"""
164 Plot the solution.
166 Parameters
167 ----------
168 u : dtype_u
169 Solution to be plotted
170 t : float
171 Time to display at the top of the figure
172 fig : matplotlib.pyplot.figure.Figure, optional
173 Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.
174 comp : str, optional
175 Name of the component to plot, ``'u'`` or ``'ux'``.
177 Returns
178 -------
179 None
180 """
181 fig = self.get_fig() if fig is None else fig
182 ax = fig.axes[0]
184 ax.plot(self.x, u[self.index(comp)])
186 if t is not None:
187 fig.suptitle(f't = {t:.2e}')
189 ax.set_xlabel(r'$x$')
190 ax.set_ylabel(r'$u$')
193class Burgers2D(GenericSpectralLinear):
194 """
195 2D viscous Burgers equation, periodic in x (FFT) and Dirichlet in z (Chebychev), diffusion implicit, advection
196 explicit.
198 See https://en.wikipedia.org/wiki/Burgers'_equation for the equation that is solved.
199 This implementation is discretized with FFTs in x and Chebychev in z.
201 Parameters:
202 nx (int): Spatial resolution in x direction
203 nz (int): Spatial resolution in z direction
204 epsilon (float): viscosity
205 fux (int): Frequency of the initial conditions in x-direction
206 fuz (int): Frequency of the initial conditions in z-direction
207 mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
209 The boundary conditions are fixed: u = 0 at z = -1 and z = 1, v = -1 at z = -1 and v = 1 at z = 1.
210 """
212 dtype_u = mesh
213 dtype_f = imex_mesh
215 def __init__(self, nx=64, nz=64, epsilon=0.1, fux=2, fuz=1, mode='T2U', **kwargs):
216 """
217 Constructor. `kwargs` are forwarded to parent class constructor.
219 Args:
220 nx (int): Spatial resolution in x direction
221 nz (int): Spatial resolution in z direction
222 epsilon (float): viscosity
223 fux (int): Frequency of the initial conditions in x-direction
224 fuz (int): Frequency of the initial conditions in z-direction
225 mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
226 """
227 self._makeAttributeAndRegister('nx', 'nz', 'epsilon', 'fux', 'fuz', 'mode', localVars=locals(), readOnly=True)
229 bases = [
230 {'base': 'fft', 'N': nx},
231 {'base': 'cheby', 'N': nz},
232 ]
233 components = ['u', 'v', 'ux', 'uz', 'vx', 'vz']
234 super().__init__(bases=bases, components=components, spectral_space=False, **kwargs)
236 self.X, self.Z = self.get_grid()
238 # prepare matrices
239 Dx = self.get_differentiation_matrix(axes=(0,))
240 Dz = self.get_differentiation_matrix(axes=(1,))
241 I = self.get_Id()
243 T2U = self.get_basis_change_matrix(axes=(1,), conv=mode)
245 self.Dx = Dx
246 self.Dz = Dz
248 # construct linear operator
249 L_lhs = {
250 'u': {'ux': -epsilon * T2U @ Dx, 'uz': -epsilon * T2U @ Dz},
251 'v': {'vx': -epsilon * T2U @ Dx, 'vz': -epsilon * T2U @ Dz},
252 'ux': {'u': -T2U @ Dx, 'ux': T2U @ I},
253 'uz': {'u': -T2U @ Dz, 'uz': T2U @ I},
254 'vx': {'v': -T2U @ Dx, 'vx': T2U @ I},
255 'vz': {'v': -T2U @ Dz, 'vz': T2U @ I},
256 }
257 self.setup_L(L_lhs)
259 # construct mass matrix
260 M_lhs = {
261 'u': {'u': T2U @ I},
262 'v': {'v': T2U @ I},
263 }
264 self.setup_M(M_lhs)
266 # boundary conditions
267 self.BCtop = 1
268 self.BCbottom = -self.BCtop
269 self.BCtopu = 0
270 self.add_BC(component='v', equation='v', axis=1, v=self.BCtop, x=1, kind='Dirichlet')
271 self.add_BC(component='v', equation='vz', axis=1, v=self.BCbottom, x=-1, kind='Dirichlet')
272 self.add_BC(component='u', equation='uz', axis=1, v=self.BCtopu, x=1, kind='Dirichlet')
273 self.add_BC(component='u', equation='u', axis=1, v=self.BCtopu, x=-1, kind='Dirichlet')
274 self.setup_BCs()
276 def u_exact(self, t=0, *args, noise_level=0, **kwargs):
277 r"""
278 Routine to return the initial conditions, which are only implemented for :math:`t=0`.
280 ``u`` is :math:`\cos(f_{ux} x) \sin(\pi f_{uz} z)` plus ``BCtopu``, ``v`` is linear in :math:`z` between
281 ``BCbottom`` and ``BCtop``, and the derivative components are set accordingly. Seeded normally distributed
282 noise, multiplied by :math:`(z-1)(z+1)`, can be added to ``v``; its derivatives are not updated for the noise.
284 Parameters
285 ----------
286 t : float, optional
287 Time, has to be ``0``.
288 *args, **kwargs
289 Not used.
290 noise_level : float, optional
291 Amplitude of the noise added to ``v``.
293 Returns
294 -------
295 me : dtype_u
296 The initial conditions in physical space.
297 """
298 me = self.u_init
300 iu, iv, iux, iuz, ivx, ivz = self.index(self.components)
301 if t == 0:
302 me[iu] = self.xp.cos(self.X * self.fux) * self.xp.sin(self.Z * np.pi * self.fuz) + self.BCtopu
303 me[iux] = -self.xp.sin(self.X * self.fux) * self.fux * self.xp.sin(self.Z * np.pi * self.fuz)
304 me[iuz] = self.xp.cos(self.X * self.fux) * self.xp.cos(self.Z * np.pi * self.fuz) * np.pi * self.fuz
306 me[iv] = (self.BCtop + self.BCbottom) / 2 + (self.BCtop - self.BCbottom) / 2 * self.Z
307 me[ivz][:] = (self.BCtop - self.BCbottom) / 2
309 # add noise
310 rng = self.xp.random.default_rng(seed=99)
311 me[iv].real += rng.normal(size=me[iv].shape) * (self.Z - 1) * (self.Z + 1) * noise_level
313 else:
314 raise NotImplementedError
316 return me
318 def eval_f(self, u, *args, **kwargs):
319 r"""
320 Routine to evaluate the right-hand side, split into diffusion and advection.
322 The implicit part is the diffusion :math:`\varepsilon (\partial_x w_x + \partial_z w_z)` of both velocity
323 components :math:`w=u,v`, computed by spectral differentiation of the auxiliary derivative components. The
324 explicit part is the advection :math:`-(u \partial_x w + v \partial_z w)`, computed pointwise. The algebraic
325 derivative equations get zero in both parts.
327 Parameters
328 ----------
329 u : dtype_u
330 Current values of the numerical solution, in physical space.
331 *args, **kwargs
332 Not used, the right-hand side does not depend on time.
334 Returns
335 -------
336 f : dtype_f
337 The right-hand side, with ``impl`` and ``expl`` parts.
338 """
339 f = self.f_init
340 iu, iv, iux, iuz, ivx, ivz = self.index(self.components)
342 u_hat = self.transform(u)
343 f_hat = self.u_init_forward
344 f_hat[iu] = self.epsilon * ((self.Dx @ u_hat[iux].flatten() + self.Dz @ u_hat[iuz].flatten())).reshape(
345 u_hat[iux].shape
346 )
347 f_hat[iv] = self.epsilon * ((self.Dx @ u_hat[ivx].flatten() + self.Dz @ u_hat[ivz].flatten())).reshape(
348 u_hat[iux].shape
349 )
350 f.impl[...] = self.itransform(f_hat).real
352 f.expl[iu] = -(u[iu] * u[iux] + u[iv] * u[iuz])
353 f.expl[iv] = -(u[iu] * u[ivx] + u[iv] * u[ivz])
354 return f
356 def compute_vorticity(self, u):
357 r"""
358 Compute the vorticity from the velocity components by spectral differentiation, as :math:`\partial_x v +
359 \partial_z u`.
361 Parameters
362 ----------
363 u : dtype_u
364 Solution in physical space.
366 Returns
367 -------
368 vorticity : np.2darray
369 The vorticity in physical space.
370 """
371 me = self.u_init_forward
373 u_hat = self.transform(u)
374 iu, iv = self.index(['u', 'v'])
376 me[iu] = (self.Dx @ u_hat[iv].flatten() - self.Dz @ u_hat[iu].flatten()).reshape(u[iu].shape)
377 return self.itransform(me)[iu].real
379 def get_fig(self): # pragma: no cover
380 """
381 Get a figure suitable to plot the solution of this problem
383 Returns
384 -------
385 self.fig : matplotlib.pyplot.figure.Figure
386 """
387 import matplotlib.pyplot as plt
388 from mpl_toolkits.axes_grid1 import make_axes_locatable
390 self.fig, axs = plt.subplots(3, 1, sharex=True, sharey=True, figsize=((8, 7)), constrained_layout=True)
391 self.cax = []
392 divider = make_axes_locatable(axs[0])
393 self.cax += [divider.append_axes('right', size='3%', pad=0.03)]
394 divider2 = make_axes_locatable(axs[1])
395 self.cax += [divider2.append_axes('right', size='3%', pad=0.03)]
396 divider3 = make_axes_locatable(axs[2])
397 self.cax += [divider3.append_axes('right', size='3%', pad=0.03)]
398 return self.fig
400 def plot(self, u, t=None, fig=None, vmin=None, vmax=None): # pragma: no cover
401 r"""
402 Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
404 Parameters
405 ----------
406 u : dtype_u
407 Solution to be plotted
408 t : float
409 Time to display at the top of the figure
410 fig : matplotlib.pyplot.figure.Figure
411 Figure with the correct structure
412 vmin : float, optional
413 Lower limit of the colour scale for the velocity components.
414 vmax : float, optional
415 Upper limit of the colour scale for the velocity components.
417 Returns
418 -------
419 None
420 """
421 fig = self.get_fig() if fig is None else fig
422 axs = fig.axes
424 iu, iv = self.index(['u', 'v'])
426 imU = axs[0].pcolormesh(self.X, self.Z, u[iu].real, vmin=vmin, vmax=vmax)
427 imV = axs[1].pcolormesh(self.X, self.Z, u[iv].real, vmin=vmin, vmax=vmax)
428 imVort = axs[2].pcolormesh(self.X, self.Z, self.compute_vorticity(u).real)
430 for i, label in zip([0, 1, 2], [r'$u$', '$v$', 'vorticity'], strict=True):
431 axs[i].set_aspect(1)
432 axs[i].set_title(label)
434 if t is not None:
435 fig.suptitle(f't = {t:.2e}')
436 axs[-1].set_xlabel(r'$x$')
437 axs[-1].set_ylabel(r'$z$')
438 fig.colorbar(imU, self.cax[0])
439 fig.colorbar(imV, self.cax[1])
440 fig.colorbar(imVort, self.cax[2])