te1.f90

这是一个计算PI的并行MPI程序 该程序已经调试通过了 可以直接计算 欢迎交流工作经验

program ex0101
implicit none

integer,parameter::imax=100                              !FDTD吸收边界
integer,parameter::imin=-100
integer,parameter::jmax=100
integer,parameter::jmin=-100

integer,parameter::itmax=80                              !连接边界
integer,parameter::itmin=-80
integer,parameter::jtmax=80
integer,parameter::jtmin=-80

real::ex(imin:imax,jmin:jmax-1)=0.0            !TE波分量的采样点的设置
real::ey(imin:imax-1,jmin:jmax)=0.0
real::hz(imin:imax,jmin:jmax)=0.0

real::ecb(-80:80,1:4)=0.                       !连接边界处E入射波分量
real::hcb(-80:80,1:4)=0.                       !连接边界处H入射波分量

real::ein(-710:710)=0.                         !入射波E采样点 
real::hin(-710:710)=0.                         !入射波H采样点
 
real::hac(1:4,1)=0.0                           !吸收边界条件角点临时变量


real,parameter::pi=3.14159265
real,parameter::mu0=pi*4.e-7
real,parameter::eps0=8.85e-12
real z                                         !自由空间阻抗
real,parameter::wavelength=632.8e-9            !波长
!real,parameter::k=2*pi/wavelength              !波矢
real,parameter::v=3.0e+8                       !自由空间电磁波传播速度
real,parameter::omiga=2*pi*v/wavelength        !角频率
real,parameter::phi=pi*(50.)/180.               !入射角度
real,parameter::istart=-710
real,parameter::iend=710 
real,parameter::sqrt2=1.41421356

real cp
real sp

real,parameter::dx=10.0e-9                     !空间步长
real,parameter::dy=10.0e-9
real,parameter::dd=10.0e-9                     
real,parameter::dt=dd/(2*3.e+8)                !时间步长
real::t1,t2                                    !临时存储变量
integer i,j,n,ii,iii,nn
integer::k
integer,parameter::timestop=2000               !总时间步数
real::media(1:4)                               !介质参数数组
character*1 outflag                            !输出标志

real ca(-100:100,-100:100)                                                                                                                                                                                                                                                                                                                                                                                                    
real cb(-100:100,-100:100)
real cr(-100:100,-100:100)
real cq(-100:100,-100:100)

real::hz1(imin:imax,jmin:jmax)=0.0             !前一时间步电场值
real::hz2(imin:imax,jmin:jmax)=0.0
real::hz3(imin:imax,jmin:jmax)=0.0
real::hz4(imin:imax,jmin:jmax)=0.0

z=sqrt(mu0/eps0)
cp=cos(phi)
sp=sin(phi)

media(1)=1.                                   !自由空间参数设置
media(2)=1.
media(3)=0.
media(4)=0.



!***************************上半空间****************************

do i=-100,100,1                                
do j=-100,100,1
ca(i,j)=(1.-media(3)*dt/(2*media(1)*eps0))/&
                      (1.+media(3)*dt/(2*media(1)*eps0))
cb(i,j)=(dt/(media(1)*eps0))/(1.+media(3)*dt/(2*media(1)*eps0))
cr(i,j)=(1.-media(4)*dt/(2*media(2)*mu0))/&
                      (1.+media(4)*dt/(2*media(2)*mu0))
cq(i,j)=(dt/(media(2)*mu0))/(1.+media(4)*dt/(2*media(2)*mu0))
end do
end do

outflag='i'

!**************************************************************************

n=0
do while(n<=timestop)             !主循环
n=n+1
if(mod(n,100)==0)then
write(*,*) n
end if

do i=istart+1,iend
ein(i)=ein(i)+0.5*(hin(i-1)-hin(i))
end do

do i=istart,iend-1
hin(i)=hin(i)-0.5*(ein(i+1)-ein(i))
end do

do i=istart,iend                    !产生入射波
hin(i)=sin(omiga*n*dt)            
end do

!* * * * * * * * * * * *连接边界入射场分量 * * * * * * * * * * * * * * 

do i=itmin,itmax
t1=float(i)*cp+float(jtmin)*sp
ii=int(t1)
if(t1>=0.) then
iii=ii+1
t1=float(iii)-t1
else
iii=ii-1
t1=t1-float(iii)
end if
hcb(i,1)=t1*(hin(ii)-hin(iii))+hin(iii)      !底部连接边界入射磁场分量      
                                                                        
t1=float(i)*cp+float(jtmax)*sp
ii=int(t1)
if(t1>=0.) then
iii=ii+1
t1=float(iii)-t1
else
iii=ii-1
t1=t1-float(iii)
end if
hcb(i,3)=t1*(hin(ii)-hin(iii))+hin(iii)    
end do                                       !顶部连接边界入射磁场分量

do j=jtmin,jtmax
t1=float(j)*sp+float(itmin)*cp
ii=int(t1)
if(t1>=0.) then
iii=ii+1
t1=float(iii)-t1
else
iii=ii-1
t1=t1-float(iii)
end if
hcb(j,4)=t1*(hin(ii)-hin(iii))+hin(iii)      !左边连接边界入射磁场分量

t1=float(j)*sp+float(itmax)*cp
ii=int(t1)
if(t1>=0.) then
iii=ii+1
t1=float(iii)-t1
else
iii=ii-1
t1=t1-float(iii)
end if
hcb(j,2)=t1*(hin(ii)-hin(iii))+hin(iii)     
end do                                       !右边连接边界入射磁场分量

do i=itmin,itmax
t1=float(i)*cp+(float(jtmin)-0.5)*sp
if(t1>.5) then 
ii=int(t1-0.5)
iii=ii+1
t1=float(iii)+.5-t1
else
ii=int(t1+.5)
iii=ii-1
t1=t1-.5-float(iii)
end if
ecb(i,1)=t1*(ein(ii)-ein(iii))+ein(iii)
ecb(i,1)=ecb(i,1)*sp                         !底部连接边界入射磁场分量(ex)

t1=float(i)*cp+(float(jtmax)+.5)*sp
if(t1>.5) then
ii=int(t1-.5)
iii=ii+1
t1=float(iii)+.5-t1
else
ii=int(t1+.5)
iii=ii-1
t1=t1-.5-float(iii)
end if
ecb(i,3)=t1*(ein(ii)-ein(iii))+ein(iii)
ecb(i,3)=ecb(i,3)*sp                        !顶部连接边界入射磁场分量(ex)
end do

do j=jtmin,jtmax
t1=float(j)*sp+(float(itmin)-.5)*cp
if(t1>.5) then
ii=int(t1-.5)
iii=ii+1
t1=float(iii)+.5-t1
else
ii=int(t1+.5)
iii=ii-1
t1=t1-.5-float(iii)
end if
ecb(j,4)=t1*(ein(ii)-ein(iii))+ein(iii)
ecb(j,4)=-ecb(j,4)*cp                      !左边连接边界入射磁场分量(ey)

t1=float(j)*sp+(float(itmax)+.5)*cp
if(t1>.5) then
ii=int(t1-.5)
iii=ii+1
t1=float(iii)+.5-t1
else
ii=int(t1+.5)
iii=ii-1
t1=t1-.5-float(iii)
end if
ecb(j,2)=t1*(ein(ii)-ein(iii))+ein(iii)
ecb(j,2)=-ecb(j,2)*cp                      !右边连接边界入射磁场分量(ey)
end do

!* * * * * * * * * * * * * HZ的FDTD迭代 * * * * * * * * * * * * * *

do i=imin+1,imax-1
do j=jmin+1,jmax-1
hz(i,j)=cr(i,j)*hz(i,j)-cq(i,j)*((ey(i,j)-ey(i-1,j))/dd-(ex(i,j)-ex(i,j-1))/dd)
end do
end do

!* * * * * * * * * * 连接边界引入入射波Hz分量 * * * * * * * * * * * *

do i=itmin+1,itmax-1
hz(i,jtmax)=cr(i,jtmax)*hz(i,jtmax)-cq(i,jtmax)*(ey(i,jtmax)-ey(i-1,jtmax)-ex(i,jtmax)+ex(i,jtmax-1)-ecb(i,3))/dd     !上边界
hz(i,jtmin)=cr(i,jtmin)*hz(i,jtmin)-cq(i,jtmin)*(ey(i,jtmin)-ey(i-1,jtmin)-ex(i,jtmin)+ex(i,jtmin-1)+ecb(i,1))/dx     !下边界
end do

do j=jtmin+1,jtmax-1
hz(itmax,j)=cr(itmax,j)*hz(itmax,j)-cq(itmax,j)*(ey(itmax,j)-ey(itmax-1,j)+ecb(j,2)-ex(itmax,j)+ex(itmax,j-1))/dy     !右边界
hz(itmin,j)=cr(itmin,j)*hz(itmin,j)-cq(itmin,j)*(ey(itmin,j)-ey(itmin-1,j)-ecb(j,4)-ex(itmin,j)+ex(itmin,j-1))/dy     !左边界
end do

!* * * * * * * * * * 连接边界角点引入入射波 * * * * * * * * * * * * *

!t1=hcb(jtmax,2)-hcb(itmax,3)

!ez(itmax,jtmax)=ez(itmax,jtmax)+.5*t1

hz(itmin,jtmin)=cr(itmin,jtmin)*hz(itmin,jtmin)-cq(itmin,jtmin)*(ey(itmin,jtmin)-ey(itmin-1,jtmin)-ex(itmin,jtmin)+ex(itmin,jtmin-1))/dx&
                -cq(itmin,jtmin)*(ecb(itmin,1)-ecb(jtmin,4))/dx

!t1=-hcb(jtmin,4)+hcb(itmin,1)
!ez(itmin,jtmin)=ez(itmin,jtmin)+.5*t1

hz(itmax,jtmax)=cr(itmax,jtmax)*hz(itmax,jtmax)-cq(itmax,jtmax)*(ey(itmax,jtmax)-ey(itmax-1,jtmax)-ex(itmax,jtmax)+ex(itmax,jtmax-1))/dx&
                -cq(itmax,jtmax)*(ecb(jtmax,2)-ecb(itmax,3))/dx

!t1=hcb(jtmin,2)+hcb(itmax,1)
!ez(itmax,jtmin)=ez(itmax,jtmin)+.5*t1

hz(itmin,jtmax)=cr(itmin,jtmax)*hz(itmin,jtmax)-cq(itmin,jtmax)*(ey(itmin,jtmax)-ey(itmin-1,jtmax)-ex(itmin,jtmax)+ex(itmin,jtmax-1))/dy&
                -cq(itmin,jtmax)*(-ecb(itmin,3)-ecb(jtmax,4))/dy

!t1=-hcb(itmin,3)-hcb(jtmax,4)
!ez(itmin,jtmax)=ez(itmin,jtmax)+.5*t1

hz(itmax,jtmin)=cr(itmax,jtmin)*hz(itmax,jtmin)-cq(itmax,jtmin)*(ey(itmax,jtmin)-ey(itmax-1,jtmin)-ex(itmax,jtmin)+ex(itmax,jtmin-1))/dy&
                -cq(itmax,jtmin)*(ecb(itmax,1)+ecb(jtmin,2))/dy

!* * * * * * * * * * 吸收边界条件Hz * * * * * * * * * * * *

do i=imin+1,imax-1                         
hz(i,jmin)=hz1(i,jmin+1)+(v*dt-dd)/(v*dt+dd)*(hz(i,jmin+1)-hz(i,jmin))-v*v*eps0*dt/(2*(v*dt+dd))&
           *(ey(i,jmin+1)-ey(i-1,jmin+1)+ey(i,jmin)-ey(i-1,jmin))           !下边界
hz1(i,jmin+1)=hz(i,jmin+1)
hz(i,jmax)=hz3(i,jmax-1)+(v*dt-dd)/(v*dt+dd)*(hz(i,jmax-1)-hz(i,jmax))-v*v*eps0*dt/(2*(v*dt+dd))&
           *(ey(i,jmax-1)-ey(i-1,jmax-1)+ey(i,jmax)-ey(i-1,jmax))           !上边界
hz3(i,jmax-1)=hz(i,jmax-1)
end do

do j=jmin+1,jmax-1                          
hz(imin,j)=hz4(imin+1,j)+(v*dt-dd)/(v*dt+dd)*(hz(imin+1,j)-hz(imin,j))+v*v*eps0*dt/(2*(v*dt+dd))&
           *(ex(imin,j)-ex(imin,j-1)+ex(imin+1,j)-ex(imin+1,j-1))           !左边界
hz4(imin+1,j)=hz(imin+1,j)
hz(imax,j)=hz2(imax-1,j)+(v*dt-dd)/(v*dt+dd)*(hz(imax-1,j)-hz(imax,j))+v*v*eps0*dt/(2*(v*dt+dd))&
           *(ex(imax,j)-ex(imax,j-1)+ex(imax-1,j)-ex(imax-1,j-1))           !右边界
hz2(imax-1,j)=hz(imax-1,j)
end do

!* * * * * * * * * * 吸收边界角点Hz * * * * * * * * * * 

hz(imin,jmin)=hac(4,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imin+1,jmin+1)-hz(imin,jmin))
hac(4,1)=hz(imin+1,jmin+1)

hz(imax,jmax)=hac(2,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imax-1,jmax-1)-hz(imax,jmax))
hac(2,1)=hz(imax-1,jmax-1)

hz(imin,jmax)=hac(3,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imin+1,jmax-1)-hz(imin,jmax))
hac(3,1)=hz(imin+1,jmax-1)

hz(imax,jmin)=hac(1,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imax-1,jmin+1)-hz(imax,jmin))
hac(1,1)=hz(imax-1,jmin+1)

!***************************Ex,Ey的FDTD迭代**************************

do i=imin,imax
do j=jmin,jmax-1
ex(i,j)=ca(i,j)*ex(i,j)+cb(i,j)*(hz(i,j+1)-hz(i,j))/dy
end do
end do

do i=imin,imax-1
do j=jmin,jmax
ey(i,j)=ca(i,j)*ey(i,j)-cb(i,j)*(hz(i+1,j)-hz(i,j))/dx
end do
end do

!* * * * * * * * * *　连接边界引入入射波Ex,Ey分量 * * * * * * * * * *

do i=itmin,itmax
ex(i,jtmin-1)=ca(i,jtmin-1)*ex(i,jtmin-1)+cb(i,jtmin-1)*(hz(i,jtmin)-hz(i,jtmin-1)-hcb(i,1))/dy
ex(i,jtmax)=ca(i,jtmax)*ex(i,jtmax)+cb(i,jtmax)*(hz(i,jtmax+1)-hz(i,jtmax)+hcb(i,3))/dy
end do
 
do j=jtmin,jtmax
ey(itmin-1,j)=ca(itmin-1,j)*ey(itmin-1,j)-cb(itmin-1,j)*(hz(itmin,j)-hz(itmin-1,j)-hcb(j,4))/dx
ey(itmax,j)=ca(itmax,j)*ey(itmax,j)-cb(itmax,j)*(hz(itmax+1,j)-hz(itmax,j)+hcb(j,2))/dx
end do

end do                                      !时间步主循环结束



open(1,file='ex.txt',access='sequential')
open(2,file='ey.txt',access='sequential')

do i=imin,imax,1
do j=jmin,jmax-1,1
write(1,*)ex(i,j)
end do
end do

do i=imin,imax-1,1
do j=jmin,jmax,1
write(2,*)ey(i,j)
end do
end do

close(1)
close(2)

stop
end