?? fdtd3d_cpml.f90
字號:
!***********************************************************************
! 3-D FDTD code with CPML absorbing boundary conditions
!***********************************************************************
!
! Program author: Jamesina J. Simpson, Ph.D. Graduate Student
! National Science Foundation Fellow
! NU Computational Electromagnetics Laboratory
! Director: Allen Taflove
! Department of Electrical Engineering and
! Computer Science
! Northwestern University
! 2145 Sheridan Road
! Evanston, IL 60208
! j-simpson@northwestern.edu
!
! Copyright 2005
!
! This FORTRAN 90 code implements the finite-difference time-domain
! solution of Maxwell's curl equations over a three-dimensional
! Cartesian space lattice. The grid is terminated by CPML absorbing
! boundary conditions. However, in a straightforward manner it can
! be altered to have PEC planes at any of the outer grid boundaries.
! Also, the number of grid cells as well as the thickness of
! the PML in each Cartesian direction can be varied independently.
!
! For illustrative purposes, the code supplied here models the PML
! numerical experiment described in Section 7.11.2. Three output
! files are generated having the following data:
!
! (1) Ez at the source for each time step;
! (2) Ey at the probe point for each time step; and
! (3) The plane of Ez values 1 mm above the source in the k-
! direction recorded at the time step, "record_grid". This
! data can be viewed in Matlab using the following commands:
!
! >> load view_grid.dat
! >> for j = 1:Jmax
! >> image(1:Imax,j) = view_grid(1+(j-1)*Imax:j*Imax);
! >> end
! >> pcolor(image'); shading flat
!
! The relative error (equation 7.135) graphed in Fig. 7.6 on
! page 318 can be reproduced by comparing the output file
! "probe.dat" with that of a much larger benchmark grid having
! Imax, Jmax, and Kmax increased to the values mentioned in the text.
!
! This code has been tested using the Intel Fortran Compiler 8.0 for
! Linux. The executable can be created by typing
! "ifort fdtd3D_CPML.f90" at the command prompt.
!
! This program is not guaranteed to be free of defects or bugs.
! Please report any bugs that may exist to:
! j-simpson@northwestern.edu
!
!
!***********************************************************************
PROGRAM fdtd3D_CPML
IMPLICIT NONE
! ..................................
! Input Fundamental Constants (MKS units)
REAL, PARAMETER :: &
pi = 3.14159265358979, C = 2.99792458E8, &
muO = 4.0 * pi * 1.0E-7, epsO = 1.0/(C*C*muO)
! ..................................
! Specify Material Relative Permittivity and Conductivity
REAL, PARAMETER:: &
epsR = 1.0, sigM1 = 0.0 ! free space
! ..................................
! Specify Number of Time Steps and Grid Size Parameters
INTEGER, PARAMETER :: &
nmax = 2100, & ! total number of time steps
Imax = 51, Jmax = 126, Kmax = 26 ! grid size corresponding to the
! number of Ez field components
! ..................................
! Specify Grid Cell Size in Each Direction and Calculate the
! Resulting Courant-Stable Time Step
REAL, PARAMETER :: &
dx = 1.0E-3, dy = 1.0E-3, dz = 1.0E-3, & ! cell size in each direction
dt = 0.99 / (C*(1.0/dx**2+1.0/dy**2+1.0/dz**2)**0.5)
! time step increment
! ..................................
! Specify the Impulsive Source (See Equation 7.134)
REAL, PARAMETER :: &
tw = 53.0E-12, tO = 4.0*tw
! ..................................
! Specify the Time Step at which the Grid is Recorded for Visualization
INTEGER, PARAMETER :: &
record_grid = 300
! ..................................
! Specify the PEC Plate Boundaries and the Source/Recording Points
INTEGER, PARAMETER :: &
istart = (Imax-1)/2-11, iend = istart+24, jstart = istart, &
jend = jstart + 99, kstart = Kmax/2, kend = kstart, &
isource = istart, jsource = jstart, ksource = kend+1, &
irecv1 = iend, jrecv1 = jend+1, krecv1 = kend ! Ey at probe point
! ..................................
! Specify the CPML Thickness in Each Direction (Value of Zero
! Corresponds to No PML, and the Grid is Terminated with a PEC)
INTEGER, PARAMETER :: &
! PML thickness in each direction
nxPML_1 = 11, nxPML_2 = nxPML_1, nyPML_1 = nxPML_1, &
nyPML_2 = nxPML_1, nzPML_1 = nxPML_1, nzPML_2 = nxPML_2
! ..................................
! Specify the CPML Order and Other Parameters
INTEGER, PARAMETER :: &
m = 3, ma = 1
REAL, PARAMETER :: &
sig_x_max = 0.75 * (0.8*(m+1)/(dx*(muO/epsO*epsR)**0.5)), &
sig_y_max = 0.75 * (0.8*(m+1)/(dy*(muO/epsO*epsR)**0.5)), &
sig_z_max = 0.75 * (0.8*(m+1)/(dz*(muO/epsO*epsR)**0.5)), &
alpha_x_max = 0.24, &
alpha_y_max = alpha_x_max, alpha_z_max = alpha_x_max, &
kappa_x_max = 15.0, &
kappa_y_max = kappa_x_max, kappa_z_max = kappa_x_max
INTEGER :: &
i,j,ii,jj,k,kk,n
REAL :: &
source, P1, P2
! TM components
REAL,DIMENSION(Imax, Jmax, Kmax) :: &
Ez, CA, CB, sig, eps
REAL,DIMENSION(Imax-1, Jmax, Kmax) :: &
Hy
REAL,DIMENSION(Imax,Jmax-1, Kmax) :: &
Hx
REAL :: &
DA, DB
! TE components
REAL,DIMENSION(Imax-1, Jmax-1, Kmax-1) :: &
Hz
REAL,DIMENSION(Imax-1, Jmax, Kmax-1) :: &
Ex
REAL,DIMENSION(Imax,Jmax-1, Kmax-1) :: &
Ey
! PML
REAL ,DIMENSION(nxPML_1,Jmax,Kmax) :: &
psi_Ezx_1
REAL ,DIMENSION(nxPML_2,Jmax,Kmax) :: &
psi_Ezx_2
REAL ,DIMENSION(nxPML_1-1,Jmax,Kmax) :: &
psi_Hyx_1
REAL ,DIMENSION(nxPML_2-1,Jmax,Kmax) :: &
psi_Hyx_2
REAL ,DIMENSION(Imax,nyPML_1,Kmax) :: &
psi_Ezy_1
REAL ,DIMENSION(Imax,nyPML_2,Kmax) :: &
psi_Ezy_2
REAL ,DIMENSION(Imax,nyPML_1-1,Kmax) :: &
psi_Hxy_1
REAL ,DIMENSION(Imax,nyPML_2-1,Kmax) :: &
psi_Hxy_2
REAL ,DIMENSION(Imax,Jmax-1,nzPML_1-1) :: &
psi_Hxz_1
REAL ,DIMENSION(Imax,Jmax-1,nzPML_2-1) :: &
psi_Hxz_2
REAL ,DIMENSION(Imax-1,Jmax,nzPML_1-1) :: &
psi_Hyz_1
REAL ,DIMENSION(Imax-1,Jmax,nzPML_2-1) :: &
psi_Hyz_2
REAL ,DIMENSION(Imax-1,Jmax,nzPML_1) :: &
psi_Exz_1
REAL ,DIMENSION(Imax-1,Jmax,nzPML_2) :: &
psi_Exz_2
REAL ,DIMENSION(Imax,Jmax-1,nzPML_1) :: &
psi_Eyz_1
REAL ,DIMENSION(Imax,Jmax-1,nzPML_2) :: &
psi_Eyz_2
REAL ,DIMENSION(nxPML_1-1,Jmax-1,Kmax-1) :: &
psi_Hzx_1
REAL ,DIMENSION(nxPML_1,Jmax-1,Kmax-1) :: &
psi_Eyx_1
REAL ,DIMENSION(nxPML_2-1,Jmax-1,Kmax-1) :: &
psi_Hzx_2
REAL ,DIMENSION(nxPML_2,Jmax-1,Kmax-1) :: &
psi_Eyx_2
REAL ,DIMENSION(Imax-1,nyPML_1-1,Kmax-1) :: &
psi_Hzy_1
REAL ,DIMENSION(Imax-1,nyPML_1,Kmax-1) :: &
psi_Exy_1
REAL ,DIMENSION(Imax-1,nyPML_2-1,Kmax-1) :: &
psi_Hzy_2
REAL ,DIMENSION(Imax-1,nyPML_2,Kmax-1) :: &
psi_Exy_2
REAL ,DIMENSION(nxPML_1) :: &
be_x_1, ce_x_1, alphae_x_PML_1, sige_x_PML_1, kappae_x_PML_1
REAL ,DIMENSION(nxPML_1-1) :: &
bh_x_1, ch_x_1, alphah_x_PML_1, sigh_x_PML_1, kappah_x_PML_1
REAL ,DIMENSION(nxPML_2) :: &
be_x_2, ce_x_2, alphae_x_PML_2, sige_x_PML_2, kappae_x_PML_2
REAL ,DIMENSION(nxPML_2-1) :: &
bh_x_2, ch_x_2, alphah_x_PML_2, sigh_x_PML_2, kappah_x_PML_2
REAL ,DIMENSION(nyPML_1) :: &
be_y_1, ce_y_1, alphae_y_PML_1, sige_y_PML_1, kappae_y_PML_1
REAL ,DIMENSION(nyPML_1-1) :: &
bh_y_1, ch_y_1, alphah_y_PML_1, sigh_y_PML_1, kappah_y_PML_1
REAL ,DIMENSION(nyPML_2) :: &
be_y_2, ce_y_2, alphae_y_PML_2, sige_y_PML_2, kappae_y_PML_2
REAL ,DIMENSION(nyPML_2-1) :: &
bh_y_2, ch_y_2, alphah_y_PML_2, sigh_y_PML_2, kappah_y_PML_2
REAL ,DIMENSION(nzPML_1) :: &
be_z_1, ce_z_1, alphae_z_PML_1, sige_z_PML_1, kappae_z_PML_1
REAL ,DIMENSION(nzPML_1-1) :: &
bh_z_1, ch_z_1, alphah_z_PML_1, sigh_z_PML_1, kappah_z_PML_1
REAL ,DIMENSION(nzPML_2) :: &
be_z_2, ce_z_2, alphae_z_PML_2, sige_z_PML_2, kappae_z_PML_2
REAL ,DIMENSION(nzPML_2-1) :: &
bh_z_2, ch_z_2, alphah_z_PML_2, sigh_z_PML_2, kappah_z_PML_2
! denominators for the update equations
REAL,DIMENSION(Imax-1) :: &
den_ex, den_hx
REAL,DIMENSION(Jmax-1) :: &
den_ey, den_hy
REAL,DIMENSION(Kmax-1) :: &
den_ez, den_hz
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! OPEN OUTPUT FILES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
OPEN (UNIT = 30, FILE = "source.dat")
OPEN (UNIT = 31, FILE = "probe.dat")
OPEN (UNIT = 33, FILE = "view_grid.dat")
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! INITIALIZE VARIABLES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P1 = 0.0
P2 = 0.0
Ez(:,:,:) = 0.0
Hz(:,:,:) = 0.0
Ex(:,:,:) = 0.0
Hx(:,:,:) = 0.0
Ey(:,:,:) = 0.0
Hy(:,:,:) = 0.0
sig(:,:,:) = sigM1
eps(:,:,:) = epsR*epsO
psi_Exy_1(:,:,:) = 0.0
psi_Exy_2(:,:,:) = 0.0
psi_Exz_1(:,:,:) = 0.0
psi_Exz_2(:,:,:) = 0.0
psi_Eyx_1(:,:,:) = 0.0
psi_Eyx_2(:,:,:) = 0.0
psi_Eyz_1(:,:,:) = 0.0
psi_Eyz_2(:,:,:) = 0.0
psi_Ezy_1(:,:,:) = 0.0
psi_Ezy_2(:,:,:) = 0.0
psi_Ezx_1(:,:,:) = 0.0
psi_Ezx_2(:,:,:) = 0.0
psi_Hxy_1(:,:,:) = 0.0
psi_Hxy_2(:,:,:) = 0.0
psi_Hxz_1(:,:,:) = 0.0
psi_Hxz_2(:,:,:) = 0.0
psi_Hyx_1(:,:,:) = 0.0
psi_Hyx_2(:,:,:) = 0.0
psi_Hyz_1(:,:,:) = 0.0
psi_Hyz_2(:,:,:) = 0.0
psi_Hzy_1(:,:,:) = 0.0
psi_Hzy_2(:,:,:) = 0.0
psi_Hzx_1(:,:,:) = 0.0
psi_Hzx_2(:,:,:) = 0.0
write(*,*)"Imax: ", Imax
write(*,*)"Jmax: ", Jmax
write(*,*)"Kmax: ", Kmax
write(*,*)"dt: ", dt
write(*,*)"nmax: ", nmax
write(*,*)"max time: ", nmax*dt
write(*,*)"record grid after ", record_grid, "dt"
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! SET CPML PARAMETERS IN EACH DIRECTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DO i = 1,nxPML_1
sige_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i) / (nxPML_1 - 1.0) )**m
alphae_x_PML_1(i) = alpha_x_max*((i-1.0)/(nxPML_1-1.0))**ma
kappae_x_PML_1(i) = 1.0+(kappa_x_max-1.0)* &
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -