toyboxfdtdbezhig.c

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

	  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 sinusoidal pulse 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.  
      // The pulse is one complete wavelength of a shifted sinusoid that 
      //     ranges from zero to one in intensity.  So the stimulus varies 
      //     sinusoidally from zero to one to zero again and then terminates
      //     --the x=0 face becomes a PEC thereafter.  
      //
      // compute the current stimulus value:
#ifdef SINUSOIDAL_PULSE_STIMULUS
      if (currentSimulatedTime <= 1.0/FREQUENCY)
	{
	  stimulus = .5*(1.0 + sin(omega*currentSimulatedTime - M_PI*.5));
	}
      else
	{
	  stimulus = 0.0;
	}
#else
      stimulus = sin(omega*currentSimulatedTime);
#endif
      // set all vectors on the x=0 face to the value of stimulus:
      for (i=0; i<(1); i++)
	{ 
	  for(j=0; j<(ny+1); j++)
	    {
	      for(k=0; k<nz; k++)
		{
		  ez[i][j][k] = stimulus;
		}
	    }
	}
#ifdef POINT_OUTPUT
      // write the current stimulus value to bezhigStimulus.dat:
      fprintf(stimulusFilePointer, "%lg %lg\n",
	      currentSimulatedTime, stimulus);
      fflush(stimulusFilePointer);
#endif
	
      /////////////////////////////////////////////////////////////////////////
      // 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:

#ifdef PMC
      // Compute the PMC symmetry boundaries:
      // j = 0 face, j = ny face
      for(i=0; i<nx; i++)
	{	
	  for(k=0; k<(nz+1); k++)
	    {
	      ex[i][0][k] = ex[i][ny-1][k];
	      ex[i][ny][k] = ex[i][1][k];
	    }
	}
      for(i=0; i<(nx+1); i++)
	{	
	  for(k=0; k<nz; k++)
	    {
	      ez[i][0][k] = ez[i][ny-1][k];
	      ez[i][ny][k] = ez[i][1][k];
	    }
	}
#endif

      // Compute the PEC boundaries:
      //
      // OK, so I'm yanking your chain on this one.  The PEC condition is 
      // enforced by setting the tangential E field components on the PEC
      // 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);
  
#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:
  
  // 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:
      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

  /////////////////////////////////////////////////////////////////////////////
  // close the various files for this simulation:
  fclose(vizFilePointer);
#ifdef POINT_OUTPUT
  fclose(stimulusFilePointer);
  fclose(centerFilePointer);
#endif

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

  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 simulationMin and simulationMax:
  fprintf(stdout, "\n");
  fprintf(stdout, "\n");
  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 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");

}// end main