toyfdtd2.c

builds an alternate memory allocation scheme into ToyFDTD1. Contributed by John Schneider, it guaran

  // 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
  //     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:

  // compute offset constants used in macros:
  ioff_ex = (ny+1)*(nz+1);
  ioff_ey = ny*(nz+1);
  ioff_ez = (ny+1)*nz;
  // allocate arrays:
  ex = calloc(nx*(ny+1)*(nz+1),sizeof(double));
  ey = calloc((nx+1)*ny*(nz+1),sizeof(double));
  ez = calloc((nx+1)*(ny+1)*nz,sizeof(double));
  // check if allocation successful
  if (ex == NULL || ey == NULL || ez == NULL) {
    fprintf(stderr,"E-field allocation failed.  Terminating...\n");
    exit(-1);
  }
  // keep rough track of allocated memory:
  allocatedBytes += ( (nx)*(ny+1)*(nz+1) * sizeof(double));
  allocatedBytes += ( (nx+1)*(ny)*(nz+1) * sizeof(double));
  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.  

  // compute offset constants used in macros:
  ioff_hx = ny*nz;
  ioff_hy = (ny-1)*nz;
  ioff_hz = ny*(nz-1);
  // allocate arrays:
  hx = calloc((nx-1)*ny*nz,sizeof(double));
  hy = calloc(nx*(ny-1)*nz,sizeof(double));
  hz = calloc(nx*ny*(nz-1),sizeof(double));
  // check if allocation successful 
  if (hx == NULL || hy == NULL || hz == NULL) {
    fprintf(stderr,"E-field allocation failed.  Terminating...\n");
    exit(-1);
  }
  // keep rough track of allocated memory:
  allocatedBytes += ( (nx-1)*(ny)*(nz) * sizeof(double));
  allocatedBytes += ( (nx)*(ny-1)*(nz) * sizeof(double));
  allocatedBytes += ( (nx)*(ny)*(nz-1) * sizeof(double));


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

  // print out some identifying information 
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "ToyFDTD2 version 1.0\n");
  fprintf(stdout, "Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, ");
  fprintf(stdout, "Matthew O'Keefe\n");
  fprintf(stdout, "Copyright (C) 1999 John Schneider\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "ToyFDTD2 is free software published under the terms\n"); 
  fprintf(stdout, "of the GNU General Public License as published by the\n"); 
  fprintf(stdout, "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 chengGbry.cmap -s %dx%dx%d *.bob\n", nx+1, ny+1, nz); 
  fprintf(stdout, "\n"); 
  // print out a viz command line: 
  fprintf(stdout, "viz ToyFDTD2c.viz\n"); 
  fprintf(stdout, "\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_MODULUS = %d\n", PLOT_MODULUS); 
  fprintf(stdout, "rectangular waveguide\n"); 
  fprintf(stdout, "Stimulus = %lg Hertz ", FREQUENCY); 
  fprintf(stdout, "continuous plane wave\n"); 
  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"); 
  fprintf(stdout, "Autoscaling every timestep\n"); 
  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("ToyFDTD2c.viz", "w")) == NULL)
    {
      fprintf(stderr, "Difficulty opening ToyFDTD2c.viz");
      perror(" ");
    }
  fprintf(vizFilePointer, "#Viz V1.0\n");
  fprintf(vizFilePointer, "time: %lg %lg\n", currentSimulatedTime, dt);
  fprintf(vizFilePointer, "color: chengGbry.cmap\n");
  fprintf(vizFilePointer, "\n");


  /////////////////////////////////////////////////////////////////////////////
  // 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);

      // 3D data output every PLOT_MODULUS 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.  :>
 
      if ( (iteration % PLOT_MODULUS) == 0)
	{// bob output section

	  // create the filename for this iteration, 
	  //     which includes the iteration number:
	  sprintf(filename, "c_%06d.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 c_%06d.bob", iteration);
	      perror(" ");
	    }

	  // 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;

	  // 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;
	  // 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 254 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++)
		    {
		      putc((int)(127.0 + 
				 scalingValue*Ez(i,j,k)), openFilePointer);
		    }
		}
	    }
	  
	  // close the output file for this iteration:
	  fclose(openFilePointer);

	  // write the dimensions and name of the output file for this