Demo - Rayleigh Benard

Demo - Rayleigh Benard#

Mikael Mortensen (email: mikaem@math.uio.no), Department of Mathematics, University of Oslo.

Date: November 21, 2019

Summary. Rayleigh-Benard convection arise due to temperature gradients in a fluid. The governing equations are Navier-Stokes coupled (through buoyancy) with an additional temperature equation derived from the first law of thermodynamics, using a linear correlation between density and temperature.

This is a demonstration of how the Python module shenfun can be used to solve for these Rayleigh-Benard cells in a 2D channel with two walls of different temperature in one direction, and periodicity in the other direction. The solver described runs with MPI without any further considerations required from the user. Note that there is also a more physically realistic 3D solver. The solver described in this demo has been implemented in a class in the RayleighBenard2D.py module in the demo folder of shenfun. Below is an example solution, which has been run at a very high Rayleigh number (Ra).

Figure 1: Temperature fluctuations in the Rayleigh Benard flow. The top and bottom walls are kept at different temperatures and this sets up the Rayleigh-Benard convection. The simulation is run at Ra =1,000,000, Pr =0.7 with 256 and 512 quadrature points in x and y-directions, respectively.

The Rayleigh Bénard equations#

The governing equations solved in domain $Ω = (- 1, 1) \times [0, 2 π)$ are

\begin{matrix} (1) & \frac{\partial u}{\partial t} + (u \cdot \nabla) u = - \nabla p + \sqrt{\frac{P r}{R a}} \nabla^{2} u + T i, \end{matrix}

\begin{matrix} (2) & \frac{\partial T}{\partial t} + u \cdot \nabla T = \frac{1}{\sqrt{R a P r}} \nabla^{2} T, \end{matrix}

\begin{matrix} (3) & \nabla \cdot u = 0, \end{matrix}

where $u (x, y, t) (= u i + v j)$ is the velocity vector, $p (x, y, t)$ is pressure, $T (x, y, t)$ is the temperature, and $i$ and $j$ are the unity vectors for the $x$ and $y$ -directions, respectively.

The equations are complemented with boundary conditions $u (\pm 1, y, t) = (0, 0), u (x, 2 π, t) = u (x, 0, t), T (- 1, y, t) = 1, T (1, y, t) = 0, T (x, 2 π, t) = T (x, 0, t)$ . Note that these equations have been non-dimensionalized according to [pandey18], using dimensionless Rayleigh number $R a = g α Δ T h^{3} / (ν κ)$ and Prandtl number $P r = ν / κ$ . Here $g i$ is the vector accelleration of gravity, $Δ T$ is the temperature difference between the top and bottom walls, $h$ is the hight of the channel in $x$ -direction, $ν$ is the dynamic viscosity coefficient, $κ$ is the heat transfer coefficient and $α$ is the thermal expansion coefficient. Note that the governing equations have been non-dimensionalized using the free-fall velocityscale $U = \sqrt{g α Δ T h}$ . See [pandey18] for more details.

The governing equations contain a non-trivial coupling between velocity, pressure and temperature. This coupling can be simplified by eliminating the pressure from the equation for the wall-normal velocity component $u$ . We accomplish this by taking the Laplace of the momentum equation in wall normal direction, using the pressure from the divergence of the momentum equation $\nabla^{2} p = - \nabla \cdot H + \partial T / \partial x$ , where $H = (H_{x}, H_{y}) = (u \cdot \nabla) u$

\begin{matrix} (4) & \frac{\partial \nabla^{2} u}{\partial t} = \frac{\partial^{2} H_{y}}{\partial x \partial y} - \frac{\partial^{2} H_{x}}{\partial y \partial y} + \sqrt{\frac{P r}{R a}} \nabla^{4} u + \frac{\partial^{2} T}{\partial y^{2}} . \end{matrix}

This equation is solved with $u (\pm 1, y, t) = \partial u / \partial x (\pm 1, y, t) = 0$ , where the latter follows from the divergence constraint. In summary, we now seem to have the following equations to solve:

\begin{matrix} (5) & \frac{\partial \nabla^{2} u}{\partial t} = \frac{\partial^{2} H_{y}}{\partial x \partial y} - \frac{\partial^{2} H_{x}}{\partial y \partial y} + \sqrt{\frac{P r}{R a}} \nabla^{4} u + \frac{\partial^{2} T}{\partial y^{2}}, \end{matrix}

\begin{matrix} (6) & \frac{\partial v}{\partial t} + H_{y} = - \frac{\partial p}{\partial y} + \sqrt{\frac{P r}{R a}} \nabla^{2} v, \end{matrix}

\begin{matrix} (7) & \frac{\partial T}{\partial t} + u \cdot \nabla T = \frac{1}{\sqrt{R a P r}} \nabla^{2} T, \end{matrix}

\begin{matrix} (8) & \nabla \cdot u = 0 . \end{matrix}

However, we note that Eqs. (5) and (7) and (8) do not depend on pressure, and, apparently, on each time step we can solve (5) for $u$ , then (8) for $v$ and finally (7) for $T$ . So what do we need (6) for? It appears to have become redundant from the elimination of the pressure from Eq. (5). It turns out that this is actually almost completely true, but (5), (7) and (8) can only provide closure for all but one of the Fourier coefficients. To see this it is necessary to introduce some discretization and basis functions that will be used to solve the problem. To this end we use $P_{N}$ , which is the set of all real polynomials of degree less than or equal to N and introduce the following finite-dimensional approximation spaces

\begin{matrix} (9) & V_{N}^{B} (x) = {v \in P_{N} | v (\pm 1) = v ´ (\pm 1) = 0}, \end{matrix}

\begin{matrix} (10) & V_{N}^{D} (x) = {v \in P_{N} | v (\pm 1) = 0}, \end{matrix}

\begin{matrix} (11) & V_{N}^{T} (x) = {v \in P_{N} | v (- 1) = 0, v (1) = 1}, \end{matrix}

\begin{matrix} (12) & V_{N}^{W} (x) = {v \in P_{N}}, \end{matrix}

\begin{matrix} (13) & V_{M}^{F} (y) = {\exp (ı l y) | l \in [- M / 2, - M / 2 + 1, \dots M / 2 - 1]} . \end{matrix}

Here $dim (V_{N}^{B}) = N - 4, dim (V_{N}^{D}) = dim (V_{N}^{W}) = dim (V_{N}^{T}) = N - 2$ and $dim (V_{M}^{F}) = M$ . We note that $V_{N}^{B}, V_{N}^{D}, V_{N}^{W}$ and $V_{N}^{T}$ can be used to approximate $u, v, T$ and $p$ , respectively, in the $x$ -direction. Also note that for $V_{M}^{F}$ it is assumed that $M$ is an even number.

We can now choose basis functions for the spaces, using Shen’s composite bases for either Legendre or Chebyshev polynomials. For the Fourier space the basis functions are already given. We leave the actual choice of basis as an implementation option for later. For now we use $ϕ^{B} (x), ϕ^{D} (x), ϕ^{W}$ and $ϕ^{T} (x)$ as common notation for basis functions in spaces $V_{N}^{B}, V_{N}^{D}, V_{N}^{W}$ and $V_{N}^{T}$ , respectively.

To get the required approximation spaces for the entire domain we use tensor products of the one-dimensional spaces in (9)-(13)

\begin{matrix} (14) & W_{B F} = V_{N}^{B} \otimes V_{M}^{F}, \end{matrix}

\begin{matrix} (15) & W_{D F} = V_{N}^{D} \otimes V_{M}^{F}, \end{matrix}

\begin{matrix} (16) & W_{T F} = V_{N}^{T} \otimes V_{M}^{F}, \end{matrix}

\begin{matrix} (17) & W_{W F} = V_{N}^{W} \otimes V_{M}^{F} . \end{matrix}

Space $W_{B F}$ has 2D tensor product basis functions $ϕ_{k}^{B} (x) \exp (ı l y)$ and similar for the others. All in all we get the following approximations for the unknowns

\begin{matrix} (18) & u_{N} (x, y, t) = \sum_{k \in k_{B}} \sum_{l \in l} {\hat{u}}_{k l} (t) ϕ_{k}^{B} (x) \exp (ı l y), \end{matrix}

\begin{matrix} (19) & v_{N} (x, y, t) = \sum_{k \in k_{D}} \sum_{l \in l} {\hat{v}}_{k l} (t) ϕ_{k}^{D} (x) \exp (ı l y), \end{matrix}

\begin{matrix} (20) & p_{N} (x, y, t) = \sum_{k \in k_{W}} \sum_{l \in l} {\hat{p}}_{k l} (t) ϕ_{k}^{W} (x) \exp (ı l y), \end{matrix}

\begin{matrix} (21) & T_{N} (x, y, t) = \sum_{k \in k_{T}} \sum_{l \in l} {\hat{T}}_{k l} (t) ϕ_{k}^{T} (x) \exp (ı l y), \end{matrix}

where $k_{x} = {0, 1, \dots dim (V_{N}^{x}) - 1}, for x \in (B, D, W, T)$ and $l = {- M / 2, - M / 2 + 1, \dots, M / 2 - 1}$ . Note that since the problem is defined in real space we will have Hermitian symmetry. This means that ${\hat{u}}_{k, l} = {\overset{―}{\hat{u}}}_{k, - l}$ , with an overbar being a complex conjugate, and similar for ${\hat{v}}_{k l}, {\hat{p}}_{k l}$ and ${\hat{T}}_{k l}$ . For this reason we can get away with solving for only the positive $l$ ’s, as long as we remember that the sum in the end goes over both positive and negative $l^{'} s$ . This is actually automatically taken care of by the FFT provider and is not much of an additional complexity in the implementation. So from now on $l = {0, 1, \dots, M / 2}$ .

We can now take a look at why Eq. (6) is needed. If we first solve (5) for ${\hat{u}}_{k l} (t), (k, l) \in k_{B} \times l$ , then we can use (8) to solve for ${\hat{v}}_{k l} (t)$ . But here there is a problem. We can see this by creating the variational form required to solve (8) by the spectral Galerkin method. To this end make $v = v_{N}$ in (8) a trial function, use $u = u_{N}$ a known function and take the weighted inner product over the domain using test function $q \in W_{D F}$

\begin{matrix} (22) & {⟨ \frac{\partial u_{N}}{\partial x} + \frac{\partial v_{N}}{\partial y}, q ⟩}_{w} = 0. \end{matrix}

Here we are using the inner product notation

\begin{matrix} (23) & {⟨ a, b ⟩}_{w} = \int_{- 1}^{1} \int_{0}^{2 π} a \overset{―}{b} d x_{w} d y_{w} (\approx \sum_{i} \sum_{j} a (x_{i}, y_{j}) \overset{―}{b} (x_{i}, y_{j}) w_{i} w_{j}), \end{matrix}

where the exact form of the weighted scalar product depends on the chosen basis; Legendre has $d x_{w} = d x$ , Chebyshev $d x_{w} = d x / \sqrt{1 - x^{2}}$ and Fourier $d y_{w} = d y / 2 / π$ . The bases also have associated quadrature weights ${w_{i}}_{i = 0}^{N - 1}$ and ${w_{j}}_{j = 0}^{M - 1}$ that are used to approximate the integrals.

Inserting now for the known $u_{N}$ , the unknown $v_{N}$ , and $q = ϕ_{m}^{D} (x) \exp (ı n y)$ the continuity equation becomes

\begin{array}{r} \int_{- 1}^{1} \int_{0}^{2 π} \frac{\partial}{\partial x} (\sum_{k \in k_{B}} \sum_{l \in l} {\hat{u}}_{k l} (t) ϕ_{k}^{B} (x) \exp (ı l y)) ϕ_{m}^{D} (x) \exp (- ı n y) d x_{w} d y_{w} + \\ \int_{- 1}^{1} \int_{0}^{2 π} \frac{\partial}{\partial y} (\sum_{k \in k_{D}} \sum_{l \in l} {\hat{v}}_{k l} (t) ϕ_{k}^{D} (x) \exp (ı l y)) ϕ_{m}^{D} (x) \exp (- ı n y) d x_{w} d y_{w} = 0. \end{array}

The $x$ and $y$ domains are separable, so we can rewrite as

\begin{array}{r} \sum_{k \in k_{B}} \sum_{l \in l} \int_{- 1}^{1} \frac{\partial ϕ_{k}^{B} (x)}{\partial x} ϕ_{m}^{D} (x) d x_{w} \int_{0}^{2 π} \exp (ı l y) \exp (- ı n y) d y_{w} {\hat{u}}_{k l} + \\ \sum_{k \in k_{D}} \sum_{l \in l} \int_{- 1}^{1} ϕ_{k}^{D} (x) ϕ_{m}^{D} (x) d x_{w} \int_{0}^{2 π} \frac{\partial \exp (ı l y)}{\partial y} \exp (- ı n y) d y_{w} {\hat{v}}_{k l} = 0. \end{array}

Now perform some exact manipulations in the Fourier direction and introduce the 1D inner product notation for the $x$ -direction

\begin{matrix} (24) & {(a, b)}_{w} = \int_{- 1}^{1} a (x) b (x) d x_{w} (\approx \sum_{j = 0}^{N - 1} a (x_{j}) b (x_{j}) w_{j}) . \end{matrix}

By also simplifying the notation using summation of repeated indices, we get the following equation

\begin{matrix} (25) & δ_{l n} {(\frac{\partial ϕ_{k}^{B}}{\partial x}, ϕ_{m}^{D})}_{w} {\hat{u}}_{k l} + ı l δ_{l n} {(ϕ_{k}^{D}, ϕ_{m}^{D})}_{w} {\hat{v}}_{k l} = 0. \end{matrix}

Now $l$ must equal $n$ and we can simplify some more

\begin{matrix} (26) & {(\frac{\partial ϕ_{k}^{B}}{\partial x}, ϕ_{m}^{D})}_{w} {\hat{u}}_{k l} + ı l {(ϕ_{k}^{D}, ϕ_{m}^{D})}_{w} {\hat{v}}_{k l} = 0. \end{matrix}

We see that this equation can be solved for ${\hat{v}}_{k l} for (k, l) \in k_{D} \times [1, 2, \dots, M / 2]$ , but try with $l = 0$ and you hit division by zero, which obviously is not allowed. And this is the reason why Eq. (6) is still needed, to solve for ${\hat{v}}_{k, 0}$ ! Fortunately, since $\exp (ı 0 y) = 1$ , the pressure derivative $\frac{\partial p}{\partial y} = 0$ , and as such the pressure is still not required. When used only for Fourier coefficient 0, Eq. (6) becomes

\begin{matrix} (27) & \frac{\partial v}{\partial t} + N_{y} = \sqrt{\frac{P r}{R a}} \frac{\partial^{2} v}{\partial x^{2}} . \end{matrix}

There is still one more revelation to be made from Eq. (26). When $l = 0$ we get

\begin{matrix} (28) & {(\frac{\partial ϕ_{k}^{B}}{\partial x}, ϕ_{m}^{D})}_{w} {\hat{u}}_{k, 0} = 0, \end{matrix}

which is trivially satisfied if ${\hat{u}}_{k, 0} = 0$ for $k \in k_{B}$ . Bottom line is that we only need to solve Eq. (5) for $l \in l / {0}$ , whereas we can use directly ${\hat{u}}_{k, 0} = 0 for k \in k_{B}$ .

To sum up, with the solution known at $t = t - Δ t$ , we solve

Equation	For unknown	With indices
(5)	${\hat{u}}_{k l} (t)$	$(k, l) \in k_{B} \times l / {0}$
(8)	${\hat{v}}_{k l} (t)$	$(k, l) \in k_{D} \times l / {0}$
(27)	${\hat{v}}_{k l} (t)$	$(k, l) \in k_{D} \times {0}$
(7)	${\hat{T}}_{k l} (t)$	$(k, l) \in k_{T} \times l$

Temporal discretization#

The governing equations are integrated in time using any one of the time steppers available in shenfun. There are several possible IMEX Runge Kutta methods, see integrators.py. The time steppers are used for any generic equation

\begin{matrix} (29) & \frac{\partial ψ}{\partial t} = N + L ψ, \end{matrix}

where $N$ and $L$ represents the nonlinear and linear contributions, respectively. The timesteppers are provided with $ψ, L$ and $N$ , and possibly some tailored linear algebra solvers, and solvers are then further assembled under the hood.

All the timesteppers split one time step into one or several stages. The classes IMEXRK222, IMEXRK3 and IMEXRK443 have 2, 3 and 4 steps, respectively, and the cost is proportional.

Implementation#

To get started we need instances of the approximation spaces discussed in Eqs. (9) - (17). When the spaces are created we also need to specify the family and the dimension of each space. Here we simply choose Chebyshev and Fourier with 100 and 256 quadrature points in $x$ and $y$ -directions, respectively. We could replace ‘Chebyshev’ by ‘Legendre’, but the former is known to be faster due to the existence of fast transforms.

from shenfun import *

N, M = 64, 128
family = 'Chebyshev'
VB = FunctionSpace(N, family, bc=(0, 0, 0, 0))
VD = FunctionSpace(N, family, bc=(0, 0))
VW = FunctionSpace(N, family)
VT = FunctionSpace(N, family, bc=(1, 0))
VF = FunctionSpace(M, 'F', dtype='d')

And then we create tensor product spaces by combining these bases (like in Eqs. (14)-(17)).

W_BF = TensorProductSpace(comm, (VB, VF))    # Wall-normal velocity
W_DF = TensorProductSpace(comm, (VD, VF))    # Streamwise velocity
W_WF = TensorProductSpace(comm, (VW, VF))    # No bc
W_TF = TensorProductSpace(comm, (VT, VF))    # Temperature
BD = VectorSpace([W_BF, W_DF])   # Velocity vector
DD = VectorSpace([W_DF, W_DF])   # Convection vector
W_DFp = W_DF.get_dealiased(padding_factor=1.5)
BDp = BD.get_dealiased(padding_factor=1.5)

Here the VectorSpae create mixed tensor product spaces by the Cartesian products BD = W_BF $\times$ W_DF and DD = W_DF $\times$ W_DF. These mixed space will be used to hold the velocity and convection vectors, but we will not solve the equations in a coupled manner and continue in the segregated approach outlined above.

We also need containers for the computed solutions. These are created as

u_  = Function(BD)     # Velocity vector, two components
T_  = Function(W_TF)   # Temperature
H_  = Function(DD)     # Convection vector
uT_ = Function(BD)     # u times T

Wall-normal velocity equation#

We implement Eq. (5) using a generic time stepper. To this end we first need to declare some test- and trial functions, as well as some model constants and the length of the time step. We store the PDE time stepper in a dictionary called pdes just for convenience:

# Specify viscosity and time step size using dimensionless Ra and Pr
Ra = 1000000
Pr = 0.7
nu = np.sqrt(Pr/Ra)
kappa = 1./np.sqrt(Pr*Ra)
dt = 0.025

# Choose timestepper and create instance of class

PDE = IMEXRK3 # IMEX222, IMEXRK443

v = TestFunction(W_BF) # The space we're solving for u in

pdes = {
    'u': PDE(v,                              # test function
             div(grad(u_[0])),               # u
             lambda f: nu*div(grad(f)),      # linear operator on u
             [Dx(Dx(H_[1], 0, 1), 1, 1)-Dx(H_[0], 1, 2), Dx(T_, 1, 2)],
             dt=dt,
             solver=chebyshev.la.Biharmonic if family == 'Chebyshev' else la.SolverGeneric1ND,
             latex=r"\frac{\partial \nabla^2 u}{\partial t} = \nu \nabla^4 u + \frac{\partial^2 N_y}{\partial x \partial y} - \frac{\partial^2 N_x}{\partial y^2}")
}
pdes['u'].assemble()

Notice the one-to-one resemblance with (5).

The right hand side depends on the convection vector $H$ , which can be computed in many different ways. Here we will make use of the identity $(u \cdot \nabla) u = - u \times (\nabla \times u) + 0.5 \nabla u \cdot u$ , where $0.5 \nabla u \cdot u$ can be added to the eliminated pressure and as such be neglected. Compute $H = - u \times (\nabla \times u)$ by first evaluating the velocity and the curl in real space. The curl is obtained by projection of $\nabla \times u$ to the no-boundary-condition space W_TF, followed by a backward transform to real space. The velocity is simply transformed backwards with padding.

# Get a mask for setting Nyquist frequency to zero
mask = W_DF.get_mask_nyquist()
Curl = Project(curl(u_), W_WF) # Instance used to compute curl

def compute_convection(u, H):
    up = u.backward(padding_factor=1.5).v
    curl = Curl().backward(padding_factor=1.5)
    H[0] = W_DFp.forward(-curl*up[1])
    H[1] = W_DFp.forward(curl*up[0])
    H.mask_nyquist(mask)
    return H

Note that the convection has a homogeneous Dirichlet boundary condition in the non-periodic direction.

Streamwise velocity#

The streamwise velocity is computed using Eq. (26) and (27). The first part is done fastest by projecting $f = \frac{\partial u}{\partial x}$ to the same Dirichlet space W_DF used by $v$ . This is most efficiently done by creating a class to do it

f = dudx = Project(Dx(u_[0], 0, 1), W_DF)

Since $f$ now is in the same space as the streamwise velocity, Eq. (26) simplifies to

\begin{matrix} (30) & ı l {(ϕ_{k}^{D}, ϕ_{m}^{D})}_{w} {\hat{v}}_{k l} = - {(ϕ_{k}^{D}, ϕ_{m}^{D})}_{w} {\hat{f}}_{k l}, \end{matrix}

\begin{matrix} (31) & {\hat{v}}_{k l} = \frac{ı {\hat{f}}_{k l}}{l} . \end{matrix}

We thus compute ${\hat{v}}_{k l}$ for all $k$ and $l > 0$ as

K = W_BF.local_wavenumbers(scaled=True)
K[1][0, 0] = 1 # to avoid division by zero. This component is computed later anyway.
u_[1] = 1j*dudx()/K[1]

which leaves only ${\hat{v}}_{k 0}$ . For this we use (27) and get

v00 = Function(VD)
v0 = TestFunction(VD)
h1 = Function(VD) # convection equal to H_[1, :, 0]
pdes1d = {
    'v0': PDE(v0,
              v00,
              lambda f: nu*div(grad(f)),
              -Expr(h1),
              dt=dt,
              solver=chebyshev.la.Helmholtz if family == 'Chebyshev' else la.Solver,
              latex=r"\frac{\partial v}{\partial t} = \nu \frac{\partial^2 v}{\partial x^2} - N_y "),
}
pdes1d['v0'].assemble()

A function that computes v for an integration stage rk is

def compute_v(rk):
    v00[:] = u_[1, :, 0].real
    h1[:] = H_[1, :, 0].real
    u_[1] = 1j*dudx()/K[1]
    pdes1d['v0'].compute_rhs(rk)
    u_[1, :, 0] = pdes1d['v0'].solve_step(rk)

Temperature#

The temperature equation (2) is implemented using a Helmholtz solver. The main difficulty with the temperature is the non-homogeneous boundary condition that requires special attention. A non-zero Dirichlet boundary condition is implemented by adding two basis functions to the basis of the function space

\begin{matrix} (32) & ϕ_{N - 2}^{D} = 0.5 (1 + x), \end{matrix}

\begin{matrix} (33) & ϕ_{N - 1}^{D} = 0.5 (1 - x), \end{matrix}

with the approximation now becoming

\begin{matrix} (34) & T_{N} (x, y, t) = \sum_{k = 0}^{N - 1} \sum_{l \in l} {\hat{T}}_{k l} ϕ_{k}^{D} (x) \exp (ı l y), \end{matrix}

\begin{matrix} (35) & = \sum_{k = 0}^{N - 3} \sum_{l \in l} {\hat{T}}_{k l} ϕ_{k}^{D} (x) \exp (ı l y) + \sum_{k = N - 2}^{N - 1} \sum_{l \in l} {\hat{T}}_{k l} ϕ_{k}^{D} (x) \exp (ı l y) . \end{matrix}

The boundary condition requires

\begin{matrix} (36) & T_{N} (1, y, t) = \sum_{k = N - 2}^{N - 1} \sum_{l \in l} {\hat{T}}_{k l} ϕ_{k}^{D} (1) \exp (ı l y), \end{matrix}

\begin{matrix} (37) & = \sum_{l \in l} {\hat{T}}_{N - 2, l} \exp (ı l y), \end{matrix}

and

\begin{matrix} (38) & T_{N} (- 1, y, t) = \sum_{k = N - 2}^{N - 1} \sum_{l \in l} {\hat{T}}_{k l} ϕ_{k}^{D} (- 1) \exp (ı l y), \end{matrix}

\begin{matrix} (39) & = \sum_{l \in l} {\hat{T}}_{N - 1, l} \exp (ı l y) . \end{matrix}

We find ${\hat{T}}_{N - 2, l}$ and ${\hat{T}}_{N - 1, l}$ using orthogonality. Multiply (37) and (39) by $\exp (- ı m y)$ and integrate over the domain $[0, 2 π]$ . We get

\begin{matrix} (40) & {\hat{T}}_{N - 2, l} = \int_{0}^{2 π} T_{N} (1, y, t) \exp (- ı l y) d y, \end{matrix}

\begin{matrix} (41) & {\hat{T}}_{N - 1, l} = \int_{0}^{2 π} T_{N} (- 1, y, t) \exp (- ı l y) d y . \end{matrix}

Using this approach it is easy to see that any inhomogeneous function $T_{N} (\pm 1, y, t)$ of $y$ and $t$ can be used for the boundary condition, and not just a constant. However, we will not get into this here. And luckily for us all this complexity with boundary conditions will be taken care of under the hood by shenfun.

A time stepper for the temperature equation is implemented as

uT_ = Function(BD)
q = TestFunction(W_TF)
pdes['T'] = PDE(q,
                T_,
                lambda f: kappa*div(grad(f)),
                -div(uT_),
                dt=dt,
                solver=chebyshev.la.Helmholtz if family == 'Chebyshev' else la.SolverGeneric1ND,
                latex=r"\frac{\partial T}{\partial t} = \kappa \nabla^2 T - \nabla \cdot \vec{u}T")
pdes['T'].assemble()

The uT_ term is computed with dealiasing as

def compute_uT(u_, T_, uT_):
    up = u_.backward(padding_factor=1.5)
    Tp = T_.backward(padding_factor=1.5)
    uT_ = BDp.forward(up*Tp, uT_)
    return uT_

Finally all that is left is to initialize the solution and integrate it forward in time.

%matplotlib inline

import matplotlib.pyplot as plt
# initialization
T_b = Array(W_TF)
X = W_TF.local_mesh(True)
#T_b[:] = 0.5*(1-X[0]) + 0.001*np.random.randn(*T_b.shape)*(1-X[0])*(1+X[0])
T_b[:] = 0.5*(1-X[0]+0.25*np.sin(np.pi*X[0]))+0.001*np.random.randn(*T_b.shape)*(1-X[0])*(1+X[0])
T_ = T_b.forward(T_)
T_.mask_nyquist(mask)

t = 0
tstep = 0
end_time = 15
while t < end_time-1e-8:
    for rk in range(PDE.steps()):
        compute_convection(u_, H_)
        compute_uT(u_, T_, uT_)
        pdes['u'].compute_rhs(rk)
        pdes['T'].compute_rhs(rk)
        pdes['u'].solve_step(rk)
        compute_v(rk)
        pdes['T'].solve_step(rk)
    t += dt
    tstep += 1
plt.contourf(X[1], X[0], T_.backward(), 100)
plt.show()

A complete solver implemented in a solver class can be found in RayleighBenard2D.py. Note that in the final solver it is also possible to use a $(y, t)$ -dependent boundary condition for the hot wall. And the solver can also be configured to store intermediate results to an HDF5 format that later can be visualized in, e.g., Paraview. The movie in the beginning of this demo has been created in Paraview.

A. Pandey, J. D. Scheel and J. Schumacher. Turbulent Superstructures in Rayleigh-B’enard Convection, Nature Communications, 9(1), pp. 2118, doi: 10.1038/s41467-018-04478-0, 2018.