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-14 08:44:43-- 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 1%[ ] 968.00K 4.66MB/s
TestFiles.zip 20%[===> ] 10.60M 26.2MB/s
TestFiles.zip 71%[=============> ] 36.41M 60.0MB/s
TestFiles.zip 100%[===================>] 50.58M 71.3MB/s in 0.7s
2024-11-14 08:44:44 (71.3 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 zerosnp.ones([dimensions])
- same thing asnp.zeros
but with all ones.np.empty([dimensoons])
- same thing asnp,zeros
but fills array with garbage. Very fast compared to the above.np.roll(x,shift,axis)
- shift arrayx
byshift
alongaxis
.
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 HERE
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
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
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
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{}\)
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 densitiesneq=Macro2Equilibrium(rho,u)
Calculate
nout = n * (1-omega) + omega * neq
a. \(\omega\) is the viscosity parameter. We will use \(\omega =1.9572953736654806\)We moved fluid where the obstacle is, so we need to undo this. Reset all
nout
where the obstacle is to be the same asn
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 vL
s.
🦉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())
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 runhow often you want to
record
the velocity and densitythe microscopic densities
n
to start running fromas well as the initial microscopic densities
n_init
for the boundary conditions (often equal ton.copy()
)and the obstacle boolean array.
Each step consists of:
Adjust boundary conditions
Collide the microscopic densities
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