三维fdtd+upml的python代码.txt

free software based on toyFDTD

scalingValue*min),max,int(127.0 + scalingValue*max)) # scale each Ez value in 
the mesh to the range of integers # from zero through 254 and write them to the 
output # file for this iteration: for k in range(0,nz+1): for j in 
range(0,ny+1): for i in range(0,nx+1): openFilePointer.write(chr(int(127.0 + 
scalingValue*Ez[i,j,k]))) # close the output file for this iteration: 
openFilePointer.close() # end bob output section 
####################################/ # Compute the stimulus: stimulus = 
sin(omega*currentSimulatedTime) # 
stimulus=exp(-pow(currentSimulatedTime/1e-11,2.0)/2.0) Hz[15,15,15]=stimulus # 
stimulate H in z direction, a mag dipole # Enforce PEC boundary Note: may or may 
not be needed ## Ey[0,:,:]=Ez[0,:,:]=0 ## Ex[:,0,:]=Ez[:,0,:]=0 ## 
Ex[:,:,0]=Ey[:,:,0]=0 ## Ey[nx,:,:]=Ez[nx,:,:]=0 ## Ex[:,ny,:]=Ez[:,ny,:]=0 ## 
Ex[:,:,nz]=Ey[:,:,nz]=0 ####################################/ # Update the 
interior of the mesh: # all vector H vector components # 
dstore[0:nx,0:ny,0:nz]=Bx[0:nx,0:ny,0:nz] Bx[0:nx,0:ny,0:nz] = 
C1z[0:nx,0:ny,0:nz]* Bx[0:nx,0:ny,0:nz] + C2z[0:nx,0:ny,0:nz] * ( 
(Ey[0:nx,0:ny,1:(nz+1)]-Ey[0:nx,0:ny,0:nz] ) / dz - 
(Ez[0:nx,1:(ny+1),0:nz]-Ez[0:nx,0:ny,0:nz]) / dy ) Hx[0:nx,0:ny,0:nz]= 
C3y[0:nx,0:ny,0:nz] * Hx[0:nx,0:ny,0:nz] + C4hy[0:nx,0:ny,0:nz] * ( 
C5x[0:nx,0:ny,0:nz] * Bx[0:nx,0:ny,0:nz] - C6x[0:nx,0:ny,0:nz] * 
dstore[0:nx,0:ny,0:nz]) dstore[0:nx,0:ny,0:nz]=By[0:nx,0:ny,0:nz] 
By[0:nx,0:ny,0:nz] = C1x[0:nx,0:ny,0:nz]* By[0:nx,0:ny,0:nz] + 
C2x[0:nx,0:ny,0:nz] * ( (Ez[1:(nx+1),0:ny,0:nz]-Ez[0:nx,0:ny,0:nz]) / dx - 
(Ex[0:nx,0:ny,1:(nz+1)]-Ex[0:nx,0:ny,0:nz]) / dz ) Hy[0:nx,0:ny,0:nz]= 
C3z[0:nx,0:ny,0:nz] * Hy[0:nx,0:ny,0:nz] + C4hz[0:nx,0:ny,0:nz] * ( 
C5y[0:nx,0:ny,0:nz] * By[0:nx,0:ny,0:nz] - C6y[0:nx,0:ny,0:nz] * 
dstore[0:nx,0:ny,0:nz]) dstore[0:nx,0:ny,0:nz]=Bz[0:nx,0:ny,0:nz] 
Bz[0:nx,0:ny,0:nz] = C1y[0:nx,0:ny,0:nz]* Bz[0:nx,0:ny,0:nz] + 
C2y[0:nx,0:ny,0:nz] * ( (Ex[0:nx,1:(ny+1),0:nz]-Ex[0:nx,0:ny,0:nz] ) / dy - 
(Ey[1:(nx+1),0:ny,0:nz]-Ey[0:nx,0:ny,0:nz] ) / dx ) Hz[0:nx,0:ny,0:nz]= 
C3x[0:nx,0:ny,0:nz] * Hz[0:nx,0:ny,0:nz] + C4hx[0:nx,0:ny,0:nz] * ( 
C5z[0:nx,0:ny,0:nz] * Bz[0:nx,0:ny,0:nz] - C6z[0:nx,0:ny,0:nz] * 
dstore[0:nx,0:ny,0:nz]) # Update the E field vector components. # The values on 
the faces of the mesh are not updated here thEy # are handled by the boundary 
condition computation # (and stimulus computation). 
dstore[1:nx,1:ny,1:nz]=Dx[1:nx,1:ny,1:nz] Dx[1:nx,1:ny,1:nz] = 
C1z[1:nx,1:ny,1:nz]* Dx[1:nx,1:ny,1:nz] + C2z[1:nx,1:ny,1:nz] * ( 
(Hz[1:nx,1:ny,1:nz]-Hz[1:nx,0:(ny-1),1:nz] ) / dy - 
(Hy[1:nx,1:ny,1:nz]-Hy[1:nx,1:ny,0:(nz-1)]) / dz ) Ex[1:nx,1:ny,1:nz]= 
C3y[1:nx,1:ny,1:nz] * Ex[1:nx,1:ny,1:nz] + C4ey[1:nx,1:ny,1:nz] * ( 
C5x[1:nx,1:ny,1:nz] * Dx[1:nx,1:ny,1:nz] - C6x[1:nx,1:ny,1:nz] * 
dstore[1:nx,1:ny,1:nz]) dstore[1:nx,1:ny,1:nz]=Dy[1:nx,1:ny,1:nz] 
Dy[1:nx,1:ny,1:nz] = C1x[1:nx,1:ny,1:nz]* Dy[1:nx,1:ny,1:nz] + 
C2x[1:nx,1:ny,1:nz] * ( (Hx[1:nx,1:ny,1:nz]-Hx[1:nx,1:ny,0:(nz-1)] ) / dz - 
(Hz[1:nx,1:ny,1:nz]-Hz[0:(nx-1),1:ny,1:nz]) / dx) Ey[1:nx,1:ny,1:nz]= 
C3z[1:nx,1:ny,1:nz] * Ey[1:nx,1:ny,1:nz] + C4ez[1:nx,1:ny,1:nz] * ( 
C5y[1:nx,1:ny,1:nz] * Dy[1:nx,1:ny,1:nz] - C6y[1:nx,1:ny,1:nz] * 
dstore[1:nx,1:ny,1:nz]) dstore[1:nx,1:ny,1:nz]=Dz[1:nx,1:ny,1:nz] 
Dz[1:nx,1:ny,1:nz] = C1y[1:nx,1:ny,1:nz]* Dz[1:nx,1:ny,1:nz] + 
C2y[1:nx,1:ny,1:nz] * ( (Hy[1:nx,1:ny,1:nz]-Hy[0:(nx-1),1:ny,1:nz]) / dx - 
(Hx[1:nx,1:ny,1:nz]-Hx[1:nx,0:(ny-1),1:nz]) / dy ) Ez[1:nx,1:ny,1:nz]= 
C3x[1:nx,1:ny,1:nz] * Ez[1:nx,1:ny,1:nz] + C4ex[1:nx,1:ny,1:nz] * ( 
C5z[1:nx,1:ny,1:nz] * Dz[1:nx,1:ny,1:nz] - C6z[1:nx,1:ny,1:nz] * 
dstore[1:nx,1:ny,1:nz]) #correct E at boundary=0 
dstore[0,1:ny,1:nz]=Dx[0,1:ny,1:nz] Dx[0,1:ny,1:nz] = C1z[0,1:ny,1:nz]* 
Dx[0,1:ny,1:nz] + C2z[0,1:ny,1:nz] * ( (Hz[0,1:ny,1:nz]-Hz[0,0:(ny-1),1:nz] ) / 
dy - (Hy[0,1:ny,1:nz]-Hy[0,1:ny,0:(nz-1)]) / dz ) Ex[0,1:ny,1:nz]= 
C3y[0,1:ny,1:nz] * Ex[0,1:ny,1:nz] + C4ey[0,1:ny,1:nz] * ( C5x[0,1:ny,1:nz] * 
Dx[0,1:ny,1:nz] - C6x[0,1:ny,1:nz] * dstore[0,1:ny,1:nz]) 
dstore[1:nx,0,1:nz]=Dy[1:nx,0,1:nz] Dy[1:nx,0,1:nz] = C1x[1:nx,0,1:nz]* 
Dy[1:nx,0,1:nz] + C2x[1:nx,0,1:nz] * ( (Hx[1:nx,0,1:nz]-Hx[1:nx,0,0:(nz-1)] ) / 
dz - (Hz[1:nx,0,1:nz]-Hz[0:(nx-1),0,1:nz]) / dx) Ey[1:nx,0,1:nz]= 
C3z[1:nx,0,1:nz] * Ey[1:nx,0,1:nz] + C4ez[1:nx,0,1:nz] * ( C5y[1:nx,0,1:nz] * 
Dy[1:nx,0,1:nz] - C6y[1:nx,0,1:nz] * dstore[1:nx,0,1:nz]) 
dstore[1:nx,1:ny,0]=Dz[1:nx,1:ny,0] Dz[1:nx,1:ny,0] = C1y[1:nx,1:ny,0]* 
Dz[1:nx,1:ny,0] + C2y[1:nx,1:ny,0] * ( (Hy[1:nx,1:ny,0]-Hy[0:(nx-1),1:ny,0]) / 
dx - (Hx[1:nx,1:ny,0]-Hx[1:nx,0:(ny-1),0]) / dy ) Ez[1:nx,1:ny,0]= 
C3x[1:nx,1:ny,0] * Ez[1:nx,1:ny,0] + C4ex[1:nx,1:ny,0] * ( C5z[1:nx,1:ny,0] * 
Dz[1:nx,1:ny,0] - C6z[1:nx,1:ny,0] * dstore[1:nx,1:ny,0]) # end mainloop 
##################################/ # Output section: # The output routine is 
repeated one last time to write out # the last data computed. # time in 
simulated seconds that the simulation has progressed: currentSimulatedTime = 
dt*iteration # print to standard output the iteration number # and current 
simulated time: print "#%d %E sec" %(iteration,currentSimulatedTime) # 3D data 
output for the last timestep: # create the filename for this iteration, # which 
includes the iteration number: filename= "c_"+`iteration`+".bob" # open a new 
data file for this iteration: openFilePointer = open(filename, "wb") # find the 
max and min values to be output this timestep: min = FLT_MAX max = -FLT_MAX 
##for i in range(0,nx+1): ## ## for j in range(0,ny+1): ## ## for k in 
range(0,nz): ## ## if (Ez[i,j,k] < min): min = Ez[i,j,k] ## if (Ez[i,j,k] > 
max): max = Ez[i,j,k] min=sort(resize(Ez,(Ezsize,)))[0] 
max=sort(resize(Ez,(Ezsize,)))[-1] # update the tracking variables for minimum 
and maximum # values for the entire simulation: if (min < simulationMin): 
simulationMin = min if (max > simulationMax): simulationMax = max # set norm to 
be max or min, whichever is greater in magnitude: if (abs(max) > abs(min)): 
norm= abs(max) else: norm=abs(min) if (norm == 0.0): # if everything is zero, 
give norm a tiny value # to avoid division by zero: norm = DBL_EPSILON 
scalingValue = 127.0/norm # write to standard output the minimum and maximum 
values # from this iteration and the minimum and maximum values # that will be 
written to the bob file this iteration: print "%E %d < Ez BoB < %E %d " 
%(min,int(127.0 + scalingValue*min),max,int(127.0 + scalingValue*max)) # scale 
each Ez value in the mesh to the range of integers # from zero through 254 and 
write them to the output # file for this iteration: for k in range(0,nz+1): for 
j in range(0,ny+1): for i in range(0,nx+1): openFilePointer.write(chr(int(127.0 
+ scalingValue*Ez[i,j,k]))) # close the output file for this iteration: 
openFilePointer.close() ######################################/ # write some 
progress notes to standard output: print "\n" print "\n" # print out how much 
memory has been allocated: print "Dynamically allocated %d bytes\n" 
%allocatedBytes print "\n" # print out some simulation parameters: print 
"PLOT_MODULUS = %d\n" %PLOT_MODULUS print "rectangular waveguide\n" print 
"Stimulus = %E Hertz " %FREQUENCY print "continuous plane wave\n" print "\n" 
print "Meshing parameters:\n" print "%dx%dx%d cells\n" %(nx,ny,nz) print "dx=%f, 
dy=%f, dz=%f meters\n" %(dx,dy,dz) print "%f x %f x %f meter cubed simulation 
region\n" %(GUIDE_LENGTH,GUIDE_HEIGHT, GUIDE_WIDTH) print "\n" print "Time 
simulated was %E seconds %d timesteps\n" %(totalSimulatedTime,MAXIMUM_ITERATION) 
print "\n" print "3D output scaling parameters:\n" print "Autoscaling every 
timestep\n" print "\n" # print out simulationMin and simulationMax: print 
"Minimum output value was %E \n" %simulationMin print "Maximum output value was 
%E \n" %simulationMax print "\n" print "\n" # print out a bob command line, 
including the dimensions # of the output files: print "bob -cmap chengGbry.cmap 
-s %dx%dx%d *.bob\n" %(nx+1,ny+1,nz+1) print "\n" # print out a viz command 
line: print "viz ToyFDTD1c.viz\n" print "\n" print "\n" print "\n" print "\n" # 
end main