fast method with numpy for 2D Heat equation

Question

I'm looking for a method for solve the 2D heat equation with python. I have already implemented the finite difference method but is slow motion (to make 100,000 simulations takes 30 minutes). The idea is to create a code in which the end can write,

for t in TIME:  
    DeltaU=f(U)
    U=U+DeltaU*DeltaT
    save(U)

How can I do that?

In the first form of my code, I used the 2D method of finite difference, my grill is 5000x250 (x, y). Now I would like to decrease the speed of computing and the idea is to find

DeltaU = f(u)

where U is a heat function. For implementation I used this source http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/ for 2D case, but the run time is more expensive for my necessity. Are there some methods to do this?

Maybe I must to work with the matrix

A=1/dx^2 (2 -1  0  0 ... 0 
         -1  2 -1  0 ... 0
         0  -1  2 -1 ... 0
         .            .
         .            .
         .            .
         0 ...        -1 2)

but how to make this in the 2D problem? How to inserting Boundary conditions in A? This is the code for the finite difference that I used:

Lx=5000   # physical length x vector in micron
Ly=250   # physical length y vector in micron
Nx = 100     # number of point of mesh along x direction
Ny = 50     # number of point of mesh along y direction
a = 0.001 # diffusion coefficent
dx = 1/Nx
dy = 1/Ny
dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2)) # it is 0.04
x = linspace(0.1,Lx, Nx)[np.newaxis] # vector to create mesh
y = linspace(0.1,Ly, Ny)[np.newaxis] # vector to create mesh
I=sqrt(x*y.T) #initial data for heat equation
u=np.ones(([Nx,Ny])) # u is the matrix referred to heat function
steps=100000
for m in range (0,steps):
        du=np.zeros(([Nx,Ny]))

        for i in range (1,Nx-1):

            for j in range(1,Ny-1):

               dux = ( u[i+1,j] - 2*u[i,j] + u[i-1, j] ) / dx**2
               duy = ( u[i,j+1] - 2*u[i,j] + u[i, j-1] ) / dy**2            
               du[i,j] = dt*a*(dux+duy)




    # Boundary Conditions
    t1=(u[:,0]+u[:,1])/2
    u[:,0]=t1
    u[:,1]=t1
    t2=(u[0,:]+u[1,:])/2
    u[0,:]=t2
    u[1,:]=t2
    t3=(u[-1,:]+u[-2,:])/2
    u[-1,:]=t3
    u[-2,:]=t3
    u[:,-1]=1


    filename1='data_{:08d}.txt'

    if m%100==0:
        np.savetxt(filename1.format(m),u,delimiter='\t' )

For elaborate 100000 steps the run time is about 30 minutes. I would to optimize this code (with the idea presented in the initial lines) to have a run time about 5/10 minutes or minus. How can I do it?

Answer 1

Have you considered paralellizing your code or using GPU acceleration.

It would help if you ran your code the python profiler (cProfile) so that you can figure out where you bottleneck in runtime is. I'm assuming it's in solving the matrix equation you get to which can be easily sped up by the methods I listed above.