toyfdtd1.c

toyFDTD一个简单的三维FDTD的程序用Fortran编的

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

  // 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");
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "ToyFDTD1 version 1.03\n");
  fprintf(stdout, "Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, ");
  fprintf(stdout, "Matthew O'Keefe\n");
  fprintf(stdout, "\n");
  fprintf(stdout, "ToyFDTD1 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 ToyFDTD1c.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("ToyFDTD1c.viz", "w")) == NULL)
    {
      fprintf(stderr, "Difficulty opening ToyFDTD1c.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++)