FDS Source: update to TEST_FLUX_LIMITER_FACE_CORRECTION to account for variable molecular weights

rmcdermo · rmcdermo · commit f2a8b5f36ae0 · 2024-12-06T12:00:36.000-05:00
diff --git a/Source/mass.f90 b/Source/mass.f90
@@ -20,15 +20,16 @@ MODULE MASS
 SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
 USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE
+USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT
 
 INTEGER, INTENT(IN) :: NM
-REAL(EB) :: TNOW
+REAL(EB) :: TNOW,MW_F,MW_G,ZZ_GET(1:N_TRACKED_SPECIES)
 REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,Z_TEMP,F_TEMP
 REAL(EB), PARAMETER :: DUMMY=0._EB
 INTEGER  :: I,J,K,N,IOR,IW,IIG,JJG,KKG,II,JJ,KK,IC
 REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP
-REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P
+REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW
 TYPE(WALL_TYPE), POINTER :: WC
 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
@@ -181,11 +182,25 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
    ! Repeat the above for DENSITY
 
-   CALL GET_SCALAR_FACE_VALUE(UU,RHOP,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER)
-   CALL GET_SCALAR_FACE_VALUE(VV,RHOP,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER)
-   CALL GET_SCALAR_FACE_VALUE(WW,RHOP,FZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER)
+   RHO_RMW=>WORK1
 
-   !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP)
+   !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC)
+   DO K=0,KBP1
+      DO J=0,JBP1
+         DO I=0,IBP1
+            ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES)
+            CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G)
+            RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G
+         ENDDO
+      ENDDO
+   ENDDO
+   !$OMP END PARALLEL DO
+
+   CALL GET_SCALAR_FACE_VALUE(UU,RHO_RMW,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER)
+   CALL GET_SCALAR_FACE_VALUE(VV,RHO_RMW,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER)
+   CALL GET_SCALAR_FACE_VALUE(WW,RHO_RMW,FZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER)
+
+   !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP,ZZ_GET,MW_F)
    WALL_LOOP_3: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
       WC=>WALL(IW)
       IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_3
@@ -211,13 +226,15 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
             CASE(-3); FZ(IIG,JJG,KKG,0)   = 0._EB
          END SELECT
       ELSE
+         ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
+         CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_F)
          SELECT CASE(IOR)
-            CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F
-            CASE(-1); FX(IIG,JJG,KKG,0)   = B1%RHO_F
-            CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F
-            CASE(-2); FY(IIG,JJG,KKG,0)   = B1%RHO_F
-            CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F
-            CASE(-3); FZ(IIG,JJG,KKG,0)   = B1%RHO_F
+            CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F/MW_F
+            CASE(-1); FX(IIG,JJG,KKG,0)   = B1%RHO_F/MW_F
+            CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F/MW_F
+            CASE(-2); FY(IIG,JJG,KKG,0)   = B1%RHO_F/MW_F
+            CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F/MW_F
+            CASE(-3); FZ(IIG,JJG,KKG,0)   = B1%RHO_F/MW_F
          END SELECT
       ENDIF
 
@@ -231,7 +248,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
                ! ///   II   ///  II+1  |  II+2  | ...
                !                       ^ WALL_INDEX(II+1,+1)
                IF ((UU(II+1,JJ,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(+1)>0)) THEN
-                  Z_TEMP(0:3,1,1) = (/RHOP(II+1,JJ,KK),RHOP(II+1:II+2,JJ,KK),DUMMY/)
+                  Z_TEMP(0:3,1,1) = (/RHO_RMW(II+1,JJ,KK),RHO_RMW(II+1:II+2,JJ,KK),DUMMY/)
                   U_TEMP(1,1,1) = UU(II+1,JJ,KK)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
                   FX(II+1,JJ,KK,0) = F_TEMP(1,1,1)
@@ -241,35 +258,35 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
                ! ... |  II-2  |  II-1  ///   II   ///
                !              ^ WALL_INDEX(II-1,-1)
                IF ((UU(II-2,JJ,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II-1,JJ,KK))%WALL_INDEX(-1)>0)) THEN
-                  Z_TEMP(0:3,1,1) = (/DUMMY,RHOP(II-2:II-1,JJ,KK),RHOP(II-1,JJ,KK)/)
+                  Z_TEMP(0:3,1,1) = (/DUMMY,RHO_RMW(II-2:II-1,JJ,KK),RHO_RMW(II-1,JJ,KK)/)
                   U_TEMP(1,1,1) = UU(II-2,JJ,KK)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
                   FX(II-2,JJ,KK,0) = F_TEMP(1,1,1)
                ENDIF
             CASE( 2) OFF_WALL_SELECT_3
                IF ((VV(II,JJ+1,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(+2)>0)) THEN
-                  Z_TEMP(1,0:3,1) = (/RHOP(II,JJ+1,KK),RHOP(II,JJ+1:JJ+2,KK),DUMMY/)
+                  Z_TEMP(1,0:3,1) = (/RHO_RMW(II,JJ+1,KK),RHO_RMW(II,JJ+1:JJ+2,KK),DUMMY/)
                   U_TEMP(1,1,1) = VV(II,JJ+1,KK)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
                   FY(II,JJ+1,KK,0) = F_TEMP(1,1,1)
                ENDIF
             CASE(-2) OFF_WALL_SELECT_3
                IF ((VV(II,JJ-2,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ-1,KK))%WALL_INDEX(-2)>0)) THEN
-                  Z_TEMP(1,0:3,1) = (/DUMMY,RHOP(II,JJ-2:JJ-1,KK),RHOP(II,JJ-1,KK)/)
+                  Z_TEMP(1,0:3,1) = (/DUMMY,RHO_RMW(II,JJ-2:JJ-1,KK),RHO_RMW(II,JJ-1,KK)/)
                   U_TEMP(1,1,1) = VV(II,JJ-2,KK)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
                   FY(II,JJ-2,KK,0) = F_TEMP(1,1,1)
                ENDIF
             CASE( 3) OFF_WALL_SELECT_3
                IF ((WW(II,JJ,KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(+3)>0)) THEN
-                  Z_TEMP(1,1,0:3) = (/RHOP(II,JJ,KK+1),RHOP(II,JJ,KK+1:KK+2),DUMMY/)
+                  Z_TEMP(1,1,0:3) = (/RHO_RMW(II,JJ,KK+1),RHO_RMW(II,JJ,KK+1:KK+2),DUMMY/)
                   U_TEMP(1,1,1) = WW(II,JJ,KK+1)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
                   FZ(II,JJ,KK+1,0) = F_TEMP(1,1,1)
                ENDIF
             CASE(-3) OFF_WALL_SELECT_3
                IF ((WW(II,JJ,KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK-1))%WALL_INDEX(-3)>0)) THEN
-                  Z_TEMP(1,1,0:3) = (/DUMMY,RHOP(II,JJ,KK-2:KK-1),RHOP(II,JJ,KK-1)/)
+                  Z_TEMP(1,1,0:3) = (/DUMMY,RHO_RMW(II,JJ,KK-2:KK-1),RHO_RMW(II,JJ,KK-1)/)
                   U_TEMP(1,1,1) = WW(II,JJ,KK-2)
                   CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
                   FZ(II,JJ,KK-2,0) = F_TEMP(1,1,1)
@@ -281,20 +298,30 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
    ENDDO WALL_LOOP_3
    !$OMP END PARALLEL DO
 
-   ! Now correct the face value of (RHO*ZZ) such that SUM(RHO*ZZ)_FACE = RHO_FACE
+   ! Now correct the max face value of (RHO*ZZ) such that SUM(RHO*ZZ/MW)_FACE = RHO_FACE/MW_FACE
+   ! (necessary condition to preserve isothermal flow)
 
-   !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
+   !$OMP PARALLEL DO PRIVATE(N,MW_G) SCHEDULE(STATIC)
    DO K=0,KBAR
       DO J=0,JBAR
          DO I=0,IBAR
             N=MAXLOC(FX(I,J,K,1:N_TRACKED_SPECIES),1)
-            FX(I,J,K,N) = MAX( 0._EB, FX(I,J,K,0) - SUM(FX(I,J,K,1:(N-1))) - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+            MW_G = SPECIES_MIXTURE(N)%MW
+            FX(I,J,K,N) = MW_G*MAX( 0._EB, FX(I,J,K,0) &
+                                           - SUM(FX(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
+                                           - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )
 
             N=MAXLOC(FY(I,J,K,1:N_TRACKED_SPECIES),1)
-            FY(I,J,K,N) = MAX( 0._EB, FY(I,J,K,0) - SUM(FY(I,J,K,1:(N-1))) - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+            MW_G = SPECIES_MIXTURE(N)%MW
+            FY(I,J,K,N) = MW_G*MAX( 0._EB, FY(I,J,K,0) &
+                                           - SUM(FY(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
+                                           - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )
 
             N=MAXLOC(FZ(I,J,K,1:N_TRACKED_SPECIES),1)
-            FZ(I,J,K,N) = MAX( 0._EB, FZ(I,J,K,0) - SUM(FZ(I,J,K,1:(N-1))) - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+            MW_G = SPECIES_MIXTURE(N)%MW
+            FZ(I,J,K,N) = MW_G*MAX( 0._EB, FZ(I,J,K,0) &
+                                           - SUM(FZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
+                                           - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )
          ENDDO
       ENDDO
    ENDDO
