[MUSIC] Hi, and welcome again to the fifth week of our class, simulation and modelling of natural processes. In the previous module, we have seen a short introduction of the lattice Boltzmann method and we have seen how it relates to the lattice gas automator we have seen in the previous lesson. We will now go into the details of the lattice Boltzmann method and develop a corresponding Python code. Let's go back to the first step of a lattice Boltzmann duration, which is a collision step. The collision step takes pre-collision populations, which we post in green, to post-collision population which we colored in red. We call the populations, we label populations with a variable which we call f. The incoming populations are the populations in pre collision state, the green population, we call them f in. In post collision state we call them the outgoing population, and name them f out. Those, f in and f out, have an index i, which labels the nine variables which we have on our every lattice cell. So there are nine variables fi in and nine variables fi out. As a reminder, there are nine because you have eight neighbors in a focal and diagonal directions, plus one rest population, which does not move during the propagation step. On this scheme, we are showing you the name of every individual variable on the left in the pre-collision state, you see that f0 to f3 are on the left hand side, f0 to f2 on the left hand side, f3, f4 and f5 are in the middle, and f6, f7, and f8 in are on the right hand side. During collision, they interact and they move to the opposite location. So after collision, f0 through f2 on the right-hand side, and f7, f6 through f8 are on the left-hand side. When you write the program for this, or in when you write out Python code, we are today going to allocate two sets of variables. One set of variable for the incoming populations, and one set of variables for the outgoing ones. So the green populations will be in their own matrix and the red populations will be in their own matrix as well. This is not necessary but it makes the code a bit easier. So this is what we are going to do today for didactic reasons. In the Python code, which you have on the righthand side, we first define the size of the system as the variables nx and ny, which are the number of cells in every space direction. Then we allocate a matrix f in, which is just a nine by nx by ny matrix, and f out just the same. Both of them are allocated with the NumPy command zeroes which initializes them to zero value. Now, how do we calculate macroscopic variables in lattice Boltzmann? We calculate them exactly the same, of course, as in the lattice gas automaton. The density are simply the sum of all the populations on one lattice node. The density, like the populations, is space and time dependent. Every cell has it's own density. To calculate the density you go in the cell, you take the nine populations and compute the sum. This gives you the macroscopic particle density on this location in space. If you do this in the Python code, you will have to calculate the density Density on every point of the space matrix on which we do our simulation, which means we start by allocating a matrix for density, again with the zeros command, which is an nx by ny matrix. This time we lost the third dimension because there is only one density per cell and not nine variables, as we have for the populations. Then we do a loop over all space directions, a loop for the ix variable going over all values of nx, and the iy variable going over all values in the y direction. We initialize the density on every cell to be zero first, then loop over the nine populations which we have on this cell, and sum them up to store them again in the density. You remember, from our introductory class to Python and high performance computing, that there are two ways of doing things in Python. Either we do this in a natural Python notation with loops, as I did it here on the left-hand side of the screen, or we do it with the NumPy array-based syntax. You can skip the loops of space dimensions by using an array based syntax which implies these loops. By saying simply that the role, which you don't need to allocate in this case, is equal to the sum of the populations. The loop over the space directions, in x and y direction, is implied here. You just need to tell NumPy which of the dimension of the populations the sum needs to be carried out. It will carried out over the dimension which has nine elements. It's the first of the three dimensions, because the populations are at nine times nx times ny matrix. We say this by stating that the parameter axis to the command sum is equal to zero, zero for the first space direction. In most cases when you do fluid dynamics, we are interested in the pressure. Now remember that lattice gas automata, they really model a gas. The Lattice Boltzmann method is the same. It is derived from lattice gas automata. So, from first principle, it is a maze to remodel a gas not a liquid. Although, we also use it afterwards to model a liquid. You do remember that [INAUDIBLE] equations are the same, on the first attempt at least, are the same equations you would use to remodel a liquid or a gas. Now in a gas, in an ideal gas, pressure is proportional to density. If the gas is isothermal then this proportionality constant is equal to the speed of sound inside this gas squared. So pressure is equal to speed of sound squared and density. Or p equals cs square row. As it is noted in the first equation, what's the speed of sound in our system? There are some lettuce Boltzmann methods in which you can trigger and change the speed of sound as you wish. But in the D2Q9 model, the speed of sound is a constant which you cannot play with. The speed of sound squared is equal to one third times this discrete parameter delta x squared over delta t squared. This is the first of the set of equations which I will just write down without actually proving them, because the purpose of this lesson is to develop, practically speaking, together with you a first lattice Boltzmann code. We won't have the time to develop all the equations theoretically. We'll just motivate them. If you're interested in more information, you will find easily further literature on the internet. What are the units of pressure? We are used to NSI units to have pressure as Pascal. But if you take the pressure as we defined it there, you take the density and multiply it by delta x squared with delta t squared. If delta x is meters and delta t is seconds, this gives you the units of meter square per second square. In order to get Pascal, you need to multiply this by the density of the material you are working with. This is a constant material property. In a gas, this could be one kilogram per cubed meter. Or in a liquid it could be 1,000 kilograms per cube meter. Multiply your meter square or your second square by kilogram, by cube meter, and you get the units of Pascal, kilogram per meter second square. The velocity is just the same, calculated, as we calculated it, in lattice gas automaton. To calculate a velocity, we need to introduce a new set of constants, which are vectors. We call them the lattice velocities. We label them v0 through v8. The lattice vectors are vectors which connect one of the cells in our simulation to one of its eight neighbors or even nine neighbors, where the zero velocity corresponds to the rest population. For example, on the left-hand side of the screen you have the nine lattice velocities which we are going to use for the code we develop today. The first one, v0, points to the upper right of the screen. And it has the components (1,1). v1 has the components (1,0), v2 (1,-1), and so on. If you have a lattice cell, which is placed at position x, you can find all of its eight neighbors by calculating x + delta t times one of the lattice velocities. The macroscopic velocity, the velocity of the flow in one of the cells, is nothing else than the sum over all populations just like the cellular density. But this time the sum is weighted. It is weighted by the value of the lattice velocities. They take the sum over the nine populations, multiply them by the corresponding lattice velocity, and what you get is density times velocity. If you are interested in just velocity, you divide the result by the local density, calculate it as we did on the previous slide. Let's do this with the Python code. First, we find the discrete velocities. Discrete velocities are a two dimensional array because we have nine velocities each of which has two components. Then let's calculate the velocity first again without using the array based on python syntax just with traditional python code. So we are sure to exactly understand what's going on. First, we allocate memory for the velocities. It is a matrix, an x by a y with two components in each cell. So, a three-dimensional two by nx or ny matrix. We have the two [INAUDIBLE] base directions, ix and iy. Start by initializing the two components of the velocity in each cell to zero, then loop over to nine populations and add the value of the populations times the component of the corresponding lattice velocity to the macroscopic velocity. We end up at the end of this loop with the right value of the macroscopic velocity in each cell. With the array-based NumPy syntax, everything, again, is much easier to write. Not necessarily much easier to understand, but at least shorter to write. And as it was mentioned in the introduction to Python and high performance computing, it is much faster true that the right way for formulating a code with Python which will run fast. In this case, we still have the loop over the nine populations, but we skip the time loop. We loop over the nine populations and add implicitly over all cells of the matrix, the value of the populations times the corresponding discrete velocity, lattice velocity, and add it up to get the microscopic velocity. In the end, we have to divide by density, because there are some [INAUDIBLE] density times velocity are not the single velocity. At this point, I probably got you completely confused, because I used two kinds of symbols which looks very similar, but have completely different meaning. On the left hand side of the screen we have the lattice velocities, and on the right hand side of the screen we have the populations. Both of which are plotted through arrows. But they don't have the same meaning. The lattice velocities are really vectors, vectors in the mathematical sense of the term. They are constants and do not change in time. On the right hand side, the populations of scalar values. We just plot them as arrows to remember that they will go, during the streaming step, they will move in one direction to a neighboring node. But the arrow does not mean that these are vectors, these are scalars, and their value is indicated by the thickness of the arrow. Also, the populations are not constant. They are the variables we are solving for in this numerical scheme. And they change in space and time. This is the end of our module on macroscopic variables and the lattice Boltzmann method. Stay tuned for the next module. [MUSIC]
