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