Answer (start)

Fluid Dynamics#

  • Author:

  • Date:

  • Time spent on this assignment:

import numpy as np
import matplotlib.pyplot as plt
import math
from matplotlib.animation import FuncAnimation
import itertools
from IPython.display import HTML
import pickle
def resetMe(keepList=[]):
    ll=%who_ls
    keepList=keepList+['resetMe','np','plt','math','FuncAnimation',
                       'HTML','itertools','pickle','testFunc']
    for iiii in keepList:
        if iiii in ll:
            ll.remove(iiii)
    for iiii in ll:
        jjjj="^"+iiii+"$"
        %reset_selective -f {jjjj}
    ll=%who_ls
    return

def testFunc(func,inFiles,outFiles):
    inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
    outputs = [pickle.load(open(f,"rb")) for f in outFiles]
    result  = func(*inputs)
    allGood = True
    if not isinstance(result, tuple): result = (result,)
    for i in range(len(outputs)):
        if np.max(np.abs(result[i]-outputs[i])) > 1e-14:
            print("Failed test for",outFiles[i],i,np.max(np.abs(result[i]-outputs[i])))
            allGood = False
    if allGood: print("Test Passed!")
    else:       print("Test Failed :(")

Download the test files here, or download from the website and upload & unzip them to wherever this notebook file is

!wget https://courses.physics.illinois.edu/phys246/fa2020/code/TestFiles.zip && unzip TestFiles.zip
--2024-11-07 16:40:55--  https://courses.physics.illinois.edu/phys246/fa2020/code/TestFiles.zip
Resolving courses.physics.illinois.edu (courses.physics.illinois.edu)... 
130.126.151.14
Connecting to courses.physics.illinois.edu (courses.physics.illinois.edu)|130.126.151.14|:443... 
connected.
HTTP request sent, awaiting response... 
200 OK
Length: 53041362 (51M) [application/x-zip-compressed]
Saving to: ‘TestFiles.zip’


TestFiles.zip         0%[                    ]       0  --.-KB/s               
TestFiles.zip         0%[                    ] 152.00K   717KB/s               
TestFiles.zip         1%[                    ] 824.00K  1.90MB/s               
TestFiles.zip         5%[>                   ]   2.84M  4.47MB/s               
TestFiles.zip        18%[==>                 ]   9.38M  11.1MB/s               
TestFiles.zip        41%[=======>            ]  20.96M  19.8MB/s               
TestFiles.zip        64%[===========>        ]  32.57M  25.7MB/s               
TestFiles.zip        87%[================>   ]  44.38M  30.0MB/s               
TestFiles.zip       100%[===================>]  50.58M  31.8MB/s    in 1.6s    

2024-11-07 16:40:57 (31.8 MB/s) - ‘TestFiles.zip’ saved [53041362/53041362]

Archive:  TestFiles.zip
  inflating: Microscopic_after_move.dat  
  inflating: Microscopic.dat         
  inflating: rho_after_Micro2Macro.dat  
  inflating: Microscopic_after_bounce.dat  
  inflating: Microscopic2.dat        
  inflating: Microscopic_after_moveDensity.dat  
  inflating: Microscopic_after_boundary.dat  
  inflating: u_after_Micro2Macro.dat  
  inflating: Microscopic_after_collision.dat  
  inflating: Microscopic_after_equilibrium.dat  

Python Warmup#

Useful functions:#

Before we begin, we want to highlight some functions that you may find useful in this assignment. The last one we’ll explore in more depth as well.

  • np.zeros([dimensions]) - create an array filled with all zeros. *Ex: n = np.zeros([9,10,10]) creates a (9,10,10) array of all zeros

  • np.ones([dimensions]) - same thing as np.zeros but with all ones.

  • np.empty([dimensoons]) - same thing as np,zeros but fills array with garbage. Very fast compared to the above.

  • np.roll(x,shift,axis) - shift array x by shift along axis.

Warmup - using np.roll#

Something we’ll see later in this assignment is the need to shift elements right,left,up,down, and diagonally. Let’s make a little test problem to illustrate how this works.

# Run me!
x = np.arange(15).reshape(3,5)
print(x)
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]]

Now we’ll want to shift every element in x to the right, where we wrap everything around the edge. Rather than figuring this out with a for loop, we’ll employ np.roll to do it fast!

print(x,'\n')
print(np.roll(x,(0,1),(0,1)))
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]] 

[[ 4  0  1  2  3]
 [ 9  5  6  7  8]
 [14 10 11 12 13]]

Now let’s disect that command. The first (0,1) means don’t shift the “y” direction, and to shift the “x” direction by +1. The next command (0,1) means to apply the shift to the two axes.
Now let’s go left:

print(x,'\n')
print(np.roll(x,(0,-1),(0,1)))
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]] 

[[ 1  2  3  4  0]
 [ 6  7  8  9  5]
 [11 12 13 14 10]]

🦉 Now we’ll try to implement a section of the Cha-Cha Slide(Lyrics)

♫Now it's time to get funky
To the right now, to the left
Take it back now y'all♪

Make your x array go right, then left, then “down”. Your output should be:

[[ 4  0  1  2  3]
 [ 9  5  6  7  8]
 [14 10 11 12 13]]
[[ 1  2  3  4  0]
 [ 6  7  8  9  5]
 [11 12 13 14 10]]
[[ 5  6  7  8  9]
 [10 11 12 13 14]
 [ 0  1  2  3  4]]

(Notice how down isn’t what you may expect…)

#!Start
print(np.roll(x,(0,1),(0,1)))
print(np.roll(x,(0,-1),(0,1)))
print(np.roll(x,(-1,0),(0,1)))
#!Stop
[[ 4  0  1  2  3]
 [ 9  5  6  7  8]
 [14 10 11 12 13]]
[[ 1  2  3  4  0]
 [ 6  7  8  9  5]
 [11 12 13 14 10]]
[[ 5  6  7  8  9]
 [10 11 12 13 14]
 [ 0  1  2  3  4]]

🦉In a special version of the song, the next line goes

♫ Take it diagonally up and to the right now y'all

Implement this array move, you should get:

[[14 10 11 12 13]
 [ 4  0  1  2  3]
 [ 9  5  6  7  8]]
print(np.roll(x,(1,1),(0,1))) #!#
[[14 10 11 12 13]
 [ 4  0  1  2  3]
 [ 9  5  6  7  8]]

Now you’re a ~~Cha-Cha Slide~~ np.roll master!

Exercise 1: Fluid Dynamics#

  • List of collaborators:

  • References you used in developing your code:

In this assignment you will produce a fluid dynamics simulations. We are going to develop a number of functions which you will store in the cell below (as well as some global variables that we will use throughout this exercise).

The global variables are

nx=520; ny=180 #size of our simulation
Answer (start)
### ANSWER HERE
Answer (end)

If you have a function that takes in an array n, always make a new variable to return:

def doMath(n):
    # nout = n.copy() #notice the .copy(), so that they aren't the same object
    nout = 2*n+5
    return nout

a. Making some obstacles#

To simulate the real world on the computer, we’ll break the tunnel into \(n_x\times n_y\) squares which we’ll call voxels.

Your simulation will be doing fluid dynamics in a \(n_x \times n_y\) “tunnel” which has some objects within it. These objects will be represented by a boolean \(n_x \times n_y\) size array obstacle=np.empty((nx,ny),dtype='bool') which should be True where the obstacle is and False otherwise. For example,

obstacle[:,:]=False
obstacle[50:70,50:70]=True

would set up a square \(20 \times 20\) object inside our tunnel.

🦉Add a GenerateCylinderObstacle() function to the “function-cell” above which generates a object which is a cylinder located at \((n_x/4,n_y/2)\) with a radius of \(n_y/9\)

If you do plt.matshow(GenerateCylinderObstacle()) it should look like this

Answer (start)
Do not use the word object as a variable!
Notice how object is green in a cell - this means it is a special keyword.

Test:#

## Call your function here to see if if works. 
plt.matshow(GenerateCylinderObstacle().transpose())
plt.title("Cool Name plot")
#plt.savefig("images/object.png")
plt.show()
            
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 2
      1 ## Call your function here to see if if works. 
----> 2 plt.matshow(GenerateCylinderObstacle().transpose())
      3 plt.title("Cool Name plot")
      4 #plt.savefig("images/object.png")

NameError: name 'GenerateCylinderObstacle' is not defined

b. Microscopic velocities \(v_i\)#

To simulate the real world on the computer, we’ll break the tunnel into \(n_x\times n_y\) squares which we’ll call voxels.

The key quantity in your simulation is nine microscopic degrees of freedom

  • \(n_k(i,j)\) where \(0\leq k \leq 8\) and \((i,j)\) are over the \(n_x \times n_y\) voxels of your simulation

which are going to be nine different densities (for each voxel \((i,j)\)) which correspond to the density of a fluid at \((i,j)\) moving in nine different directions (or velocities \(v_k\)):

  • \(n_0\): stationary fluid moving at velocity \(v_0=(0,0)\)

  • \(n_1\): fluid moving up moving at velocity \(v_1=(0,1)\)

  • \(n_2\): fluid moving down at velocity \(v_2=(0,-1)\)

  • \(n_3\): fluid moving right at velocity \(v_3=(1,0)\)

  • \(n_4\): fluid moving left at velocity \(v_4=(-1,0)\)

  • \(n_5\): fluid moving left-down at velocity \(v_5=(-1,-1)\)

  • \(n_6\): fluid moving left-up at velocity \(v_6=(-1,1)\)

  • \(n_7\): fluid moving right-down at velocity \(v_7=(1,-1)\)

  • \(n_8\): fluid moving right-up at velocity \(v_8=(1,1)\)

It will be useful to have access to these microscopic velocities as global variables. 🦉Go ahead and add to the top of your “function-cell” the various velocities:

v=np.zeros((9,2),dtype='int')
v[0,:]=[0,0]
v[1,:]=[0,1]
...

Plot these (below) as

for i in range(0,9):
    plt.arrow(0,0,v[i,0],v[i,1],head_width=0.05, head_length=0.1, fc='k', ec='k')
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.show()

and make sure it looks like

Answer (start)

Test:#

## Plot your microscopic velocities

for i in range(0,9):
    plt.arrow(0,0,v[i,0],v[i,1],head_width=0.05, head_length=0.1, fc='k', ec='k')
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
#plt.savefig("images/densityDirection.png")
plt.show()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/4062129112.py in <module>
      2 
      3 for i in range(0,9):
----> 4     plt.arrow(0,0,v[i,0],v[i,1],head_width=0.05, head_length=0.1, fc='k', ec='k')
      5 plt.xlim(-1.2,1.2)
      6 plt.ylim(-1.2,1.2)

NameError: name 'v' is not defined

We will represent the microscopic densities by a variable

  • n=np.zeros(9,nx,ny)

The high level algorithm for fluid dynamics is very simple. After setting up the initial conditions, we loop many times doing

  • Adjust boundary conditions

  • Collide the microscopic densities

  • Move the microscopic densities

c. Computing macroscopic quantities from the microscopic density#

Given the microscopic densities, there are two macroscopic quanties:

  • the macroscopic density \(\rho(i,j)=\sum_k n_k\) (size: \(n_x \times n_y\))

  • the macroscopic velocity \(\vec{u}(i,j) \equiv (u_x(i,j),u_y(i,j))\) (size: \(2 \times n_x \times n_y\)) where

    • \(u_x(i,j) = 1/\rho \sum_k v_{k,x} n_k(i,j)\)

    • \(u_y(i,j) = 1/\rho \sum_k v_{k,y} n_k(i,j)\)

which you will compute using the function (you will write) (rho,u) = Micro2Macro(n)

🦉Go ahead and write this function and add it to your ‘function-cell’

You can test it using these lines of code:

testFunc(Micro2Macro,["Microscopic.dat"],["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"])

and see if it print two zeros (which is a success!)

Tip: Because u is \((2,n_x,n_y)\), we can get ux by doing u[0] (which is then a \((n_x,n_y)\) array) and uy by doing u[1].

optional - np.einsum#

Figuring out the summations effeciently can be a pain, but we can use np.einsum help out here. Figure out how to rewrite the above problem with einstein summation notation. To use the function, if we have $\( A_{ijk}\cdot B_{dk} = C_{ijd}, \)$ we’d use then call:

C = np.einsum("ijk,dk->ijd",A,B)

Test:#

#!Start
#import pickle
#np.random.seed(30)
#microscopicDensity=np.random.random((9,nx,ny))
#pickle.dump( microscopicDensity, open( "Microscopic.dat", "wb" ) )
#(rho,u)=Micro2Macro(microscopicDensity)
#pickle.dump( rho, open( "rho_after_Micro2Macro.dat", "wb" ) )
#pickle.dump( u, open( "u_after_Micro2Macro.dat", "wb" ) )
#!Stop
testFunc(Micro2Macro,["Microscopic.dat"],["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/2093376147.py in <module>
      8 #pickle.dump( u, open( "u_after_Micro2Macro.dat", "wb" ) )
      9 #!Stop
---> 10 testFunc(Micro2Macro,["Microscopic.dat"],["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"])

NameError: name 'Micro2Macro' is not defined

d. Getting the equilibrium microscopic densities.#

Given the macroscopic densities \(\rho\) and \(\vec{u}\), there is a equilibrium microscopic densities \(n_{eq}\) (size: \(9 \times n_x \times n_y\)) which you will compute using n_eq=Macro2Equilibrium(rho,u)

The relevant formula for this is

\[ n^{eq}(v_k) \equiv n_k^{eq} = \omega_k \rho \left(1 + 3 \vec{v_k} \cdot \vec{u} + \frac{9}{2}(\vec{v_k} \cdot \vec{u})^2- \frac{3}{2}(\vec{u}\cdot \vec{u}) \right) \tag{A} \]

where

  • \(\omega_{0}=4/9\)

  • \(\omega_{1-4}=1/9\)

  • \(\omega_{5-8}=1/36\)

Given some microscopic densities \(n\) you should be able then to figure out the equilibrium microscopic densities \(n_{eq}\).

🦉Write the functions Macro2Equilibrium adding them to your “function-cell”

You can test it by checking if

testFunc(Macro2Equilibrium,["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"],["Microscopic_after_equilibrium.dat"])

is equal to zero.

Test:#

#pickle.dump( microscopicDensity, open( "Microscopic_after_equilibrium.dat", "wb" ) ) #!#
testFunc(Macro2Equilibrium,["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"],["Microscopic_after_equilibrium.dat"])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/2887498537.py in <module>
      1 #pickle.dump( microscopicDensity, open( "Microscopic_after_equilibrium.dat", "wb" ) ) #!#
----> 2 testFunc(Macro2Equilibrium,["rho_after_Micro2Macro.dat","u_after_Micro2Macro.dat"],["Microscopic_after_equilibrium.dat"])

NameError: name 'Macro2Equilibrium' is not defined

d. Implementing collision#

Once you can compute these quantities, the next step is to implement the collision step.

🦉Write a function Collision(n,obstacle) which returns the density n after collision. Put it in your function-cell above.

When we collide the microscopic densities, we get new microscopic densities given as

  • nout = n * (1-omega) + omega * neq

Here are the steps for the function: \(\phantom{}\)

  1. Calculate neq
    a. take the microscopic densities \(\rightarrow\) compute the macroscopic density and velocity (rho,u) = Micro2Macro(n)
    b. take the macroscopic density and velocity \(\rightarrow\) compute the equilibrium microscopic densities neq=Macro2Equilibrium(rho,u)

  2. Calculate nout = n * (1-omega) + omega * neq
    a. \(\omega\) is the viscosity parameter. We will use \(\omega =1.9572953736654806\)

  3. We moved fluid where the obstacle is, so we need to undo this. Reset all nout where the obstacle is to be the same as n

Test with

testFunc(lambda x: Collision(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_collision.dat"])

Test:#

#pickle.dump( Collision(microscopicDensity,GenerateCylinderObstacle()), open( "Microscopic_after_collision.dat", "wb" ) ) #!#

testFunc(lambda x: Collision(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_collision.dat"])
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/727106097.py in <module>
      1 #pickle.dump( Collision(microscopicDensity,GenerateCylinderObstacle()), open( "Microscopic_after_collision.dat", "wb" ) ) #!#
      2 
----> 3 testFunc(lambda x: Collision(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_collision.dat"])

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in testFunc(func, inFiles, outFiles)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in <listcomp>(.0)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

FileNotFoundError: [Errno 2] No such file or directory: 'Microscopic.dat'

e. Moving#

We’ll now move on to getting the fluid to move. This will be a function called, surprise, Move(n,obstacle). But first we need to helper functions:

  • Bounce(n,obstacle)

    • bounce the velocities - that is every velocity where the obstacle is gets reversed

  • MoveDensity(n)

    • Move all your densities in the direction of the velocity (hint the velocity tells you where to move it)

    • Assume periodic boundary conditions here

Thus a Move(n,obstacle) call is just a Bounce then MoveDensity. 🦉Write this function and test!

Test:

testFunc(lambda x: Bounce(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_bounce.dat"])

testFunc(MoveDensity,["Microscopic.dat"],["Microscopic_after_moveDensity.dat"])

testFunc(lambda x: Move(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_move.dat"])

Tests:#

#!Start
#pickle.dump( MoveDensity(microscopicDensity), open( "Microscopic_after_moveDensity.dat", "wb" ) )
##pickle.dump( MoveDensity(microscopicDensity,GenerateCylinderObstacle()), open( "Microscopic_after_collision.dat", "wb" ) )
##print(np.max(np.abs(Collision(microscopicDensity,GenerateCylinderObstacle())-pickle.load( open( "Microscopic_after_collision.dat", "rb" ) ))))
#microscopicDensity=pickle.load( open( "Microscopic.dat", "rb" ) )
#pickle.dump( Bounce(microscopicDensity,GenerateCylinderObstacle()), open( "Microscopic_after_bounce.dat", "wb" ))
#microscopicDensity=pickle.load( open( "Microscopic.dat", "rb" ) )
#pickle.dump( Move(microscopicDensity,GenerateCylinderObstacle()), open( "Microscopic_after_move.dat", "wb" ) )
#!Stop

testFunc(lambda x: Bounce(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_bounce.dat"])
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/1549029420.py in <module>
      9 #!Stop
     10 
---> 11 testFunc(lambda x: Bounce(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_bounce.dat"])

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in testFunc(func, inFiles, outFiles)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in <listcomp>(.0)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

FileNotFoundError: [Errno 2] No such file or directory: 'Microscopic.dat'
testFunc(MoveDensity,["Microscopic.dat"],["Microscopic_after_moveDensity.dat"])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3243438587.py in <module>
----> 1 testFunc(MoveDensity,["Microscopic.dat"],["Microscopic_after_moveDensity.dat"])

NameError: name 'MoveDensity' is not defined
testFunc(lambda x: Move(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_move.dat"])
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/1149084007.py in <module>
----> 1 testFunc(lambda x: Move(x,GenerateCylinderObstacle()),["Microscopic.dat"],["Microscopic_after_move.dat"])

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in testFunc(func, inFiles, outFiles)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/3739739110.py in <listcomp>(.0)
     20 
     21 def testFunc(func,inFiles,outFiles):
---> 22     inputs  = [pickle.load(open(f,"rb")) for f in inFiles]
     23     outputs = [pickle.load(open(f,"rb")) for f in outFiles]
     24     result  = func(*inputs)

FileNotFoundError: [Errno 2] No such file or directory: 'Microscopic.dat'

f. Boundary conditions#

In the y-direction, we are just going to assume that the boundary conditions are periodic. For the x-direction, we are going to assign certain boundary conditions on the left and right.

Let’s call this function FixBoundary(n,n_init) which applies these boundary conditions and returns the new density n.

Left: we are going to assume that there is a flow - this means that the microscopic densities are the same at each step. Therefore, what we should do is simply replace the current densities on the left-most row with the initial microscopic densities (n_init).

Right: we are going to assume that the gradient is zero - i.e. the important physics has dissapeared by this point. To do this, we will set

# loop over all microscopic velocities vL that are going left 
   nout[vL,-1,:] = nout[vL,-2,:] 

Hint: there are 3 vLs.

🦉Put all of this together and run the next test!

Testing:

testFunc(FixBoundary,["Microscopic.dat","Microscopic2.dat"], ["Microscopic_after_boundary.dat"])

Test:#

#!Start
#microscopicDensity2=np.random.random((9,nx,ny))
#pickle.dump(microscopicDensity2,open( "Microscopic2.dat","wb"))
#pickle.dump( FixBoundary(microscopicDensity,microscopicDensity2), open( "Microscopic_after_boundary.dat", "wb" ) )
#!Stop
testFunc(FixBoundary,["Microscopic.dat","Microscopic2.dat"], ["Microscopic_after_boundary.dat"])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/239426050.py in <module>
      4 #pickle.dump( FixBoundary(microscopicDensity,microscopicDensity2), open( "Microscopic_after_boundary.dat", "wb" ) )
      5 #!Stop
----> 6 testFunc(FixBoundary,["Microscopic.dat","Microscopic2.dat"], ["Microscopic_after_boundary.dat"])

NameError: name 'FixBoundary' is not defined

g. Setting up the initial conditions#

We need to generate an initial conditions. The way that we are going to do that is

  • pick some macroscopic density \(\rho\) of size \((n_x,n_y)\) (uniformly equal to 1.0) and a

  • macroscopic velocity \(\vec{u}(i,j)\) (sizes \((2,n_x,n_y)\))

    • 0.04*(1.0+1e-4*np.sin(y/(ny)*2*np.pi)) in the x-direction (aka 0). This introduces a very tiny anisotropy to the system.

    • zero in the y-direction (aka 1)

  • compute the equilibrium density \(n_{eq}\) associated with the initial density and velocities.
    Use (and return) this as the initial conditions

🦉Write a function Setup() which returns the initial microscopic densities. Remember to save this as you will need it to apply your boundary conditions.

I get the following when I visualize it as plt.matshow(Setup()[3].transpose())

Answer (start)

A quick test here is just to check that the following:

print("All good?",(np.abs(np.max(Setup())-0.443377991)<1e-5) and ((np.abs(np.min(Setup())-0.0245774711)<1e-5)))

Test:#

print("All good?",(np.abs(np.max(Setup())-0.443377991)<1e-5) and ((np.abs(np.min(Setup())-0.0245774711)<1e-5)))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/4075721444.py in <module>
----> 1 print("All good?",(np.abs(np.max(Setup())-0.443377991)<1e-5) and ((np.abs(np.min(Setup())-0.0245774711)<1e-5)))

NameError: name 'Setup' is not defined

h. Putting it all together#

🦉Write a function Run(steps,record,n,n_init,obstacle) (up in your function-cell) which takes

  • the number of steps steps you want to run

  • how often you want to record the velocity and density

  • the microscopic densities n to start running from

  • as well as the initial microscopic densities n_init for the boundary conditions (often equal to n.copy())

  • and the obstacle boolean array.

Each step consists of:

  1. Adjust boundary conditions

  2. Collide the microscopic densities

  3. Move the microscopic densities

It should return

  • a list of macroscopic densities rhos

  • a list of macroscopic velocities us

    • This will be u2=\(\sqrt{\vec{u}\cdot\vec{u}}\)

    • To help us later, store these as u2.transpose()

  • the last microscopic density n

Now (below) call your Setup and Run functions. Then plot the final configuration (something like plt.matshow(us[-1])).

🦉Run for 2001 steps recording every 100 steps.

Answer:#

#!Start
fin=Setup()
feq_init=fin.copy()
obstacle=GenerateCylinderObstacle()
#!Stop
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/1772239070.py in <module>
      1 #!Start
----> 2 fin=Setup()
      3 feq_init=fin.copy()
      4 obstacle=GenerateCylinderObstacle()
      5 #!Stop

NameError: name 'Setup' is not defined
#!Start
fins=rhos=us= [None] * 21
fins[0]=fin.copy()
for i in range(0,20):
    %time (rhos[i],us[i],fins[i+1])=Run(2001,100,fins[i],feq_init,obstacle)
#!Stop
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/104084777.py in <module>
      1 #!Start
      2 fins=rhos=us= [None] * 21
----> 3 fins[0]=fin.copy()
      4 for i in range(0,20):
      5     get_ipython().run_line_magic('time', '(rhos[i],us[i],fins[i+1])=Run(2000,100,fins[i],feq_init,obstacle)')

NameError: name 'fin' is not defined

i. Animation#

🦉Now write a function AnimateMe(us_flat,vMax) which is going to take a list of velocities us_flat (and maximum for the vMax) and going to return an animation which is going to be then produced by calling

anim = AnimateMe(us,vMax)
HTML(anim.to_jshtml())

This function will look like:

def AnimateMe(us_flat,vMax):
    fig, ax = plt.subplots()
    cax = ax.imshow(us_flat[1],cmap=plt.cm.Reds,vmin=0,vmax=vMax)
    plt.close(fig)
    def animate(i):
         cax.set_array(us_flat[i])

    anim = FuncAnimation(fig, animate, interval=100, frames=len(us_flat))
    return anim

Answer:#

#!Start
us_flat = list(itertools.chain.from_iterable(us[:-1]))
#us_flat = list(itertools.chain.from_iterable(us[0:17]))
anim=AnimateMe(us_flat,0.07)
HTML(anim.to_jshtml())
#!Stop
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/1102925857.py in <module>
      1 #!Start
----> 2 us_flat = list(itertools.chain.from_iterable(us[:-1]))
      3 #us_flat = list(itertools.chain.from_iterable(us[0:17]))
      4 anim=AnimateMe(us_flat,0.07)
      5 HTML(anim.to_jshtml())

TypeError: 'NoneType' object is not iterable

j. Speed Test (EC)#

5 Points

Our code for for %time Run(2000,3000,Setup(),Setup(),GenerateCylinderObstacle()) takes between 45 and 60 seconds. We’ll give 5 points if you can make your code faster than 40 seconds runtime (on the server). Warning: we aren’t sure this is possible!

Answer:#

Exercise 2: Walls#

  • List of collaborators:

  • References you used in developing your code:

In this exercise we are going to use all the same code but generate a different object. Here generate an object which consists of two walls:

  • one spans \(50 \leq x \leq 60\) and \(n_y/4 \leq y \leq 3n_y/4\)

  • the other spans \(200 \leq x \leq 210\) and \(n_y/4 \leq y \leq 3n_y/4\)

🦉Run your code again with this new obstacle and generate a new animation.

Answer:#

#!Start
fin=Setup()
feq_init=fin.copy()

obstacle=GenerateWallObstacle()
fins=rhos=us= [None] * 21
fins[0]=fin.copy()
for i in range(0,20):
    %time (rhos[i],us[i],fins[i+1])=Run(2000,100,fins[i],feq_init,obstacle)
#!Stop
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/165298796.py in <module>
      1 #!Start
----> 2 fin=Setup()
      3 feq_init=fin.copy()
      4 
      5 obstacle=GenerateWallObstacle()

NameError: name 'Setup' is not defined
#!Start
us_flat = list(itertools.chain.from_iterable(us[0:-1]))
anim=AnimateMe(us_flat,0.2)
HTML(anim.to_jshtml())
#!Stop
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/5z/1xxxzf610_7gxv7scl0v99240000gn/T/ipykernel_1449/129420736.py in <module>
      1 #!Start
----> 2 us_flat = list(itertools.chain.from_iterable(us[0:-1]))
      3 anim=AnimateMe(us_flat,0.2)
      4 HTML(anim.to_jshtml())
      5 #!Stop

TypeError: 'NoneType' object is not iterable

Exercise 3: Something new… (EC)#

(Extra credit: 5 points)

  • List of collaborators:

  • References you used in developing your code:

🦉Come up with something new and cool to do with this animation code. Maybe try out some other interesting obstacles are add some friction to the walls or such.

Answer:#

Acknowledgement: This assignment originally inspired from code from flowkit.com

  • Bryan Clark (original)

Copyright: 2021