toyfdtd1.f90

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

                  write(unit=69,rec=record) char(int(127.0d0+scalingValue*ez(i,j,k)))
                  record = record + 1
               end do
            end do
         end do
       
	 !! close the output file for this iteration:
         close(unit=69) 
	 !! write the dimensions and name of the output file for this 
	 !!     iteration to the viz file:
         write(unit=33,fmt=334) nx+1,ny+1,nz,filename
334      format(I4.4,"x",I4.4,"x",I4.4," ",A14)
	 !! write identification of the corners of the mesh and the max
	 !!     and min values for this iteration to the viz file:
         write(unit=33,fmt=335) dx*nx,dy*ny,dz*nz,min,max
335      format("bbox: 0.0 0.0 0.0 ",F7.4,F7.4,F7.4," ",F8.4," ",F8.4)

      end if !! 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)
      do i=0,0,1
         do j=0,ny,1
            do k=0,(nz-1),1
               ez(i,j,k) = stimulus
            end do
         end do
      end do
    
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! 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:
      do i=0,(nx-2),1
         do j=0,(ny-1),1
            do k=0,(nz-1),1
               hx(i,j,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))
            end do
         end do
      end do
    
      !! Update the hy values:
      do i=0,(nx-1),1
         do j=0,(ny-2),1
            do k=0,(nz-1),1
               hy(i,j,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))
            end do
         end do
      end do
    
      !! Update the hz values:
      do i=0,(nx-1),1
         do j=0,(ny-1),1
            do k=0,(nz-2),1
               hz(i,j,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))
            end do
         end do
      end do
    
      !! 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:
      do i=0,(nx-1),1
         do j=1,(ny-1),1
            do k=1,(nz-1),1
               ex(i,j,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))
            end do
         end do
      end do
    
      !! Update the ey values:
      do i=1,(nx-1),1
         do j=0,(ny-1),1
            do k=1,(nz-1),1
               ey(i,j,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))
            end do
         end do
      end do
    
      !! Update the ez values:
      do i=1,(nx-1),1
         do j=1,(ny-1),1
            do k=0,(nz-1),1
               ez(i,j,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))
            end do
         end do
      end do
    
      print*," "
    
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! 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 do !! 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*iteration                       
  !! print to standard output the iteration number 
  !!     and current simulated time:
  write(*,372) iteration, currentSimulatedTime

  !! 3D data output for the last timestep: 
  !! create the filename for this iteration, 
  !!     which includes the iteration number:

  !! 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 min and max values for this iteration
  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.0) then
     norm = DBL_EPSILON
  end if
  scalingValue = 127.0d0/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)
  
  !! 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
           write(unit=69,rec=record) char(int(127.0d0+scalingValue*ez(i,j,k)))
           record = record + 1
        end do
     end do
  end do
  
  !! close the output file for this iteration:
  close(unit=69) 
  !! write the dimensions and name of the output file for this 
  !!     iteration to the viz file:
  write(unit=33,fmt=334) 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:
  write(unit=33,fmt=335) dx*nx,dy*ny,dz*nz,min,max

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! close the viz file for this simulation:
  close(unit=33)

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

  print*," "
  print*," "
  !! print out how much memory has been allocated: 
  print*,"Dynamically allocated ", allocatedBytes, " bytes"
  !! 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 was ",E12.5E2," seconds, ",I5.5," timesteps")', &
       & totalSimulatedTime, MAXIMUM_ITERATION
  print*," "
  print*,"3D output scaling parameters:"
  print*,"Autoscaling every timestep"
  print*," "
  print*," "
  !! print out simulationMin and simulationMax: 
  print'(" Minimum output value was : ", f10.5)', simulationMin
  print'(" Maximum output value was : ", f10.5)', simulationMax
  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*," "
  print*," "

end program ToyFDTD1 !! end main