toyfdtd1.c

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

		    {
		      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 
	  //     iteration to the viz file:
	  fprintf(vizFilePointer, "%dx%dx%d %s\n", nx+1, ny+1, nz, filename);
	  // write identification of the corners of the mesh and the max
	  //     and min values for this iteration to the viz file:
	  fprintf(vizFilePointer, "bbox: 0.0 0.0 0.0 %lg %lg %lg %lg %lg\n",
		  dx*(double)nx, dy*(double)ny, dz*(double)nz, min, max);

	}// end bob output section

      /////////////////////////////////////////////////////////////////////////
      // Compute the stimulus: a plane wave emanates from the x=0 face:
      //     The length of the guide lies in the x-direction, the width of the 
      //     guide lies in the y-direction, and the height of the guide lies
      //     in the z-direction.  So the guide is sourced by all the ez 
      //     components on the stimulus face.  

      stimulus = sin(omega*currentSimulatedTime);
      for (i=0; i<(1); i++)
	{ 
	  for(j=0; j<(ny+1); j++)
	    {
	      for(k=0; k<nz; k++)
		{
		  ez[i][j][k] = stimulus;
		}
	    }
	}

      /////////////////////////////////////////////////////////////////////////
      // Update the interior of the mesh:
      //    all vector components except those on the faces of the mesh.  
      //
      // Update all the H field vector components within the mesh:
      //     Since all H vectors are internal, all H values are updated here.
      //     Note that the normal H vectors on the faces of the mesh are not 
      //     computed here, and in fact were never allocated -- the normal 
      //     H components on the faces of the mesh are never used to update 
      //     any other value, so they are left out of the memory allocation  
      //     entirely.   
      
      // Update the hx values:
      for(i=0; i<(nx-1); i++)
	{  
	  for(j=0; j<(ny); j++)
	    {
	      for(k=0; k<(nz); k++)
		{
		  hx[i][j][k] += (dtmudz*(ey[i+1][j][k+1] - ey[i+1][j][k]) - 
		                  dtmudy*(ez[i+1][j+1][k] - ez[i+1][j][k]));
		}
	    }
	}

      // Update the hy values:
      for(i=0; i<(nx); i++)
	{  
	  for(j=0; j<(ny-1); j++)
	    {
	      for(k=0; k<(nz); k++)
		{
		  hy[i][j][k] +=  (dtmudx*(ez[i+1][j+1][k] - ez[i][j+1][k]) -
      		                   dtmudz*(ex[i][j+1][k+1] - ex[i][j+1][k]));
		}
	    }
	}

      // Update the hz values:
      for(i=0; i<(nx); i++)
	{  
	  for(j=0; j<(ny); j++)
	    {
	      for(k=0; k<(nz-1); k++)
		{
		  hz[i][j][k] +=  (dtmudy*(ex[i][j+1][k+1] - ex[i][j][k+1]) -
      		                   dtmudx*(ey[i+1][j][k+1] - ey[i][j][k+1]));
		}
	    }
	}

      // Update the E field vector components.  
      // The values on the faces of the mesh are not updated here; they 
      //      are handled by the boundary condition computation 
      //      (and stimulus computation).  

      // Update the ex values:
      for(i=0; i<(nx); i++)
	{  
	  for(j=1; j<(ny); j++)
	    {
	      for(k=1; k<(nz); k++)
		{
		  ex[i][j][k] += (dtepsdy*(hz[i][j][k-1] - hz[i][j-1][k-1]) -
				  dtepsdz*(hy[i][j-1][k] - hy[i][j-1][k-1]));
		}
	    }
	}      

      // Update the ey values:
      for(i=1; i<(nx); i++)
	{  
	  for(j=0; j<(ny); j++)
	    {
	      for(k=1; k<(nz); k++)
		{
		  ey[i][j][k] += (dtepsdz*(hx[i-1][j][k] - hx[i-1][j][k-1]) -
				  dtepsdx*(hz[i][j][k-1] - hz[i-1][j][k-1]));
		}
	    }
	}

      // Update the ez values:
      for(i=1; i<(nx); i++)
	{  
	  for(j=1; j<(ny); j++)
	    {
	      for(k=0; k<(nz); k++)
		{
		  ez[i][j][k] += (dtepsdx*(hy[i][j-1][k] - hy[i-1][j-1][k]) -
				  dtepsdy*(hx[i-1][j][k] - hx[i-1][j-1][k]));
		}
	    }
	}
      
      fprintf(stdout, "\n");

      /////////////////////////////////////////////////////////////////////////
      // Compute the boundary conditions:

      // OK, so I'm yanking your chain on this one.  The PEC condition is 
      // enforced by setting the tangential E field components on all the
      // faces of the mesh to zero every timestep (except the stimulus 
      // face).  But the lazy/efficient way out is to initialize those 
      // vectors to zero and never compute them again, which is exactly 
      // what happens in this code.  
      
    }// end mainloop

  ///////////////////////////////////////////////////////////////////// 
  // Output section:  
  // The output routine is repeated one last time to write out  
  // the last data computed.   
  
  // 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 for the last timestep: 
  // 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 
  //     iteration to the viz file:
  fprintf(vizFilePointer, "%dx%dx%d %s\n", nx+1, ny+1, nz, filename);
  // write identification of the corners of the mesh and the max
  //     and min values for this iteration to the viz file:
  fprintf(vizFilePointer, "bbox: 0.0 0.0 0.0 %lg %lg %lg %lg %lg\n",
	  dx*(double)nx, dy*(double)ny, dz*(double)nz, min, max);

  /////////////////////////////////////////////////////////////////////////////
  // close the viz file for this simulation:
  fclose(vizFilePointer);

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

  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 was %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"); 
  // print out simulationMin and simulationMax: 
  fprintf(stdout, "Minimum output value was %lg \n", simulationMin); 
  fprintf(stdout, "Maximum output value was %lg \n", simulationMax); 
  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"); 
  fprintf(stdout, "\n"); 
  fprintf(stdout, "\n"); 

}// end main