myfdtd_1d_parallel.f90

一个Fortran的1D-FDTD并行程序

program main

    include 'mpif.h'
    
    real,parameter::c=3.0e8
    real,parameter::dx=0.005
    real,parameter::N=1.0
    integer,parameter::Nx0=int(N/dx)
    real,parameter::dt=(dx/(2.0*c))
    real,parameter::T=1.0e-9
    integer,parameter::Nt=320
    
    real,allocatable::Ey(:)
    
    real,allocatable::node_Ey(:)
    real,allocatable::node_Hz(:)
    
    real node_Ey_right
    real node_Hz_left
    
    real node_Eyoldleft
    real node_Eyoldright
    
    integer source_id,ksn
    
    integer,allocatable::displs(:),recvcounts(:)
    integer num_Ey
    
    integer i,nstep
    integer nx,nxtemp
    logical flag
    
    real w
    real t0
    real ks
    real pulse
    
    integer processorID,ieer
    integer numberofprocessors
    integer status(MPI_Status_size)
    integer irequest
    integer,parameter::main_id=0
    integer,parameter::tagnx=1,tagEy=2,tagHz=3
    integer,parameter::SendType=2
    
    call MPI_Init(ieer)
    call MPI_Comm_size(MPI_Comm_world,numberofprocessors,ieer);
    call MPI_Comm_rank(MPI_Comm_world,processorID,ieer);
    
    w=12.0
    t0=40.0
    ks=Nx0/2
    
!*******************************Main Process********************************
    
    if(processorID==main_id)then
        write(*,*)'Number of time steps=',Nt
        write(*,*)'The problem will be calculated by',numberofprocessors,'processors'
        
        nx=int(Nx0/numberofprocessors)
        write(*,*)'The main processor',processorID,'nx=',nx
        
        if(ks>=nx)then
            flag=.true.
        else
            source_id=main_id
            ksn=ks
            flag=.false.
        end if
        
        allocate(displs(0:numberofprocessors-1))
        allocate(recvcounts(0:numberofprocessors-1))
        displs(:)=0
        recvcounts(:)=0
        recvcounts(0)=nx
        
        do i=1,numberofprocessors-1
            nxtemp=int(Nx0*(i+1)/numberofprocessors)-int(Nx0*i/numberofprocessors)
            
            displs(i)=int(Nx0*i/numberofprocessors)
            recvcounts(i)=nxtemp
           
    !***************************the number of Ey is one more than the number of Hz in the last process******************************************       
            if(i==numberofprocessors-1)then
                recvcounts(numberofprocessors-1)=nxtemp+1
            end if
      
            call MPI_Send(nxtemp,1,MPI_integer,i,tagnx,MPI_Comm_world,ieer)
            
            if(flag)then
                if(ks<int(Nx0*(i+1)/numberofprocessors))then
                    source_id=i
                    ksn=ks-int(Nx0*i/numberofprocessors)
                    flag=.false.
                end if
            end if
        end do
        
        write(*,*)'The source will be calculated in processor',source_id
        write(*,*)'The source position in this processor is',ksn

!*************************************slave process************************************************

    else
        call MPI_Recv(nx,1,MPI_integer,main_id,tagnx,MPI_Comm_world,status,ieer)
        
        write(*,*)'The processor',processorID,':nx=',nx
    end if
    
    
    call MPI_BCAST(source_id,1,MPI_integer,main_id,MPI_Comm_world,ieer)
    call MPI_BCAST(ksn,1,MPI_integer,main_id,MPI_Comm_world,ieer)
    
    
!***********************************Start FDTD***************************************************

    allocate(Ey(0:Nx0))
    allocate(node_Ey(0:nx))
    allocate(node_Hz(0:nx-1))
    
    Ey(:)=0.0d0
    node_Ey(:)=0.0d0
    node_Hz(:)=0.0d0
    node_Ey_right=0.0d0
    node_Hz_left=0.0d0
    node_Eyoldleft=0.0d0
    node_Eyoldright=0.0d0
    
    do nstep=0,Nt
    
!****************************************************H FDTD********************************************

        if(processorID==numberofprocessors-1)then
            do i=0,nx-1
                node_Hz(i)=node_Hz(i)+0.5*(node_Ey(i)-node_Ey(i+1))
            end do
        else
            do i=0,(nx-1)-1
                node_Hz(i)=node_Hz(i)+0.5*(node_Ey(i)-node_Ey(i+1))
            end do
            node_Hz(nx-1)=node_Hz(nx-1)+0.5*(node_Ey(nx-1)-node_Ey_right)
        end if
        
        
        if(processorID==main_id)then
            call MPI_Send(node_Hz(nx-1),1,MPI_Real,1,tagHz,MPI_Comm_world,ieer)
        else if(processorID==numberofprocessors-1)then
            call MPI_Recv(node_Hz_left,1,MPI_Real,numberofprocessors-2,tagHz,MPI_Comm_world,status,ieer)
        else
            call MPI_Send(node_Hz(nx-1),1,MPI_Real,processorID+1,tagHz,MPI_Comm_world,ieer)
            call MPI_Recv(node_Hz_left,1,MPI_Real,processorID-1,tagHz,MPI_Comm_world,status,ieer)
        end if
        
!****************************************************E FDTD********************************************
        
        if(processorID==main_id)then
            node_Eyoldleft=node_Ey(1)
            do i=1,nx-1
                node_Ey(i)=node_Ey(i)+0.5*(node_Hz(i-1)-node_Hz(i))
            end do
            
            node_Ey(0)=node_Eyoldleft-(node_Ey(1)-node_Ey(0))/3.0
        else if(processorID==numberofprocessors-1)then
            node_Eyoldright=node_Ey(nx-1)
            node_Ey(0)=node_Ey(0)+0.5*(node_Hz_left-node_Hz(0))
            do i=1,nx-1
                node_Ey(i)=node_Ey(i)+0.5*(node_Hz(i-1)-node_Hz(i))
            end do
            node_Ey(nx)=node_Eyoldright-(node_Ey(nx-1)-node_Ey(nx))/3.0
        else
            node_Ey(0)=node_Ey(0)+0.5*(node_Hz_left-node_Hz(0))
            do i=1,nx-1
                node_Ey(i)=node_Ey(i)+0.5*(node_Hz(i-1)-node_Hz(i))
            end do
        end if
        
!****************************************Pulse*************************************************

        if(processorID==source_id)then
            pulse=exp(-0.5*((t0-nstep)**2/(w*w)))
            node_Ey(ksn)=pulse;
        end if
        
        if(processorID==numberofprocessors-1)then
            call MPI_Send(node_Ey(0),1,MPI_Real,numberofprocessors-2,tagEy,MPI_Comm_world,ieer)
        else if(processorID==main_id)then
            call MPI_Recv(node_Ey_right,1,MPI_Real,1,tagEy,MPI_Comm_world,status,ieer)
        else
            call MPI_Send(node_Ey(0),1,MPI_Real,processorID-1,tagEy,MPI_Comm_world,ieer)
            call MPI_Recv(node_Ey_right,1,MPI_Real,processorID+1,tagEy,MPI_Comm_world,status,ieer)
        end if
        
        
        if(processorID==main_id) write(*,*)nstep
    end do
    
    
!*****************************************Output*************************************************    
  
    call MPI_BARRIER(MPI_Comm_world,ieer)
    num_Ey=nx
    if(processorID==numberofprocessors-1)num_Ey=nx+1
    call MPI_Gatherv(node_Ey,num_Ey,MPI_Real,Ey,recvcounts,displs,MPI_Real,main_id,MPI_Comm_world,ieer)
    
    if(processorID==main_id)then
        open(1,file='Ey.dat')
        do i=0,Nx0
            write(1,*)i,Ey(i)
        end do
        close(1)
    end if
    
    if(processorID==main_id)write(*,*)'Finish'
    call MPI_Finalize(ieer)
        
end program