CODEWITHSUNDEEP
pythonpython-3.xnumpyheat
Riccardo De Nigris
Riccardo De Nigris

Reputation: 41

fast method with numpy for 2D Heat equation

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?

Upvotes: 4

Views: 11430

Answers (3)

pyano
pyano

Reputation: 1988

There are some simple but tremendous improvements possible.

Just by introducing Dxx, Dyy = 1/(dx*dx), 1/(dy*dy) the runtime drops 25%. By using slices and avoiding for-loops, the code is now 400 times faster.

import numpy as np

def f_for(u):
    for m in range(0, steps):
        du = np.zeros_like(u)
        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)
    
    return du


def f_slice(u):
    du = np.zeros_like(u)
    Dxx, Dyy = 1/dx**2, 1/dy**2
    i = slice(1, Nx-1)
    iw = slice(0, Nx-2)
    ie = slice(2, Nx)
    j = slice(1, Ny-1)
    js = slice(0, Ny-2)
    jn = slice(2, Ny)
    for m in range(0, steps):
        dux = Dxx * (u[ie, j] - 2*u[i, j] + u[iw, j])
        duy = Dyy * (u[i, jn] - 2*u[i, j] + u[i, js])
        du[i, j] = dt*a*(dux + duy)
    
    return du


Nx = 100        # number of mesh points in the x-direction
Ny = 50         # number of mesh points in the y-direction
a = 0.001       # diffusion coefficent
dx = 1/Nx
dy = 1/Ny
dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2))
steps = 10000
U = np.ones((Nx, Ny))

%timeit f_for(U)
%timeit f_slice(U)

Upvotes: 2

Ozgur Taner
Ozgur Taner

Reputation: 1

I might be wrong but in your code in the loop you create for the time steps, "m in range(steps)" in the one line below you continue with; Du =np.zeros(----). Not an expert of python but this may be resulting in creating a sparse matrix in the number of steps, in this case 100k times.

Upvotes: 0

Intrepid Traveller
Intrepid Traveller

Reputation: 95

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.

Upvotes: 1

Related Questions