toyfdtd1.f90

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

  else
    dx = dz
  end if
  !! 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.0d0/(LIGHT_SPEED*sqrt(1.0d0/(dx*dx) + 1.0d0/(dy*dy) + 1.0d0/(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.
  !! 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:
  allocate(ex(0:nx-1,0:ny,0:nz))
  do i=0,(nx-1),1
     do j=0,ny,1
        do k=0,nz,1
           ex(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = allocatedBytes + ( (nx)*(ny+1)*(nz+1) * SIZEOF_DOUBLE)

  !! allocate the array of ey components:
  allocate(ey(0:nx,0:ny-1,0:nz))
  do i=0,nx,1
     do j=0,(ny-1),1
        do k=0,nz,1
           ey(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = allocatedBytes + ( (nx+1)*(ny)*(nz+1) * SIZEOF_DOUBLE)

  !! allocate the array of ez components:
  allocate(ez(0:nx,0:ny,0:nz-1))
  do i=0,nx,1
     do j=0,ny,1
        do k=0,(nz-1),1
           ez(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = 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:
  allocate(hx(0:nx-2,0:ny-1,0:nz-1))
  do i=0,(nx-2),1
     do j=0,(ny-1),1
        do k=0,(nz-1),1
           hx(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = allocatedBytes + ( (nx-1)*(ny)*(nz) * SIZEOF_DOUBLE)

  !! allocate the array of hy components:
  allocate(hy(0:nx-1,0:ny-2,0:nz-1))
  do i=0,(nx-1),1
     do j=0,(ny-2),1
        do k=0,(nz-1),1
           hy(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = allocatedBytes + ( (nx)*(ny-1)*(nz) * SIZEOF_DOUBLE)

  !! allocate the array of hz components:
  allocate(hz(0:nx-1,0:ny-1,0:nz-2))
  do i=0,(nx-1),1
     do j=0,(ny-1),1
        do k=0,(nz-2),1
           hz(i,j,k) = 0.0d0
        end do
     end do
  end do
  allocatedBytes = allocatedBytes + ( (nx)*(ny)*(nz-1) * SIZEOF_DOUBLE)

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! write some progress notes to standard output

  !! print out some identifying information 
  print*," "
  print*," "
  print*," "
  print*," "
  print*,"ToyFDTD1 version 1.03 (F90)"
  print*,"Copyright (C) 1998, 1999 Laurie E. Miller, Paul Hayes, "
  print*,"Matthew O'Keefe, Max Smirnoff, Matt Rundquist"
  print*," "
  print*,"ToyFDTD1 is free software published under the terms" 
  print*,"of the GNU General Public License as published by the" 
  print*,"Free Software Foundation. " 
  print*," "
  print*," "
  !! print out a bob command line, 
  !! including the dimensions of the output files:
  print'(" bob -cmap chengGbry.cmap -s ",I5.5,"x",I5.5,"x",I5.5," *.bob")',nx+1,ny+1,nz
  print*," "
  !! print out a viz command line:
  print*,"viz ToyFDTD1f90.viz"
  print*," "
  print*," "
  !! print out how much memory has been allocated: 
  print*,"Dynamically allocated ", allocatedBytes, " bytes"
  print*," "
  !! print out some simulation parameters: 
  print*,"PLOT_MODULUS = ", PLOT_MODULUS
  print*,"rectangular waveguide"
  print*,"Stimulus = ", FREQUENCY, " Hertz continuous plane wave"
  print*," "
  print*,"Meshing parameters:"
  print '(I4,"x",I3,"x",I3," cells")', nx, ny, nz
  print '(" dx=",F11.8," dy=",F11.8," dz=",F11.8," meters")', dx, dy, dz
  print '(F11.8, " x ", F11.8, " x ", F11.8, " meter^3 simulation region")', &
       & GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda
  print*," "
  print '(" Time simulated will be ",E12.5E2," seconds, ",I5.5," timesteps")', &
       & totalSimulatedTime, MAXIMUM_ITERATION
  print*," "
  print*,"3D output scaling parameters:"
  print*,"Autoscaling every timestep"
  print*," "
  print*," "
  !! print out some info on each timestep: 
  print*,"Following is the iteration number and current "
  print*,"simulated time for each timestep/iteration of "
  print*,"the simulation.  For each timestep that data is "
  print*,"output to file, the maximum and minimum data "
  print*,"values are printed here with the maximum and "
  print*,"minimum scaled values in parentheses. "
  print*," "

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! open and start writing the .viz file
  !!    this file will be handy to feed parameters to viz if desired

  open(unit=33,file="ToyFDTD1f90.viz",iostat=ios)
  !! if the file fails to open, print an error message to 
  !!     standard output:
  if(ios /= 0) then
    print*,"Difficulty opening ToyFDTD1f90.viz"
    STOP
  end if
  write(unit=33, fmt='("#Viz V1.0")')
  write(unit=33, fmt='("time: ",F5.3,ES12.5E2)') currentSimulatedTime,dt
  write(unit=33, fmt='("color: chengGbry.cmap")')
  write(unit=33, fmt='("")')                                     

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! main loop:

  do iteration=0,(MAXIMUM_ITERATION-1),1
     
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! Output section: 

      !! time in simulated seconds that the simulation has progressed:
      currentSimulatedTime = dt*iteration
      !! print to standard output the iteration number 
      !!     and current simulated time:
      write(*,372) iteration, currentSimulatedTime
372   format ("#",i5," ", ES14.5E2, "sec", $)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! 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(mod(iteration,PLOT_MODULUS) == 0) then

         !! Construct the filename:
         !! The numerical part is obtained by successively
         !!      dividing the current iteration by ten, and the
         !!      remainders of these divisions are the digits.
         !!      For example, 12345 would be converted as follows:
         !!      12345/10 = 1234+5/10 (5 is the first digit in the number)
         !!      1234/10 = 123 + 4/10 (4 is the next digit)
         !!      ....and so on..this is performed six times to create
         !!      a filename out of any 6-digit number.
         temp = iteration
         filename(1:4) = "f90_"
         filename(10:10) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(9:9) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(8:8) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(7:7) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(6:6) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(5:5) = numbers(mod(temp,10)+1:mod(temp,10)+1)
         temp = int(temp/10)
         filename(11:14) = ".bob"

         !! open a new data file for this iteration:
         open(unit=69,file=filename,access="direct",iostat=ios,recl=1)
         !! if the file fails to open, print an error message to 
	 !!     standard output:
         if(ios /= 0) then
           print*,"Difficulty opening ",filename
           STOP
         end if

         !! find the max and min values to be output this timestep:  
         min = FLT_MAX
         max = -FLT_MAX
         do i=0,(nx),1
            do j=0,(ny),1
               do k=0,(nz-1),1
                  if(ez(i,j,k) < min) then
                     min = ez(i,j,k)
                  end if
                  if(ez(i,j,k) > max) then
                     max = ez(i,j,k)
                  end if
               end do
            end do
         end do

         !!  update the tracking variables for minimum and maximum 
         !!  values for the entire simulation:
         if(min < simulationMin) then
            simulationMin = min
         end if
         if(max > simulationMax) then
            simulationMax = max
         end if
       
         !! set norm to be max or min, whichever is greater in magnitude:
         if(abs(max) > abs(min)) then
            norm = abs(max)
         else
            norm = abs(min)
         end if

         !! if everything is zero, give norm a tiny value 
	 !!     to avoid division by zero:
         if(norm == 0.0d0) then
            norm = DBL_EPSILON
         end if
         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:
         write(*, 373) min, int(127.0d0 + scalingValue * min), max, int(127.0d0 + scalingValue * max)
373      format(" ", f9.5, " (", i3, ") < ez BoB < ", f9.5, " (", i3, ")", $)

         !! 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:
         record = 1
         do k=0,(nz-1),1
            do j=0,(ny),1
               do i=0,(nx),1