三维fdtd+upml的python代码.txt

free software based on toyFDTD

#!/usr/bin/env python # # # # UPML FDTD by Rob Lytle Magnetic Dipole Simulation 
v-0.2-beta # # # # This work is partially based on that of the ToyFDTD authors: 
# Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe # # 
Copyright (C) 2001 Rob Lytle, for Numerical Python Port and UPML additions # 
This program is free software you can redistribute it and/or # modify it under 
the terms of the GNU General Public License # as published by the Free Software 
Foundation either version 2 # of the License, or any later version, with the 
following conditions # attached in addition to any and all conditions of the GNU 
# General Public License: # # This program is distributed in the hope that it 
will be useful, # but WITHOUT ANY WARRANTY without even the implied warranty # 
of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. # See the GNU General 
Public License for more details. # # You should have received a copy of the GNU 
General Public License # along with this program if not, write to the Free 
Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA # 02111-1307 
USA # Boundaries: 4 layers of UPML- Unaxial Perfectly Matched Layer # Stimulus: 
a single cell by sinusoidal H field # 3D output: The electric field intensity 
vector components in the direction # of the height of the guide (Ez) are output 
to file every # PLOT_MODULUS timesteps (and for the last timestep), scaled to # 
the range of integers from zero through 254. (The colormap # included in the tar 
file assigns rgb values to the range zero through # 255.) Scaling is performed 
on each timestep individually, since # it is not known in advance what the 
maximum and minimum # values will be for the entire simulation. The integer 
value 127 is # held to be equal to the data zero for every timestep. This method 
# of autoscaling every timestep can be very helpful in a simulation # where the 
intensities are sometimes strong and sometimes faint, # since it will highlight 
the presence and structure of faint signals # when stronger signals have left 
the mesh. # Each timestep has it's own output file. This data output file format 
# can be used in several visualization tools, such as animabob and viz. # Other 
output: Notes on the progress of the simulation are written to standard # output 
as the program runs. # Some terminology used here: # # This code implements a 
Cartesian mesh with space differentials # of dx, dy, dz. # This means that a 
point in the mesh has neighboring points dx meters # away in the direction of 
the positive and negative x-axis, # and neighboring points dy meters away in the 
directions # of the +- y-axis, and neighboring points dz meters away # in the 
directions of the +- z-axis, # The mesh has nx cells in the x direction, ny in 
the y direction, # and nz in the z direction. # Ex, Ey, and Ez refer to the 
arrays of electric field intensity vectors # -- for Example, Ex is a 
3-dimensional array consisting of the # x component of the E field intensity 
vector for every point in the # mesh. Ex[i,j,k] refers to the x component of the 
E field intensity # vector at point [i,j,k]. # Hx, Hy, and Hz refer to the 
arrays of magnetic field intensity vectors. # # dt is the time differential -- 
the length of each timestep in seconds. # # bob is a file format that stands for 
"brick of bytes", meaning a string # of bytes that can be interpreted as a 
3-dimensional array of byte # values (integers from zero through 255). animabob 
is a free # visualization tool that displays and animates a sequence of bricks # 
of bytes. For more information on animabob or to download a copy, # see the 
ToyFDTD website at http:#www.borg.umn.edu/toyfdtd/ # You need Numeric Python 
available on the links at www.members.home.net/europax from Numeric import * # 
program control constants MAXIMUM_ITERATION= 1000 # total number of timesteps to 
be computed PLOT_MODULUS= 5 # The program will output 3D data every PLOT_MODULUS 
timesteps, # Except for the last iteration computed, which is always # output. 
So if MAXIMUM_ITERATION is not an integer # multiple of PLOT_MODULUS, the last 
timestep output will # come after a shorter interval than that separating # 
previous outputs. FREQUENCY= 10e9 # frequency of the stimulus in Hertz 
GUIDE_WIDTH= 0.02 # meters GUIDE_HEIGHT= 0.02 # meters GUIDE_LENGTH= 0.02 # 
meters NUMCELLS=1500.0 #number per meter Determined by experiment. 
##CELLS_PER_WAVELENGTH=25 ## ##LENGTH_IN_WAVELENGTHS=5 # pHysical constants 
LIGHT_SPEED= 299792458.0 # speed of light in a vacuum in meters/second 
LIGHT_SPEED_SQUARED= 89875517873681764.0 # m^2/s^2 MU_0= 
1.2566370614359172953850573533e-6 # permeability of free space in henry/meter 
EPSILON_0= 8.85418781762038985053656e-12 # permittivity of free space in 
farad/meter # The value used for pi is M_PI as found in /usr/include/math.h on 
SGI # IRIX 6.2. Other such constants used here are DBL_EPSILON and FLT_MAX 
FLT_MAX=3.40282347e38 DBL_EPSILON=2.2204460492503131e-16 M_PI=pi # main 
######################################/ # variable declarations #int i,j,k # 
indices of the 3D array of cells #int nx, ny, nz # total number of cells along 
the x, y, and z axes, respectively allocatedBytes = 0 # a counter to track 
number of bytes allocated iteration = 0 # counter to track how many timesteps 
have been computed stimulus = 0.0 # value of the stimulus at a given timestep 
currentSimulatedTime = 0.0 # time in simulated seconds that the simulation has 
progressed totalSimulatedTime = 0.0 # time in seconds that will be simulated by 
the program # omega # angular frequency in radians/second # lambdaa # wavelength 
of the stimulus in meters # dx, dy, dz # space differentials (or dimensions of a 
single cell) in meters # dt # time differential (how much time between 
timesteps) in seconds # Ex, Ey, Ez # arrays of Ex, Ey, and Ez values # Hx, Hy, 
Hz # arrays of Hx, Hy, and Hz values # bob output routine variables: # filename 
# filename variable for 3D bob files simulationMin = FLT_MAX # tracks minimum 
value output by the entire simulation simulationMax = -FLT_MAX # tracks maximum 
value output by the entire simulation # min, max # these track minimum and 
maximum values output in one timestep # norm # norm is set each iteration to be 
max or min, whichever is # greater in magnitude # scalingValue # multiplier used 
in output scaling, calculated every timestep #openFilePointer # handle to 3D bob 
output file ######################################/ # setting up the problem to 
be modeled # Rectangular box # This is a just a convenient boundary to form the 
UPML # # Choosing nx, ny, and nz: # There should be at least 20 cells per 
wavelength in each direction, # The number of cells along the width of the guide 
and the width of # those cells should fit the guide width Exactly, so that ny*dy 
# = GUIDE_WIDTH meters. # The same should be true for nz*dz = GUIDE_HEIGHT 
meters. # dx is chosen to be dy or dz -- whichever is smaller # nx is chosen 
similarly # # # dt is chosen for Courant stability the timestep must be kept 
small # enough so that the plane wave only travels one cell length # (one dx) in 
a single timestep. Otherwise FDTD cannot keep up # with the signal propagation, 
since FDTD computes a cell only from # it's immediate neighbors. # wavelength in 
meters: lambdaa = LIGHT_SPEED/FREQUENCY # angular frequency in radians/second: 
omega = 2.0*M_PI*FREQUENCY nx= int(NUMCELLS*GUIDE_LENGTH) ny= 
int(NUMCELLS*GUIDE_HEIGHT) nz= int(NUMCELLS*GUIDE_WIDTH) dx= GUIDE_LENGTH/nx dy= 
GUIDE_HEIGHT/ny dz= GUIDE_WIDTH/nz if dx if ( (iteration % PLOT_MODULUS) == 0): 
# bob output section # create the filename for this iteration, # which includes 
the iteration number: filename="c_%04d.bob" % iteration # 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 
Ezsize=(nx+1)*(ny+1)*(nz+1) ## 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] #Replace above code by better Numpy version 
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 +