toyboxfdtdbezhig.c

This simple simulation of a pulse traveling down a parallel-plate guide makes a handy test code for

  //     increment ny and recalculate dy:
  while(dy >= lambda/25.0)
    {
      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/25.0)
    {
      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/(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);
  

#ifdef NO_AUTOSCALING
  /////////////////////////////////////////////////////////////////////////////
  // set the output scaling value for the main output files:
  //   -- This value is the same for every timestep when no autoscaling is 
  //      used, and can be set outside the main loop.  If autoscaling is used,
  //      scalingValue must be set within the main loop every output timestep.
  min = SCALING_MIN;
  max = SCALING_MAX;
  norm = (fabs(max) > fabs(min)) ? 
    fabs(max) : fabs(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;
#endif

  /////////////////////////////////////////////////////////////////////////////
  // 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
  //     and 4 in the y-direction.  That face also has 12 ey components,
  //     4 in the x-direction and 3 in the y-direction.  

  // Allocate memory for the E field arrays:

  // allocate the array of ex components:
  ex = (double ***)malloc((nx)*sizeof(double **));
  for(i=0; i<(nx); i++)
    {  
      ex[i] = (double **)malloc((ny+1)*sizeof(double *));
      for(j=0; j<(ny+1); j++)
	{
	  ex[i][j] = (double *)malloc((nz+1)*sizeof(double));
	  for(k=0; k<(nz+1); k++)
	    {
	      ex[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx)*(ny+1)*(nz+1) * sizeof(double));
  
  // allocate the array of ey components:
  ey = (double ***)malloc((nx+1)*sizeof(double **));
  for(i=0; i<(nx+1); i++)
    {  
      ey[i] = (double **)malloc((ny)*sizeof(double *));
      for(j=0; j<(ny); j++)
	{
	  ey[i][j] = (double *)malloc((nz+1)*sizeof(double));
	  for(k=0; k<(nz+1); k++)
	    {
	      ey[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx+1)*(ny)*(nz+1) * sizeof(double));

  // allocate the array of ez components:
  ez = (double ***)malloc((nx+1)*sizeof(double **));
  for(i=0; i<(nx+1); i++)
    {  
      ez[i] = (double **)malloc((ny+1)*sizeof(double *));
      for(j=0; j<(ny+1); j++)
	{
	  ez[i][j] = (double *)malloc((nz)*sizeof(double));
	  for(k=0; k<(nz); k++)
	    {
	      ez[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx+1)*(ny+1)*(nz) * sizeof(double));

  // Allocate the H field arrays:
  // Since the H arrays are staggered half a step off 
  //     from the E arrays in every direction, the H 
  //     arrays are one cell smaller in the x, y, and z 
  //     directions than the corresponding E arrays. 
  // By this arrangement, the outer faces of the mesh
  //     consist of E components only, and the H 
  //     components lie only in the interior of the mesh.  

  // allocate the array of hx components:
  hx = (double ***)malloc((nx-1)*sizeof(double **));
  for(i=0; i<(nx-1); i++)
    {  
      hx[i] = (double **)malloc((ny)*sizeof(double *));
      for(j=0; j<(ny); j++)
	{
	  hx[i][j] = (double *)malloc((nz)*sizeof(double));
	  for(k=0; k<(nz); k++)
	    {
	      hx[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx-1)*(ny)*(nz) * sizeof(double));
  
  // allocate the array of hy components:
  hy = (double ***)malloc((nx)*sizeof(double **));
  for(i=0; i<(nx); i++)
    {  
      hy[i] = (double **)malloc((ny-1)*sizeof(double *));
      for(j=0; j<(ny-1); j++)
	{
	  hy[i][j] = (double *)malloc((nz)*sizeof(double));
	  for(k=0; k<(nz); k++)
	    {
	      hy[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx)*(ny-1)*(nz) * sizeof(double));
  
  // allocate the array of hz components:
  hz = (double ***)malloc((nx)*sizeof(double **));
  for(i=0; i<(nx); i++)
    {  
      hz[i] = (double **)malloc((ny)*sizeof(double *));
      for(j=0; j<(ny); j++)
	{
	  hz[i][j] = (double *)malloc((nz-1)*sizeof(double));
	  for(k=0; k<(nz-1); k++)
	    {
	      hz[i][j][k] = 0.0;
	    }
	}
    }
  allocatedBytes += ( (nx)*(ny)*(nz-1) * sizeof(double));

  /////////////////////////////////////////////////////////////////////////////
  // write some progress notes to standard output

  // print out some identifying information 
  fprintf(stdout, "\n\nToyBoxFDTDbezhig version 1.0\n");
  fprintf(stdout, "Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, ");
  fprintf(stdout, "Matthew O'Keefe\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "ToyBoxFDTDbezhig is free software published under the\n"); 
  fprintf(stdout, "terms of the GNU General Public License as published by\n"); 
  fprintf(stdout, "the Free Software Foundation.\n");  
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  // print out a bob command line, 
  //     including the dimensions of the output files:
  fprintf(stdout, "bob -cmap chengGbry2.cmap -s %dx%dx%d *.bob\n", 
	  nx+1, ny+1, nz);
  fprintf(stdout, "\n");
  // print out a viz command line:
  fprintf(stdout, "viz ToyBoxFDTDbezhig.viz\n");
  fprintf(stdout, "\n");
  // print out how much memory has been allocated: 
  fprintf(stdout, "Dynamically allocated %d bytes\n", allocatedBytes);
  fprintf(stdout, "\n");
  // print out some simulation parameters:
  fprintf(stdout, "PLOT_MODULO = %d\n", PLOT_MODULO);
#ifdef PMC
  fprintf(stdout, "parallel plate waveguide\n");
#else
  fprintf(stdout, "rectangular waveguide\n");
#endif
  fprintf(stdout, "Stimulus = %lg Hz ", FREQUENCY);
#ifdef SINUSOIDAL_PULSE_STIMULUS
  fprintf(stdout, "sinusoidal pulse\n");
#else
  fprintf(stdout, "continuous plane wave\n");
#endif
  fprintf(stdout, "\n");
  fprintf(stdout, "Meshing parameters:\n");
  fprintf(stdout, "%dx%dx%d cells\n", nx, ny, nz);
  fprintf(stdout, "dx=%lg, dy=%lg, dz=%lg meters\n", dx, dy, dz);
  fprintf(stdout, "%lg x %lg x %lg meter^3 simulation region\n", 
	  GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda);
  fprintf(stdout, "\n");
  fprintf(stdout, "Time simulated will be %lg seconds, %d timesteps\n", 
	  totalSimulatedTime, MAXIMUM_ITERATION);
  fprintf(stdout, "\n");
  fprintf(stdout, "3D output scaling parameters:\n");
#ifndef NO_AUTOSCALING
  fprintf(stdout, "Autoscaling every timestep\n");
#endif
#ifdef NO_AUTOSCALING
  fprintf(stdout, "Global scaling over all timesteps\n");
  fprintf(stdout, "SCALING_MAX = %lg\n", SCALING_MAX);
  fprintf(stdout, "SCALING_MIN = %lg\n", SCALING_MIN);
#ifdef IDIOTPROOF_SCALING
  fprintf(stdout, "IDIOTPROOF_SCALING: test for overflow and underflow:\n");
  fprintf(stdout, "OVERFLOW_DEFAULT = %d\n", OVERFLOW_DEFAULT);
  fprintf(stdout, "UNDERFLOW_DEFAULT = %d\n", UNDERFLOW_DEFAULT);
#endif
#endif
  fprintf(stdout, "\n");
#ifdef POINT_OUTPUT
  fprintf(stdout, "Stimulus value for every timestep written to file\n");
  fprintf(stdout, "Ez at center of mesh for every timestep written to file\n");
#endif
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  // print out some info on each timestep:
  fprintf(stdout, "Following is the iteration number and current\n");
  fprintf(stdout, "simulated time for each timestep/iteration of\n");
  fprintf(stdout, "the simulation.  For each timestep that 3D data is\n");
  fprintf(stdout, "output to file, the maximum and minimum data\n");
  fprintf(stdout, "values are printed here with the maximum and\n");
  fprintf(stdout, "minimum scaled values in parentheses.\n");
  fprintf(stdout, "\n");

  /////////////////////////////////////////////////////////////////////////////
  // open and start writing the .viz file:
  //    this file will be handy to feed parameters to viz if desired
  while ((vizFilePointer = fopen("ToyBoxFDTDbezhig.viz", "w")) == NULL)
    {
      fprintf(stderr, "Difficulty opening ToyBoxFDTDbezhig.viz");
      perror(" ");
    }
  fprintf(vizFilePointer, "#Viz V1.0\n");
  fprintf(vizFilePointer, "time: %lg %lg\n", currentSimulatedTime, dt);
  fprintf(vizFilePointer, "color: chengGbry2.cmap\n");
  fprintf(vizFilePointer, "\n");

#ifdef POINT_OUTPUT
  /////////////////////////////////////////////////////////////////////////////
  // open the stimulus data and center data output files:
  centerFilePointer = fopen("bezhigCenter.dat", "w");
  stimulusFilePointer = fopen("bezhigStimulus.dat", "w");
#endif

  /////////////////////////////////////////////////////////////////////////////
  // main loop:

  for(iteration = 0; iteration < MAXIMUM_ITERATION; iteration++)
    {// mainloop

      /////////////////////////////////////////////////////////////////////////
      // Output section:

      // time in simulated seconds that the simulation has progressed:
      currentSimulatedTime = dt*(double)iteration;  
      // print to standard output the iteration number 
      //     and current simulated time:
      fprintf(stdout, "#%d %lgsec", iteration, currentSimulatedTime);

#ifdef POINT_OUTPUT
      // write values to the center data file:
      fprintf(centerFilePointer, "%lg %lg\n", currentSimulatedTime, 
	      ez[nx/2][ny/2][nz/2]);
      fflush(centerFilePointer);
#endif

      // 3D data output every PLOT_MODULO timesteps:
      //     The first time through the main loop all the data written to 
      //     file will be zeros.  If anything is nonzero, there's a bug.  :>

      // write out a bob file every PLOT_MODULO timesteps:
      if ( (iteration % PLOT_MODULO) == 0)
	{// bob output section

	  // print to standard output the filename for this iteration, 
	  //     which includes the iteration number:
	  sprintf(filename, "eh%05d.bob", iteration);
	  // open a new data file for this iteration:
	  while ((openFilePointer = fopen(filename, "wb")) == NULL)
	    {// if the file fails to open, print an error message to 
	     //     standard output:
	      fprintf(stderr, "Difficulty opening eh%05d.bob", iteration);
	      perror(" ");
	    }

#ifndef NO_MIN_MAX_CALC
	  // find the max and min values to be output this timestep:  
	  min = FLT_MAX;
	  max = -FLT_MAX;
	  for(i=0; i<(nx+1); i++)
	    { 
	      for(j=0; j<(ny+1); j++)
		{
		  for(k=0; k<(nz); k++)
		    {
		      if (ez[i][j][k] < min) min = ez[i][j][k];
		      if (ez[i][j][k] > max) max = ez[i][j][k];
		    }
		}
	    }

	  // update the tracking variables for minimum and maximum 
	  //      values for the entire simulation:
	  if (min < simulationMin) simulationMin = min;
	  if (max > simulationMax) simulationMax = max;
#endif

#ifndef NO_AUTOSCALING
	  // set the output scaling value for this timestep:
	  // set norm to be max or min, whichever is greater in magnitude:
	  norm = (fabs(max) > fabs(min)) ? fabs(max) : fabs(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;
#endif
	  // 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:
	  fprintf(stdout, "\t%lg(%d) < ez BoB < %lg(%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 255 and write them to the output 
	  //     file for this iteration:
	  for(k=0; k<(nz); k++)
	    { 
	      for(j=0; j<(ny+1); j++)
		{
		  for(i=0; i<(nx+1); i++)
		    {
#ifdef NO_AUTOSCALING
#ifdef IDIOTPROOF_SCALING
		      nextOut = (int)(127.0 + scalingValue*ez[i][j][k]);
		      if (nextOut < 0) nextOut = UNDERFLOW_DEFAULT;
		      else if (nextOut > 255) nextOut = OVERFLOW_DEFAULT;
		      putc(nextOut, openFilePointer);
#else 
		      putc((int)(127.0 + 
				 scalingValue*ez[i][j][k]), openFilePointer);
#endif
#endif

#ifndef NO_AUTOSCALING
		      putc((int)(127.0 + 
				 scalingValue*ez[i][j][k]), openFilePointer);
#endif
		    }
		}
	    }

	  // close the output file for this iteration: