pres.f90

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

MODULE PRES
 
! Find the perturbation pressure by solving Poisson's Equation
 
USE PRECISION_PARAMETERS
USE MESH_POINTERS

IMPLICIT NONE

PRIVATE
CHARACTER(255), PARAMETER :: presid='$Id: pres.f90 567 2007-09-11 20:48:44Z drjfloyd $'
CHARACTER(255), PARAMETER :: presrev='$Revision: 567 $'
CHARACTER(255), PARAMETER :: presdate='$Date: 2007-09-11 16:48:44 -0400 (Tue, 11 Sep 2007) $'
PUBLIC PRESSURE_SOLVER,COMPUTE_A_B,UPDATE_PRESSURE,COMPUTE_CORRECTION_PRESSURE,COMPUTE_C,GET_REV_PRES
 
CONTAINS
 
SUBROUTINE PRESSURE_SOLVER(NM)
USE POIS, ONLY: H3CZSS,H2CZSS,H2CYSS,H3CSSS
USE COMP_FUNCTIONS, ONLY: SECOND
USE GLOBAL_CONSTANTS
 
INTEGER, INTENT(IN) :: NM
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW
REAL(EB), POINTER, DIMENSION(:) :: UWP
INTEGER :: I,J,K,IW,IOR,BC_TYPE
REAL(EB) :: TRM1,TRM2,TRM3,TRM4,RES,LHSS,RHSS,HH, DWDT,DVDT,DUDT,HQ2,RFODT,U2,V2,W2,HFAC,H0RR(6),TNOW
LOGICAL :: GET_H
 
IF (SOLID_PHASE_ONLY) RETURN

TNOW=SECOND()
CALL POINT_TO_MESH(NM)
 
IF (PREDICTOR) THEN
   UU => U
   VV => V
   WW => W
   UWP=> UWS
ELSE
   UU => US
   VV => VS
   WW => WS
   UWP=> UW
ENDIF
 
RFODT = RELAXATION_FACTOR/DT
HFAC  = 1._EB-RHOA/RHO_AVG
H0RR  = 0.
IF (U0>=0._EB) H0RR(1) = H0*RHOA/RHO_AVG
IF (U0<=0._EB) H0RR(2) = H0*RHOA/RHO_AVG
IF (V0>=0._EB) H0RR(3) = H0*RHOA/RHO_AVG
IF (V0<=0._EB) H0RR(4) = H0*RHOA/RHO_AVG
IF (W0>=0._EB) H0RR(5) = H0*RHOA/RHO_AVG
IF (W0<=0._EB) H0RR(6) = H0*RHOA/RHO_AVG
 
! Apply pressure boundary conditions at external cells.
! If Neumann, BXS, BXF, etc., contain dH/dx(x=XS), dH/dx(x=XF), etc.
! If Dirichlet, BXS, BXF, etc., contain H(x=XS), H(x=XF), etc.
! LBC, MBC and NBC are codes used be Poisson solver to denote type
! of boundary condition at x, y and z boundaries. See Crayfishpak
! manual for details.
 
WALL_CELL_LOOP: DO IW=1,NEWC
 
   I   = IJKW(1,IW)
   J   = IJKW(2,IW)
   K   = IJKW(3,IW)
   IOR = IJKW(4,IW)
   SELECT CASE(IOR)
      CASE( 1)
         IF (LBC==3 .OR. LBC==4 .OR. LBC==6) BC_TYPE = NEUMANN
         IF (LBC==1 .OR. LBC==2 .OR. LBC==5) BC_TYPE = DIRICHLET
      CASE(-1)
         IF (LBC==2 .OR. LBC==3 .OR. LBC==6) BC_TYPE = NEUMANN
         IF (LBC==1 .OR. LBC==4 .OR. LBC==5) BC_TYPE = DIRICHLET
      CASE( 2)
         IF (MBC==3 .OR. MBC==4) BC_TYPE = NEUMANN
         IF (MBC==1 .OR. MBC==2) BC_TYPE = DIRICHLET
      CASE(-2)
         IF (MBC==3 .OR. MBC==2) BC_TYPE = NEUMANN
         IF (MBC==1 .OR. MBC==4) BC_TYPE = DIRICHLET
      CASE( 3)
         IF (NBC==3 .OR. NBC==4) BC_TYPE = NEUMANN
         IF (NBC==1 .OR. NBC==2) BC_TYPE = DIRICHLET
      CASE(-3)
         IF (NBC==3 .OR. NBC==2) BC_TYPE = NEUMANN
         IF (NBC==1 .OR. NBC==4) BC_TYPE = DIRICHLET
   END SELECT
      
   IF_NEUMANN: IF (BC_TYPE==NEUMANN) THEN
      SELECT CASE(IOR)
         CASE( 1)
            BXS(J,K) = HX(0)   *(-FVX(0,J,K)    + DUWDT(IW))
         CASE(-1)
            BXF(J,K) = HX(IBP1)*(-FVX(IBAR,J,K) - DUWDT(IW))
         CASE( 2)
            BYS(I,K) = HY(0)   *(-FVY(I,0,K)    + DUWDT(IW))
         CASE(-2)
            BYF(I,K) = HY(JBP1)*(-FVY(I,JBAR,K) - DUWDT(IW))
         CASE( 3)
            BZS(I,J) = HZ(0)   *(-FVZ(I,J,0)    + DUWDT(IW))
         CASE(-3)
            BZF(I,J) = HZ(KBP1)*(-FVZ(I,J,KBAR) - DUWDT(IW))
      END SELECT
   ENDIF IF_NEUMANN
 
   IF_DIRICHLET: IF (BC_TYPE==DIRICHLET) THEN
 
      NOT_OPEN: IF (BOUNDARY_TYPE(IW)/=OPEN_BOUNDARY) THEN
 
         SELECT CASE(IOR)
            CASE( 1)
               DUDT = -RFODT*(UU(0,J,K)   +UWP(IW))
               HH = H(1,J,K)
            CASE(-1) 
               DUDT = -RFODT*(UU(IBAR,J,K)-UWP(IW))
               HH = H(IBAR,J,K)
            CASE( 2)
               DVDT = -RFODT*(VV(I,0,K)   +UWP(IW)) 
               HH = H(I,1,K)
            CASE(-2) 
               DVDT = -RFODT*(VV(I,JBAR,K)-UWP(IW))
               HH = H(I,JBAR,K)
            CASE( 3)
               DWDT = -RFODT*(WW(I,J,0)   +UWP(IW))
               HH = H(I,J,1)
            CASE(-3) 
               DWDT = -RFODT*(WW(I,J,KBAR)-UWP(IW))
               HH = H(I,J,KBAR)
         END SELECT
 
         GET_H = .FALSE.
 
         IF (BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) THEN
            SELECT CASE(IOR)
               CASE( 1)
                  IF (UU(0,J,K)   >0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
               CASE(-1)
                  IF (UU(IBAR,J,K)<0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
               CASE( 2)
                  IF (VV(I,0,K)   >0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
               CASE(-2)
                  IF (VV(I,JBAR,K)<0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
               CASE( 3)
                  IF (WW(I,J,0)   >0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
               CASE(-3)
                  IF (WW(I,J,KBAR)<0._EB .AND. UWP(IW)<0._EB) GET_H=.TRUE.
            END SELECT
            IF (GET_H) HH = OMESH(IJKW(9,IW))%H(IJKW(10,IW),IJKW(11,IW),IJKW(12,IW))
         ENDIF
 
         SELECT CASE(IOR)
            CASE( 1)
               BXS(J,K) = HH + 0.5_EB*DX(0)   *(DUDT+FVX(0,J,K))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BXS(J,K) = HH 
            CASE(-1) 
               BXF(J,K) = HH - 0.5_EB*DX(IBP1)*(DUDT+FVX(IBAR,J,K))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BXF(J,K) = HH
            CASE( 2) 
               BYS(I,K) = HH + 0.5_EB*DY(0)   *(DVDT+FVY(I,0,K))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BYS(I,K) = HH
            CASE(-2) 
               BYF(I,K) = HH - 0.5_EB*DY(JBP1)*(DVDT+FVY(I,JBAR,K))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BYF(I,K) = HH
            CASE( 3) 
               BZS(I,J) = HH + 0.5_EB*DZ(0)   *(DWDT+FVZ(I,J,0))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BZS(I,J) = HH
            CASE(-3) 
               BZF(I,J) = HH - 0.5_EB*DZ(KBP1)*(DWDT+FVZ(I,J,KBAR))
               IF (GET_H.AND.BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) BZF(I,J) = HH
         END SELECT
 
      ENDIF NOT_OPEN
 
      OPEN: IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN
 
         SELECT CASE(IOR)
            CASE( 1)
               U2  = UU(0,J,K)**2
               V2  = .25_EB*(VV(1,J,K)+VV(1,J-1,K))**2
               W2  = .25_EB*(WW(1,J,K)+WW(1,J,K-1))**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (UU(0,J,K)<0._EB) THEN
                  BXS(J,K) = HQ2
               ELSE
                  BXS(J,K) = H0RR(1) + HQ2*HFAC
               ENDIF
            CASE(-1)
               U2  = UU(IBAR,J,K)**2
               V2  = .25_EB*(VV(IBAR,J,K)+VV(IBAR,J-1,K))**2
               W2  = .25_EB*(WW(IBAR,J,K)+WW(IBAR,J,K-1))**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (UU(IBAR,J,K)>0._EB) THEN
                  BXF(J,K) = HQ2
               ELSE
                  BXF(J,K) = H0RR(2) + HQ2*HFAC
               ENDIF
            CASE( 2)
               U2  = .25_EB*(UU(I,1,K)+UU(I-1,1,K))**2
               V2  = VV(I,0,K)**2
               W2  = .25_EB*(WW(I,1,K)+WW(I,1,K-1))**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (VV(I,0,K)<0._EB) THEN
                  BYS(I,K) = HQ2
               ELSE
                  BYS(I,K) = H0RR(3) + HQ2*HFAC
               ENDIF
            CASE(-2)
               U2  = .25_EB*(UU(I,JBAR,K)+UU(I-1,JBAR,K))**2
               V2  = VV(I,JBAR,K)**2
               W2  = .25_EB*(WW(I,JBAR,K)+WW(I,JBAR,K-1))**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (VV(I,JBAR,K)>0._EB) THEN
                  BYF(I,K) = HQ2
               ELSE
                  BYF(I,K) = H0RR(4) + HQ2*HFAC
               ENDIF
            CASE( 3)
               U2  = .25_EB*(UU(I,J,1)+UU(I-1,J,1))**2
               V2  = .25_EB*(VV(I,J,1)+VV(I,J-1,1))**2
               W2  = WW(I,J,0)**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (WW(I,J,0)<0._EB) THEN
                  BZS(I,J) = HQ2
               ELSE
                  BZS(I,J) = H0RR(5) + HQ2*HFAC
               ENDIF
            CASE(-3)
               U2  = .25_EB*(UU(I,J,KBAR)+UU(I-1,J,KBAR))**2
               V2  = .25_EB*(VV(I,J,KBAR)+VV(I,J-1,KBAR))**2
               W2  = WW(I,J,KBAR)**2
               HQ2 = MIN(5000._EB,0.5_EB*(U2+V2+W2))
               IF (WW(I,J,KBAR)>0._EB) THEN
                  BZF(I,J) = HQ2
               ELSE
                  BZF(I,J) = H0RR(6) + HQ2*HFAC
               ENDIF
         END SELECT
    
      ENDIF OPEN
 
   ENDIF IF_DIRICHLET
 
ENDDO WALL_CELL_LOOP
! Compute the RHS of the Poisson equation
 
IF (CYLINDRICAL) THEN
   DO K=1,KBAR
      DO I=1,IBAR
         TRM1 = (R(I-1)*FVX(I-1,1,K)-R(I)*FVX(I,1,K))*RDX(I)*RRN(I)
         TRM3 = (FVZ(I,1,K-1)-FVZ(I,1,K))*RDZ(K)
         TRM4 = -DDDT(I,1,K)
         PRHS(I,1,K) = TRM1 + TRM3 + TRM4
      ENDDO
   ENDDO
ENDIF
 
IF ( (IPS<=1 .OR. IPS==4 .OR. IPS==7) .AND. .NOT.CYLINDRICAL) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            TRM1 = (FVX(I-1,J,K)-FVX(I,J,K))*RDX(I)
            TRM2 = (FVY(I,J-1,K)-FVY(I,J,K))*RDY(J) 
            TRM3 = (FVZ(I,J,K-1)-FVZ(I,J,K))*RDZ(K)
            TRM4 = -DDDT(I,J,K)
            PRHS(I,J,K) = TRM1 + TRM2 + TRM3 + TRM4
         ENDDO
      ENDDO
   ENDDO
END IF
 
IF (IPS==2) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            TRM1 = (FVX(I-1,J,K)-FVX(I,J,K))*RDX(I)
            TRM2 = (FVY(I,J-1,K)-FVY(I,J,K))*RDY(J)
            TRM3 = (FVZ(I,J,K-1)-FVZ(I,J,K))*RDZ(K)
            TRM4 = -DDDT(I,J,K)
            PRHS(J,I,K) = TRM1 + TRM2 + TRM3 + TRM4
         ENDDO
      ENDDO
   ENDDO
   BZST = TRANSPOSE(BZS)
   BZFT = TRANSPOSE(BZF)
ENDIF
 
IF (IPS==3 .OR. IPS==6) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            TRM1 = (FVX(I-1,J,K)-FVX(I,J,K))*RDX(I)
            TRM2 = (FVY(I,J-1,K)-FVY(I,J,K))*RDY(J)
            TRM3 = (FVZ(I,J,K-1)-FVZ(I,J,K))*RDZ(K)
            TRM4 = -DDDT(I,J,K)
            PRHS(K,J,I) = TRM1 + TRM2 + TRM3 + TRM4
         ENDDO
      ENDDO
   ENDDO
   BXST = TRANSPOSE(BXS)
   BXFT = TRANSPOSE(BXF)
   BYST = TRANSPOSE(BYS)
   BYFT = TRANSPOSE(BYF)
   BZST = TRANSPOSE(BZS)
   BZFT = TRANSPOSE(BZF)
END IF
 
IF (IPS==5) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            TRM1 = (FVX(I-1,J,K)-FVX(I,J,K))*RDX(I)
            TRM2 = (FVY(I,J-1,K)-FVY(I,J,K))*RDY(J)
            TRM3 = (FVZ(I,J,K-1)-FVZ(I,J,K))*RDZ(K)
            TRM4 = -DDDT(I,J,K)
            PRHS(I,K,J) = TRM1 + TRM2 + TRM3 + TRM4
         ENDDO
      ENDDO
   ENDDO
   BXST = TRANSPOSE(BXS)
   BXFT = TRANSPOSE(BXF)
END IF

! Poisson solve with stretching the 1st coordinate
 
IF (IPS<=1) THEN
   IF (.NOT.TWO_D) CALL H3CZSS(BXS,BXF,BYS,BYF,BZS,BZF,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HX)
   IF (TWO_D .AND. .NOT. CYLINDRICAL) CALL H2CZSS(BXS,BXF,BZS,BZF,ITRN,PRHS,POIS_PTB,SAVE,WORK,HX)
   IF (TWO_D .AND. CYLINDRICAL) CALL H2CYSS(BXS,BXF,BZS,BZF,ITRN,PRHS,POIS_PTB,SAVE,WORK)
ENDIF
 
IF (IPS==2) CALL H3CZSS(BYS,BYF,BXS,BXF,BZST,BZFT,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HY)
 
IF (IPS==3) THEN
   IF (.NOT.TWO_D) CALL H3CZSS(BZST,BZFT,BYST,BYFT,BXST,BXFT,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HZ)
   IF (TWO_D)CALL H2CZSS(BZS,BZF,BXS,BXF,ITRN,PRHS,POIS_PTB,SAVE,WORK,HZ)
ENDIF
 
! Poisson solve with stretching the 1st and 2nd coordinate
 
IF (IPS==4) CALL H3CSSS(BXS,BXF,BYS,BYF,BZS,BZF,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HX,HY)
 
IF (IPS==5) THEN
   IF (.NOT.TWO_D) CALL H3CSSS(BXST,BXFT,BZS,BZF,BYS,BYF,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HX,HZ)
   IF (TWO_D) CALL H2CZSS(BZS,BZF,BXS,BXF,ITRN,PRHS,POIS_PTB,SAVE,WORK,HZ)
ENDIF
 
IF (IPS==6) CALL H3CSSS(BZST,BZFT,BYST,BYFT,BXST,BXFT,ITRN,JTRN,PRHS,POIS_PTB,SAVE,WORK,HZ,HY)
 
IF (IPS==7) CALL H2CZSS(BXS,BXF,BYS,BYF,ITRN,PRHS,POIS_PTB,SAVE,WORK,HX)
 
IF (IPS<=1 .OR. IPS==4 .OR. IPS==7) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            H(I,J,K) = PRHS(I,J,K)
         ENDDO
      ENDDO
   ENDDO
ENDIF
 
IF (IPS==2) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            H(I,J,K) = PRHS(J,I,K)
         ENDDO
      ENDDO
   ENDDO
ENDIF
 
IF (IPS==3 .OR. IPS==6) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            H(I,J,K) = PRHS(K,J,I)
         ENDDO
      ENDDO
   ENDDO
ENDIF
 
IF (IPS==5) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            H(I,J,K) = PRHS(I,K,J)
         ENDDO
      ENDDO
   ENDDO
END IF
 
! Apply boundary conditions to H
 
DO K=1,KBAR
   DO J=1,JBAR
      IF (LBC==3 .OR. LBC==4)             H(0,J,K)    = H(1,J,K)    - DXI*BXS(J,K)
      IF (LBC==3 .OR. LBC==2 .OR. LBC==6) H(IBP1,J,K) = H(IBAR,J,K) + DXI*BXF(J,K)
      IF (LBC==1 .OR. LBC==2)             H(0,J,K)    =-H(1,J,K)    + 2._EB*BXS(J,K)
      IF (LBC==1 .OR. LBC==4 .OR. LBC==5) H(IBP1,J,K) =-H(IBAR,J,K) + 2._EB*BXF(J,K)