# Modelling of partcle flow under electric field

This is my first modeling exercise on the particle flow. Recently, I read a paper https://acp.copernicus.org/articles/20/3181/2020/. After reading the paper, thought of modeling it in Python. It seems like that the model came out well but I am not entirely sure. Since I don't know how to fit the parameters alpha and beta to the experimental curve. I directly used those values. I am uncertain of the way I used the call function. Is there a better way to apply it?

"""
Created on Wed Sep  9 20:08:25 2020

@author: vishal
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

class ParticleFlow:

def __init__(self, d, ec, el, rp, rf, g, na):
self.d = d
self.ec = ec
self.el = el
self.rp = rp
self.rf = rf
self.g = g
self.na = na

def __call__(self, v0, t):

def cd():
re = (self.rf*np.linalg.norm(v0, axis = 1)*self.d)/self.na
t = np.zeros(nop)
for i, j in enumerate(re):
t[i] = (24/j)*(1 + 0.15*pow(j,0.681)) + 0.407/(1+(8710/j))
return t

def fg():
fgc = np.zeros(v0.shape)
fgc[:,1] = (np.pi*pow(self.d,3)*self.rp*self.g)/6
return fgc

def fe():
fec = np.zeros(v0.shape)
fec[:,1] = np.pi*pow(self.d,2)*self.ec*self.el
return fec

def fd():
fdc = np.zeros(v0.shape)
fdc [:,0] = -(1/8)*(cd()*self.rf*np.pi*pow(self.d,2)*pow(v0[:,0],2))*np.reciprocal(np.linalg.norm(v0, axis = 1))
fdc [:,1] = -(1/8)*(cd()*self.rf*np.pi*pow(self.d,2)*pow(v0[:,1],2))*np.reciprocal(np.linalg.norm(v0, axis = 1))
return (fdc*v0)

a = np.zeros(v0.shape)
a[:,0] = (6/(np.pi*pow(self.d,3)*rp))*(fg()[:,0] + fe()[:,0] + fd()[:,0])
a[:,1] = (6/(np.pi*pow(self.d,3)*rp))*(fg()[:,1] + fe()[:,1] + fd()[:,1])

return a

class ExplicitEuler:
"""
Euler scheme for the numerical resolution of
a differentiel equation.
"""
def __init__(self,f):
self.f = f

def iterate(self,v0,t,dt):
return v0+dt*self.f(v0,t)

class Integrator:
"""
Integration of a differential equation between tMin and tMax
with N discretization steps and x0 as an initial condition
"""
def __init__(self,method,v0,tMin,tMax,N):
self.v0   = v0
self.tMin = tMin
self.tMax = tMax
self.dt   = (tMax - tMin)/(N-1)

self.f = method

def getIntegrationTime(self):
return np.arange(self.tMin,self.tMax+self.dt,self.dt)

def integrate(self):
v = np.array(self.v0)
for t in np.arange(self.tMin,self.tMax,self.dt):
v = np.append( v, self.f.iterate(v[-nop:],t,self.dt),axis=0)
return v

'''
rp: density of particle (quartz)
rf: density of fluid (air)
da: average diameter of particle
sd: SD for the particle size distribution
g: acceleration due to gravity
ec: electric charge
alpha: SD in electric charge
el: electric field
nop: # of particles

'''

rp = 2600
rf = 1.18
da = 132*1e-06
sd = 45*1e-06
g = 9.8
beta = 11000*1e-12
alpha = 1.4*1e-06
vh = 4.95
dim = 2
na = 1.5*1e-05
tmin = 1e-06
tmax = 18e-03
nop = 100

d0 = np.random.normal(da, sd, nop*10)

d0 = np.random.choice(d0[d0 > 0], nop) #to avoid any negative value of diameter

ec0 = np.random.normal(0, alpha, nop)

# a0 = np.zeros((n, dim))

v0 = np.ones(nop)*vh
#v1 = np.array([norm(0, beta/pow(i,2), i) for i in d0])
v1 = np.array([np.random.normal(0, beta/(i**2), 1) for i in d0])[:,0]
#v1 = np.random.normal(0, beta/(d0**2), nop)

vint = np.dstack((v0, v1))[0]

eul = Integrator(ExplicitEuler(ParticleFlow(d0, ec0, 1e05, rp, rf, g, na)),vint,tmin,tmax,2000)

temp = eul.integrate()

'''
Intial positions of particles set to zero
and subsequent positions were calculated by multiplying with velocity.
'''

ri = np.zeros(vint.shape)
r = ri
k = 0

dt = (tmax - tmin)/(2000-1)

lv = len(np.arange(tmin, tmax, dt))

for i in range(lv):
ri = ri + temp[i*nop:(i+1)*nop]*dt
r = np.append(r, ri, axis =0)

#plt.scatter(r[-100:,0], r[-100:,1], c=np.sign(ec0), cmap="bwr", s = np.abs(ec0)*4e07)
#plt.scatter(d0, v1, c = np.sign(ec0), cmap = 'bwr', s = np.abs(ec0)*4e07)
#plt.xlim(d0.min()-1e-05, d0.max()+1e-05)
#plt.ylim(v1.min()-1, v1.max()+1)
#plt.xlabel('y', fontsize = 20)
#plt.ylabel('x', fontsize = 20)
##plt.savefig('position_1.png', dpi = 600, bbox_inches = 'tight')

# animate

fig = plt.figure()
ax = plt.axes(xlim=(r[:,0].min(), r[:,0].max()), ylim=(r[:,1].min(), r[:,1].max()))
particles, = ax.plot([], [], 'o', ms = 2)

plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.rc('xtick', labelsize=20)
plt.rc('ytick', labelsize=20)
#plt.xlabel(r'x',fontsize=20)
#plt.ylabel(r'$y$', fontsize=20)

def init(): # initialize animation
particles.set_data([], [])
return particles,

def animate(i): # define animation Euler
particles.set_data(r[i*nop:(i+1)*nop,0], r[i*nop:(i+1)*nop,1])
return particles,

anim=animation.FuncAnimation(fig, animate, init_func=init, frames = 2000, interval=2,blit=True,repeat=False)
#anim.save('Particle_flow_1000.mp4',fps=200,dpi=400)

'''
scatter plot to give different colors
based on the polarity of the charge
'''

# scatter

def PlotScatter(i, r):
plt.figure()
plt.ioff()
plt.axes(xlim=(r[:,0].min(), r[:,0].max()), ylim=(r[:,1].min(), r[:,1].max()))
plt.scatter(r[i*nop:(i+1)*nop, 0], r[i*nop:(i+1)*nop,1], c=np.sign(ec0),
cmap="bwr", s = np.abs(ec0)*4e07)
plt.savefig('{}'.format(i), dpi = 200, bbox_inches = 'tight')
plt.clf()
plt.cla()
plt.close()

i = 0
while i <= lv:
PlotScatter(i, r)
i = i + 50


## 1 Answer

TLDR: Vectorize and eliminate redundant calculations. Rewriting most of your code gave a 400-500x speed up

Also, use functions and name your variables in ways that people can read

Long version:

1. Legibility

The code in ParticleFlow, ExplicitEuler, and Integrate was well organize but slow and hard to follow. The variable names were cryptic. They get repeated in different functions with different meanings. Method parameters and global variables are intermixed randomly. And one case a variable is created, explicitly passed through 4 functions, and then thrown away without ever being using.

If there hadn't been an academic paper for reference I don't think I could have figured out what was happening beyond something is getting integrated.

Then the second half of the code was hard to follow because organizing functions and comments suddenly vanish. For example:

for i in range(lv):
ri = ri + temp[i*nop:(i+1)*nop]*dt
r = np.append(r, ri, axis =0)


Is there without any explanation of what it does or why. What is temp and why is it getting sliced between i*nop and (i+1)*nop? In the rewrite I tried to give examples of better naming and structure so I'll stop harping on legibility now and move on the the interesting stuff of

1. Performance

Once I read through the paper the code was based on on I could follow the gist of what it calculates and it found it does so pretty inefficiently. Running on my laptop it takes about 3.4 seconds to use Euler's Method to approximate the solution for 100 particles. The rewrite below does the exact same thing in 0.007 seconds.

The major performance issues are that you repeat calculations based on static variables hundreds/thousands of times and you repeat them in the slowest way possible so your method has to slog through the same number crunching again and again and again.

Look at cd():

   def cd():
re = (self.rf*np.linalg.norm(v0, axis = 1)*self.d)/self.na
t = np.zeros(nop)
for i, j in enumerate(re):
t[i] = (24/j)*(1 + 0.15*pow(j,0.681)) + 0.407/(1+(8710/j))
return t


The beauty of numpy is that it lets you performed vectorized calculations entire arrays at once. That for loop isn't doing anything that couldn't be written in one line. And once you start using vectorized results you can strip out the overhead of creating empty N-length vectors (or rather nop-length, an ambiguously named global variable defined halfway down) and filling it element-by-element every iteration. Vectorize whatever you can. It'll give roughly an order of magnitude speedup for free.

def fe():
fec = np.zeros(v0.shape)
fec[:,1] = np.pi*pow(self.d,2)*self.ec*self.el
return fec


Is better but you're still creating things inefficiently. np.zeros(v0.shape) could be np.zeros_like(v0). You're only using values in fec[:,1] but are creating a 2D array to hold them. But, the thing that really kills efficiency: You go through the whole process of creating an unnecessarily large array, putting some values in it, and then returning it 4,000 separate times for the same result.

The differential equation you're approximating a solution to only depends on the variable V. Everything in fg is a constant defined at run time. That whole function can be rewritten as a single line that gets run once on initialization and saved:

self.fe = np.pi*pow(self.p.d, 2)*self.p.ec*self.env.el


Cutting out out almost a million operations (2D array x 100 elements x 4000 calls = 800,000 allocations) just for the empty array. If you increase the number of particles it saves even more. I can see how you copied the equations from the article but pay attention to what changes and what does to see what you can calculate once and save.

This next example I'm not sure sped anything up much but highlights a way of thinking that can help save time in the long run. You wrote ExplicitEuler as a class:

class ExplicitEuler:
"""
Euler scheme for the numerical resolution of
a differentiel equation.
"""
def __init__(self,f):
self.f = f

def iterate(self,v0,t,dt):
return v0+dt*self.f(v0,t)


When you look at the code that runs it though it becomes apparent that again, nothing it does ever changes. There's no actual state being used. When you start looking for bits and pieces like that you can realize that all you actually want from it is an interface to define a specific function composition (a function using a function). You don't need a class with initialization and methods to do that. You can use a simple function:

def ExplicitEuler(func):
def f(v, dt):
return v+dt*func(v)
return f


Same thing with Integrator. It gets passed some variables once and returns a single result. You don't need a class for it. You can even change them both to functions still keep your moneyshot line:

eul = Integrator(ExplicitEuler(ParticleFlow(d0, ec0, 1e05, rp, rf, g, na)),vint,tmin,tmax,2000)
temp = eul.integrate()


Becomes

vv = Integrate(ExplicitEuler(ParticleFlow(p, e)), p.v0, tmin, tmax, nop)


with just a bit of refactoring.

You can read through my full rewrite below. There are a few odds and ends I won't go into detail on (grouping variables as tuples, adding some functions to spaghetti code sections, adding human readable names and comments, etc).

But again, the three big things to pay attention to are: Vectorization, Redudant/Expensive Calculation, and Legibility

P.S. There are two lines you might notice I cut completely. When I read the article it stated: > In the horizontal direction, we do not consider forces on the particles and we assumed that the particles travel at the velocity of the wind. < It wasn't entirely clear in your code but you seem to partially calculate horizontal drag on the particles

# author: cmelgreen
# based on vishal's codereview question at
# https://codereview.stackexchange.com/questions/249239/modelling-of-partcle-flow-under-electric-field

import numpy as np
import scipy.stats as stats
from collections import namedtuple
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import time

# Setup Particle Arrays and Environemntal Consts
def NewParticles(mean, sd, alpha, beta, n):
# Truncate Norm Dist at 0
diameter = stats.truncnorm((0 - mean) / sd, np.inf, loc=mean, scale=sd).rvs(n)
charge = np.random.normal(0, alpha, n)
velocity_0 = np.random.normal(size=n)*(beta/(diameter**2))

return namedtuple('Particles', 'd ec v0')(diameter, charge, velocity_0)

def NewEnvironment(rf, na, rp, g, el, vh):
return namedtuple('Environment', 'rf na rp g vh el')(rf, na, rp, g, vh, el)

# Setup DiffEq for Particle motion in given Environment
# Vectorize as much as possible for performance
class ParticleFlow:
def __init__(self, particles, environment):
self.p = particles
self.env = environment

# Create const horizontal velocity array for vectorization
self.vh = np.zeros_like(particles.d)*self.env.vh

# Forces independent of velocity are constant and only need calcuated once
self.fg = (np.pi*pow(self.p.d, 3)*self.env.rp*self.env.g)/6
self.fe = np.pi*pow(self.p.d, 2)*self.p.ec*self.env.el

def __call__(self, v):
# Drag Coefficient
def cd(v):
re = (self.env.rf *
np.linalg.norm([self.vh, v])*self.p.d)/self.env.na
return (24/re)*(1 + 0.15*pow(re, 0.681)) + 0.407/(1+(8710/re))

# Drag
def fd(v):
return (1/8)*(cd(v)*self.env.rf*np.pi*pow(self.p.d, 2)*pow(v, 2))*np.reciprocal(np.linalg.norm([self.vh, v]))

# Return Total Forces
return (6/(np.pi*pow(self.p.d, 3)*self.env.rp))*(self.fg + self.fe + fd(v))

def ExplicitEuler(func):
# Wrap ParticleFlow for Integration
def f(v, dt):
return v+dt*func(v)
return f

def Integrate(method, v0, tMin, tMax, n):
dt = (tMax - tMin)/(n-1)
v = np.zeros((int((tMax - tMin) / dt), n))
p = np.zeros((2, *v.shape))
v[0] = v0

# Calculate Integral
for i in range(1, len(v)):
v[i] = method(v[i - 1], dt)

return v

# Used solved velocities to create array of Particle positions over time
def PositionOverTime(vv, vh, tMin, tMax):
p = np.zeros((*vv.shape, 2))
dt = (tMax - tMin) / (len(vv) - 1)

for i in range(1, len(p)):
p[i, :, 0] = p[i-1, :, 0] + vh*dt
p[i, :, 1] = p[i-1, :, 1] + vv[i-1]*dt

return p

def GeneratePoints():
# Bunch of constants that could be passed in as tuples but hardcoded for example
(rf, na, rp, g, el, vh) = (1.18, 1.5*1e-05, 2600, 9.8, 1e05, 4.95)
(da, sd, alpha, beta, nop) = (132 * 1e-06, 45 * 1e-06, 1.4 * 1e-06, 11000 * 1e-12, 100)
(tmin, tmax) = (1e-06, 18e-03)

# Time Setup and Euler Method
start = time.time()

# Create Environment and Particle Arrays
e = NewEnvironment(rf, na, rp, g, el, vh)
p = NewParticles(da, sd, alpha, beta, nop)

# Run Euler Method
vv = Integrate(ExplicitEuler(ParticleFlow(p, e)), p.v0, tmin, tmax, nop)

end = time.time()
print("Ran in: ", end - start)

return PositionOverTime(vv, np.ones(nop) * e.vh, tmin, tmax), p.d, p.ec

# Plot the points based on diameter and charge
def Plot(p, d, c):
frames = 99
interval = 65
scale = 1.7e6
alpha = .67

def init():
return scatter,

def animate(i):
X = np.c_[p[i, :, 0], p[i, :, 1]]
scatter.set_offsets(X)
return scatter,

# Scale marker size array once instead of recalculaing for each frame
s = (d * scale)

# Setup figure and axes
fig = plt.figure()
ax = plt.axes(
xlim=(p[0, :, 0].min(), p[-1, :, 0].max()),
ylim=(p[-1, :, 1].min(), p[-1, :, 1].max())
)
scatter = ax.scatter(p[1, :, 0], p[1, :, 1], c=c, cmap='bwr', s=s, alpha=alpha)

# Start animation and display
_ = animation.FuncAnimation(fig, animate, init_func=init, frames=frames, interval=interval, blit=True, repeat=False)
plt.show()

# Run
if __name__ == "__main__":
Plot(*GeneratePoints())


• Thanks a lot @Coupcoup. Yes, I saw that the fg and fe do not change over time... so need to calculate them repeatedly... Sep 13 '20 at 7:53
• I used temp to store velocity at various dt steps... and used those to calculate position r. The variable naming was really terrible at many places and I invoked variables like k = 0 but didn't use them... I will keep in mind those things.... and also give proper explanations for the terms... Sep 13 '20 at 7:58
• youtube.com/watch?v=FfKRoBn5ISM - This simulation was produced using the code that I provided... you can see clearly that particles disperse both in horizontal and vertical direction...The colors represent polarity and size the amount of electric charge, respectively. Sep 13 '20 at 8:45
• And one last thing I forgot to mention was generating the particle distributions! If you check out NewParticles() you'll see how to sample directly from a truncated normal distribution instead of the hack-y sampling and throwing out negatives. I also worked using beta directly into a vectorized sampling function. When you have problems like that generally you can sample from a standard normal dist and then simply multiply the vector by whatever standard deviation you want to reshape it and then add/subtract the mean i.e. doing the opposite of normalization Sep 13 '20 at 16:34
• After adding horizontal component of drag force. This is the final result youtu.be/0KVvG16CK1k Total time taken was 0.014 sec. Now, I am thinking of adding collision between particles and use momentum conservation for that. Thank you, again. Sep 13 '20 at 20:57