toyfdtd1 c source code listing.htm

toyFDTD GNU

#define MAXIMUM_ITERATION 1000
          // total number of timesteps to be computed
#define 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.  
#define FREQUENCY 10.0e9     
          // frequency of the stimulus in Hertz
#define GUIDE_WIDTH 0.0229
          // meters
#define GUIDE_HEIGHT 0.0102
          // meters
#define LENGTH_IN_WAVELENGTHS 5.0
          // length of the waveguide in wavelengths of the stimulus wave
#define CELLS_PER_WAVELENGTH 25.0
          // minimum number of grid cells per wavelength in the x, y, and
          //     z directions

// physical constants
#define LIGHT_SPEED 299792458.0       
          // speed of light in a vacuum in meters/second
#define LIGHT_SPEED_SQUARED 89875517873681764.0        
          // m^2/s^2
#define MU_0 1.25663706143591729538505735331180115367886775975
      00423283899778369231265625144835994512139301368468271e-6
          // permeability of free space in henry/meter
#define EPSILON_0 8.854187817620389850536563031710750260608370
      1665994498081024171524053950954599821142852891607182008932e-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

main()  
{// 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
  int allocatedBytes = 0;          
           // a counter to track number of bytes allocated

  int iteration = 0;          
           // counter to track how many timesteps have been computed
  double stimulus = 0.0;       
           // value of the stimulus at a given timestep
  double currentSimulatedTime = 0.0;
           // time in simulated seconds that the simulation has progressed
  double totalSimulatedTime = 0.0;       
           // time in seconds that will be simulated by the program
  double omega;       
           // angular frequency in radians/second
  double lambda;       
           // wavelength of the stimulus in meters
  double dx, dy, dz; 
           // space differentials (or dimensions of a single cell) in meters
  double dt;
           // time differential (how much time between timesteps) in seconds
  double dtmudx, dtepsdx;       
           // physical constants used in the field update equations 
  double dtmudy, dtepsdy;       
           // physical constants used in the field update equations 
  double dtmudz, dtepsdz;       
           // physical constants used in the field update equations 

  double ***ex, ***ey, ***ez;  
           // pointers to the arrays of ex, ey, and ez values
  double ***hx, ***hy, ***hz;
           // pointers to the arrays of hx, hy, and hz values

  // bob output routine variables:
  char filename[1024];
           // filename variable for 3D bob files
  double simulationMin = FLT_MAX;
           // tracks minimum value output by the entire simulation
  double simulationMax = -FLT_MAX;
           // tracks maximum value output by the entire simulation
  double min, max;
           // these track minimum and maximum values output in one timestep
  double norm;
           // norm is set each iteration to be max or min, whichever is 
           //     greater in magnitude
  double scalingValue;
           // multiplier used in output scaling, calculated every timestep
  FILE *openFilePointer;
           // pointer to 3D bob output file
  FILE *vizFilePointer;
           // pointer to viz file


  /////////////////////////////////////////////////////////////////////////////
  // setting up the problem to be modeled
  //
  // David K. Cheng, Field and Wave Electromagnetics, 2nd ed., 
  //     pages 554-555. 
  // Rectangular waveguide, interior width = 2.29cm, interior height = 1.02cm.
  // This is a WG-16 waveguide useful for X-band applications.
  //
  // Choosing nx, ny, and nz:
  // There should be at least 20 cells per wavelength in each direction, 
  //     but we'll go with 25 so the animation will look prettier.   
  //     (CELLS_PER_WAVELENGTH was set to 25.0 in the global 
  //     constants at the beginning of the code.)
  // 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 to make the guide LENGTH_IN_WAVELENGTHS 
  //     wavelengths long.  
  // 
  // 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:
  lambda = LIGHT_SPEED/FREQUENCY; 
  // angular frequency in radians/second:
  omega = 2.0*M_PI*FREQUENCY; 

  // set ny and dy:
  // start with a small ny:
  ny = 3;  
  // calculate dy from the guide width and ny:
  dy = GUIDE_WIDTH/ny;
  // until dy is less than a twenty-fifth of a wavelength,
  //     increment ny and recalculate dy:
  while(dy >= lambda/CELLS_PER_WAVELENGTH)
    {
      ny++;
      dy = GUIDE_WIDTH/ny;
    }

  // set nz and dz:
  // start with a small nz:
  nz = 3;  
  // calculate dz from the guide height and nz:
  dz = GUIDE_HEIGHT/nz;
  // until dz is less than a twenty-fifth of a wavelength,
  //     increment nz and recalculate dz:
  while(dz >= lambda/CELLS_PER_WAVELENGTH)
    {
      nz++;
      dz = GUIDE_HEIGHT/nz;
    }

  // set dx, nx, and dt:
  // set dx equal to dy or dz, whichever is smaller:
  dx = (dy < dz) ? dy : dz;
  // choose nx to make the guide LENGTH_IN_WAVELENGTHS 
  //     wavelengths long:  
  nx = (int)(LENGTH_IN_WAVELENGTHS*lambda/dx);
  // chose dt for Courant stability:
  dt = 1.0/(LIGHT_SPEED*sqrt(1.0/(dx*dx) + 1.0/(dy*dy) + 1.0/(dz*dz)));
  // time in seconds that will be simulated by the program:
  totalSimulatedTime = MAXIMUM_ITERATION*dt;

  // constants used in the field update equations:
  dtmudx = dt/(MU_0*dx);
  dtepsdx = dt/(EPSILON_0*dx);
  dtmudy = dt/(MU_0*dy);
  dtepsdy = dt/(EPSILON_0*dy);
  dtmudz = dt/(MU_0*dz);
  dtepsdz = dt/(EPSILON_0*dz);
  

  /////////////////////////////////////////////////////////////////////////////
  // memory allocation for the FDTD mesh:
  // There is a separate array for each of the six vector components, 
  //      ex, ey, ez, hx, hy, and hz.
  // The mesh is set up so that tangential E vectors form the outer faces of 
  //     the simulation volume.  There are nx*ny*nz cells in the mesh, but 
  //     there are nx*(ny+1)*(nz+1) ex component vectors in the mesh.
  //     There are (nx+1)*ny*(nz+1) ey component vectors in the mesh.
  //     There are (nx+1)*(ny+1)*nz ez component vectors in the mesh.
  // If you draw out a 2-dimensional slice of the mesh, you'll see why
  //     this is.  For example if you have a 3x3x3 cell mesh, and you 
  //     draw the E field components on the z=0 face, you'll see that
  //     the face has 12 ex component vectors, 3 in the x-direction