Merge pull request #15586 from drjfloyd/master

FDS Source: Revert Flux MW

Author: drjfloyd

Source/cons.f90

@@ -278,7 +278,6 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE.             !< Flag for tracking residence time of spreading fire front over a surface
 LOGICAL :: STORE_LS_SPREAD_RATE=.FALSE.             !< Flag for outputting local level set spread rate magnitude
 LOGICAL :: TEST_NEW_CHAR_MODEL=.FALSE.              !< Flag to envoke new char model
-LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.TRUE.       !< Flag for MW correction ensure consistent equation of state at face
 
 INTEGER, ALLOCATABLE, DIMENSION(:) :: CHANGE_TIME_STEP_INDEX      !< Flag to indicate if a mesh needs to change time step
 INTEGER, ALLOCATABLE, DIMENSION(:) :: SETUP_PRESSURE_ZONES_INDEX  !< Flag to indicate if a mesh needs to keep searching for ZONEs

Source/data.f90

@@ -1360,6 +1360,42 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(553)%UNITS  = 'm/s'
 OUTPUT_QUANTITY(553)%SHORT_NAME  = 'V_LS'
 
+OUTPUT_QUANTITY(560)%NAME  = 'BFX'
+OUTPUT_QUANTITY(560)%UNITS  = ''
+OUTPUT_QUANTITY(560)%SHORT_NAME  = 'bfx'
+OUTPUT_QUANTITY(560)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(560)%IOR = 1
+
+OUTPUT_QUANTITY(561)%NAME  = 'BFY'
+OUTPUT_QUANTITY(561)%UNITS  = ''
+OUTPUT_QUANTITY(561)%SHORT_NAME  = 'bfy'
+OUTPUT_QUANTITY(561)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(561)%IOR = 2
+
+OUTPUT_QUANTITY(562)%NAME  = 'BFZ'
+OUTPUT_QUANTITY(562)%UNITS  = ''
+OUTPUT_QUANTITY(562)%SHORT_NAME  = 'bfz'
+OUTPUT_QUANTITY(562)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(562)%IOR = 3
+
+OUTPUT_QUANTITY(563)%NAME  = 'BFX MINUS'
+OUTPUT_QUANTITY(563)%UNITS  = ''
+OUTPUT_QUANTITY(563)%SHORT_NAME  = 'bfx-'
+OUTPUT_QUANTITY(563)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(563)%IOR = 1
+
+OUTPUT_QUANTITY(564)%NAME  = 'BFY MINUS'
+OUTPUT_QUANTITY(564)%UNITS  = ''
+OUTPUT_QUANTITY(564)%SHORT_NAME  = 'bfy-'
+OUTPUT_QUANTITY(564)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(564)%IOR = 2
+
+OUTPUT_QUANTITY(565)%NAME  = 'BFZ MINUS'
+OUTPUT_QUANTITY(565)%UNITS  = ''
+OUTPUT_QUANTITY(565)%SHORT_NAME  = 'bfz-'
+OUTPUT_QUANTITY(565)%CELL_POSITION = CELL_FACE
+OUTPUT_QUANTITY(565)%IOR = 3
+
 ! Boundary Quantities (Negative indices)
 
 OUTPUT_QUANTITY(-1)%NAME = 'RADIATIVE HEAT FLUX'

Source/divg.f90

@@ -8583,7 +8583,20 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
    CASE(553) ! V_LS
       GAS_PHASE_OUTPUT_RES = V_LS(II,JJ)
 
-END SELECT IND_SELECT
+   CASE(560) ! BFX
+      GAS_PHASE_OUTPUT_RES = BFX(II,JJ,KK)
+   CASE(561) ! BFY
+      GAS_PHASE_OUTPUT_RES = BFY(II,JJ,KK)
+   CASE(562) ! BFZ
+      GAS_PHASE_OUTPUT_RES = BFZ(II,JJ,KK)
+   CASE(563) ! BFX MINUS
+      GAS_PHASE_OUTPUT_RES = BFX(MAX(0,II-1),JJ,KK)
+   CASE(564) ! BFY MINUS
+      GAS_PHASE_OUTPUT_RES = BFY(II,MAX(0,JJ-1),KK)
+   CASE(565) ! BFZ MINUS
+      GAS_PHASE_OUTPUT_RES = BFZ(II,JJ,MAX(0,KK-1))
+
+ END SELECT IND_SELECT
 
 ! Fill GAS_PHASE_OUTPUT for CUT_CELLs.
 ! Some variables have already been filled in fire.f90

Source/dump.f90

@@ -1451,6 +1451,159 @@ END FUNCTION MP5
 END SUBROUTINE GET_SCALAR_FACE_VALUE
 
 
+SUBROUTINE GET_SCALAR_FACE_COEF(A,U,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER)
+
+USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER
+
+REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:)
+REAL(EB), INTENT(INOUT) :: BF(0:,0:,0:)
+INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR
+REAL(EB) :: R,B,DU_UP,DU_LOC
+INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2
+
+SELECT CASE(LIMITER)
+   CASE(GODUNOV_LIMITER,CENTRAL_LIMITER)
+      BF=0._EB
+      RETURN
+END SELECT
+
+SELECT CASE(IOR)
+   CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0
+   CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0
+   CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2
+END SELECT
+
+!$OMP PARALLEL DEFAULT(NONE) IF(I2>1) &
+!$OMP& SHARED(U,A,BF,I1,I2,J1,J2,K1,K2,IP1,JP1,KP1,IM1,JM1,KM1,IP2,JP2,KP2,LIMITER) &
+!$OMP& PRIVATE(I,J,K,DU_UP,DU_LOC,B,R)
+
+SELECT CASE (LIMITER)
+
+   CASE (SUPERBEE_LIMITER)
+      !$OMP DO COLLAPSE(3) SCHEDULE(STATIC)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+               DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K)
+               IF (A(I,J,K) > 0._EB) THEN
+                  DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1)
+               ELSE
+                  DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1)
+               ENDIF
+               R = 0._EB ; B = 0._EB
+               IF (ABS(DU_LOC) > TWO_EPSILON_EB) R = DU_UP/DU_LOC
+               IF (R > TWO_EPSILON_EB) B = MAX(0._EB, MIN(2._EB*R,1._EB), MIN(R,2._EB))
+               BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K)))
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+
+   CASE (CHARM_LIMITER)
+      !$OMP DO COLLAPSE(3) SCHEDULE(STATIC)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+               DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K)
+               IF (A(I,J,K) > 0._EB) THEN
+                  DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1)
+               ELSE
+                  DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1)
+               ENDIF
+               R = 0._EB ; B = 0._EB
+               IF (ABS(DU_UP) > TWO_EPSILON_EB) R = DU_LOC/DU_UP
+               IF (R > TWO_EPSILON_EB) B = R*(3._EB*R+1._EB)/((R+1._EB)**2)
+               BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K)))
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+
+END SELECT
+!$OMP END PARALLEL
+
+END SUBROUTINE GET_SCALAR_FACE_COEF
+
+
+SUBROUTINE GET_SCALAR_FACE_VALUE_NEW(A,U,F,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER)
+
+USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER
+
+REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:),BF(0:,0:,0:)
+REAL(EB), INTENT(OUT) :: F(0:,0:,0:)
+INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR
+REAL(EB) :: DU_UP,DU_LOC
+INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2
+
+SELECT CASE(IOR)
+   CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0
+   CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0
+   CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2
+END SELECT
+
+!$OMP PARALLEL IF(I2>1)
+SELECT CASE(LIMITER)
+   CASE(CENTRAL_LIMITER) ! central differencing
+      !$OMP DO SCHEDULE(STATIC)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+                  F(I,J,K) = 0.5_EB*(U(I,J,K) + U(I+IP1,J+JP1,K+KP1))
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+   CASE(GODUNOV_LIMITER) ! first-order upwinding
+      !$OMP DO SCHEDULE(STATIC)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+               IF (A(I,J,K)>0._EB) THEN
+                  F(I,J,K) = U(I,J,K)
+               ELSE
+                  F(I,J,K) = U(I+IP1,J+JP1,K+KP1)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+   CASE(SUPERBEE_LIMITER) ! SUPERBEE, Roe (1986)
+      !$OMP DO SCHEDULE(STATIC) PRIVATE(DU_LOC)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+               DU_LOC = U(I+IP1,J+JP1,K+KP1)-U(I,J,K)
+               IF (A(I,J,K)>0._EB) THEN
+                  F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_LOC
+               ELSE
+                  F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_LOC
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+   CASE(CHARM_LIMITER) ! CHARM
+      !$OMP DO SCHEDULE(STATIC) PRIVATE(DU_UP)
+      DO K=K1,K2
+         DO J=J1,J2
+            DO I=I1,I2
+               IF (A(I,J,K)>0._EB) THEN
+                  DU_UP  = U(I,J,K)-U(I+IM1,J+JM1,K+KM1)
+                  F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_UP
+               ELSE
+                  DU_UP  = U(I+IP2,J+JP2,K+KP2)-U(I+IP1,J+JP1,K+KP1)
+                  F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_UP
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !$OMP END DO
+END SELECT
+!$OMP END PARALLEL
+
+END SUBROUTINE GET_SCALAR_FACE_VALUE_NEW
+
+
 !> \brief Random fluctuation, Theta'(t+dt) = R^2*Theta'(t) + Normal(0,sqrt(1-R^2)*SIGMA) ; R = exp(-dt/TAU)
 !> \param SIGMA Standard deviation of time series
 !> \param TAU Time scale or period of the time function (s)

Source/func.f90

@@ -499,7 +499,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 USE RADCONS, ONLY: UIIDIM
 USE CONTROL_VARIABLES
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
-INTEGER :: N,I,J,K,IZERO,N_LOWER_SCALARS,NS,IW
+INTEGER :: N,I,J,K,IZERO,NS,IW
 REAL(EB), INTENT(IN) :: DT
 INTEGER, INTENT(IN) :: NM
 REAL(EB) :: MU_N,CS,DELTA
@@ -595,20 +595,6 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 
 ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','RSUM',IZERO) ; M%RSUM=RSUM0
 
-! Allocate scalar face values
-
-! Required for cell face density correction for multicomponent mixtures
-
-IF (FLUX_LIMITER_MW_CORRECTION) THEN
-   N_LOWER_SCALARS=0
-ELSE
-   N_LOWER_SCALARS=1
-ENDIF
-
-ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB
-ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB
-ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB
-
 ! Allocate storage for scalar total fluxes
 
 IF (STORE_SPECIES_FLUX) THEN
@@ -696,9 +682,9 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 ALLOCATE(M%WORK9(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK9',IZERO)
 
 ALLOCATE(M%IWORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)                                 ; CALL ChkMemErr('INIT','IWORK1',IZERO)
-ALLOCATE(M%SWORK1(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK1',IZERO)
-ALLOCATE(M%SWORK2(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK2',IZERO)
-ALLOCATE(M%SWORK3(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK3',IZERO)
+ALLOCATE(M%SWORK1(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK1',IZERO)
+ALLOCATE(M%SWORK2(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK2',IZERO)
+ALLOCATE(M%SWORK3(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK3',IZERO)
 ALLOCATE(M%SWORK4(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)               ; CALL ChkMemErr('INIT','SWORK4',IZERO)
 M%IWORK1=0
 M%SWORK1=0._EB
@@ -833,6 +819,23 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
    M%ZZS(:,:,:,N) = INITIAL_UNMIXED_FRACTION
 ENDDO
 
+! Allocate scalar face values
+
+ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB
+ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB
+ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB
+DO N=1,N_TRACKED_SPECIES
+   M%FX(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0
+   M%FY(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0
+   M%FZ(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0
+ENDDO
+
+! Allocate flux limiter coefficients
+
+ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB
+ALLOCATE( M%BFY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFY',IZERO) ; M%BFY=0._EB
+ALLOCATE( M%BFZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFZ',IZERO) ; M%BFZ=0._EB
+
 ! Allocate and Initialize Mesh-Dependent Radiation Arrays
 
 M%QR  = 0._EB