Locally-Hosted Collector/ diffraction of plane wave on thin perfectly conducting half-plane

diffraction of plane wave on thin perfectly conducting half-plane

  1. diffraction of plane wave on thin perfectly conducting half-plane
    1. The program
    2. See also

The program

! Copyright (C) 2010 Vsevolod Scherbinin
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
! 02110-1301 USA


!===================================================================
!   MODULE FOR FIELD COMPUTING
!===================================================================
Module Vars
    !Fract      Number of fraction for numerical integration
    !SinTh, CosTh   Incident angle sine and cosine
    !Kappa      Parameter of contour shifting
    !Xi     Dimensionless X coordinate
    !Zeta       Dimensionless Z coordinate
    !Kappa      Integration parameter
    !Hlim       Integration limit
    Integer:: Fract
    Double Precision:: SinTh, CosTh, Mult, Kappa, Xi, Zeta, Hlim

Contains
!===================================================================
!   BASIC SUBROUTINE FOR FIELD AMPLITUDES COMPUTING
!===================================================================
Subroutine HalfPlaneField(cHy, cSf)
    Use Integrals
    Implicit None
    !Sx, Sz     Poynting vector components
    Double Precision:: Sx, Sz
    !cHy        Magnetic field
    !cSf        Poynting vwctor
    !Ecur       Electric field
    !Iw     Incident wave
    Double Complex:: cHy, cSf, Intg, Iw, Ecur
    !VecIntg        Integration result
    Double Complex, Dimension(3):: VecIntg
    Interface 
        Function HfPn(X, N)
            Integer:: N
            Double Precision:: X 
            Double Complex, Dimension(N):: HfPn
        End Function HfPn
    End interface
    Iw = Exp((0._8, -6.28318530718_8)*(Xi*SinTh + Zeta*CosTh))
    Call IntGauss16(-Hlim, Hlim, Fract, VecIntg, HfPn, 3)
    cHy = Iw - Dcmplx(0._8, Sign(Mult,Zeta))*VecIntg(1)
    Ecur = CosTh*Iw + Dcmplx(0._8, Mult)*VecIntg(2)
    Sz = Real(Ecur*Conjg(cHy))
    Ecur = -SinTh*Iw + Dcmplx(0._8, Sign(Mult,Zeta))*VecIntg(3)
    Sx = Real(Ecur*Conjg(cHy))
    cSf = -Dcmplx(Sx, Sz)
End Subroutine HalfPlaneField

End Module Vars

!===================================================================
!   ROUTINE FOR INTEGRAND COMPUTING
!===================================================================
Function HfPn(X, N)
  Use Vars
  Implicit None
  !N        Result array size
  Integer:: N
  !X        Integration variable
  Double Precision:: X
  !HfPn     Integration result
  !Alph     Wavenumber by X
  !W        Wavenumber by Z
  !Exx      Exponent in integrand numerator
  !Den      Integrand denominator
  Double Complex:: Alph, W, Exx, Den
  Double Complex, Dimension(N):: HfPn
  Alph = Dcmplx(X, kappa)
  W = Sqrt((1._8,0._8) - Alph*Alph)
  W = Dcmplx(Real(W), Abs(Aimag(W)))
  Den = (Alph - SinTh)*Sqrt((1._8,0._8) - Alph)
  Exx = Exp((0._8,6.28318530718_8)*(Abs(Zeta)*W - Xi*Alph))
  HfPn(1) = Exx/Den
  HfPn(2) = W*HfPn(1)
  HfPn(3) = Alph*HfPn(1)
End Function HfPn


!===================================================================
!   MAIN PROGRAM
!===================================================================
Program HalfPlane
    Use Vars
    Use Service
    Implicit None
    Include 'mpif.h'
    !===========================================================
    !   Variables for MPI managing
    !===========================================================
    !MpiNum     Number of processes
    !MpiRank    Process rank (identifier)
    !MpiTag     Message tag (number)
    !MpiErr     Error status
    Integer:: MpiNum, MpiRank, MpiTag, MpiErr, CheckPoint = 1
    Integer:: STATUS(MPI_STATUS_SIZE)
    !MpiId      Identifier
    Integer::   MpiId
    !Xbeg, Xend Begin and end of band for computing
    !Xstp       Step for X range splitting
    !Xrng       Quantity of X points for local process
    Integer:: Xbeg, Xend, Xstp, Xrng, ArIndx, tSndBuf, tRcvBuf
    Integer:: K0, KA, IntPrm(3), Shp(2), Sshp(1)
    Double Precision:: DprPrm(7)
    !===========================================================
    !   Input and output variables
    !===========================================================
    !Xin, Zin   Numbers of points for field computing
    !Theta      Incident angle (degrees) counted out z axis
    !Xmax, Zmax Maximum values of X and Z in wavelength
    Integer:: Xin, Zin
    Double Precision:: Theta, Xmax, Zmax
    !===========================================================
    !   Internal variables
    !===========================================================
    !IOS1       Input/output status
    !Xcnt, Zcnt Counters for coordinates
    Integer:: IOS1, Xcnt, Zcnt
    !Phas       Phase of magnetic field
    !Xis, Zis   Step by X and Z variable
    Double Precision:: Phas, Xis, Zis
    !Hycurr, Sfcurr Current field and Poynting vector
    Double Complex:: Hycurr, Sfcurr
    !Xdts, Zdts X and Z coordinates
    !Hy, Sf     Arrays for field and flow
    Double Precision, Allocatable, Dimension(:):: Xdts, Zdts
    Double Complex, Allocatable, Dimension(:,:):: Hy, Sf
    Double Complex, Allocatable, Dimension(:):: locHy, locSf
    !===========================================================
    !   String variables
    !===========================================================
    !Date       System date
    !Stime, Ftime   Start and finish time of job
    !Prefix     Output filenames prefix
    !argbuffer  Buffer for command line options
    !configfile File of initial data
    !basename   Base of output files names
    Character(len=8):: Date
    Character(len=10):: STime, Ftime, Prefix
    Character(len=40):: argbuffer, configfile, basename
    Character(len=40):: outfilea, outfilep, outfiles
    !===========================================================
    !   Namelists and interfaces
    !===========================================================
    Namelist /params/ Theta, Xin, Zin, Xmax, Zmax
    !===========================================================
    !   COMMAND LINE ARGUMENTS PROCESSING
    !===========================================================
    Call Getarg(1,argbuffer)
    Read(argbuffer,*,IOSTAT=IOS1) configfile
    If (IOS1.LT.0) configfile='halfplane.ini'
    !===========================================================
    !   FILES OPENING AND VARIABLES INITIALIZING
    !===========================================================
    Call MPI_INIT(MpiErr)
    Call MPI_COMM_SIZE(MPI_COMM_WORLD, MpiNum, MpiErr)
    Call MPI_COMM_RANK(MPI_COMM_WORLD, MpiRank, MpiErr)
!===================================================================
!   PARAMETERS READING, INITIALIZING (BY 0 PROCESS)
!===================================================================
3   Continue
    If (MpiRank.EQ.0) Then
        Xin = 601
        Zin = 600
        Xmax = 4._8
        Zmax = 4._8
        Open(Unit=3, File='current.log')
        Open(Unit=1, File=Trim(configfile))
2       Read(1,'(A)',END=1) basename
        If ((Scan(Trim(basename),'#').EQ.1).OR.((Scan(Trim(basename),';').EQ.1))) Then
          Write(*,*) 'Comment found. Program will skip the string.'
          Goto 2
        EndIf
        If (Trim(basename).EQ.'stop') Goto 1
        Call Date_And_Time(Date,STime)
        Write(*,*) Trim(basename),' starts at ',Stime
        Read(1,NML=params)
        !===================================================
        !Data for all processes preparing
        !===================================================
        !Higher integration limit
        Hlim = 4._8
        !Number of fraction for numerical integrating
        Fract = Int(2._8*Hlim)*5
        !Incident angle processing
        Theta = Theta*0.0174532925199_8 !Incident angle to radian
        SinTh = Sin(Theta)
        CosTh = Cos(Theta)
        !Multiplyer 1D-1 cannot be minimized, cause pole set on integration contour
        Kappa = (SinTh-1._8)*1D-1
        !Zin must be positive, even and not less than 2
        Zin = 2*Int(Abs(Zin)/2); If (Zin.LE.2) Zin = 2
        !Xin must positive, odd and not less than 3
        Xin = 2*Int(Abs(Xin)/2)+1; If (Xin.LE.3) Xin = 3
        !Xmax and Zmax must bee positive
        Xmax = Abs(Xmax)
        Zmax = Abs(Zmax)
        !Normalizing multiplier
        Mult = CosTh/(6.28318530718_8*Sqrt(1._8+SinTh))
        !===================================================
        !Local data for 0 process
        !===================================================
        !File names preparing. Prefix needs to determine area for plotting.
        !Files for ampl[itude], phas[e] and [power ]flow.
        Prefix='m'//Trim(NumToString(Int(Max(Xmax,Zmax))))//'_'
        outfilea=Trim(Prefix)//Trim(basename)//'_ampl.dat'
        outfilep=Trim(Prefix)//Trim(basename)//'_phas.dat'
        outfiles=Trim(Prefix)//Trim(basename)//'_flow.dat'
        !===================================================
        !   Preparing broadcasting
        !===================================================
        IntPrm(1) = Fract; IntPrm(2) = Zin; IntPrm(3) = Xin
        DprPrm(1) = Xmax; DprPrm(2) = Zmax; DprPrm(3) = Hlim
        DprPrm(4) = SinTh; DprPrm(5) = CosTh
        DprPrm(6) = Kappa; DprPrm(7) = Mult
        !===================================================
        !Array for field making
        !===================================================
        Allocate(Hy(Xin, Zin))
        Allocate(Sf(Xin, Zin))
    EndIf 
!===================================================================
!   PARAMETERS BROADCASTING (EXCHANGE BETWEEN ALL PROCESSES)
!===================================================================
    Call MPI_BCAST(IntPrm, 3, MPI_Integer, 0, MPI_COMM_WORLD, MpiErr)
    If ((Sum(IntPrm).EQ.0).AND.(MpiRank.NE.0)) Goto 4
    Call MPI_BCAST(DprPrm, 7, MPI_Double_Precision, 0, MPI_COMM_WORLD, MpiErr)
    Fract = IntPrm(1); Zin = IntPrm(2); Xin = IntPrm(3)
    Xmax = DprPrm(1); Zmax = DprPrm(2); Hlim = DprPrm(3); SinTh = DprPrm(4)
    CosTh = DprPrm(5); Kappa = DprPrm(6); Mult = DprPrm(7)
!===================================================================
!   ARRAYS FOR DOTS AND LOCAL ARRAYS FOR FIELDS PREPARING (FOR ALL PROCESSES)
!===================================================================
    !Limitation for X coordinate computing
    KA = Xin/MpiNum
    K0 = Xin-KA*(MpiNum-1)
    If (MpiRank.Eq.0) Then
        Xbeg = 1
        Xend = K0
    Else
        Xbeg = K0 + KA*(MpiRank-1) + 1
        Xend = K0 + KA*MpiRank
    EndIf
    !===================================================
    !   Buffer arrays making
    !===================================================
    Allocate(locHy(KA*Zin))
    Allocate(locSf(KA*Zin))
    Allocate(Xdts(Xin))
    Allocate(Zdts(Zin))
    ! ARRAYS OF DOTS FILLING
    ! May be distributed between processes but I think that time for exchanging
    ! same to computing time. This computation carried out by every process.
    Xis = Xmax/Dble((Xin-1)/2)
    Do Xcnt = 1, Xin
        Xdts(Xcnt) = - Xmax + (Xcnt-1)*Xis
    EndDo
    Zis = Zmax/Dble(Zin/2)
    Do Zcnt = 1, Zin/2
        Zdts(Zcnt) = - Zmax + (Zcnt-1)*Zis
    EndDo
    Do Zcnt = Zin/2+1, Zin
        Zdts(Zcnt) = - Zmax + Zcnt*Zis
    EndDo
!===================================================================
!   COMPUTING (FOR ALL PROCESSES)
!===================================================================
    !Computing itself
    If (MpiRank.Eq.0) Then
        Do Xcnt = 1, K0
        Xi = Xdts(Xcnt)
        Do Zcnt = 1, Zin
            Zeta = Zdts(Zcnt)
            Call HalfPlaneField(Hy(Xcnt, Zcnt),Sf(Xcnt, Zcnt))
        EndDo
        EndDo
    Else
        Do Xcnt = 1, KA
        Xi = Xdts(Xcnt + Xbeg -1)
        Do Zcnt = 1, Zin
            Zeta = Zdts(Zcnt)
            ArIndx = (Xcnt-1)*Zin + Zcnt
            Call HalfPlaneField(locHy(ArIndx), locSf(ArIndx))
        EndDo
        EndDo
    EndIf
    Shp(2) = KA
    Shp(1) = Zin
    Call MPI_BARRIER(MPI_COMM_WORLD, MpiErr)
!===================================================================
!   DATA GATHERING (FOR ALL PROCESSES)
!===================================================================
    If (MpiRank.EQ.0) Then
        Do MpiId = 1, MpiNum-1
            Xbeg = K0 + KA*(MpiId-1) + 1
            Xend = K0 + KA*MpiId
            Call MPI_RECV(locHy, KA*Zin, MPI_Double_Complex, MpiId, 1, MPI_COMM_WORLD, Status, MpiErr)
            Hy(Xbeg:Xend,:) = Transpose(Reshape(locHy, Shp))
            Call MPI_RECV(locSf, KA*Zin, MPI_Double_Complex, MpiId, 2, MPI_COMM_WORLD, Status, MpiErr)
            Sf(Xbeg:Xend,:) = Transpose(Reshape(locSf, Shp))
        EndDo
    Else
        Call MPI_SSEND(locHy, KA*Zin, MPI_Double_Complex, 0, 1, MPI_COMM_WORLD, MpiErr)
        Call MPI_SSEND(locSf, KA*Zin, MPI_Double_Complex, 0, 2, MPI_COMM_WORLD, MpiErr)
    EndIf
!===================================================================
!   DATA OUTPUT (FOR MAIN PROCESSES WITH 0 RANK)
!===================================================================
    If (MpiRank.EQ.0) Then
    !====================================
    !Field output starts
!       Write(*,*) Trim(outfilea), Trim(outfilep), Trim(outfiles)
        Open(Unit=2, file=Trim(outfilea))
        Open(Unit=4, file=Trim(outfilep))
        Open(Unit=5, file=Trim(outfiles))
        Do Xcnt = 1, Xin
            Do Zcnt = 1, Zin
            Write(2,'(2(F6.2,1X),F8.4)') Xdts(Xin+1-Xcnt), Zdts(Zcnt), Abs(Hy(Xcnt,Zcnt))
            Phas = PolarAngle(Real(Hy(Xcnt,Zcnt)), Aimag(Hy(Xcnt,Zcnt)))*57.2957795131
            Write(4,'(2(F6.2,1X),F8.4)') Xdts(Xin+1-Xcnt), Zdts(Zcnt), Phas
            Sf(Xcnt,Zcnt) = 0.05_8*Sf(Xcnt,Zcnt)/Abs(Sf(Xcnt,Zcnt)) !Poyting vector normalizing
            Write(5,'(2(F6.2,1X),2F9.4)') Xdts(Xin+1-Xcnt), Zdts(Zcnt), Real(Sf(Xcnt,Zcnt)), Aimag(Sf(Xcnt,Zcnt))
            EndDo
            Write(2,*) ''
            Write(4,*) ''
            Write(5,*) ''
        EndDo
        Write(*,*) 'Output finished '
        Close(5)
        Close(4)
        Close(2)
        Write(*,*) 'Files closing'
    !==============================================================================
    !Finalizing
    !==============================================================================
    !Array clearing
        If (Allocated(Hy)) Deallocate(Hy)
        If (Allocated(Sf)) Deallocate(Sf)
        If (Allocated(locHy)) Deallocate(locHy)
        If (Allocated(locSf)) Deallocate(locSf)
        If (Allocated(Xdts)) Deallocate(Xdts)
        If (Allocated(Zdts)) Deallocate(Zdts)
        Call Date_And_Time(Date,FTime)
        Write(*,*) 'Task ',Trim(basename),' has been computed during ',TimeInterval(Stime,Ftime)
        Write(3,*) 'Task ',Trim(basename),' has been computed during ',TimeInterval(Stime,Ftime)
        Goto 2
1       Close(1)
        Close(3)
    Else
        If (Allocated(locHy)) Deallocate(locHy)
        If (Allocated(locSf)) Deallocate(locSf)
        If (Allocated(Xdts)) Deallocate(Xdts)
        If (Allocated(Zdts)) Deallocate(Zdts)
        Goto 3
    EndIf
    IntPrm = 0
    Call MPI_BCAST(IntPrm, 3, MPI_Integer, 0, MPI_COMM_WORLD, MpiErr)
4   Call MPI_FINALIZE(MpiErr)
    Stop
End Program

See also

A non-vectorized, non-multiprocessing simple C version.