fdtd，pml.txt

ＦＤＴＤ !-- Fortran code for FDTD with Berenger PMLs, version 1.0, May 1999 !-- by Jos Bergervoet.

program test_fdtd

  !-- Fortran code for FDTD with Berenger PMLs, version 1.0, May 1999
  !-- by Jos Bergervoet.
  !-- Plot field and/or Poynting vector S around radiating linear dipole

    implicit none
    integer, parameter :: wp = selected_real_kind(5) ! Choose working precision
                                                     ! 5 or 15 (single/double)
    real(wp),parameter :: c = 299792458,  &
                          nu = 2          ! nu=order of Berenger increase

    real(wp), allocatable, dimension(:)        :: sxe,sye,sze, sxh,syh,szh
    real(wp), allocatable, dimension(:,:,:,:)  :: exs,eys,ezs, hxs,hys,hzs
    character(len=20)  :: name
    integer            :: nx,ny,nz, nb, nt, ns, n
    real(wp)           :: Pi, mu0, Z0, eps0,  &
                          fun, g, t, period, smax, F_in, dx, dy, dz, dt

!    g(t) = sin(t)**6                      ! smooth, 5 times differentiable
!    g(t) = sin(2*t) * sin(t)**5           ! `double-delta'
!    fun(t) = g(min(t/period,1.0_wp)*Pi)   ! scaled g(t)
    fun(t) = sin(2*Pi*t/period)            ! CW

    print "(1x,a)",  &
      "Give:  nx,ny,nz, nt, nb, smax, period, n_s",  &
      "For instance: 60 60 30  310  10 4e8  110 30",  &
      "",  &
      "nb and smax are Berenger parameters, n_s = nr. of segm. in antenna",  &
      ""

    read *, nx,ny,nz, nt, nb,smax, period, ns

    call setup()
    
    n = 0
    Fdtd_loop: do
        n=n+1

        !-- divide source signal over both E-field species
        F_in = fun(n+0.0_wp) 
        exs(nx/2-1:nx/2,ny/2,nz/2,:) = - F_in/2

        if (mod(n, max(nint(period/30),1))==0) then
          write(name,*) n
          name = trim(adjustl(name)) // ".gif"
          !-- Plot S and some field component
          call plot_log("S"//name, compute_S_z(k_z=nz/2+ns/10))
          call plot("E"//name, exs(:,:,nz/2+ns/10,1)+exs(:,:,nz/2+ns/10,2) )
        end if

        if (n>nt) exit Fdtd_loop

        !-- update the fields including the Berenger PMLs
        call Berenger(exs,eys,ezs, hxs,hys,hzs)

        !-- apply PEC condition
        exs((nx-ns)/2:(nx+ns)/2-1, ny/2, nz/2, :) = 0
        
    end do Fdtd_loop

    stop


contains

  subroutine setup()
    real(wp) :: rho, x
    integer  :: nmax, i
    rho(x,nmax) = max(-x,x-nmax,0.0_wp)/max(nb,1)

    Pi  = acos(-1d0)
    mu0 = 4d-7*Pi
    Z0  = mu0*c
    eps0= 1/(Z0*c)

    dx = 1; dy = 1; dz = 1
    
    dt = 1 / ( c*sqrt(1/dx**2+1/dy**2+1/dz**2) )         ! acc. to Courant

    allocate( exs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &   ! All arrays include
              eys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &   ! the Berenger PML
              ezs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &   ! region (nb cells)
              hxs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &   ! All fields have two
              hys(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2),  &   ! species (last index)
              hzs(-nb:nx+nb,-nb:ny+nb,-nb:nz+nb,2) )
    allocate( sxe(-nb:nx+nb),  &
              sye(-nb:ny+nb),  &                         ! these are the sigma
              sze(-nb:nz+nb),  &                         ! arrays (scaled by
              sxh(-nb:nx+nb),  &                         ! eps0 and mu0)
              syh(-nb:ny+nb),  &
              szh(-nb:nz+nb) )

    exs = 0; eys = 0; ezs = 0; hxs = 0; hys = 0; hzs = 0

    sxe = (/( smax * rho(i+1.0_wp,nx)**nu, i=-nb,nx+nb )/) ! 1.0 is clumsy ?
    sxh = (/( smax * rho(i+0.5_wp,nx)**nu, i=-nb,nx+nb )/)
    sye = (/( smax * rho(i+1.0_wp,ny)**nu, i=-nb,ny+nb )/)
    syh = (/( smax * rho(i+0.5_wp,ny)**nu, i=-nb,ny+nb )/)
    sze = (/( smax * rho(i+1.0_wp,nz)**nu, i=-nb,nz+nb )/)
    szh = (/( smax * rho(i+0.5_wp,nz)**nu, i=-nb,nz+nb )/)

    write(*,*) "R_normal = ", exp(-2.0_wp/(nu+1) * smax*nb/c)

  end subroutine
        

  function compute_S(k_z)  result(S)
    integer, intent(in)  :: k_z
    real(wp), dimension(-nb+1:nx+nb-1, -nb+1:ny+nb-1)  :: S
    real(wp), dimension(-nb:nx+nb, -nb:ny+nb, 3)       :: E, H
    real(wp), dimension(-nb:nx+nb, -nb:ny+nb)          :: temp
    integer  :: i
    
    do i=1,3
      E(:,:,1) = exs(:,:,k_z,1) + exs(:,:,k_z,2)  ! problem: E and H are not
      E(:,:,2) = eys(:,:,k_z,1) + eys(:,:,k_z,2)  !   available for exactly the
      E(:,:,3) = ezs(:,:,k_z,1) + ezs(:,:,k_z,2)  !   same positions
      H(:,:,1) = hxs(:,:,k_z,1) + hxs(:,:,k_z,2)
      H(:,:,2) = hys(:,:,k_z,1) + hys(:,:,k_z,2)
      H(:,:,3) = hzs(:,:,k_z,1) + hzs(:,:,k_z,2)
    end do
    
!-- Use:  S^2 = (E x H)^2 = E^2 H^2 - (E dot H)^2
    temp = sqrt( sum(E*E, dim=3)*sum(H*H, dim=3) - sum(E*H, dim=3)**2 )

    temp( (nx-ns)/2:(nx+ns)/2-1, ny/2 ) = 0    ! Just to show the dipole-wire
    S = temp(-nb+1:nx+nb-1, -nb+1:ny+nb-1)
    
  end function compute_S
        

  function compute_S_z(k_z)  result(S_z)
    integer, intent(in)  :: k_z
    real(wp), dimension(-nb+1:nx+nb-1, -nb+1:ny+nb-1)  :: S_z
    real(wp), dimension(-nb:nx+nb, -nb:ny+nb, 3)       :: E, H
    real(wp), dimension(-nb:nx+nb, -nb:ny+nb)          :: temp
    integer  :: i
    
    do i=1,3
      E(:,:,1) = exs(:,:,k_z,1) + exs(:,:,k_z,2)  ! problem: E and H are not
      E(:,:,2) = eys(:,:,k_z,1) + eys(:,:,k_z,2)  !   available for exactly the
      E(:,:,3) = ezs(:,:,k_z,1) + ezs(:,:,k_z,2)  !   same positions
      H(:,:,1) = hxs(:,:,k_z,1) + hxs(:,:,k_z,2)
      H(:,:,2) = hys(:,:,k_z,1) + hys(:,:,k_z,2)
      H(:,:,3) = hzs(:,:,k_z,1) + hzs(:,:,k_z,2)
    end do
    
!-- Use:  S_z = E_x H_y - E_y H_x
    temp = E(:,:,1)*H(:,:,2) - E(:,:,2)*H(:,:,1) 

    temp( (nx-ns)/2:(nx+ns)/2-1, ny/2 ) = 0    ! Just to show the dipole-wire
    S_z = temp(-nb+1:nx+nb-1, -nb+1:ny+nb-1)
    
  end function compute_S_z
  
  
  subroutine plot_log(Name, v)
  !-- logscale picture, peak-hold range between calls
    use gif_util
    character(len=*) :: Name
    real(wp)         :: range=0, ndecades=3, v(:,:)
    integer          :: i, Pixel(size(v,1), size(v,2))
    type(color)      :: ColMap(256)

    ColMap = (/  ( color(2*i-2,0,0), i=1,128 ),    &        ! black-to-red
                 ( color(255,5*i/2,0), i=1,101 ),  &        ! red-to-yellow
                 ( color(255,255,9*i), i=1,27 )   /)        ! yellow-to-white

    range = max(range, maxval(abs(v)), 1e-30_wp )
    
    Pixel = nint( max(0.0,1+log10(abs(v)/range)/ndecades) * 255 )

    write(*,*) "Max(|",trim(Name),"|) =", maxval(abs(v)) 
    call writegif (Name, Pixel, ColMap)
  end subroutine plot_log
  

  subroutine plot(Name, v)
  !-- linear scale picture, peak-hold range between calls
    use gif_util
    character(len=*) :: Name
    real(wp)         :: range=0, v(:,:)
    integer          :: i, Pixel(size(v,1), size(v,2))
    type(color)      :: ColMap(256)

    ColMap = (/( color(i-1,i-1,256-i), i=1,256 )/)

    range = max(range, maxval(abs(v)), 1e-30_wp )
    
    Pixel = nint( (v+range)/range/2 * 255 )

    write(*,*) "Max(|",trim(Name),"|) =", maxval(abs(v)) 
    call writegif (Name, Pixel, ColMap)
  end subroutine plot
  

  subroutine Berenger(Ex,Ey,Ez,Hx,Hy,Hz)
  !-- Updates all fields, except the outer PEC faces of the domain.
  !-- Uses internal fields and E-fields on faces (not H on faces)

    real(wp), dimension(0:,0:,0:,:) :: Ex,Ey,Ez,Hx,Hy,Hz

    call add_diff_y(syh, -mu0,  Hx(1:,0:,0:,1), Ez(1:,0:,0:,:) )
    call add_diff_z(szh,  mu0,  Hx(1:,0:,0:,2), Ey(1:,0:,0:,:) )

    call add_diff_z(szh, -mu0,  Hy(0:,1:,0:,1), Ex(0:,1:,0:,:) )
    call add_diff_x(sxh,  mu0,  Hy(0:,1:,0:,2), Ez(0:,1:,0:,:) )
 
    call add_diff_x(sxh, -mu0,  Hz(0:,0:,1:,1), Ey(0:,0:,1:,:) )
    call add_diff_y(syh,  mu0,  Hz(0:,0:,1:,2), Ex(0:,0:,1:,:) )

    call add_diff_y(sye,  eps0, Ex(0:,1:,1:,1), Hz(0:,0:,1:,:) )
    call add_diff_z(sze, -eps0, Ex(0:,1:,1:,2), Hy(0:,1:,0:,:) )
 
    call add_diff_z(sze,  eps0, Ey(1:,0:,1:,1), Hx(1:,0:,0:,:) )
    call add_diff_x(sxe, -eps0, Ey(1:,0:,1:,2), Hz(0:,0:,1:,:) )
 
    call add_diff_x(sxe,  eps0, Ez(1:,1:,0:,1), Hy(0:,1:,0:,:) )
    call add_diff_y(sye, -eps0, Ez(1:,1:,0:,2), Hx(1:,0:,0:,:) )
 
  end subroutine


  subroutine add_diff_x(s, pfac, v, w)
    integer  :: i,j,k
    real(wp) :: v(:,:,:), w(:,:,:,:)
    real(wp) :: pfac, alpha(size(v,1)), beta(size(v,1)), s(*)

  !-- Using exponential stepping, as described by Katz, Thiele en Tavlove
  !-- in IEEE microw., Vol 4, Nr. 8, p268(3), 1994

    call coeff(pfac, dx, s, alpha, beta, size(v,1))
    do k=1,size(v,3)-1
      do j=1,size(v,2)-1
        do i=1,size(v,1)-1
          v(i,j,k) = alpha(i) * v(i,j,k) + beta(i) *  &
              ( w(i+1,j,k,1) + w(i+1,j,k,2)               &
              - w(i,  j,k,1) - w(i,  j,k,2) )
        end do
      end do
    end do

  end subroutine

  subroutine add_diff_y(s, pfac, v, w)
    integer  :: i,j,k
    real(wp) :: v(:,:,:), w(:,:,:,:)
    real(wp) :: pfac, alpha(size(v,2)), beta(size(v,2)), s(*)

    call coeff(pfac, dy, s, alpha, beta, size(v,2))
    do k=1,size(v,3)-1
      do j=1,size(v,2)-1
        do i=1,size(v,1)-1
          v(i,j,k) = alpha(j) * v(i,j,k) + beta(j) *  &
              ( w(i,j+1,k,1) + w(i,j+1,k,2)               &
              - w(i,j,  k,1) - w(i,j,  k,2) )
        end do
      end do
    end do

  end subroutine

  subroutine add_diff_z(s, pfac, v, w)
    integer  :: i,j,k
    real(wp) :: v(:,:,:), w(:,:,:,:)
    real(wp) :: pfac, alpha(size(v,3)), beta(size(v,3)), s(*)

    call coeff(pfac, dz, s, alpha, beta, size(v,3))
    do k=1,size(v,3)-1
      do j=1,size(v,2)-1
        do i=1,size(v,1)-1
          v(i,j,k) = alpha(k) * v(i,j,k) + beta(k) *  &
              ( w(i,j,k+1,1) + w(i,j,k+1,2)               &
              - w(i,j,k  ,1) - w(i,j,k  ,2) )
        end do
      end do
    end do

  end subroutine

  subroutine coeff(pfac, d, s, alpha, beta, n)
    integer  :: n, i
    real(wp) :: pfac, d, s(n), alpha(n), beta(n)
    real*8   :: t  ! double precision here, whatever the working precision

    do i=1,n
      t = exp(dble(-s(i)*dt))
      alpha(i) = t
      if(t .gt. 1d0 - 1d-7) then
        beta(i) = dt/(pfac*d)
      else
        beta(i) = (1-t)/(pfac*s(i)*d)
      end if
    end do

  end subroutine

end program test_fdtd