main_mpi.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

      CALL MPI_RECV(MESHES(NOM)%J_VN,1,MPI_INTEGER, NOM-1,1,MPI_COMM_WORLD,STATUS,IERR)
      CALL MPI_RECV(MESHES(NOM)%K_VN,1,MPI_INTEGER, NOM-1,1,MPI_COMM_WORLD,STATUS,IERR)
      CALL MPI_RECV(MESHES(NOM)%NLP,1,MPI_INTEGER, NOM-1,1,MPI_COMM_WORLD,STATUS,IERR)
   ENDDO
 
ENDIF
 
TUSED(11,:) = TUSED(11,:) + SECOND() - TNOW
END SUBROUTINE EXCHANGE_DIAGNOSTICS


SUBROUTINE CORRECT_PRESSURE(CODE) !!!! Experimental Code, DO NOT USE

REAL(EB), DIMENSION(NMESHES*4,NMESHES*4) :: A,A_LOC
REAL(EB), DIMENSION(NMESHES*4) :: B,B_LOC
TYPE (MESH_TYPE), POINTER :: M
TYPE (OMESH_TYPE), POINTER :: OM,M3
INTEGER :: IERROR,NM,SNODE,RNODE,LL,II,JJ,KK,CODE,I,J,K,IOR
LOGICAL :: FLAG

! Post Receives for arrays containing Pressure boundary info

N_PREQ = 0
NM = MYID+1

OTHER_MESH_LOOP: DO NOM=1,NMESHES

   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)==0) CYCLE OTHER_MESH_LOOP
   IF (.NOT.ACTIVE_MESH(NOM)) CYCLE OTHER_MESH_LOOP

   IF (DEBUG) THEN
      WRITE(0,*) NM,' posting pressure receives from ',NOM,' code=',code
      IF (CODE==0) WRITE(0,'(A,I2,A,I2,A,I5)') 'NIC(',NM,',',NOM,')=',NIC(NM,NOM)
   ENDIF

   M =>MESHES(NM)
   M3=>MESHES(NM)%OMESH(NOM)

   RNODE = NOM-1
   TAG   = TAGS(NM,NOM,CODE)

   N_PREQ = N_PREQ+1
   CALL MPI_IRECV(M3%RPKG3(1),NIC(NM,NOM)*6+1, MPI_DOUBLE_PRECISION,RNODE,TAG,MPI_COMM_WORLD, PREQ(N_PREQ),IERR)

ENDDO OTHER_MESH_LOOP

! Exchange boundary values of H and FVX, FVY, FVZ

SEND_OTHER_MESH_LOOP: DO NOM=1,NMESHES

   IF (NIC(NOM,NM)==0 .AND. NIC(NM,NOM)==0) CYCLE SEND_OTHER_MESH_LOOP

   M =>MESHES(NM)
   M3=>MESHES(NM)%OMESH(NOM)
   M4=>MESHES(NOM)

   SNODE = NOM-1
   RNODE = NM-1

   TAG = TAGS(NM,NOM,CODE)

   LL = 0
   IWW = 0
   PACK_SPKG3: DO IW=1,M4%NEWC
      IF (M3%IJKW(9,IW)/=NM .OR. M3%BOUNDARY_TYPE(IW)/=INTERPOLATED_BOUNDARY) CYCLE PACK_SPKG3
      IWW = IWW + 1
      II = M3%IJKW(10,IW)
      JJ = M3%IJKW(11,IW)
      KK = M3%IJKW(12,IW)
      IOR = M3%IJKW(4,IW)
      M3%SPKG3(LL+1) = REAL(IW,EB)
      M3%SPKG3(LL+2) = M%H(II,JJ,KK)
      SELECT CASE(IOR)
         CASE(-1) 
            M3%SPKG3(LL+3) = M%H(II-1,JJ,KK)
            M3%SPKG3(LL+4) = M%FVX(II-1,JJ,KK)
         CASE( 1) 
            M3%SPKG3(LL+3) = M%H(II+1,JJ,KK)
            M3%SPKG3(LL+4) = M%FVX(II,JJ,KK)
         CASE(-2) 
            M3%SPKG3(LL+3) = M%H(II,JJ-1,KK)
            M3%SPKG3(LL+5) = M%FVY(II,JJ-1,KK)
         CASE( 2) 
            M3%SPKG3(LL+3) = M%H(II,JJ+1,KK)
            M3%SPKG3(LL+5) = M%FVY(II,JJ,KK)
         CASE(-3) 
            M3%SPKG3(LL+3) = M%H(II,JJ,KK-1)
            M3%SPKG3(LL+6) = M%FVZ(II,JJ,KK-1)
         CASE( 3) 
            M3%SPKG3(LL+3) = M%H(II,JJ,KK+1)
            M3%SPKG3(LL+6) = M%FVZ(II,JJ,KK)
      END SELECT
      LL = LL+6
   ENDDO PACK_SPKG3
   M3%SPKG3(IWW*6+1) = -999.0_EB
   N_PREQ=N_PREQ+1
   CALL MPI_ISEND(M3%SPKG3(1),IWW*6+1, MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD, PREQ(N_PREQ),IERR)

   IF (DEBUG) THEN
      WRITE(0,*) NM,' sending P data to ',NOM, ' tag=',TAGS(NM,NOM,CODE),' PREQ=',PREQ(N_PREQ)
   ENDIF

ENDDO SEND_OTHER_MESH_LOOP

CALL MPI_WAITALL(N_PREQ,PREQ(1:N_PREQ),ARRAY_OF_STATUSES2,IERR)

NOM = MYID+1

RECEIVE_OTHER_MESH_LOOP: DO NM=1,NMESHES

IF (NIC(NOM,NM)==0 .AND. NIC(NM,NOM)==0)  CYCLE RECEIVE_OTHER_MESH_LOOP

M =>MESHES(NM)
M2=>MESHES(NOM)%OMESH(NM)
M4=>MESHES(NOM)

SNODE = NOM-1
RNODE = NM-1
TAG = TAGS(NM,NOM,CODE)

LL = 0
UNPACK_RPKG3: DO
   IW = NINT(M2%RPKG3(LL+1))
   IF (IW==-999) EXIT UNPACK_RPKG3
   II = M4%IJKW(10,IW)
   JJ = M4%IJKW(11,IW)
   KK = M4%IJKW(12,IW)
   IOR = M4%IJKW(4,IW)
   M2%H(II,JJ,KK)                = M2%RPKG3(LL+2)
   SELECT CASE(IOR)
      CASE(-1) 
         M2%H(II-1,JJ,KK)   = M2%RPKG3(LL+3)
         M2%FVX(II-1,JJ,KK) = M2%RPKG3(LL+4)
      CASE( 1) 
         M2%H(II+1,JJ,KK)   = M2%RPKG3(LL+3)
         M2%FVX(II,JJ,KK)   = M2%RPKG3(LL+4)
      CASE(-2) 
         M2%H(II,JJ-1,KK)   = M2%RPKG3(LL+3)
         M2%FVY(II,JJ-1,KK) = M2%RPKG3(LL+5)
      CASE( 2) 
         M2%H(II,JJ+1,KK)   = M2%RPKG3(LL+3)
         M2%FVY(II,JJ,KK)   = M2%RPKG3(LL+5)
      CASE(-3) 
         M2%H(II,JJ,KK-1)   = M2%RPKG3(LL+3)
         M2%FVZ(II,JJ,KK-1) = M2%RPKG3(LL+6)
      CASE( 3) 
         M2%H(II,JJ,KK+1)   = M2%RPKG3(LL+3)
         M2%FVZ(II,JJ,KK)   = M2%RPKG3(LL+6)
   END SELECT
   LL = LL+6
ENDDO UNPACK_RPKG3

ENDDO RECEIVE_OTHER_MESH_LOOP

! Construct "correction matrix" A

A = 0.
B = 0.
A_LOC = 0.
B_LOC = 0.

MESH_LOOP_1: DO NM=MYID+1,NMESHES,NUMPROCS
   CALL COMPUTE_A_B(A_LOC,B_LOC,NM)
ENDDO MESH_LOOP_1

CALL MPI_ALLREDUCE(A_LOC(1,1),A(1,1),(NMESHES*4)**2, MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
CALL MPI_ALLREDUCE(B_LOC(1),B(1),(NMESHES*4), MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)

IF (MYID==0) THEN
   DO I=1,4*NMESHES
      WRITE(0,'(16F7.1,3X,F7.1)') (A(I,J),J=1,4*NMESHES),B(I)
   ENDDO
ENDIF

CALL GAUSSJ(A,NMESHES*4,NMESHES*4,B,1,1,IERROR)

IF (MYID==0) THEN
   IF (IERROR>0) WRITE(0,*) ' IERROR= ',IERROR
   WRITE(0,*)
   DO I=1,4*NMESHES
      WRITE(0,'(F12.5)') B(I)
   ENDDO
   WRITE(0,*)
ENDIF

MESH_LOOP_2: DO NM=MYID+1,NMESHES,NUMPROCS
   CALL UPDATE_PRESSURE(B,NM)
ENDDO MESH_LOOP_2

END SUBROUTINE CORRECT_PRESSURE


SUBROUTINE DUMP_GLOBAL_OUTPUTS(T)

! Dump HRR data to CHID_hrr.csv, MASS data to CHID_mass.csv, DEVICE data to _devc.csv

REAL(EB) :: T
INTEGER :: N

! Dump out HRR info  after first "gathering" data to node 0

IF_DUMP_HRR: IF (T>=HRR_CLOCK) THEN
   CALL MPI_ALLGATHER(HRR_COUNT(MYID+1),1,MPI_DOUBLE_PRECISION, HRR_COUNT,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
   IF (MINVAL(HRR_COUNT,MASK=.NOT.EVACUATION_ONLY)>0._EB) THEN
      CALL MPI_GATHER(HRR_SUM(MYID+1), 1, MPI_DOUBLE_PRECISION, HRR_SUM, 1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHER(RHRR_SUM(MYID+1),1, MPI_DOUBLE_PRECISION, RHRR_SUM,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHER(CHRR_SUM(MYID+1),1, MPI_DOUBLE_PRECISION, CHRR_SUM,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHER(FHRR_SUM(MYID+1),1, MPI_DOUBLE_PRECISION, FHRR_SUM,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHER(MLR_SUM(MYID+1), 1, MPI_DOUBLE_PRECISION, MLR_SUM, 1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      IF (MYID==0) CALL DUMP_HRR(T)
      HRR_CLOCK = HRR_CLOCK + DT_HRR
      HRR_SUM   = 0._EB
      RHRR_SUM  = 0._EB
      CHRR_SUM  = 0._EB
      FHRR_SUM  = 0._EB
      MLR_SUM   = 0._EB
      HRR_COUNT = 0._EB
   ENDIF
ENDIF IF_DUMP_HRR

!     Dump out Evac info: EVAC_TODO next lines should be
!     made to work also in the parallel code.
!
!      IF (T>=EVAC_CLOCK .AND. ANY(EVACUATION_GRID)) THEN
!        CALL DUMP_EVAC_CSV(T)
!        EVAC_CLOCK = EVAC_CLOCK + DTHRR
!      ENDIF

! Dump out Mass info after first "gathering" data to node 0

IF_DUMP_MASS: IF (T>=MINT_CLOCK) THEN
   CALL MPI_ALLGATHER(MINT_COUNT(MYID+1),1,MPI_DOUBLE_PRECISION, MINT_COUNT,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
   IF (MINVAL(MINT_COUNT,MASK=.NOT.EVACUATION_ONLY)>0.) THEN
      CALL MPI_GATHER(MINT_SUM(0,MYID+1),21,MPI_DOUBLE_PRECISION, MINT_SUM,21,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD, IERR)
      IF (MYID==0) CALL DUMP_MASS(T)
      MINT_CLOCK    = MINT_CLOCK + DT_MASS
      MINT_SUM   = 0._EB
      MINT_COUNT = 0._EB
   ENDIF
ENDIF IF_DUMP_MASS

! Exchange DEVICE parameters among meshes and dump out DEVICE info after first "gathering" data to node 0
 
IF_DUMP_DEVC: IF (T>=DEVC_CLOCK .AND. N_DEVC>0) THEN
  
   ! Exchange the CURRENT_STATE of each DEViCe

   STATE_LOC(1:N_DEVC) = .FALSE.  ! _LOC is a temporary array that holds the STATE value for the devices on each node
   DO N=1,N_DEVC
      IF (DEVICE(N)%MESH==MYID+1) STATE_LOC(N) = DEVICE(N)%CURRENT_STATE 
   ENDDO
   CALL MPI_ALLREDUCE(STATE_LOC(1),STATE_GLB(1),N_DEVC,MPI_LOGICAL,MPI_LXOR,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%CURRENT_STATE = STATE_GLB(1:N_DEVC)

   ! Exchange the INSTANT_VALUE of each DEViCe

   TC_LOC(1:N_DEVC) = 0._EB 
   DO N=1,N_DEVC
      IF (DEVICE(N)%MESH==MYID+1) TC_LOC(N) = DEVICE(N)%INSTANT_VALUE
   ENDDO
   CALL MPI_ALLREDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%INSTANT_VALUE = TC_GLB(1:N_DEVC)

   ! Exchange the INSTANT_VALUE of each DEViCe

   TC_LOC(1:N_DEVC) = 0._EB
   DO N=1,N_DEVC
      IF (DEVICE(N)%MESH==MYID+1) TC_LOC(N) = DEVICE(N)%T_CHANGE
   ENDDO
   CALL MPI_ALLREDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%T_CHANGE = TC_GLB(1:N_DEVC)

   ! Exchange the current COUNT of each DEViCe

   COUNT_LOC(1:N_DEVC) = DEVICE(1:N_DEVC)%COUNT
   CALL MPI_ALLREDUCE(COUNT_LOC(1),COUNT_GLB(1),N_DEVC,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR)

   ! Get the current VALUEs of all DEViCes into DEVICE(:)%VALUE on node 0

   IF (MINVAL(COUNT_GLB)>0) THEN
      TC_LOC(1:N_DEVC) = DEVICE(1:N_DEVC)%VALUE 
      CALL MPI_REDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,IERR)
      IF (MYID==0) THEN
         DEVICE(1:N_DEVC)%VALUE = TC_GLB(1:N_DEVC)
         DEVICE(1:N_DEVC)%COUNT = COUNT_GLB(1:N_DEVC)
         CALL DUMP_DEVICES(T)
      ENDIF
      DEVC_CLOCK = DEVC_CLOCK + DT_DEVC
      DEVICE(1:N_DEVC)%VALUE = 0._EB
      DEVICE(1:N_DEVC)%COUNT = 0
   ENDIF
ENDIF IF_DUMP_DEVC

! Dump CONTROL info. No gathering required as CONTROL is updated on all meshes

IF (T>=CTRL_CLOCK .AND. N_CTRL>0) THEN
   IF (MYID==0) CALL DUMP_CONTROLS(T)
   CTRL_CLOCK = CTRL_CLOCK + DT_CTRL
ENDIF

END SUBROUTINE DUMP_GLOBAL_OUTPUTS

SUBROUTINE GET_REVISION_NUMBER(REV_NUMBER,REV_DATE)
USE isodefs, ONLY : GET_REV_smvv
USE POIS, ONLY : GET_REV_pois
USE COMP_FUNCTIONS, ONLY : GET_REV_func
USE MESH_POINTERS, ONLY : GET_REV_mesh
USE RADCALV, ONLY : GET_REV_irad
USE DCDFLIB, ONLY : GET_REV_ieva
INTEGER,INTENT(INOUT) :: REV_NUMBER
CHARACTER(255),INTENT(INOUT) :: REV_DATE
INTEGER :: MODULE_REV
CHARACTER(255) :: MODULE_DATE

CALL GET_REV_cons(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_ctrl(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_devc(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_divg(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_dump(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
!EVAC:CALL GET_REV_evac(MODULE_REV,MODULE_DATE)
!EVAC:IF (MODULE_REV > REV_NUMBER) THEN
!EVAC:   REV_NUMBER = MODULE_REV
!EVAC:   WRITE(REV_DATE,'(A)') MODULE_DATE
!EVAC:ENDIF
CALL GET_REV_fire(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_func(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_ieva(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_init(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_irad(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_mass(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_mesh(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_part(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_pois(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_prec(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_pres(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_radi(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_read(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_smvv(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_type(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_velo(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_wall(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF

END SUBROUTINE GET_REVISION_NUMBER

END PROGRAM FDS