main_mpi.f90

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

PROGRAM FDS  

! Fire Dynamics Simulator, Main Program, Multiple CPU version.

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
USE TRAN
USE DUMP
USE READ_INPUT
USE INIT
USE DIVG
USE PRES
USE MASS
USE PART
USE VELO
USE RAD
USE MEMORY_FUNCTIONS
USE COMP_FUNCTIONS, ONLY : SECOND, WALL_CLOCK_TIME, SHUTDOWN
USE MATH_FUNCTIONS, ONLY : GAUSSJ
USE DEVICE_VARIABLES
USE WALL_ROUTINES
USE FIRE
USE RADCONS
USE CONTROL_FUNCTIONS
!EVAC:USE EVAC

IMPLICIT NONE

! Miscellaneous declarations

CHARACTER(255), PARAMETER :: mainmpiid='$Id: main_mpi.f90 719 2007-10-01 17:09:23Z mcgratta $'
CHARACTER(255), PARAMETER :: mainmpirev='$Revision: 719 $'
CHARACTER(255), PARAMETER :: mainmpidate='$Date: 2007-10-01 13:09:23 -0400 (Mon, 01 Oct 2007) $'
LOGICAL  :: EX,DIAGNOSTICS,EXCHANGE_RADIATION=.TRUE.
INTEGER  :: LO10,NM,IZERO,DATE_TIME(8),NN,REVISION_NUMBER
CHARACTER(10) :: BIG_BEN(3)
CHARACTER(255) :: REVISION_DATE
REAL(EB) :: T_MAX,T_MIN
REAL(EB), ALLOCATABLE, DIMENSION(:) ::  T,TC_GLB,TC_LOC,DT_SYNC, DTNEXT_SYNC,DSUM_ALL,PSUM_ALL,USUM_ALL
INTEGER, ALLOCATABLE, DIMENSION(:) ::  MESH_STOP_STATUS,COUNT_GLB,COUNT_LOC
LOGICAL, ALLOCATABLE, DIMENSION(:) ::  ACTIVE_MESH,STATE_GLB,STATE_LOC
INTEGER NOM,IWW,IW
INTEGER, PARAMETER :: N_DROP_ADOPT_MAX=10000
TYPE (MESH_TYPE), POINTER :: M,M4
TYPE (OMESH_TYPE), POINTER :: M2,M3
 
! MPI stuff

! INCLUDE '/usr/local/include/mpif.h'  ! Uncomment this line only if the compiler needs help finding mpif.h
INCLUDE 'mpif.h'
INTEGER :: N,MYID=0,NUMPROCS=1,I,IERR,STATUS(MPI_STATUS_SIZE)
INTEGER :: RNODE,BUFFER_SIZE,TAG,PNAMELEN
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: TAGS
INTEGER, ALLOCATABLE, DIMENSION(:) :: REQ,PREQ
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ARRAY_OF_STATUSES
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ARRAY_OF_STATUSES2
INTEGER :: N_REQ,N_PREQ
CHARACTER(MPI_MAX_PROCESSOR_NAME) PNAME
 
! Initialize MPI (First executable lines of code)
 
CALL MPI_INIT(IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS, IERR)
CALL MPI_GET_PROCESSOR_NAME(PNAME, PNAMELEN, IERR)
 
WRITE(LU_ERR,'(A,I2,A,I2,A,A)') 'Mesh ',MYID+1,' of ', NUMPROCS,' is alive on ',PNAME(1:PNAMELEN)

! Start wall clock timing

WALL_CLOCK_START = WALL_CLOCK_TIME()
 
! Assign a compilation date (All Nodes)
WRITE(REVISION_DATE,'(A)') mainmpirev(INDEX(mainmpirev,':')+1:LEN_TRIM(mainmpirev)-2)
READ (REVISION_DATE,'(I5)') REVISION_NUMBER
WRITE(REVISION_DATE,'(A)') mainmpidate
CALL GET_REVISION_NUMBER(REVISION_NUMBER,REVISION_DATE)
SVN_REVISION_NUMBER = REVISION_NUMBER
WRITE(COMPILE_DATE,'(A)') REVISION_DATE(INDEX(REVISION_DATE,'(')+1:INDEX(REVISION_DATE,')')-1)
WRITE(VERSION_STRING,'(A)') '5.0.0'
VERSION_NUMBER = 5.0  ! Just use to indicate the major version
PARALLEL       = .TRUE.
 
! Read input from CHID.data file (All Nodes)
 
CALL READ_DATA(MYID)
 
IF (NMESHES/=NUMPROCS) CALL SHUTDOWN('ERROR: Number of meshes not equal to '// 'number of threads')
 
! Read input for EVACUATION routine
 
! IF (ANY(EVACUATION_GRID)) CALL READ_EVAC
 
! Open and write to Smokeview file (Master Node Only)
 
CALL ASSIGN_FILE_NAMES
IF (MYID==0) CALL WRITE_SMOKEVIEW_FILE
 
! Stop all the processes if this is just a set-up run
 
IF (SET_UP) CALL SHUTDOWN('Stop FDS, Set-up only')
 
! Set up Time arrays (All Nodes)
 
ALLOCATE(ACTIVE_MESH(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','ACTIVE_MESH',IZERO)
ALLOCATE(T(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','T',IZERO)
ALLOCATE(DT_SYNC(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','DT_SYNC',IZERO)
ALLOCATE(DTNEXT_SYNC(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','DTNEXT_SYNC',IZERO)
ALLOCATE(MESH_STOP_STATUS(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','MESH_STOP_STATUS',IZERO)

! Set up dummy arrays to hold various arrays that must be exchanged among meshes

ALLOCATE(COUNT_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','COUNT_LOC',IZERO) 
ALLOCATE(COUNT_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','COUNT_GLB',IZERO) 
ALLOCATE(STATE_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','STATE_GLB',IZERO) 
ALLOCATE(STATE_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','STATE_LOC',IZERO) 
ALLOCATE(TC_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TC_GLB',IZERO)
ALLOCATE(TC_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TC_LOC',IZERO)

! Allocate a few arrays needed to exchange divergence and pressure info among meshes

IF (N_ZONE > 0) THEN
   ALLOCATE(DSUM_ALL(N_ZONE),STAT=IZERO)
   ALLOCATE(PSUM_ALL(N_ZONE),STAT=IZERO)
   ALLOCATE(USUM_ALL(N_ZONE),STAT=IZERO)
ENDIF

! Start the clock

T     = T_BEGIN
MESH_STOP_STATUS = NO_STOP
 
! Create unique tags for all mesh exchanges
 
ALLOCATE(REQ(NMESHES*NMESHES*10)) 
REQ = MPI_REQUEST_NULL
ALLOCATE(PREQ(NMESHES*NMESHES*10)) 
PREQ = MPI_REQUEST_NULL
ALLOCATE(ARRAY_OF_STATUSES(MPI_STATUS_SIZE,NMESHES*NMESHES*10))
ALLOCATE(ARRAY_OF_STATUSES2(MPI_STATUS_SIZE,NMESHES*NMESHES*10))
ALLOCATE(TAGS(NMESHES,NMESHES,0:2))
TAG = 0
DO NM=1,NMESHES
   DO NOM=NM,NMESHES
      TAG = TAG+1
      TAGS(NM,NOM,0) = TAG
      TAGS(NOM,NM,0) = TAG
   ENDDO
ENDDO
TAGS(:,:,1) = TAGS(:,:,0) + 1000
TAGS(:,:,2) = TAGS(:,:,0) + 2000
 
! Initialize global parameters (All Nodes)
 
CALL INITIALIZE_GLOBAL_VARIABLES
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize radiation (All Nodes)
 
IF (RADIATION) CALL INIT_RADIATION
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Allocate and initialize mesh-specific variables
 
DO NM=MYID+1,NMESHES,NUMPROCS
   CALL INITIALIZE_MESH_VARIABLES(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Allocate and initialize mesh variable exchange arrays
 
DO NM=MYID+1,NMESHES,NUMPROCS
CALL INITIALIZE_MESH_EXCHANGE(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
CALL MPI_ALLGATHER(I_MIN(1,MYID+1),NMESHES,MPI_INTEGER,I_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(I_MAX(1,MYID+1),NMESHES,MPI_INTEGER,I_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(J_MIN(1,MYID+1),NMESHES,MPI_INTEGER,J_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(J_MAX(1,MYID+1),NMESHES,MPI_INTEGER,J_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(K_MIN(1,MYID+1),NMESHES,MPI_INTEGER,K_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(K_MAX(1,MYID+1),NMESHES,MPI_INTEGER,K_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHER(NIC(1,MYID+1),  NMESHES,MPI_INTEGER,NIC,   NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)
 
I_MIN = TRANSPOSE(I_MIN)
I_MAX = TRANSPOSE(I_MAX)
J_MIN = TRANSPOSE(J_MIN)
J_MAX = TRANSPOSE(J_MAX)
K_MIN = TRANSPOSE(K_MIN)
K_MAX = TRANSPOSE(K_MAX)
NIC   = TRANSPOSE(NIC)
 
DO NM=MYID+1,NMESHES,NUMPROCS
   CALL DOUBLE_CHECK(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Potentially read data from a previous calculation 
 
DO NM=MYID+1,NMESHES,NUMPROCS
   IF (RESTART) CALL READ_RESTART(T(NM),NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize output files containing global data (Master Node Only)
 
IF (MYID==0) CALL INITIALIZE_GLOBAL_DUMPS
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize output files that are mesh-specific
 
DO NM=MYID+1,NMESHES,NUMPROCS
   CALL INITIALIZE_MESH_DUMPS(NM)
   CALL INITIALIZE_DROPLETS(NM)
   CALL INITIALIZE_TREES(NM)
!     IF (ANY(EVACUATION_GRID)) CALL INITIALIZE_EVACUATION(NM)
   CALL POST_RECEIVES(NM,0)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Write out character strings to .smv file
 
CALL WRITE_STRINGS
 
! Initialize Mesh Exchange Arrays (All Nodes)
 
CALL MESH_EXCHANGE(0)
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

! Make an initial dump of ambient values

DO NM=MYID+1,NMESHES,NUMPROCS
   CALL UPDATE_OUTPUTS(T(NM),NM)      
   CALL DUMP_MESH_OUTPUTS(T(NM),NM)
ENDDO

CALL MPI_ALLGATHER(T(MYID+1),1,MPI_DOUBLE_PRECISION,T,1, MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
CALL UPDATE_CONTROLS(T)
CALL DUMP_GLOBAL_OUTPUTS(T(1))
 
! Check for changes in VENT or OBSTruction control and device status at t=T_BEGIN

OBST_VENT_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS
   CALL OPEN_AND_CLOSE(T(NM),NM)  
ENDDO OBST_VENT_LOOP

!***********************************************************************************************************************************
!                                                           MAIN TIMESTEPPING LOOP
!***********************************************************************************************************************************
 
MAIN_LOOP: DO  
 
   ICYC  = ICYC + 1 
 
   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,' starts iteration',ICYC,' at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
   EXCHANGE_RADIATION = .FALSE.
   IF (RADIATION) THEN
      IF (MOD(ICYC,ANGLE_INCREMENT*TIME_STEP_INCREMENT)==0) EXCHANGE_RADIATION = .TRUE.
   ENDIF
 
   ! Synchronize clocks
   
   CALL MPI_ALLGATHER(T(MYID+1),1,MPI_DOUBLE_PRECISION,T,1, MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)

   ! Check for program stops
 
   INQUIRE(FILE=TRIM(CHID)//'.stop',EXIST=EX)
   IF (EX) MESH_STOP_STATUS = USER_STOP
   CALL MPI_ALLGATHER(MESH_STOP_STATUS(MYID+1),1,MPI_INTEGER,MESH_STOP_STATUS,1,MPI_INTEGER,MPI_COMM_WORLD,IERR)
 
   ! Figure out fastest and slowest meshes
 
   T_MAX = -1000000._EB
   T_MIN =  1000000._EB
   DO NM=1,NMESHES
      T_MIN = MIN(T(NM),T_MIN)
      T_MAX = MAX(T(NM),T_MAX)
      IF (MESH_STOP_STATUS(NM)/=NO_STOP) GLOBAL_STOP_STATUS = MESH_STOP_STATUS(NM)
   ENDDO
 
   ! Determine time step
 
   IF (SYNCHRONIZE) THEN
      DT_SYNC(MYID+1) = MESHES(MYID+1)%DTNEXT
      CALL MPI_ALLGATHER(DT_SYNC(MYID+1),1,MPI_DOUBLE_PRECISION, DT_SYNC,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
      IF (SYNC_TIME_STEP(MYID+1)) THEN
         MESHES(MYID+1)%DTNEXT = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)
         T(MYID+1) = MINVAL(T,MASK=SYNC_TIME_STEP)
         ACTIVE_MESH(MYID+1) = .TRUE.
      ELSE
         ACTIVE_MESH(MYID+1) = .FALSE.
         IF (T(MYID+1)+MESHES(MYID+1)%DTNEXT <= T_MAX)  ACTIVE_MESH(MYID+1) = .TRUE.
         IF (GLOBAL_STOP_STATUS/=NO_STOP) ACTIVE_MESH(MYID+1) = .TRUE.
      ENDIF
   ELSE
      ACTIVE_MESH = .FALSE.
      DO NM=1,NMESHES
         IF (T(NM)+MESHES(NM)%DTNEXT <= T_MAX) ACTIVE_MESH(NM) = .TRUE.
         IF (GLOBAL_STOP_STATUS/=NO_STOP) ACTIVE_MESH(NM) = .TRUE.
      ENDDO
   ENDIF
 
   ! Determine when to dump out diagnostics to the .out file
 
   DIAGNOSTICS = .FALSE.
   LO10 = LOG10(REAL(ICYC,EB))
   IF (MOD(ICYC,10**LO10)==0 .OR. MOD(ICYC,100)==0 .OR. T_MIN>=T_END .OR. GLOBAL_STOP_STATUS/=NO_STOP) DIAGNOSTICS = .TRUE.
 
   ! Give every processor the full ACTIVE_MESH array
 
   CALL MPI_ALLGATHER(ACTIVE_MESH(MYID+1), 1, MPI_LOGICAL,  ACTIVE_MESH,1, MPI_LOGICAL, MPI_COMM_WORLD, IERR)
 
   ! If no meshes are due to be updated, update them all
 
   IF (ALL(.NOT.ACTIVE_MESH)) ACTIVE_MESH = .TRUE.
 
!     Do not do EVACuation if past the max iteration criteria
!
!     IF (ANY(EVACUATION_GRID)) THEN
!       EVAC_DT = 1000000.
!       DO NM=1,NMESHES
!         IF (.NOT.EVACUATION_ONLY(NM))
!    .         EVAC_DT = MIN(EVAC_DT,MESHES(NM)%DTNEXT)
!       ENDDO
!       DO NM=1,NMESHES
!         IF (EVACUATION_ONLY(NM)) THEN
!           IF (ICYC > EVAC_TIME_ITERATIONS) THEN
!             ACTIVE_MESH(NM) = .FALSE.
!             EVAC_DT = MIN(EVAC_DT, EVAC_DT_STEADY_STATE)
!             MESHES(NM)%DT     = EVAC_DT
!             T(NM)           = T(NM) + MESHES(NM)%DT
!             IF (EVACUATION_GRID(NM) ) THEN
!               CALL EVACUATE_HUMANS(T(NM),NM)
!               IF (T(NM)>=PART_CLOCK(NM)) THEN
!                 CALL DUMP_EVAC(PART_CLOCK(NM),NM)
!                 DO
!                   PART_CLOCK(NM) = PART_CLOCK(NM) + WPAR
!                   IF (PART_CLOCK(NM)>=T(NM)) EXIT
!                 ENDDO
!               ENDIF
!             ENDIF
!           ELSE
!             ACTIVE_MESH(NM) = .TRUE.
!             EVAC_DT = MIN(EVAC_DT, EVAC_DT_FLOWFIELD)
!           ENDIF
!         ENDIF
!       ENDDO
!       DO NM=1,NMESHES
!         IF (EVACUATION_ONLY(NM)) MESHES(NM)%DTNEXT = EVAC_DT
!       ENDDO
!     ENDIF
!
   PREDICTOR = .TRUE.
   CORRECTOR = .FALSE.

   ! Diagnostic calls

   IF (DEBUG) WRITE(0,*) 'Cycle ',ICYC,' Mesh ',MYID+1, ' starting',ACTIVE_MESH(MYID+1)

   IF (MOD(ICYC,3)==0 .AND. TIMING) THEN
      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
      IF (ACTIVE_MESH(MYID+1)) WRITE(0,'(A,I2,A,I2,A,I3.3)') ' Thread ',MYID+1,' is active at ', DATE_TIME(7),'.',DATE_TIME(8)
   ENDIF
 
   ! Begin the PREDICTOR step
 
   COMPUTE_FINITE_DIFFERENCES_1: DO NM=MYID+1,NMESHES,NUMPROCS
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_1
      MESHES(NM)%DT = MESHES(NM)%DTNEXT
      NTCYC(NM)   = NTCYC(NM) + 1
 
      CALL INSERT_DROPLETS_AND_PARTICLES(T(NM),NM)
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM)
      CALL UPDATE_PARTICLES(T(NM),NM)
 
!     IF (EVACUATION_ONLY(NM)) CALL EVACUATE_HUMANS(T(NM),NM)
 
      IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) CALL MASS_FINITE_DIFFERENCES(NM)
 
   ENDDO COMPUTE_FINITE_DIFFERENCES_1
 
   CHANGE_TIME_STEP_LOOP: DO
 
      COMPUTE_DIVERGENCE_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_DIVERGENCE_LOOP
         IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
            CALL DENSITY(NM)
            CALL WALL_BC(T(NM),NM)
         ENDIF
         CALL DIVERGENCE_PART_1(T(NM),NM)
      ENDDO COMPUTE_DIVERGENCE_LOOP

      CALL EXCHANGE_DIVERGENCE_INFO

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

      IF (PRESSURE_CORRECTION) CALL CORRECT_PRESSURE(1)
 
!     IF (ANY(EVACUATION_GRID) .AND. EVACUATION_ONLY(NM)) THEN
!        PRESSURE_ITERATION_LOOP: DO N=1,EVAC_PRESSURE_ITERATIONS
!           CALL NO_FLUX
!           CALL PRESSURE_SOLVER(NM)
!        ENDDO PRESSURE_ITERATION_LOOP
!     ENDIF
 
      PREDICT_VELOCITY_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE PREDICT_VELOCITY_LOOP
         CALL VELOCITY_PREDICTOR(T(NM),NM,MESH_STOP_STATUS(NM))
      ENDDO PREDICT_VELOCITY_LOOP
 
      IF (SYNCHRONIZE) THEN
         NM = MYID+1
         CALL MPI_ALLGATHER(CHANGE_TIME_STEP(NM),1,MPI_LOGICAL, CHANGE_TIME_STEP,1,MPI_LOGICAL,MPI_COMM_WORLD,IERR)
         CALL MPI_ALLGATHER(MESH_STOP_STATUS(NM),1,MPI_INTEGER,MESH_STOP_STATUS,1, MPI_INTEGER,MPI_COMM_WORLD,IERR)