📄 myfdtd_2d_parallel.f90
字号:
Program main
use omp_lib
real Ez(-100:100,-100:100)
real Hx(-100:100,-100:100)
real Hy(-100:100,-100:100)
real Ein(-800:800)
real Hin(-800:800)
integer IncidentStart,IncidentEnd,Isource
real EBin(1:4)
real ECB(-100:100,1:4)
real HCB(-100:100,1:4)
real EAB(-100:100,1:4)
real EAC(4,2,0:1)
integer ob(-100:100,-100:100)
real FE1(0:10,0:10)
real FE2(0:10,0:10)
real FH1(0:10,0:10)
real FH2(0:10,0:10)
integer TEMflag
real FE
real FH
real MU0
real EPS0
real SMU0
real SEPS0
real PI
real OMIGA
integer WL
real Phi
integer TimeStop
real MaxX
real MinX
real MaxY
real MinY
real Step
integer XNum
integer YNum
integer Imax
integer Imin
integer Jmax
integer Jmin
integer Itmax
integer Itmin
integer Jtmax
integer Jtmin
integer Iomax
integer Iomin
integer Jomax
integer Jomin
integer Ismax
integer Ismin
integer Jsmax
integer Jsmin
real cp
real sp
character * 30 FileName
character * 1 Blank
integer FileNameLength
logical FileExist
character * 1 OUTflag
character * 1 Respond
character * 2 WaveMode
character * 30 rem
real k,Z,WaveLength
real Temp,T1,T2,T3,T4,T5,T6,EI0,ER0
complex CT1,CT2,CT3,CT4
integer I,J,N,II,III,IFlag,NN
real Media(1:4,0:10)
integer MediaNumber
complex Jz
complex Mx
complex My
complex Ctemp,Ctemp1
complex E(-800:800,4)
complex H(-800:800,4)
real StartAngle,EndAngle,DegreeStep
real rk
integer tid,imi,ima
!*************************************Parallel Initializing**************************************************************
call cpu_time(time0)
nthreads=2
call omp_set_num_threads(nthreads)
!*************************************Initializing**************************************************************
WaveLength=1
WL=40
TEMFlag=1 !1TM,-1TE
Ismin=-60
Ismax=60
Jsmin=-60
Jsmax=60
Itmin=-65
Itmax=65
Jtmin=-65
Jtmax=65
T1=Itmin*Itmin+Jtmin*Jtmin
T2=Itmax*Itmax+Jtmin*Jtmin
T3=Itmin*Itmin+Jtmax*Jtmax
T4=Itmax*Itmax+Jtmax*Jtmax
Isource=-sqrt(max(T1,T2,T3,T4))-3
Imin=-75
Imax=75
Jmin=-75
Jmax=75
Iomin=-70
Iomax=70
Jomin=-70
Jomax=70
IncidentStart=-710
IncidentEnd=710
Phi=0.0 !degree
Phi=Phi/180.0*Pi !rad
TimeStop=1200
! write(*,*) char(7)
! write(*,'(11X,23h Now initializing)')
MediaNumber=1
Media(1,1)=1.0
Media(2,1)=1.0
Media(3,1)=3.72e+7
Media(4,1)=0.0
Media(1,0)=1.0
Media(2,0)=1.0
Media(3,0)=0.0
Media(4,0)=0.0
do I=Imin,Imax
do J=Jmin,Jmax
ob(I,J)=0
end do
end do
do I=Ismin,Ismax
do J=Jsmin,Jsmax
if((I.ge.-WL).and.(I.le.WL-1).and.(J.ge.-WL).and.(J.le.WL-1))then
ob(I,J)=1
end if
end do
end do
Pi=3.14159265
MU0=Pi*0.0000004
EPS0=8.85e-12
Z=sqrt(MU0/EPS0)
k=2*Pi/WaveLength
do I=1,MediaNumber
Media(3,I)=-Media(3,I)*Z/k
Media(4,I)=-Media(4,I)/Z/k
end do
if(TEMFlag.eq.-1)then
Temp=EPS0
EPS0=MU0
MU0=Temp
Z=sqrt(MU0/EPS0)
end if
SMU0=sqrt(MU0)
SEPS0=sqrt(EPS0)
OMIGA=Pi/WL
cp=cos(Phi)
sp=sin(Phi)
do i=Imin,Imax
do j=Jmin,Jmax
Ez(i,j)=0
Hy(i,j)=0
Hx(i,j)=0
end do
end do
do i=IncidentStart,IncidentEnd
Ein(i)=0
Hin(i)=0
end do
do i=1,4
EBin(i)=0
end do
do i=Imin,Imax
do j=1,4
ECB(i,j)=0
HCB(i,j)=0
EAB(i,j)=0
end do
end do
do i=1,4
do j=1,2
do T1=0,1
EAC(i,j,T1)=0
end do
end do
end do
FE=TEMFlag*0.5
FH=TEMFlag*0.5
if(TEMFlag.eq.1)then
do I=0,MediaNumber
do J=0,MediaNumber
FE1(I,J)=((Media(1,I)+Media(1,J))+(Media(3,I)+Media(3,J))/2*Pi/WL)/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
FE2(I,J)=1.0/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
FH1(I,J)=((Media(2,I)+Media(2,J))+(Media(4,I)+Media(4,J))/2*Pi/WL)/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
FH2(I,J)=1.0/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
end do
end do
else
do I=0,MediaNumber
do J=0,MediaNumber
FH1(I,J)=((Media(1,I)+Media(1,J))+(Media(3,I)+Media(3,J))/2*Pi/WL)/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
FH2(I,J)=-1.0/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
FE1(I,J)=((Media(2,I)+Media(2,J))+(Media(4,I)+Media(4,J))/2*Pi/WL)/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
FE2(I,J)=-1.0/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
end do
end do
end if
!************************************Main Loop***************************************************************
OUTflag='I'
N=0
999 Continue
do while (N.LT.TimeStop)
N=N+1
write(*,1)N
1 format(1H+,10X,'Time step',I5,' is in processing......')
!***************************Incidentwave**********************************
do I=IncidentStart,IncidentEnd-1
Hin(i)=Hin(i)-FH*(Ein(i+1)-Ein(i))
end do
do I=IncidentStart+1,IncidentEnd-1
Ein(i)=Ein(i)+FE*(Hin(i-1)-Hin(i))
end do
Ein(Isource)=sin(OMIGA*N)
if(N.le.WL) Ein(Isource)=Ein(Isource)*0.5*(1.0-cos(OMIGA*N))
Ein(IncidentStart)=EBin(4)
EBin(4)=EBin(3)
EBin(3)=Ein(IncidentStart+1)
Ein(IncidentEnd)=EBin(2)
EBin(2)=EBin(1)
EBin(1)=Ein(IncidentEnd-1)
!***********************Connective Boundary*****************************
do I=Itmin,Itmax
T1=float(I)*cp+float(Jtmin)*sp
II=int(T1)
if(T1.GE.0.0)then
III=II+1
T1=float(III)-T1
else
III=II-1
T1=T1-float(III)
end if
ECB(I,1)=T1*(Ein(II)-Ein(III))+Ein(III) !down
T1=float(I)*cp+float(Jtmax)*sp
II=int(T1)
if(T1.GE.0.0)then
III=II+1
T1=float(III)-T1
else
III=II-1
T1=T1-float(III)
end if
ECB(I,3)=T1*(Ein(II)-Ein(III))+Ein(III) !up
end do
do J=Jtmin,Jtmax
T1=float(Itmin)*cp+float(J)*sp
II=int(T1)
if(T1.GE.0.0)then
III=II+1
T1=float(III)-T1
else
III=II-1
T1=T1-float(III)
end if
ECB(J,4)=T1*(Ein(II)-Ein(III))+Ein(III) !left
T1=float(Itmax)*cp+float(J)*sp
II=int(T1)
if(T1.GE.0.0)then
III=II+1
T1=float(III)-T1
else
III=II-1
T1=T1-float(III)
end if
ECB(J,2)=T1*(Ein(II)-Ein(III))+Ein(III) !right
end do
do I=Itmin,Itmax
T1=float(I)*cp+(float(Jtmin)-0.5)*sp
if(T1.GT.0.5)then
II=int(T1-0.5)
III=II+1
T1=float(III)+0.5-T1
else
II=int(T1+0.5)
III=II-1
T1=T1-0.5-float(III)
end if
HCB(i,1)=(T1*(Hin(II)-Hin(III))+Hin(III))*sp !down
T1=float(I)*cp+(float(Jtmax)+0.5)*sp
if(T1.GT.0.5)then
II=int(T1-0.5)
III=II+1
T1=float(III)+0.5-T1
else
II=int(T1+0.5)
III=II-1
T1=T1-0.5-float(III)
end if
HCB(i,3)=(T1*(Hin(II)-Hin(III))+Hin(III))*sp !up
end do
do J=Jtmin,Jtmax
T1=float(J)*sp+(float(Itmin)-0.5)*cp
if(T1.GT.0.5)then
II=int(T1-0.5)
III=II+1
T1=float(III)+0.5-T1
else
II=int(T1+0.5)
III=II-1
T1=T1-0.5-float(III)
end if
HCB(J,4)=-(T1*(Hin(II)-Hin(III))+Hin(III))*cp !left
T1=float(J)*sp+(float(Itmax)+0.5)*cp
if(T1.GT.0.5)then
II=int(T1-0.5)
III=II+1
T1=float(III)+0.5-T1
else
II=int(T1+0.5)
III=II-1
T1=T1-0.5-float(III)
end if
HCB(J,2)=-(T1*(Hin(II)-Hin(III))+Hin(III))*cp !right
end do
!********************Ez FDTD*****************************
!$omp parallel private(tid,imi,ima,I,J,T1)
tid=omp_get_thread_num()
imi=Imin+(Imax-Imin-1)*tid/nthreads+1
ima=Imin+(Imax-Imin-1)*(tid+1)/nthreads
do I=imi,ima
do J=Jmin+1,Jmax-1
T1=Hy(I,J)-Hy(I-1,J)-Hx(I,J)+Hx(I,J-1)
Ez(I,J)=FE1(ob(I,J),ob(I,J))*Ez(I,J)+FE2(ob(I,J),ob(I,J))*T1
end do
end do
!$omp end parallel
!********************Ez Connective Boundary*************************
do I=Itmin+1,Itmax-1
Ez(I,Jtmax)=Ez(I,Jtmax)-FE*HCB(I,3)
Ez(I,Jtmin)=Ez(I,Jtmin)+FE*HCB(I,1)
end do
do J=Jtmin+1,Jtmax-1
Ez(Itmax,J)=Ez(Itmax,J)+FE*HCB(J,2)
Ez(Itmin,J)=Ez(Itmin,J)-FE*HCB(J,4)
end do
T1=HCB(Jtmax,2)-HCB(Itmax,3)
Ez(Itmax,Jtmax)=Ez(Itmax,Jtmax)+FE*T1
T1=-HCB(Jtmin,4)+HCB(Itmin,1)
Ez(Itmin,Jtmin)=Ez(Itmin,Jtmin)+FE*T1
T1=HCB(Jtmin,2)+HCB(Itmax,1)
Ez(Itmax,Jtmin)=Ez(Itmax,Jtmin)+FE*T1
T1=-HCB(Jtmax,4)-HCB(Itmin,3)
Ez(Itmin,Jtmax)=Ez(Itmin,Jtmax)+FE*T1
!**********************Ez Absorbing Boundary*****************
do I=Imin+1,Imax-1
T1=(Hy(I,Jmin)-Hy(I-1,Jmin)+Hy(I,Jmin+1)-Hy(I-1,Jmin+1))*TEMFlag
T2=Ez(I,Jmin+1)-Ez(I,Jmin)
Ez(I,Jmin)=EAB(I,1)-0.33333333*T2
Ez(I,Jmin)=Ez(I,Jmin)+0.16666667*T1
EAB(I,1)=Ez(I,Jmin+1)
T1=(Hy(I,Jmax)-Hy(I-1,Jmax)+Hy(I,Jmax-1)-Hy(I-1,Jmax-1))*TEMFlag
T2=Ez(I,Jmax-1)-Ez(I,Jmax)
Ez(I,Jmax)=EAB(I,3)-0.33333333*T2
Ez(I,Jmax)=Ez(I,Jmax)+0.16666667*T1
EAB(I,3)=Ez(I,Jmax-1)
end do
do J=Jmin+1,Jmax-1
T1=(Hx(Imin,J)-Hx(Imin,J-1)+Hx(Imin+1,J)-Hx(Imin+1,J-1))*TEMFlag
T2=Ez(Imin+1,J)-Ez(Imin,J)
Ez(Imin,J)=EAB(J,4)-0.33333333*T2
Ez(Imin,J)=Ez(Imin,J)-0.16666667*T1
EAB(J,4)=Ez(Imin+1,J)
T1=(Hx(Imax,J)-Hx(Imax,J-1)+Hx(Imax-1,J)-Hx(Imax-1,J-1))*TEMFlag
T2=Ez(Imax-1,J)-Ez(Imax,J)
Ez(Imax,J)=EAB(J,2)-0.33333333*T2
Ez(Imax,J)=Ez(Imax,J)-0.16666667*T1
EAB(J,2)=Ez(Imax-1,J)
end do
Iflag=1-Iflag
Ez(Imin,Jmin)=0.2928932*EAC(1,1,Iflag)+0.7071068*EAC(1,2,Iflag)
EAC(1,1,Iflag)=Ez(Imin,Jmin)
EAC(1,2,Iflag)=Ez(Imin+1,Jmin+1)
Ez(Imax,Jmin)=0.2928932*EAC(2,1,Iflag)+0.7071068*EAC(2,2,Iflag)
EAC(2,1,Iflag)=Ez(Imax,Jmin)
EAC(2,2,Iflag)=Ez(Imax-1,Jmin+1)
Ez(Imin,Jmax)=0.2928932*EAC(3,1,Iflag)+0.7071068*EAC(3,2,Iflag)
EAC(3,1,Iflag)=Ez(Imin,Jmax)
EAC(3,2,Iflag)=Ez(Imin+1,Jmax-1)
Ez(Imax,Jmax)=0.2928932*EAC(4,1,Iflag)+0.7071068*EAC(4,2,Iflag)
EAC(4,1,Iflag)=Ez(Imax,Jmax)
EAC(4,2,Iflag)=Ez(Imax+1,Jmax+1)
!*******************Hx,Hy FDTD*********************
!$omp parallel private(tid,imi,ima,I,J,T1)
tid=omp_get_thread_num()
imi=Imin-1+(Imax-Imin+1)*tid/nthreads+1
ima=Imin-1+(Imax-Imin+1)*(tid+1)/nthreads
do I=imi,ima
do J=Jmin,Jmax-1
T1=Ez(I,J+1)-Ez(I,J)
Hx(I,J)=FH1(ob(i,j),ob(i,j+1))*Hx(I,j)-FH2(ob(i,j),ob(i,j+1))*T1
end do
end do
imi=Imin-1+(Imax-Imin)*tid/nthreads+1
ima=Imin-1+(Imax-Imin)*(tid+1)/nthreads
do I=imi,ima
do J=Jmin,Jmax
T1=Ez(I+1,J)-Ez(I,J)
Hy(I,J)=FH1(ob(i,j),ob(i+1,j))*Hy(I,j)+FH2(ob(i,j),ob(i+1,j))*T1
end do
end do
!$omp end parallel
!****************Hx,Hy Connective Boundary*************************
do I=Itmin,Itmax
Hx(I,Jtmin-1)=Hx(I,Jtmin-1)+FH*ECB(I,1)
Hx(I,Jtmax)=Hx(I,Jtmax)-FH*ECB(I,3)
end do
do J=Jtmin,Jtmax
Hy(Itmin-1,J)=Hy(Itmin-1,J)-FH*ECB(J,4)
Hy(Itmax,J)=Hy(Itmax,J)+FH*ECB(J,2)
end do
end do
!************************Output*************************************
if(OUTflag.eq.'I')then
! write(*,'(1H+,10X,40h Now outputing IMAGE part)')
Open(1,FILE='FDTD_2D.EI')!,ACCESS='DIRECT',RECL=4)
Open(4,FILE='FDTD_2D.HxI')!,ACCESS='DIRECT',RECL=4)
Open(5,FILE='FDTD_2D.HyI')!,ACCESS='DIRECT',RECL=4)
Open(2,FILE='FDTD_2D.BI')
EI0=Ein(0)
else
! write(*,'(1H+,10X,40h Now outputing REAL part)')
Open(1,FILE='FDTD_2D.ER')!,ACCESS='DIRECT',RECL=4)
Open(4,FILE='FDTD_2D.HxR')!,ACCESS='DIRECT',RECL=4)
Open(5,FILE='FDTD_2D.HyR')!,ACCESS='DIRECT',RECL=4)
Open(2,FILE='FDTD_2D.BR')
ER0=Ein(0)
end if
NN=0
do I=Imin,Imax
do J=Jmin,Jmax
NN=NN+1
write(1,*)Ez(I,J)
write(4,*)Hx(I,J)/Z
write(5,*)Hy(I,J)/Z
end do
end do
Close(1)
Close(4)
Close(5)
do I=Iomin,Iomax
T1=0.5*(Ez(I,Jomin)+Ez(I,Jomin-1))
write(2,*)T1,Hx(I,Jomin-1)/Z
T1=0.5*(Ez(I,Jomax)+Ez(I,Jomax+1))
write(2,*)T1,Hx(I,Jomax)/Z
end do
do J=Jomin,Jomax
T1=0.5*(Ez(Iomax,J)+Ez(Iomax+1,J))
write(2,*)T1,Hy(Iomax,J)/Z
T1=0.5*(Ez(Iomin,J)+Ez(Iomin-1,J))
write(2,*)T1,Hy(Iomin-1,J)/Z
end do
Close(2)
if (OUTflag.eq.'I')then
OUTflag='R'
TimeStop=TimeStop+WL/2
GOTO 999
end if
!*******************R+I=C**************************************
! write(*,'(1H+,10X,40h Make BND Files)')
CT1=(1.0,0.0)/CMPLX(ER0,EI0)
CT2=CEXP(CMPLX(0,-0.5*OMIGA))*CT1
Open(1,FILE='FDTD_2D.BI')
Open(2,FILE='FDTD_2D.BR')
FileName='rect.BND'
Open(3,FILE=FileName)!,ACCESS='DIRECT',RECL=8)
! NN=1
do III=1,10000
read(1,*,end=111)T1,T2
read(2,*)T3,T4
CT3=CMPLX(T3,T1)*CT1
write(3,*)CT3
CT4=CMPLX(T4,T2)*CT2
write(3,*)CT4
! NN=NN+1
end do
111 Continue
Close(1)
Close(2)
Close(3)
!************************************************RCS*****************************************
FileName='rect.BND'
Open(1,File=FileName)!,ACCESS='DIRECT',RECL=8)
rk=2*Pi/float(WL)
! NN=0
do I=Iomin,Iomax
! NN=NN+1
read(1,*)E(I,1)
! NN=NN+1
read(1,*)H(I,1)
! NN=NN+1
read(1,*)Ctemp
E(I,3)=-Ctemp
! NN=NN+1
read(1,*)Ctemp1
H(I,3)=-Ctemp1
end do
do I=Jomin,Jomax
! NN=NN+1
read(1,*)E(I,2)
! NN=NN+1
read(1,*)H(I,2)
! NN=NN+1
read(1,*)Ctemp
E(I,4)=-Ctemp
! NN=NN+1
read(1,*)Ctemp1
H(I,4)=-Ctemp1
end do
Close(1)
if(TEMFlag.eq.-1)then
do I=Iomin,Iomax
E(I,1)=E(I,1)
E(I,3)=E(I,3)
H(I,1)=-H(I,1)
H(I,3)=-H(I,3)
end do
do I=Jomin,Jomax
E(I,2)=E(I,2)
E(I,4)=E(I,4)
H(I,2)=-H(I,2)
H(I,4)=-H(I,4)
end do
end if
Open(2,FILE='rect.dat')
!*****************************RCS Input******************
StartAngle=0
EndAngle=180
DegreeStep=3
!**********************Jz,Mx,My*************************
do dg=StartAngle,EndAngle,DegreeStep
Jz=(0.0,0.0)
Mx=(0.0,0.0)
My=(0.0,0.0)
! write(*,1)dg
!2 format(1H+,10x,'Now calculate degree')
Phi=dg*Pi/180.0
do I=Iomin,Iomax
Ctemp=cmplx(0.0,(float(I)*cos(Phi)+(Jomin-0.5)*sin(Phi))*rk)
Ctemp=cexp(Ctemp)
Mx=Mx+E(I,1)*Ctemp
Jz=Jz+H(i,1)*Ctemp
Ctemp=cmplx(0.0,(float(I)*cos(Phi)+(Jomax+0.5)*sin(Phi))*rk)
Ctemp=cexp(Ctemp)
Mx=Mx+E(I,3)*Ctemp
Jz=Jz+H(i,3)*Ctemp
end do
do I=Jomin,Jomax
Ctemp=cmplx(0.0,(float(I)*sin(Phi)+(Iomax+0.5)*cos(Phi))*rk)
Ctemp=cexp(Ctemp)
My=My+E(I,2)*Ctemp
Jz=Jz+H(i,2)*Ctemp
Ctemp=cmplx(0.0,(float(I)*sin(Phi)+(Iomin-0.5)*cos(Phi))*rk)
Ctemp=cexp(Ctemp)
My=My+E(I,4)*Ctemp
Jz=Jz+H(i,4)*Ctemp
end do
!**************************RCS Calculation*************************
RCS=abs(Mx*sin(Phi)-My*cos(Phi)+Z*Jz)
RCS=0.25*rk*RCS*RCS*WaveLength/float(WL)/WaveLength
! write(*,1)dg,rcs
RCS=10.0*log10(RCS)
write(2,*)RCS
end do
Close(2)
call cpu_time(time1)
write (*,*) 'Elapsed time: ', (time1 - time0)
read(*,*)
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -