In [1]:

```
from lbmpy.session import *
from pystencils.timeloop import TimeLoop
```

In this notebook we demonstrate how to run a thermal lattice Boltzmann simulation. We use a separate set of distribution functions to solve the advection-diffusion equation for the temperature. The zeroth moment of these additional pdfs corresponds to temperature.

The thermal LB step is coupled to the normal hydrodynamic LB scheme using the Boussinesq approximation. The force on the liquid is proportional to the relative temperature. The hydrodynamic LB method computes the fluid velocity which in turn enters the thermal scheme, completing the two-way coupling.

To set this up in *lbmpy* we create first a `data handling`

object and create an array to store the temperature.

In [2]:

```
domain_size = (100, 50)
gpu = False
dh = ps.create_data_handling(domain_size)
temperature_field = dh.add_array("T", gpu=gpu)
dh.fill('T', val=1.0)
```

Next, we define how to compute the local force from the temperature field:

In [3]:

```
gravity = sp.Matrix([0, -1e-2])
force = -gravity * (temperature_field(0) - 1.0)
```

Now, we can create both LB steps.

The coupling is created by passing the following parameters to the hydrodynamic step,

`compute_velocity_in_every_step`

: usually the velocity is not computed/stored in every time step, only for output and plotting reasons. In the coupled algorithm we have to make sure that after every time step the current velocity is stored in an array, since it enters the thermal LB scheme.`force`

: we can simply pass our sympy expression for the force here, as long as the`LatticeBoltzmannStep`

operates on a data handling that stores the fields that are referenced in the expression. This is why we have to create the data handling first and pass it to both Step objects

and to the thermal step

`compute_density_in_every_step`

: density corresponds to the temperature here, which we need to be computed in every time step, since it enters the force expression for the hydrodynamic scheme`equilibrium_order`

: for the thermal LB method a first order accurate equilibrium is sufficient. This is slightly faster to compute than the normal equilibrium of order 2`velocity_input_array_name`

: the velocity entering the thermal equilibrium equation is not computed as first moment of the thermal pdfs. Instead, the hydrodynamic velocity is used here.

In [4]:

```
optimization = {'target': 'cpu' if gpu else 'cpu', 'openmp': 4}
hydro_step = LatticeBoltzmannStep(data_handling=dh, name='hydro', optimization=optimization,
relaxation_rate=1.8,
compute_velocity_in_every_step=True,
force=force)
thermal_step = LatticeBoltzmannStep(data_handling=dh, name='thermal', optimization=optimization,
relaxation_rate=1.8, density_data_name="T",
compute_density_in_every_step=True,
equilibrium_order=1,
velocity_input_array_name=hydro_step.velocity_data_name)
```

We add `NoSlip`

boundary conditions on all four walls for both schemes with the exception of the left wall of the thermal scheme. Here we set a heated wall, where the temperature is fixed to `1.01`

. This kind of Dirichlet boundary condition can be set using a `FixedDensity`

LB boundary.

In [5]:

```
add_box_boundary(hydro_step.boundary_handling)
add_box_boundary(thermal_step.boundary_handling)
thermal_step.boundary_handling.set_boundary(FixedDensity(1.01), slice_from_direction('W', dh.dim))
plt.figure(figsize=(20, 5), dpi=200)
plt.subplot(1, 2, 1)
plt.title("Hydrodynamic")
plt.boundary_handling(hydro_step.boundary_handling)
plt.subplot(1, 2, 2)
plt.title("Thermal")
plt.boundary_handling(thermal_step.boundary_handling)
```

In [6]:

```
def run(time_steps):
hydro_step.pre_run()
thermal_step.pre_run()
for t in range(time_steps):
hydro_step.time_step()
thermal_step.time_step()
hydro_step.post_run()
thermal_step.post_run()
hydro_step.time_steps_run += time_steps
thermal_step.time_steps_run += time_steps
```

In [7]:

```
run(5000)
assert dh.max('T') < 2
plt.figure(figsize=(16, 5), dpi=200)
plt.subplot(1, 2, 1)
plt.scalar_field(thermal_step.density[:, :])
plt.subplot(1, 2, 2)
plt.vector_field(hydro_step.velocity[:, :]);
```

In [8]:

```
assert np.isfinite(hydro_step.velocity[:, :].max())
assert np.isfinite(thermal_step.density[:, :].max())
```