[Tutor] numpy speed problems
Bob Gailer
bgailer at alum.rpi.edu
Fri Jun 9 20:57:32 CEST 2006
Jeff Peery wrote:
> hello, I am having some trouble with the speed of numpy. I'm crunching
> some numbers (see the attached script) and in total I have 1,000,000
> grid points over which I am integrating. I'm doing a bunch of adding,
> mulitply, divide, powers, etc, but in total there are 1,000,000 points
> to do these operations over and it really shouldn't take this long...
> as far as I know. my last simulation took about 8+ hours.
>
> What might I be doing wrong in my code to cause this to be so slow?
> big thanks!!
I don't have time to analyze all your code but I'll point out that it is
repeating some calculations many times. Example:
x_coord[ii])**2.0
y_coord[jj])**2.0
are calculated every time GetPressure is called. I suggest you calculate
these one time outside any loops and pass the result into GetPressure.
Inside GetPressure:
1j*w_mn*rho
exp(1j*k_mn*r)*dx*dy/(2.0*pi*r
are repeatedly calculated inside the for loop. Move them out of the loop.
Then look for other "loop invariants".
sin(n[mode_n]*pi*(L/2.0+x_coord[i]))*sin(m[mode_m]*pi*(W/2.0+y_coord[j]))
is another place to look.
>
> ------------------------------------------------------------------------
>
> from numpy import *
>
> """
> function definitions
> """
> def GetPressure(x_coord, y_coord, x_r, y_r, z_r, dx, dy, w_mn, rho, v_mn, k_mn, n, m):
> # intialize pressure
> p = 0.0 + 1j
>
> # sum contributions from all point sources to receiver
> for ii in range(n):
> for jj in range(m):
> r = ((x_r - x_coord[ii])**2.0 + (y_r - y_coord[jj])**2.0 + (z_r - 0.0)**2)**0.5
> p += (1j*w_mn*rho*v_mn[ii][jj])*exp(1j*k_mn*r)*dx*dy/(2.0*pi*r)
>
> p = sqrt(p*conjugate(p))
> return abs(p)
>
>
> """
> vaiables and constants
> """
>
> """problem definition parameter"""
> n = arange(1,70) #mode number in x direction
> m = arange(1,70) #mode number in y direction
> mode_n = 10 #mode number - 1
> mode_m = 10 #mode number - 1
> L = 1.2 #piston length (m)
> W = 0.6 #piston width (m)
>
> """material properties fluid"""
> c = 343.0 #speed sound in water (m/s)
> rho_a = 1.21 #density of air (kg/m^3)
>
> """piston material properties"""
> E = 7.0e10 #youngs modulus (N/m^2) (stainless steel)
> nu = 0.29 #poisson's ratio (stainless steel
> rho = 2700.0 #density of piston (stainless steel) (kg/m^3)
> t = 0.0015 #piston thickness (m)
>
> """wave speed, wave number, frequency"""
> c_l = (E/(rho*(1 - nu**2)))**0.5 #longitudinal wave speed in piston
> k_x = n*pi/W #wave number x direction
> k_y = m*pi/L #wave number y direction
> k_mn = (k_x[mode_n]**2 + k_y[mode_m]**2)**0.5 #bending wave number for n and m mode
> w_c = (c**2)/(1.8*c_l*t) #critical frequency (Hz)
> w_mn = (k_mn**2)*1.8*c_l*t/(2.0*pi)**2 #frequency for n and m (see notes 5/15/06)
> k_c = 2.0*pi*(w_c/(1.8*c_l*t))**0.5 #critical wave number
> k_a = 2.0*pi*w_mn/c #wave number in acoustic medium (air)
>
> """piston grid"""
> dx = 1.0/k_a #x direction step in space (m)
> dy = 1.0/k_a #y direction step in space (m)
>
> if dx < 0.005:
> dx = 0.005
> dy = 0.005
>
> num_y = int(L/dy) #number of nodes in y direction
> num_x = int(W/dx) #number of nodes in x direction
>
> #piston grid coordinates
> x_coord = arange(num_x, dtype=float)*W/num_x - W/2.0
> y_coord = arange(num_y, dtype=float)*L/num_y - L/2.0
>
> """field grid"""
> a = 1
> b = 50
> d = 50
> x_r = arange(a, dtype=float)*1.0/float(a) #x position of receiver (m)
> y_r = arange(b, dtype=float)*1.0/float(b) #y position of receiver (m)
> z_r = arange(d, dtype=float)*10.0/float(d) #z position of receiver (m)
>
> """acoustic variables"""
> p = 0 #amplitude of pressure at receiver (Pa)
> r = 0 #distance from origin to receiver (m)
> p_field = zeros((a,b,d), dtype=float) #pressure field (m)
>
> """calculate piston surface velocity amplitude"""
> U_mn = zeros((len(x_coord), len(y_coord)), dtype=float)
> for i in range(len(x_coord)):
> for j in range(len(y_coord)):
> #amplitude of piston surface displacement
> U_mn[i][j] = sin(n[mode_n]*pi*(L/2.0+x_coord[i]))*sin(m[mode_m]*pi*(W/2.0+y_coord[j]))
> #amplitude of piston surface velocity (m/s)
> V_mn = w_mn*U_mn
>
>
> """
> numerical integration of Raleigh's equation
> """
> for i in range(a):
> for j in range(b):
> for k in range(d):
> p_field[i][j][k] = GetPressure(x_coord, y_coord, x_r[i], y_r[j], z_r[k], dx, dy, w_mn, rho, V_mn, k_mn, num_x, num_y)
> print '%d Percent Complete'%(100.0*j/b)
>
> p_field = 20.0*log10(p_field)
>
> fileHandle = file('beam pattern.dat', "w")
> fileHandle.write('TITLE = "HW 4"\n')
> fileHandle.write('VARIABLES = "x"\n"y"\n"z"\n"pressure"\n')
> fileHandle.write('ZONE T="%s"\n' % 'pressure field')
> fileHandle.write('I=%d, J=1, ZONETYPE=Ordered\n' % (a*b*d))
> fileHandle.write('DATAPACKING=POINT DT=(DOUBLE DOUBLE DOUBLE DOUBLE)\n')
> for ii in range(a):
> for jj in range(b):
> for kk in range(d):
> fileHandle.write('%f %f %f %f\n' % (x_r[ii], y_r[jj], z_r[kk], p_field[ii][jj][kk]))
> fileHandle.close()
>
>
>
> """
> contour of surface velocity
> """
>
> fileHandle = file('mode shape.dat', "w")
> fileHandle.write('TITLE = "HW 4"\n')
> fileHandle.write('VARIABLES = "x"\n"y"\n"velocity"\n')
> fileHandle.write('ZONE T="%s"\n' % 'velocity field')
> fileHandle.write('I=%d, J=1, ZONETYPE=Ordered\n' % (num_y*num_x))
> fileHandle.write('DATAPACKING=POINT DT=(DOUBLE DOUBLE DOUBLE)\n')
> for ii in range(num_x):
> for jj in range(num_y):
> fileHandle.write('%f %f %f\n' % (x_coord[ii], y_coord[jj], V_mn[ii][jj]))
> fileHandle.close()
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Tutor maillist - Tutor at python.org
> http://mail.python.org/mailman/listinfo/tutor
>
--
Bob Gailer
510-978-4454
More information about the Tutor
mailing list