Skip to content

FDS Source: update to TEST_FLUX_LIMITER_FACE_CORRECTION to account fo… #13854

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 1 commit into from
Dec 6, 2024
Merged

FDS Source: update to TEST_FLUX_LIMITER_FACE_CORRECTION to account fo… #13854

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
73 changes: 50 additions & 23 deletions Source/mass.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
Loading