Skip to content

FDS Source: remove FLUX_LIMITER_MW_CORRECTION #14996

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rmcdermo merged 4 commits into from
Aug 20, 2025
Merged

FDS Source: remove FLUX_LIMITER_MW_CORRECTION #14996

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5f1f3d1
FDS Source: remove FLUX_LIMITER_MW_CORRECTION scheme
rmcdermo Aug 20, 2025
f863e72
FDS User Guide: update user guide to remove FLUX_LIMITER_MW_CORRECTIO…
rmcdermo Aug 20, 2025
552da90
FDS Source: fix accidentally committed paramater
rmcdermo Aug 20, 2025
7a2ea71
FDS Source: fix unused warnings
rmcdermo Aug 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Manuals/FDS_User_Guide/FDS_User_Guide.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11912,7 +11912,7 @@ \chapter{Alphabetical List of Input Parameters}
% ignorenamelistkw: /ISOF/DEBUG
% ignorenamelistkw: /MISC/PERIODIC_TEST, /MISC/POSITIVE_ERROR_TEST, /MISC/PROFILING, /MISC/RADIATION
% ignorenamelistkw: /MISC/STRATIFICATION, /MISC/SUPPRESSION, /MISC/UVW_FILE, /MISC/TENSOR_DIFFUSIVITY
% ignorenamelistkw: /MISC/CC_IBM, /MISC/CCVOL_LINK, /MISC/TEST_NEW_CHAR_MODEL, /MISC/FLUX_LIMITER_MW_CORRECTION
% ignorenamelistkw: /MISC/CC_IBM, /MISC/CCVOL_LINK, /MISC/TEST_NEW_CHAR_MODEL, /MISC/PR, /MISC/SC
% ignorenamelistkw: /PART/DEBUG, /PART/EMBER_SNAG_FACTOR
% ignorenamelistkw: /PRES/WRITE_PARCSRPCG_MATRIX
% ignorenamelistkw: /REAC/C, /REAC/H, /REAC/O, /REAC/N, /REAC/FORMULA
Expand Down
1 change: 0 additions & 1 deletion Source/cons.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,6 @@ MODULE GLOBAL_CONSTANTS
LOGICAL :: TENSOR_DIFFUSIVITY=.FALSE. !< If true, use experimental tensor diffusivity model for spec and tmp
LOGICAL :: OXPYRO_MODEL=.FALSE. !< Flag to use oxidative pyrolysis mass transfer model
LOGICAL :: OUTPUT_WALL_QUANTITIES=.FALSE. !< Flag to force call to WALL_MODEL
LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.FALSE. !< Flag for MW correction ensure consistent equation of state at face
LOGICAL :: FLUX_LIMITER_SINGLE_COEF=.FALSE. !< Flag to base flux coefficients off of a single worst-case scalar gradient
LOGICAL :: STORE_FIRE_ARRIVAL=.FALSE. !< Flag for tracking arrival of spreading fire front over a surface
LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE. !< Flag for tracking residence time of spreading fire front over a surface
Expand Down
126 changes: 2 additions & 124 deletions Source/divg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -973,13 +973,11 @@ SUBROUTINE SPECIES_ADVECTION_PART_1

USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE
USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT
! INTEGER, INTENT(IN) :: N
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),MW_G
REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP
INTEGER :: I,J,K,IW,N !,II,JJ,KK,IOR,IC,IIG,JJG,KKG
INTEGER :: I,J,K,IW,N
TYPE(WALL_TYPE), POINTER :: WC
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
Expand Down Expand Up @@ -1078,126 +1076,6 @@ SUBROUTINE SPECIES_ADVECTION_PART_1

ENDDO SPECIES_LOOP

FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN

! Repeat the above for DENSITY

RHO_RMW=>WORK_PAD

!$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC)
DO K=-1,KBP1+1
DO J=-1,JBP1+1
DO I=-1,IBP1+1
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_ZZ(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER)
CALL GET_SCALAR_FACE_VALUE(VV,RHO_RMW,FY_ZZ(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER)
CALL GET_SCALAR_FACE_VALUE(WW,RHO_RMW,FZ_ZZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER)

!$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,U_TEMP,Z_TEMP,F_TEMP,ZZ_GET)
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
BC=>BOUNDARY_COORD(WC%BC_INDEX)
B1=>BOUNDARY_PROP1(WC%B1_INDEX)

! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell

OFF_WALL_IF_3: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN

OFF_WALL_SELECT_3: SELECT CASE(BC%IOR)
CASE( 1) OFF_WALL_SELECT_3
! ghost FX/UU(II+1)
! /// II /// II+1 | II+2 | ...
! ^ WALL_INDEX(II+1,+1)
IF ((UU(BC%II+1,BC%JJ,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II+1,BC%JJ,BC%KK))%WALL_INDEX(+1)>0)) THEN
Z_TEMP(0:2,1,1) = (/RHO_Z_P(BC%II+1,BC%JJ,BC%KK),RHO_Z_P(BC%II+1:BC%II+2,BC%JJ,BC%KK)/)
U_TEMP(1,1,1) = UU(BC%II+1,BC%JJ,BC%KK)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
FX_ZZ(BC%II+1,BC%JJ,BC%KK,0) = F_TEMP(1,1,1)
ENDIF
CASE(-1) OFF_WALL_SELECT_3
! FX/UU(II-2) ghost
! ... | II-2 | II-1 /// II ///
! ^ WALL_INDEX(II-1,-1)
IF ((UU(BC%II-2,BC%JJ,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II-1,BC%JJ,BC%KK))%WALL_INDEX(-1)>0)) THEN
Z_TEMP(1:3,1,1) = (/RHO_Z_P(BC%II-2:BC%II-1,BC%JJ,BC%KK),RHO_Z_P(BC%II-1,BC%JJ,BC%KK)/)
U_TEMP(1,1,1) = UU(BC%II-2,BC%JJ,BC%KK)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
FX_ZZ(BC%II-2,BC%JJ,BC%KK,0) = F_TEMP(1,1,1)
ENDIF
CASE( 2) OFF_WALL_SELECT_3
IF ((VV(BC%II,BC%JJ+1,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ+1,BC%KK))%WALL_INDEX(+2)>0)) THEN
Z_TEMP(1,0:2,1) = (/RHO_Z_P(BC%II,BC%JJ+1,BC%KK),RHO_Z_P(BC%II,BC%JJ+1:BC%JJ+2,BC%KK)/)
U_TEMP(1,1,1) = VV(BC%II,BC%JJ+1,BC%KK)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
FY_ZZ(BC%II,BC%JJ+1,BC%KK,0) = F_TEMP(1,1,1)
ENDIF
CASE(-2) OFF_WALL_SELECT_3
IF ((VV(BC%II,BC%JJ-2,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ-1,BC%KK))%WALL_INDEX(-2)>0)) THEN
Z_TEMP(1,1:3,1) = (/RHO_Z_P(BC%II,BC%JJ-2:BC%JJ-1,BC%KK),RHO_Z_P(BC%II,BC%JJ-1,BC%KK)/)
U_TEMP(1,1,1) = VV(BC%II,BC%JJ-2,BC%KK)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
FY_ZZ(BC%II,BC%JJ-2,BC%KK,0) = F_TEMP(1,1,1)
ENDIF
CASE( 3) OFF_WALL_SELECT_3
IF ((WW(BC%II,BC%JJ,BC%KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK+1))%WALL_INDEX(+3)>0)) THEN
Z_TEMP(1,1,0:2) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK+1),RHO_Z_P(BC%II,BC%JJ,BC%KK+1:BC%KK+2)/)
U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK+1)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
FZ_ZZ(BC%II,BC%JJ,BC%KK+1,0) = F_TEMP(1,1,1)
ENDIF
CASE(-3) OFF_WALL_SELECT_3
IF ((WW(BC%II,BC%JJ,BC%KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK-1))%WALL_INDEX(-3)>0)) THEN
Z_TEMP(1,1,1:3) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK-2:BC%KK-1),RHO_Z_P(BC%II,BC%JJ,BC%KK-1)/)
U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-2)
CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
FZ_ZZ(BC%II,BC%JJ,BC%KK-2,0) = F_TEMP(1,1,1)
ENDIF
END SELECT OFF_WALL_SELECT_3

ENDIF OFF_WALL_IF_3

ENDDO WALL_LOOP_3
!$OMP END PARALLEL DO

! 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,MW_G) SCHEDULE(STATIC)
DO K=0,KBAR
DO J=0,JBAR
DO I=0,IBAR
N=MAXLOC(FX_ZZ(I,J,K,1:N_TRACKED_SPECIES),1)
MW_G = SPECIES_MIXTURE(N)%MW
FX_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FX_ZZ(I,J,K,0) &
- SUM(FX_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
- SUM(FX_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )

N=MAXLOC(FY_ZZ(I,J,K,1:N_TRACKED_SPECIES),1)
MW_G = SPECIES_MIXTURE(N)%MW
FY_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FY_ZZ(I,J,K,0) &
- SUM(FY_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
- SUM(FY_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )

N=MAXLOC(FZ_ZZ(I,J,K,1:N_TRACKED_SPECIES),1)
MW_G = SPECIES_MIXTURE(N)%MW
FZ_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FZ_ZZ(I,J,K,0) &
- SUM(FZ_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) &
- SUM(FZ_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

ENDIF FACE_CORRECTION_IF

END SUBROUTINE SPECIES_ADVECTION_PART_1


Expand Down
22 changes: 7 additions & 15 deletions Source/init.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -597,17 +597,9 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)

! 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( 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

IF (FLUX_LIMITER_SINGLE_COEF) THEN
ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB
Expand Down Expand Up @@ -702,9 +694,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
Expand Down
155 changes: 2 additions & 153 deletions Source/mass.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,14 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT

INTEGER, INTENT(IN) :: NM
REAL(EB) :: TNOW,MW_F,MW_G,ZZ_GET(1:N_TRACKED_SPECIES)
REAL(EB) :: TNOW
REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_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,RHO_RMW
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P
TYPE(WALL_TYPE), POINTER :: WC
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
Expand Down Expand Up @@ -179,157 +179,6 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)

ENDDO SPECIES_LOOP

FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN

! Repeat the above for DENSITY

RHO_RMW=>WORK_PAD

!$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC)
DO K=-1,KBP1+1
DO J=-1,JBP1+1
DO I=-1,IBP1+1
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
BC=>BOUNDARY_COORD(WC%BC_INDEX)
B1=>BOUNDARY_PROP1(WC%B1_INDEX)

II = BC%II
JJ = BC%JJ
KK = BC%KK
IIG = BC%IIG
JJG = BC%JJG
KKG = BC%KKG
IOR = BC%IOR
IC = CELL_INDEX(II,JJ,KK)

IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR) THEN
SELECT CASE(IOR)
CASE( 1); FX(IIG-1,JJG,KKG,0) = 0._EB
CASE(-1); FX(IIG,JJG,KKG,0) = 0._EB
CASE( 2); FY(IIG,JJG-1,KKG,0) = 0._EB
CASE(-2); FY(IIG,JJG,KKG,0) = 0._EB
CASE( 3); FZ(IIG,JJG,KKG-1,0) = 0._EB
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/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

! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell

OFF_WALL_IF_3: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN

OFF_WALL_SELECT_3: SELECT CASE(IOR)
CASE( 1) OFF_WALL_SELECT_3
! ghost FX/UU(II+1)
! /// 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) = (/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)
ENDIF
CASE(-1) OFF_WALL_SELECT_3
! FX/UU(II-2) ghost
! ... | 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,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) = (/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,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) = (/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,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)
ENDIF
END SELECT OFF_WALL_SELECT_3

ENDIF OFF_WALL_IF_3

ENDDO WALL_LOOP_3
!$OMP END PARALLEL DO

! 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,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)
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)
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)
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
!$OMP END PARALLEL DO

ENDIF FACE_CORRECTION_IF

T_USED(3)=T_USED(3)+CURRENT_TIME()-TNOW
END SUBROUTINE MASS_FINITE_DIFFERENCES

Expand Down
Loading
Loading