[julia-users] Re: Plasma Actuator Simulation
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,:]=1.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,:]=1.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)
>
[julia-users] Re: Plasma Actuator Simulation
In my experience it is helpful to have a Fortran version to compare speeds
against when learning Julia.
This would be easy w.r.t. your problem.
It is easy to right functional Julia code that is very slow, but looks
okay. The best way for me was to have
have a basis of comparison.
On Monday, September 28, 2015 at 9:14:32 PM UTC-5, 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,:]=1.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,:]=1.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[
[julia-users] Re: Plasma Actuator Simulation
Ok...wow...I am not surprised it is slow. Well here is what you could try.
1. Wrap everything into a function. Julia has difficulty with globals.
2. You are defining pointless functions in a loop. You can easily change
the code to exclude them
3. In your loops you are accessing arrays in row major. This will cause
alot of cache misses in a column major language like Julia
4. I don't think you need to code up your own version of the gauss seidel
method. There are plenty of built in ways to solve a system of linear
equations in Julia.
Fix those things and it should run faster
-Luke
On Monday, September 28, 2015 at 7:14:32 PM UTC-7, 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,:]=1.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,:]=1.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
Re: [julia-users] Re: Plasma Actuator Simulation
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 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,:]=1.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[
