I don't know anything about plasma physics, but this entire code is in global scope and uses a bunch of non-const global arrays. This is the very first thing the performance tips warn against:
http://julia.readthedocs.org/en/latest/manual/performance-tips/ Refactor your code into functions that take arguments and return values instead of operating on globals. On Mon, Sep 28, 2015 at 10:54 PM, John Gibson <johnfgib...@gmail.com> wrote: > Marius: I think you'd be better off in learning to write fast Julia code > if you presented an straightforward implementation of a classic PDE problem > and asked for help in optimizing that. Maybe 2D heat equation would be a > good place to start. The code you've presented is not really intelligible > without a better description of the equations, the boundary conditions, the > discretization methods, more detailed comments, explanations of magic > numbers, etc. And requiring plasma physics expertise makes your audience > here too small, probably. > > best regards, > > John > > > > On Monday, September 28, 2015 at 10:14:32 PM UTC-4, Marius wrote: >> >> 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 >> >>
