main_mpi.f90

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

         IF (ANY(MESH_STOP_STATUS/=NO_STOP)) THEN
            IF (ANY(MESH_STOP_STATUS==INSTABILITY_STOP)) GLOBAL_STOP_STATUS = INSTABILITY_STOP
            EXIT CHANGE_TIME_STEP_LOOP
         ENDIF
         IF (ANY(CHANGE_TIME_STEP)) THEN
            CHANGE_TIME_STEP = .TRUE.
            DT_SYNC(NM) = MESHES(NM)%DT
            DTNEXT_SYNC(NM) = MESHES(NM)%DTNEXT
            CALL MPI_ALLGATHER(DT_SYNC(NM),1,MPI_DOUBLE_PRECISION, DT_SYNC,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
            CALL MPI_ALLGATHER(DTNEXT_SYNC(NM),1,MPI_DOUBLE_PRECISION, DTNEXT_SYNC,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
            MESHES(NM)%DTNEXT = MINVAL(DTNEXT_SYNC,MASK=SYNC_TIME_STEP)
            MESHES(NM)%DT     = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)
         ENDIF
      ENDIF
 
      IF (MESH_STOP_STATUS(MYID+1)/=NO_STOP) THEN
         GLOBAL_STOP_STATUS = MESH_STOP_STATUS(MYID+1)
         EXIT CHANGE_TIME_STEP_LOOP
      ENDIF
 
      IF (.NOT.ANY(CHANGE_TIME_STEP)) EXIT CHANGE_TIME_STEP_LOOP
 
   ENDDO CHANGE_TIME_STEP_LOOP
 
   CHANGE_TIME_STEP = .FALSE.
 
   ! Do the first step in the Runge-Kutta update scheme for sprinklers and detectors
 
   UPDATE_TIME: DO NM=MYID+1,NMESHES,NUMPROCS
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE UPDATE_TIME
      CALL POST_RECEIVES(NM,1)
      T(NM) = T(NM) + MESHES(NM)%DT  ! Advance the time and start the CORRECTOR part of the time step
   ENDDO UPDATE_TIME
 
!====================================================================================================================
 
   ! Exchange information among meshes
 
   IF (MOD(ICYC,3)==0.AND.TIMING.AND.ACTIVE_MESH(MYID+1)) THEN
      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
      WRITE(0,'(A,I2,A,I2,A,I3.3)')  ' Thread ',MYID+1,' enters Mesh Exchange 1 at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
   CALL MESH_EXCHANGE(1)
 
   IF (MOD(ICYC,3)==0.AND.TIMING.AND.ACTIVE_MESH(MYID+1)) THEN
   CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
   WRITE(0,'(A,I2,A,I2,A,I3.3)')  ' Thread ',MYID+1,' exits  Mesh Exchange 1 at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
!+=============================================+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   CORRECTOR = .TRUE.
   PREDICTOR = .FALSE.
 
   COMPUTE_FINITE_DIFFERENCES_2: DO NM=MYID+1,NMESHES,NUMPROCS
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_2
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM)
      IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
         CALL MASS_FINITE_DIFFERENCES(NM)
         CALL DENSITY(NM)
         ! Do combustion, then apply thermal, species and density boundary conditions and solve for radiation
         IF (N_REACTIONS > 0) CALL COMBUSTION (NM)
         CALL WALL_BC(T(NM),NM)
         CALL COMPUTE_RADIATION(NM)
      ENDIF
!     IF (EVACUATION_ONLY(NM)) CALL EVACUATE_HUMANS(T(NM),NM)
      CALL UPDATE_PARTICLES(T(NM),NM)
      CALL DIVERGENCE_PART_1(T(NM),NM)
   ENDDO COMPUTE_FINITE_DIFFERENCES_2

   CALL EXCHANGE_DIVERGENCE_INFO

   COMPUTE_FINITE_DIFFERENCES_3: DO NM=MYID+1,NMESHES,NUMPROCS
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_3
      CALL DIVERGENCE_PART_2(NM)
      CALL PRESSURE_SOLVER(NM)
   ENDDO COMPUTE_FINITE_DIFFERENCES_3

   IF (PRESSURE_CORRECTION) CALL CORRECT_PRESSURE(2)
 
!     IF (ANY(EVACUATION_GRID) .AND. EVACUATION_ONLY(NM)) THEN
!        PRESSURE_ITERATION_LOOP2: DO N=1,EVAC_PRESSURE_ITERATIONS
!           CALL NO_FLUX
!           CALL PRESSURE_SOLVER(NM)
!        ENDDO PRESSURE_ITERATION_LOOP2
!     ENDIF
 
   CORRECT_VELOCITY_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE CORRECT_VELOCITY_LOOP
      CALL POST_RECEIVES(NM,2)   ! Post Receive Arrays
      CALL OPEN_AND_CLOSE(T(NM),NM)   ! Doors, windows, etc.
      CALL VELOCITY_CORRECTOR(T(NM),NM)
      CALL UPDATE_OUTPUTS(T(NM),NM)
      CALL DUMP_MESH_OUTPUTS(T(NM),NM)
      IF (DIAGNOSTICS) CALL CHECK_DIVERGENCE(NM)
   ENDDO CORRECT_VELOCITY_LOOP
 
   !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   ! Exchange information among meshes
 
   IF (MOD(ICYC,3)==0.AND.TIMING.AND.ACTIVE_MESH(MYID+1)) THEN
      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
      WRITE(0,'(A,I2,A,I2,A,I3.3)')  ' Thread ',MYID+1,' enters Mesh Exchange 2 at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
   CALL MESH_EXCHANGE(2)
   IF (MOD(ICYC,3)==0.AND.TIMING.AND.ACTIVE_MESH(MYID+1)) THEN
      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
      WRITE(0,'(A,I2,A,I2,A,I3.3)')  ' Thread ',MYID+1,' exits  Mesh Exchange 2 at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
   !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   ! Write character strings out to the .smv file
 
   IF (DIAGNOSTICS) CALL WRITE_STRINGS
 
   ! Exchange info for diagnostic print out
 
   IF (DIAGNOSTICS) CALL EXCHANGE_DIAGNOSTICS

   ! Check for dumping end of timestep outputs

   CALL DUMP_GLOBAL_OUTPUTS(T_MIN)
   CALL UPDATE_CONTROLS(T)
 
   ! Dump out diagnostics

   IF (MYID==0 .AND. DIAGNOSTICS) CALL WRITE_DIAGNOSTICS(T)
 
   ! Stop the run

   IF (T_MIN>=T_END .OR. GLOBAL_STOP_STATUS/=NO_STOP) EXIT MAIN_LOOP
 
   ! Flush Buffers (All Nodes)
 
   IF (MOD(ICYC,10)==0 .AND. FLUSH_FILE_BUFFERS) THEN
      IF (MYID==0) CALL FLUSH_GLOBAL_BUFFERS
      DO NM=MYID+1,NMESHES,NUMPROCS
         CALL FLUSH_LOCAL_BUFFERS(NM)
      ENDDO
   ENDIF

   IF (MOD(ICYC,3) ==0 .AND. TIMING) THEN
      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
      WRITE(0,'(A,I2,A,I6,A,I2,A,I3.3)')  ' Thread ',MYID+1,' ends iteration',ICYC,' at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
ENDDO MAIN_LOOP
 
!****************************************************************************************************************************
!                                                   END OF TIME STEPPING LOOP
!****************************************************************************************************************************
 
TUSED(1,MYID+1) = SECOND() - TUSED(1,MYID+1)
CALL MPI_GATHER(TUSED(1,MYID+1),N_TIMERS,MPI_DOUBLE_PRECISION, TUSED,N_TIMERS,MPI_DOUBLE_PRECISION,0, MPI_COMM_WORLD,IERR)
IF (MYID==0) CALL TIMINGS
CALL MPI_FINALIZE(IERR)
 
SELECT CASE(GLOBAL_STOP_STATUS)
   CASE(NO_STOP)
      CALL SHUTDOWN('STOP: FDS completed successfully')
   CASE(INSTABILITY_STOP) 
      CALL SHUTDOWN('STOP: Numerical Instability')
   CASE(USER_STOP) 
      CALL SHUTDOWN('STOP: FDS stopped by user')
END SELECT
 
 
CONTAINS
 
 
SUBROUTINE EXCHANGE_DIVERGENCE_INFO

! Exchange information mesh to mesh needed for divergence integrals

IF (N_ZONE > 0) THEN
   CALL MPI_ALLREDUCE(DSUM(1,MYID+1),DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   CALL MPI_ALLREDUCE(PSUM(1,MYID+1),PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   CALL MPI_ALLREDUCE(USUM(1,MYID+1),USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   DSUM(1:N_ZONE,MYID+1) = DSUM_ALL(1:N_ZONE)
   PSUM(1:N_ZONE,MYID+1) = PSUM_ALL(1:N_ZONE)
   USUM(1:N_ZONE,MYID+1) = USUM_ALL(1:N_ZONE)
ENDIF

END SUBROUTINE EXCHANGE_DIVERGENCE_INFO

 
SUBROUTINE INITIALIZE_MESH_EXCHANGE(NM)
 
! Create arrays by which info is to exchanged across meshes
 
INTEGER IMIN,IMAX,JMIN,JMAX,KMIN,KMAX,NOM,IOR,IW
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M2,M
LOGICAL FOUND
 
M=>MESHES(NM)
ALLOCATE(MESHES(NM)%OMESH(NMESHES))
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NOM==NM) CYCLE OTHER_MESH_LOOP
 
   M2=>MESHES(NOM)
   IMIN=0 
   IMAX=M2%IBP1
   JMIN=0 
   JMAX=M2%JBP1
   KMIN=0 
   KMAX=M2%KBP1
   NIC(NOM,NM) = 0
   FOUND = .FALSE.
   SEARCH_LOOP: DO IW=1,M%NEWC
      IF (M%IJKW(9,IW)/=NOM) CYCLE SEARCH_LOOP
      NIC(NOM,NM) = NIC(NOM,NM) + 1
      FOUND = .TRUE.
      IOR = M%IJKW(4,IW)
      SELECT CASE(IOR)
         CASE( 1) 
            IMIN=MAX(IMIN,M%IJKW(10,IW)-1)
         CASE(-1) 
            IMAX=MIN(IMAX,M%IJKW(10,IW))
         CASE( 2) 
            JMIN=MAX(JMIN,M%IJKW(11,IW)-1)
         CASE(-2) 
            JMAX=MIN(JMAX,M%IJKW(11,IW))
         CASE( 3) 
            KMIN=MAX(KMIN,M%IJKW(12,IW)-1)
         CASE(-3) 
            KMAX=MIN(KMAX,M%IJKW(12,IW))
      END SELECT
   ENDDO SEARCH_LOOP
 
   IF ( M2%XS>=M%XS .AND. M2%XF<=M%XF .AND. M2%YS>=M%YS .AND. M2%YF<=M%YF .AND. &
      M2%ZS>=M%ZS .AND. M2%ZF<=M%ZF ) FOUND = .TRUE.
 
   IF (.NOT.FOUND) CYCLE OTHER_MESH_LOOP
 
   I_MIN(NOM,NM) = IMIN
   I_MAX(NOM,NM) = IMAX
   J_MIN(NOM,NM) = JMIN
   J_MAX(NOM,NM) = JMAX
   K_MIN(NOM,NM) = KMIN
   K_MAX(NOM,NM) = KMAX
 
   ALLOCATE(M%OMESH(NOM)% TMP(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%TMP = TMPA
   ALLOCATE(M%OMESH(NOM)%   H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%H = 0.
   ALLOCATE(M%OMESH(NOM)%   U(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%U = U0
   ALLOCATE(M%OMESH(NOM)%   V(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%V = V0
   ALLOCATE(M%OMESH(NOM)%   W(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%W = W0
   ALLOCATE(M%OMESH(NOM)% FVX(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% FVY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% FVZ(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
 
   IF (N_SPECIES>0) THEN
      ALLOCATE(M%OMESH(NOM)%  YY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,N_SPECIES))
      ALLOCATE(M%OMESH(NOM)% YYS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,N_SPECIES))
      DO N=1,N_SPECIES
      M%OMESH(NOM)%YY(:,:,:,N)  = SPECIES(N)%YY0
      M%OMESH(NOM)%YYS(:,:,:,N) = SPECIES(N)%YY0
      ENDDO
   ENDIF
 
   ! Wall arrays
 
   ALLOCATE(M%OMESH(NOM)%IJKW(12,M2%NEWC))
   ALLOCATE(M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC))
   ALLOCATE(M%OMESH(NOM)%WALL(0:M2%NEWC))
 
   ! Particle and Droplet Orphan Arrays
 
   IF (DROPLET_FILE) THEN
      M%OMESH(NOM)%N_DROP_ORPHANS = 0
      M%OMESH(NOM)%N_DROP_ORPHANS_DIM = 1000
      ALLOCATE(M%OMESH(NOM)%DROPLET(M%OMESH(NOM)%N_DROP_ORPHANS_DIM), STAT=IZERO)
      CALL ChkMemErr('INIT','DROPLET',IZERO)
   ENDIF
 
ENDDO OTHER_MESH_LOOP
 
END SUBROUTINE INITIALIZE_MESH_EXCHANGE
 
 

SUBROUTINE DOUBLE_CHECK(NM)
 
! Double check exchange pairs
 
INTEGER NOM
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M2,M
 
M=>MESHES(NM)
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NOM==NM) CYCLE OTHER_MESH_LOOP
 
   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)>0) THEN
      M2=>MESHES(NOM)
      ALLOCATE(M%OMESH(NOM)%IJKW(12,M2%NEWC))
      ALLOCATE(M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC))
      ALLOCATE(M%OMESH(NOM)%WALL(0:M2%NEWC))
   ENDIF
 
ENDDO OTHER_MESH_LOOP
 
END SUBROUTINE DOUBLE_CHECK
 
 

SUBROUTINE POST_RECEIVES(NM,CODE)
USE RADCONS, ONLY: NRA,NSB
 
INTEGER, INTENT(IN) :: NM,CODE
 
N_REQ = 0
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)==0) CYCLE OTHER_MESH_LOOP
   IF (CODE>0 .AND. .NOT.ACTIVE_MESH(NOM)) CYCLE OTHER_MESH_LOOP
 
   IF (DEBUG) THEN
      WRITE(0,*) NM,' posting 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)
   M4=>MESHES(NOM)
   M3=>MESHES(NM)%OMESH(NOM)
 
   RNODE = NOM-1
   TAG   = TAGS(NM,NOM,CODE)
 
   INITIALIZATION_IF: IF (CODE==0) THEN
 
      IF (NIC(NM,NOM)>0) THEN
         ALLOCATE(M3%RPKG1(NIC(NM,NOM)*(3+N_SPECIES)+1))
         ALLOCATE(M3%RPKG2(NIC(NM,NOM)*(9+N_SPECIES)+1))
         IF (PRESSURE_CORRECTION) ALLOCATE(M3%RPKG3(NIC(NM,NOM)*(6      )+1))
         ALLOCATE(M3%WRPKG((NRA*NSB+1)*NIC(NM,NOM)+1))
      ENDIF
 
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%IJKW(1,1),12*M4%NEWC, MPI_INTEGER,RNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
 
      IF (NIC(NM,NOM)>0 .OR. NIC(NOM,NM)>0) THEN
         ALLOCATE(M3%R_RDBUF(13*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%R_IDBUF( 2*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%R_LDBUF(   N_DROP_ADOPT_MAX))
      ENDIF
 
   ENDIF INITIALIZATION_IF
 
   PREDICTOR: IF (CODE==1 .AND. NIC(NM,NOM)>0) THEN
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%RPKG1(1),NIC(NM,NOM)*(3+N_SPECIES)+1, MPI_DOUBLE_PRECISION,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
   ENDIF PREDICTOR
 
   CORRECTOR: IF (CODE==0 .OR. CODE==2) THEN
 
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%BOUNDARY_TYPE(0),M4%NEWC+1, MPI_INTEGER,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
 
      IF (CODE==2 .AND. NIC(NM,NOM)>0) THEN
         N_REQ=N_REQ+1
         CALL MPI_IRECV(M3%RPKG2(1),NIC(NM,NOM)*(9+N_SPECIES)+1, MPI_DOUBLE_PRECISION,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
      ENDIF
      IF (EXCHANGE_RADIATION .AND. NIC(NM,NOM)>0 .AND. CODE==2) THEN
         IWW = NIC(NM,NOM)
         N_REQ=N_REQ+1
         CALL MPI_IRECV(M3%WRPKG(1),(NRA*NSB+1)*IWW+1, MPI_DOUBLE_PRECISION,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
      ENDIF
 
   ENDIF CORRECTOR
 
   ! Droplet Orphan Numbers
 
   IF (DROPLET_FILE .AND.  (NIC(NM,NOM)>0 .OR. NIC(NOM,NM)>0)) THEN
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%N_DROP_ADOPT, 1,MPI_INTEGER,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
   ENDIF
 
   ! Droplet Buffer Arrays

   IF (DROPLET_FILE .AND. (NIC(NM,NOM)>0 .OR. NIC(NOM,NM)>0)) THEN
      BUFFER_SIZE=13*N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%R_RDBUF(1),BUFFER_SIZE, MPI_DOUBLE_PRECISION,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
      BUFFER_SIZE=2*N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%R_IDBUF(1),BUFFER_SIZE, MPI_INTEGER,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
      BUFFER_SIZE=N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%R_LDBUF(1),BUFFER_SIZE, MPI_LOGICAL,RNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
   ENDIF
 
ENDDO OTHER_MESH_LOOP
 
END SUBROUTINE POST_RECEIVES
 
 
 
SUBROUTINE MESH_EXCHANGE(CODE)
USE RADCONS, ONLY: NRA,NSB
 
! Exchange Information between Meshes
 
REAL(EB) :: TNOW
INTEGER, INTENT(IN) :: CODE
INTEGER NM,II,JJ,KK,LL,NC,N,NN,SNODE
INTEGER :: NN1,NN2
 
TNOW = SECOND()
 
! Send Information to other meshes
 
NM = MYID+1
 
SEND_OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NIC(NOM,NM)==0 .AND. NIC(NM,NOM)==0)  CYCLE SEND_OTHER_MESH_LOOP
 
   IF (CODE>0) THEN
      IF (.NOT.ACTIVE_MESH(NM) .OR. .NOT.ACTIVE_MESH(NOM))  CYCLE SEND_OTHER_MESH_LOOP
   ENDIF