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
>
>
