toyfdtd2.c

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

// ToyFDTD2, version 1.0 
//      The if-I-can-do-it-you-can-do-it FDTD! 
// Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe
// Copyright (C) 1999 John Schneider 

// This program is free software; you can redistribute it and/or 
//     modify it under the terms of the GNU General Public License 
//     as published by the Free Software Foundation; either version 2
//     of the License, or any later version, with the following conditions
//     attached in addition to any and all conditions of the GNU
//     General Public License:
//     When reporting or displaying any results or animations created
//     using this code or modification of this code, make the appropriate
//     citation referencing ToyFDTD2 by name and including the version
//     number.  
//
// This program is distributed in the hope that it will be useful,
//     but WITHOUT ANY WARRANTY; without even the implied warranty 
//     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//     See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
//     along with this program; if not, write to the Free Software
//     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  
//     02111-1307  USA

// Contacting the perpetratorss:
//
// Laurie E. Miller, Paul Hayes, Matthew O'Keefe
// Department of Electrical and Computer Engineering
//      200 Union Street S. E.
//      Minneapolis, MN 55455
//
// [e-mail]
// 
// http://www.borg.umn.edu/toyfdtd/
// http://www.borg.umn.edu/toyfdtd/ToyFDTD2.html
// http://www.toyfdtd.org/

// This code is here for everyone, but not everyone will need something 
//      so simple, and not everyone will need to read all the comments.  
// This file is over 700 lines long, but less than 400 of that is actually
//      code. 

///////////////////////////////////////////////////////////////////////////////
// ToyFDTD2 builds into ToyFDTD1 a way of allocating memory so that each 
//      field verctor array is forced to be be contiguous, which the
//      ToyFDTD1 scheme does not.  In many cases this will be faster
//      than the ToyFDTD1 scheme, but in some cases it will be slower;
//      which will be better is highly platform-dependent.  Some optimizing
//      compilers simply force the ToyFDTD1 arrays to be contiguous.  It's
//      worth testing more than one option before undertaking a large run.
//      Both schemes are designed to be easy to read.  In ToyFDTD2, macros
//      are defined for accessing the arrays to make the notation easy to 
//      read: i.e. Ez(i,j,k).  The modifications were written by 
//      John Schneider.
// This ToyFDTD2 is a stripped-down, minimalist 3D FDTD code.  It 
//      illustrates the minimum factors that must be considered to 
//      create a simple FDTD simulation.  
///////////////////////////////////////////////////////////////////////////////

// This is a very simple Yee algorithm 3D FDTD code in C implementing
//      the free space form of Maxwell's equations on a Cartesian grid.
// There are no internal materials or geometry.  
// The code as delivered simulates an idealized rectangular waveguide 
//      by treating the interior of the mesh as free space/air and enforcing
//      PEC (Perfect Electric Conductor) conditions on the faces of the mesh.
// The problem is taken from Field and Wave Electromagnetics, 2nd ed., by 
//      David K. Cheng, pages 554-555.  It is a WG-16 waveguide 
//      useful for X-band applications, interior width = 2.29cm, 
//      interior height = 1.02cm.  The frequency (10 GHz) is chosen to be 
//      in the middle of the frequency range for TE10 operation.  
// Boundaries: PEC (Perfect Electric Conductor).
// Stimulus: A simplified sinusoidal plane wave emanates from x = 0 face.
// 3D output: The electric field intensity vector components in the direction 
//      of the height of the guide (ez) are output to file every 
//      PLOT_MODULUS timesteps (and for the last timestep), scaled to 
//      the range of integers from zero through 254.  (The colormap 
//      included in the tar file assigns rgb values to the range zero through
//      255.)   Scaling is performed on each timestep individually, since 
//      it is not known in advance what the maximum and minimum 
//      values will be for the entire simulation.  The integer value 127 is 
//      held to be equal to the data zero for every timestep.  This method 
//      of autoscaling every timestep can be very helpful in a simulation 
//      where the intensities are sometimes strong and sometimes faint, 
//      since it will highlight the presence and structure of faint signals 
//      when stronger signals have left the mesh.  
//      Each timestep has it's own output file.  This data output file format
//      can be used in several visualization tools, such as animabob and viz. 
// Other output: Notes on the progress of the simulation are written to standard
//      output as the program runs.  
//      A .viz file is output to feed parameters to viz, should viz later be
//      used to view the data files.

// Some terminology used here:
//
// This code implements a Cartesian mesh with space differentials 
//     of dx, dy, dz.
// This means that a point in the mesh has neighboring points dx meters 
//     away in the direction of the positive and negative x-axis,
//     and neighboring points dy meters away in the directions 
//     of the +- y-axis, and neighboring points dz meters away 
//     in the directions of the +- z-axis,
// The mesh has nx cells in the x direction, ny in the y direction, 
//     and nz in the z direction.
// ex, ey, and ez refer to the arrays of electric field intensity vectors 
//     -- for example, ex is a 3-dimensional array consisting of the 
//     x component of the E field intensity vector for every point in the 
//     mesh.  ex[i][j][k] refers to the x component of the E field intensity 
//     vector at point [i][j][k].  
// hx, hy, and hz refer to the arrays of magnetic field intensity vectors.
//
// dt is the time differential -- the length of each timestep in seconds.
//
// bob is a file format that stands for "brick of bytes", meaning a string 
//     of bytes that can be interpreted as a 3-dimensional array of byte 
//     values (integers from zero through 255).  animabob is a free 
//     visualization tool that displays and animates a sequence of bricks 
//     of bytes.  For more information on animabob or to download a copy, 
//     see the ToyFDTD website at http://www.borg.umn.edu/toyfdtd/
//
// viz is another free visualization tool that displays and animates 
//     brick-of-byte files. For more information on viz or to download a copy, 
//     see the ToyFDTD website at http://www.borg.umn.edu/toyfdtd/

#include <math.h> 
#include <stdio.h>
#include <stdlib.h>
#include <float.h>

// program control constants
#define MAXIMUM_ITERATION 1000
          // total number of timesteps to be computed
#define PLOT_MODULUS 5
          // The program will output 3D data every PLOT_MODULUS timesteps,
          //     except for the last iteration computed, which is always
          //     output.  So if MAXIMUM_ITERATION is not an integer
          //     multiple of PLOT_MODULUS, the last timestep output will
          //     come after a shorter interval than that separating
          //     previous outputs.  
#define FREQUENCY 10.0e9     
          // frequency of the stimulus in Hertz
#define GUIDE_WIDTH 0.0229
          // meters
#define GUIDE_HEIGHT 0.0102
          // meters
#define LENGTH_IN_WAVELENGTHS 5.0
          // length of the waveguide in wavelengths of the stimulus wave
#define CELLS_PER_WAVELENGTH 25.0
          // minimum number of grid cells per wavelength in the x, y, and
          //     z directions

// physical constants
#define LIGHT_SPEED 299792458.0       
          // speed of light in a vacuum in meters/second
#define LIGHT_SPEED_SQUARED 89875517873681764.0        
          // m^2/s^2
#define MU_0 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6
          // permeability of free space in henry/meter
#define EPSILON_0 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12
          // permittivity of free space in farad/meter

// The value used for pi is M_PI as found in /usr/include/math.h on SGI 
//     IRIX 6.2.  Other such constants used here are DBL_EPSILON and FLT_MAX

main()	
{// main

  /////////////////////////////////////////////////////////////////////////////
  // variable declarations
  int i,j,k;     
           // indices of the 3D array of cells
  int ioff_ex, ioff_ey, ioff_ez, ioff_hx, ioff_hy, ioff_hz;
           // offsets used in the macros below for calculating position
           // within the data arrays
  int nx, ny, nz;         
           // total number of cells along the x, y, and z axes, respectively
  int allocatedBytes = 0;          
           // a counter to track number of bytes allocated

  int iteration = 0;          
           // counter to track how many timesteps have been computed
  double stimulus = 0.0;       
           // value of the stimulus at a given timestep
  double currentSimulatedTime = 0.0;
           // time in simulated seconds that the simulation has progressed
  double totalSimulatedTime = 0.0;       
           // time in seconds that will be simulated by the program
  double omega;       
           // angular frequency in radians/second
  double lambda;       
           // wavelength of the stimulus in meters
  double dx, dy, dz; 
           // space differentials (or dimensions of a single cell) in meters
  double dt;
           // time differential (how much time between timesteps) in seconds
  double dtmudx, dtepsdx;       
           // physical constants used in the field update equations 
  double dtmudy, dtepsdy;       
           // physical constants used in the field update equations 
  double dtmudz, dtepsdz;       
           // physical constants used in the field update equations 

  double *ex, *ey, *ez; 
           // pointers to the arrays of ex, ey, and ez values
#define Ex(I,J,K) ex[ioff_ex*(I) + (nz+1)*(J) + (k)]
#define Ey(I,J,K) ey[ioff_ey*(I) + (nz+1)*(J) + (k)]
#define Ez(I,J,K) ez[ioff_ez*(I) + nz*(J) + (k)]
           // macros to make accessing Ex, Ey, and Ez convenient
           // i.e. Ez(i,j,k)
  double *hx, *hy, *hz;
           // pointers to the arrays of hx, hy, and hz values
#define Hx(I,J,K) hx[ioff_hx*(I) + nz*(J) + (k)]
#define Hy(I,J,K) hy[ioff_hy*(I) + nz*(J) + (k)]
#define Hz(I,J,K) hz[ioff_hz*(I) + (nz-1)*(J) + (k)]
           // macros to make accessing Hx, Hy, and Hz convenient
           // i.e. Hz(i,j,k)

  // bob output routine variables:
  char filename[1024];
           // filename variable for 3D bob files
  double simulationMin = FLT_MAX;
           // tracks minimum value output by the entire simulation
  double simulationMax = -FLT_MAX;
           // tracks maximum value output by the entire simulation
  double min, max;
           // these track minimum and maximum values output in one timestep
  double norm;
           // norm is set each iteration to be max or min, whichever is 
           //     greater in magnitude
  double scalingValue;
           // multiplier used in output scaling, calculated every timestep
  FILE *openFilePointer;
           // pointer to 3D bob output file
  FILE *vizFilePointer;
           // pointer to viz file


  /////////////////////////////////////////////////////////////////////////////
  // setting up the problem to be modeled
  //
  // David K. Cheng, Field and Wave Electromagnetics, 2nd ed., 
  //     pages 554-555. 
  // Rectangular waveguide, interior width = 2.29cm, interior height = 1.02cm.
  // This is a WG-16 waveguide useful for X-band applications.
  //
  // Choosing nx, ny, and nz:
  // There should be at least 20 cells per wavelength in each direction, 
  //     but we'll go with 25 so the animation will look prettier.   
  //     (CELLS_PER_WAVELENGTH was set to 25.0 in the global 
  //     constants at the beginning of the code.)
  // The number of cells along the width of the guide and the width of 
  //     those cells should fit the guide width exactly, so that ny*dy 
  //     = GUIDE_WIDTH meters.  
  //     The same should be true for nz*dz = GUIDE_HEIGHT meters.  
  // dx is chosen to be dy or dz -- whichever is smaller
  // nx is chosen to make the guide LENGTH_IN_WAVELENGTHS 
  //     wavelengths long.  
  // 
  // dt is chosen for Courant stability; the timestep must be kept small 
  //     enough so that the plane wave only travels one cell length 
  //     (one dx) in a single timestep.  Otherwise FDTD cannot keep up 
  //     with the signal propagation, since FDTD computes a cell only from 
  //     it's immediate neighbors.  

  // wavelength in meters:
  lambda = LIGHT_SPEED/FREQUENCY; 
  // angular frequency in radians/second:
  omega = 2.0*M_PI*FREQUENCY;