Hi everybody,
I wrote a simulation code for Plasma Actuator. This is based on my previous
code written in Python and uses Suzen model combined with Navier-Stokes
equations. It's about 5 times faster than Python but still simulation takes
days to finish. I am new to Julia and I tried to follow the performance
tips but so far this is how much I could do. Any suggestions for improving
speed are highly appreciated. This is the code:
using PyPlot
using PyCall
const nx=441
const ny=441
const nt=100
const Length=0.011#11 mm
const dx=Length/(nx-1)
const dy=Length/(ny-1)
x=linspace(-Length/2,Length/2,nx)
y=linspace(-Length/2,Length/2,ny)
################################
phi1 = zeros(ny,nx) #potential
eps=zeros(ny,nx)#permitivity
roc=ones(ny,nx)#charge
ld=ones(ny,nx)#Debye length
############### All Electrodes ########################
# initial boundary conditions
phi1[221:222,49:57]=1400.0
phi1[218:219,57:97]=0.0
#phi1[221:222,97:105]=0.0
phi1[221:222,145:153]=1400.0
phi1[218:219,153:193]=0.0
#phi1[221:222,193:201]=0.0
#phi1[221:222,241:249]=0.0
phi1[218:219,249:289]=0.0
phi1[221:222,289:297]=1400.0
#phi1[221:222,337:345]=0.0
phi1[218:219,345:385]=0.0
phi1[221:222,385:393]=1400.0
eps[1:220,:]=2.7
eps[221:end,:]=1.0
eps[221:221,:]=(eps[222:222,1:end].*eps[220:220,1:end])./(((dx/(2*dx))*eps[222:222,1:end])+
(((dx/(2*dx))*eps[220:220,1:end])))
###############################################################
#convergence parameters
to1=1e-8
###Gauss-Seidel solver
max_phi1_diff=1
while max_phi1_diff > to1
phi1_old=copy(phi1)
phi1[2:nx-1,2:ny-1]=((eps[2:nx-1,2:ny-1].*((phi1[3:nx,2:ny-1]/(dx^2))+
(phi1[2:nx-1,3:ny]/(dy^2)))+((eps[1:nx-2,2:ny-1].*phi1[1:nx-2,2:ny-1])/(dx^2))+
((eps[2:nx-1,1:ny-2].*phi1[2:nx-1,1:ny-2]/(dy^2)))))./(((eps[2:nx-1,2:ny-1]+
eps[1:nx-2,2:ny-1])/(dx^2))+
((eps[2:nx-1,2:ny-1]+eps[2:nx-1,1:ny-2])/(dy^2)))
phi1[end,:]=phi1[end-1,:]
phi1[1,:]=phi1[2,:]
phi1[:,end]=phi1[:,end-1]
phi1[:,1]=phi1[:,2]
phi1[221:222,49:57]=1400.0
phi1[218:219,57:97]=0.0
#phi1[221:222,97:105]=0.0
phi1[221:222,145:153]=1400.0
phi1[218:219,153:193]=0.0
#phi1[221:222,193:201]=0.0
#phi1[221:222,241:249]=0.0
phi1[218:219,249:289]=0.0
phi1[221:222,289:297]=1400.0
#phi1[221:222,337:345]=0.0
phi1[218:219,345:385]=0.0
phi1[221:222,385:393]=1400.0
eps[1:220,:]=2.7
eps[221:end,:]=1.0
phi1_diff = phi1.-phi1_old
max_phi1_diff=maximum(abs(phi1_diff))
end
fig = plt.figure(figsize = (11,7), dpi=100)
plt.contourf(x*1000,y*1000,phi1,50)
cbar = plt.colorbar()
savefig("wflow1368AC29SeptemberElectricField.png")
################ CHARGE ####################
ld[221:end,:]=0.00017 #Debye lenght
ld[1:220,:]=10000.0#Debye lenght
roc[218:219,57:97]=0.00750
roc[218:219,153:193]=0.00750
roc[218:219,249:289]=0.00750
roc[218:219,345:385]=0.00750
eps[1:220,:]=2.7
eps[221:end,:]=1.0
eps[221:221,:]=(eps[222:222,1:end].*eps[220:220,1:end])./(((dx/(2*dx))*eps[222:222,1:end])+
(((dx/(2*dx))*eps[220:220,1:end])))
roc[:,end] =0
roc[end,:] = 0
roc[:,1] =0
roc[1,:] = 0
#convergence parameters
to1=1e-8
###Gauss-Seidel solver
max_roc_diff=1
while max_roc_diff > to1
roc_old=copy(roc)
roc[2:nx-1,2:ny-1]=((eps[2:nx-1,2:ny-1].*((roc[3:nx,2:ny-1]/(dx^2)).+(roc[2:nx-1,3:ny]/dy^2)))+
((eps[1:nx-2,2:ny-1].*roc[1:nx-2,2:ny-1])/dx^2)+
((eps[2:nx-1,1:ny-2].*roc[2:nx-1,1:ny-2])/dy^2))./(((eps[2:nx-1,2:ny-1].+eps[1:nx-2,2:ny-1])/dx^2)+
((eps[2:nx-1,2:ny-1]+eps[2:nx-1,1:ny-2])/dy^2)+(1./(ld[2:nx-1,2:ny-1].^2)))
ld[221:end,:]=0.00017 #charge
ld[1:220,:]=10000.0#Debye lenght
roc[218:219,57:97]=0.00750
roc[218:219,153:193]=0.00750
roc[218:219,249:289]=0.00750
roc[218:219,345:385]=0.00750
roc_diff = roc.-roc_old
max_roc_diff=maximum(abs(roc_diff))
end
fig = plt.figure(figsize = (11,7), dpi=100)
plt.contourf(x*1000,y*1000,roc,50)
cbar = plt.colorbar()
savefig("wflow1368AC29SeptemberCharge.png")
######### BODY FORCE #######
F1=ones(ny,nx)
Fx1=ones(ny,nx)
Fy1=ones(ny,nx)
F1[2:nx-1,2:ny-1]=roc[2:nx-1,2:ny-1].*(-(((phi1[2:nx-1,2:ny-1]-phi1[1:nx-2,2:ny-1])/dx)+((phi1[2:nx-1,2:ny-1]-phi1[2:nx-1,1:ny-2])/dy)))
F1[221:222,49:57]=0.0
F1[218:219,57:97]=0.0
F1[221:222,97:105]=0.0
F1[221:222,145:153]=0.0
F1[218:219,153:193]=0.0
F1[221:222,193:201]=0.0
F1[221:222,241:249]=0.0
F1[218:219,249:289]=0.0
F1[221:222,289:297]=0.0
F1[221:222,337:345]=0.0
F1[218:219,345:385]=0.0
F1[221:222,385:393]=0.0
F1[1:220,:]=0.0
###### Body force on x and Y ###############
Fx1[2:nx-1,2:ny-1]=roc[2:nx-1,2:ny-1].*(-(((phi1[2:nx-1,2:ny-1]-phi1[1:nx-2,2:ny-1])/dx)))
Fy1[2:nx-1,2:ny-1]=roc[2:nx-1,2:ny-1].*(-((phi1[2:nx-1,2:ny-1]-phi1[2:nx-1,1:ny-2])/dy))
fig = plt.figure(figsize = (11,7), dpi=100)
plt.contourf(x*1000,y*1000,F1,50)
cbar = plt.colorbar()
savefig("wflow1368AC29SeptemberForce.png")
#######################################################
#Navier Stokes
const dt=0.0000025;#time step
const nstep=24001;#number of steps
const mu=0.000018;#kinematic viscosity
const maxit=100;
const beta=1.2;
const h=Length/(nx-1);
u=zeros(nx,ny);
v=zeros(nx,ny);
p=zeros(nx,ny);
ut=zeros(nx,ny);
vt=zeros(nx,ny);
for iter =1:nstep
function Fxn1(iter)
(abs(sin(2*pi*20000*(iter*dt))))
end
Fxn11=Fx1.*Fxn1(iter)
function Fyn1(iter)
(abs(sin(2*pi*20000*(iter*dt))))
end
Fyn11=Fy1.*Fyn1(iter)
for i =2:nx-1
for j=2:ny-1
ut[i,j]=u[i,j]+dt*(-(0.5/h)*(u[i,j].*(u[i+1,j]-u[i-1,j])+v[i,j]*(u[i,j+1]-u[i,j-1]))+
(mu/h^2)*(u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]-4*u[i,j]))+dt.*(Fxn11[i,j])
vt[i,j]=v[i,j]+dt*(-(0.5/h)*(u[i,j].*(v[i+1,j]-v[i-1,j])+v[i,j]*(v[i,j+1]-v[i,j-1]))+
(mu/h^2)*(v[i+1,j]+v[i-1,j]+v[i,j+1]+v[i,j-1]-4*v[i,j]))+dt.*(Fyn11[i,j])
end
end
vt[2:nx-1,1]=(mu*dt/h^2)*(v[2:nx-1,3]-2*v[2:nx-1,2]);
vt[2:nx-1,ny-1]=(mu*dt/h^2)*(v[2:nx-1,ny-3]-2*v[2:nx-1,ny-2]);
ut[1,2:ny-1]=(mu*dt/h^2)*(u[3,2:ny-1]-2*u[2,2:ny-1]);
ut[nx-1,2:ny-1]=(mu*dt/h^2)*(u[nx-3,2:ny-1]-2*u[nx-2,2:ny-1]);
for it=1:maxit
p[2:nx-1,1]=p[2:nx-1,2]+(h/dt)*vt[2:nx-1,1];
p[2:nx-1,ny-1]=p[2:nx-1,ny-2]-(h/dt)*vt[2:nx-1,ny-1];
p[1,1:ny-1]=p[2,1:ny-1]+(h/dt)*ut[1,1:ny-1];
p[nx-1,2:ny-1]=p[nx-2,2:ny-1]-(h/dt)*ut[nx-1,2:ny-1];
for i =2:nx-1
for j=2:ny-1
p[i,j]=beta*0.25*(p[i+1,j]+p[i-1,j]+p[i,j+1]+p[i,j-1]-(0.5*h/dt)*(ut[i+1,j]-ut[i-1,j]+vt[i,j+1]-vt[i,j-1]))+(1-beta)*p[i,j];
end
end
p[floor(nx/2),floor(ny/2)]=0.0;
end
u[2:nx-1,2:ny-1]=ut[2:nx-1,2:ny-1]-(0.5*dt/h)*(p[3:nx,2:ny-1]-p[1:nx-2,2:ny-1]);
v[2:nx-1,2:ny-1]=vt[2:nx-1,2:ny-1]-(0.5*dt/h)*(p[2:nx-1,3:ny]-p[2:nx-1,1:ny-2]);
u[1:220,:] = 0
u[:,end] = 0
v[1:220,:] = 0
v[:,end]=0
u[221:222,49:57]=0.0
u[218:219,57:97]=0.0
u[221:222,97:105]=0.0
u[221:222,145:153]=0.0
u[218:219,153:193]=0.0
u[221:222,193:201]=0.0
u[221:222,241:249]=0.0
u[218:219,249:289]=0.0
u[221:222,289:297]=0.0
u[221:222,337:345]=0.0
u[218:219,345:385]=0.0
u[221:222,385:393]=0.0
v[221:222,49:57]=0.0
v[218:219,57:97]=0.0
v[221:222,97:105]=0.0
v[221:222,145:153]=0.0
v[218:219,153:193]=0.0
v[221:222,193:201]=0.0
v[221:222,241:249]=0.0
v[218:219,249:289]=0.0
v[221:222,289:297]=0.0
v[221:222,337:345]=0.0
v[218:219,345:385]=0.0
v[221:222,385:393]=0.0
w=sqrt(v.^2+u.^2)
if mod(iter,500)==0
grid_x = [i for i in x, j in y]
grid_y = [j for i in x, j in y]
fig = plt.figure(figsize = (11,5), dpi=100)
ax1 = axes()
#plt.quiver(x*1000,y*1000,v,u,w)
plt.quiver(grid_y[1:4:end,1:4:end]*1000,grid_x[1:4:end,1:4:end]*1000,v[1:4:end,1:4:end],u[1:4:end,1:4:end],w[1:4:end,1:4:end])
cbar = plt.colorbar()
title("Flow")
xlabel("X Axis")
ylabel("Y axis")
ax1[:set_xlim](-5.5, 5.5)
ax1[:set_ylim](-1.5, 5.5)
savefig("wflow1368AC29September="*string(iter)*".png")
end
end
