What might I be doing wrong in my code to cause this to be so slow? big thanks!!
Jeff
__________________________________________________
Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around
http://mail.yahoo.com
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@python.org http://mail.python.org/mailman/listinfo/tutor
