4

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?

Riccardo De Nigris
  • 41
  • 1
  • 1
  • 3
  • 1
    Welcome to Stack overflow! Your question is pretty broad as it stands and it's not easy to understand where the specific problem is. You've also asked [an identical question](http://http://stackoverflow.com/questions/30547023), asking the same question twice is usually a bad idea. You've got some good suggestions in the comments on that question and it's unlikely that you'll get better help unless you show in more detail your exact implementation and where you think it is too slow as a [minimal, complete and verifiable example](http://stackoverflow.com/help/mcve) – J Richard Snape May 30 '15 at 23:27
  • Sorry, but I'm new and I don't known the guide lines of this forum.I just wanted to know if there is a faster method than that of finite differences on the basis of the idea outlined above. I apologize if I've created two identical questions but I could not understand the suspension of the first application. Now what should I do with this question? – Riccardo De Nigris May 31 '15 at 08:19
  • No problem, thanks for responding. I think you could edit this into a better question for help. Include the section of code that actually performs the finite difference, the number of points you calculate at (i.e. your mesh size) and how fast it runs vs how fast you think it could / would like it to – J Richard Snape May 31 '15 at 08:31
  • Then, open another question or place a comment on this? – Riccardo De Nigris Jun 01 '15 at 08:16
  • Edit this one (you can still edit when it's on hold). Here's a suggestion. Show the exact code where the looping and finite difference logic is executed. Your version, not simply a link, put the code in the question. Lookup up python timeit module and post the results for your situation. Write why you think a matrix based solution might be quicker and where you're stuck (i.e. Why you can't just go and try it). Then it may be eligible to be reopened – J Richard Snape Jun 01 '15 at 08:28
  • now is it possible to reopen the question? – Riccardo De Nigris Jun 01 '15 at 09:28
  • It looks likely. The way the site works, members vote on reopening (look it up on help / meta parts of site for details) this Q currently has 4 reopen votes, so it's likely your edits have done enough. :) – J Richard Snape Jun 01 '15 at 10:16

3 Answers3

2

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)
jared
  • 4,165
  • 1
  • 8
  • 31
pyano
  • 1,885
  • 10
  • 28
  • I assume your naming convention of `ie`, `iw`, `jn`, and `js` refers to the cardinal directions. But `i` indexes the rows, so it would adjust "looking" up and down, which would be north and south, while `j` is "looking" left and right, which would be east and west, right? – jared Dec 29 '22 at 18:21
  • Almost :-) . `i` is linked to the west-east direction, `iw`= west of i =`i-1`. And `j` is linked to the south - north direction. This is called compass notation. In a geo/graphical view, west-east direction is left-right. BUT even then I use the mathematical matrix orientation. Try `x, y = np.meshgrid(np.arange(5), np.arange(10), indexing='ij')` to see what I mean: `i/x` is downward and `j/y` to the right now. Due to my experience this is on the long run less confusing than the geo/graphical matrix orientation with `x, y = np.meshgrid(np.arange(5), np.arange(10), indexing='xy')` . – pyano Jan 02 '23 at 19:39
  • In your solution, isn't `i` controlling the up-down traversal and `j` controlling the left-right traversal? That would imply to me that `i` is for north-south and `j` is for east-west. – jared Jan 02 '23 at 19:50
  • No, this time your interpretation does not coincidence with my code. There is an unpleasant difference between mathemacial matrix orientation and geo/graphical matrix orientation. You have to find your own way to overcome this. For that I recomend to try the meshgrid example and to print the matrixes x and y and U and to plot some functions on (x,y), `U = x*x + y*y` e.g. using plt.contour and plt.pcolormesh. – pyano Jan 02 '23 at 20:15
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.

Intrepid Traveller
  • 95
  • 6
  • the code works fine, I must make it faster and I think to use a function DeltaU=f(u) to decrease run time. Also the code should run on slow machines therefore I would like to optimize it. – Riccardo De Nigris May 30 '15 at 18:53
0

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.

Ozgur Taner
  • 1
  • 1