Skip to content

FDS Source: Rename PR variable to avoid confusion with Prandtl number #14943

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
cxp484 merged 1 commit into from
Aug 11, 2025
Merged

FDS Source: Rename PR variable to avoid confusion with Prandtl number #14943

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
52 changes: 26 additions & 26 deletions Source/chem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,15 +98,15 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
TYPE(USERDATA), INTENT(IN):: USER_DATA

REAL(EB) :: R_F,MIN_SPEC(N_TRACKED_SPECIES), KG, TMP, RHO, &
K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PR
K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PRES
INTEGER :: I,NS, ITMP
REAL(EB) :: ZZ(N_TRACKED_SPECIES), CP, HS_I, DG, TMPI
REAL(EB) :: ZETA, MIXING_FACTOR, VOL_CHANGE_TERM, SUM_OMEGA_DOT, SUM_CC, MW0, MW, SUM_H_ZZ0, EXPONENT, CONC
TYPE(REACTION_TYPE), POINTER :: RN

TMP = MAX(CVEC(N_TRACKED_SPECIES+1), MIN_CHEM_TMP)
TMPI = 1._EB/TMP
PR = CVEC(N_TRACKED_SPECIES+2) ! PA
PRES = CVEC(N_TRACKED_SPECIES+2) ! PA
RRTMP = 1._EB/(R0*TMP)
ZETA = USER_DATA%ZETA0*EXP(-TN/USER_DATA%TAU_MIX)
MIXING_FACTOR = 0._EB
Expand All @@ -118,7 +118,7 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
CALL MOLAR_CONC_TO_MASS_FRAC(CVEC(1:N_TRACKED_SPECIES), ZZ(1:N_TRACKED_SPECIES))
CALL GET_MOLECULAR_WEIGHT(ZZ(1:N_TRACKED_SPECIES),MW)
ELSE ! Enters at the first timestep when ZETA0=1.
RHO = PR*MW0/R0/TMP
RHO = PRES*MW0/R0/TMP
MW = MW0
ZZ(1:N_TRACKED_SPECIES) = USER_DATA%ZZ_0(1:N_TRACKED_SPECIES)
ENDIF
Expand Down Expand Up @@ -160,7 +160,7 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
IF (RN%N_THIRD > 0) THEN
THIRD_BODY_ENHANCEMENT = DOT_PRODUCT(CVEC(1:N_SPECIES),RN%THIRD_EFF(1:N_SPECIES))
ELSE
THIRD_BODY_ENHANCEMENT = PR*RRTMP
THIRD_BODY_ENHANCEMENT = PRES*RRTMP
ENDIF

IF (RN%REACTYPE==THREE_BODY_ARRHENIUS_TYPE) THEN
Expand Down Expand Up @@ -343,7 +343,7 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)
TYPE(USERDATA), INTENT(IN):: USER_DATA

REAL(EB) :: R_F,DCVEC1,DCVEC2, MIN_SPEC(N_TRACKED_SPECIES), KG, TMP, RHO, &
K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PR
K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PRES
REAL(EB) :: ZZ(N_TRACKED_SPECIES), CP_I(N_TRACKED_SPECIES), HS_I(N_TRACKED_SPECIES)
REAL(EB) :: DKCDTBYKC, DBIDC(N_TRACKED_SPECIES), DBIDT, CP, DCPDT, DKINFDTMPBYKINF, DTMPDT, DG, TMPI, RHOI, CPI
REAL(EB) :: ZETA, MIXING_FACTOR, SUM_OMEGA_DOT, SUM_CC, SUM_OMEGA_DOT_BY_CC, SUM_CC_I, ENRG_TERM, DUMMY1, DUMMY2, DUMMY3
Expand All @@ -354,7 +354,7 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)

TMP = MAX(CVEC(N_TRACKED_SPECIES+1), MIN_CHEM_TMP)
TMPI = 1._EB/TMP
PR = CVEC(N_TRACKED_SPECIES+2) ! PA
PRES = CVEC(N_TRACKED_SPECIES+2) ! PA
RRTMP = 1._EB/(R0*TMP)
ZETA = USER_DATA%ZETA0*EXP(-TN/USER_DATA%TAU_MIX)
MIXING_FACTOR = 0._EB
Expand All @@ -366,7 +366,7 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)
CALL MOLAR_CONC_TO_MASS_FRAC(CVEC(1:N_TRACKED_SPECIES), ZZ(1:N_TRACKED_SPECIES))
CALL GET_MOLECULAR_WEIGHT(ZZ(1:N_TRACKED_SPECIES),MW)
ELSE ! Enters at the first timestep when ZETA0=1.
RHO = PR*MW0/R0/TMP
RHO = PRES*MW0/R0/TMP
MW = MW0
ZZ(1:N_TRACKED_SPECIES) = USER_DATA%ZZ_0(1:N_TRACKED_SPECIES)
ENDIF
Expand Down Expand Up @@ -409,7 +409,7 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)
IF (RN%N_THIRD > 0) THEN
THIRD_BODY_ENHANCEMENT = DOT_PRODUCT(CVEC(1:N_SPECIES),RN%THIRD_EFF(1:N_SPECIES))
ELSE
THIRD_BODY_ENHANCEMENT = PR*RRTMP
THIRD_BODY_ENHANCEMENT = PRES*RRTMP
ENDIF
IF (RN%REACTYPE==THREE_BODY_ARRHENIUS_TYPE) THEN
R_F = R_F * THIRD_BODY_ENHANCEMENT
Expand Down Expand Up @@ -661,13 +661,13 @@ END SUBROUTINE PRINT_JMAT
!> \param RNI is the reaction index.
!> \param K0 is the low pressure rate coeff.
!> \param KINF is the high pressure rate coeff.
!> \param PR is the pressure ratio.
!> \param P_RATIO is the pressure ratio.
!> \param F is the falloff function value.
!> \param DBIDC is the derivative of modification factor w.r.t concentration (out).
!> \param DBIDT is the derivative of modification factor w.r.t temperature (out).

SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT(TMP, RNI, K0, KINF, PR, F, DBIDC, DBIDT)
REAL(EB), INTENT(IN) :: TMP, PR, K0, KINF, F
SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT(TMP, RNI, K0, KINF, P_RATIO, F, DBIDC, DBIDT)
REAL(EB), INTENT(IN) :: TMP, P_RATIO, K0, KINF, F
INTEGER, INTENT(IN) :: RNI
REAL(EB), INTENT(INOUT) :: DBIDC(N_TRACKED_SPECIES)
REAL(EB), INTENT(INOUT) :: DBIDT
Expand All @@ -682,37 +682,37 @@ SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT(TMP, RNI, K0, KINF, PR, F, DBIDC, DBIDT)
DPRDBI = -RN%THIRD_EFF(NS )*K0/KINF
DFDBI = 0._EB
IF (RN%REACTYPE==FALLOFF_TROE_TYPE) THEN
DFDBI = DDC_TROE(PR, F, DPRDBI, TMP, RNI)
DFDBI = DDC_TROE(P_RATIO, F, DPRDBI, TMP, RNI)
ENDIF
DBIDC(NS) = (DPRDBI/(PR*(1 + PR)) + DFDBI/F)
DBIDC(NS) = (DPRDBI/(P_RATIO*(1 + P_RATIO)) + DFDBI/F)
ENDDO


DPRDT = PR/TMP*( RN%N_T_LOW_PR + RN%E_LOW_PR*RRTMP - RN%N_T - RN%E*RRTMP - 1)
DPRDT = P_RATIO/TMP*( RN%N_T_LOW_PR + RN%E_LOW_PR*RRTMP - RN%N_T - RN%E*RRTMP - 1)
DFDT = 0._EB
IF (RN%REACTYPE==FALLOFF_TROE_TYPE) THEN
DFDT = DDTMP_TROE(PR, F, DPRDT, TMP, RNI)
DFDT = DDTMP_TROE(P_RATIO, F, DPRDT, TMP, RNI)
ENDIF
DBIDT = (DPRDT/(PR*(1 + PR)) + DFDT/F)
DBIDT = (DPRDT/(P_RATIO*(1 + P_RATIO)) + DFDT/F)

RETURN
END SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT

!> \brief Calculate derivative of TROE function w.r.t concentration
!> \param PR is the pressure ratio.
!> \param P_RATIO is the pressure ratio.
!> \param F is the falloff function value.
!> \param DPRDC is the derivative of TROE function w.r.t concentration.
!> \param TMP is the current temperature.
!> \param RNI is the reaction index
REAL(EB) FUNCTION DDC_TROE(PR, F, DPRDC, TMP, RNI)
REAL(EB), INTENT(IN) :: PR, F, DPRDC, TMP
REAL(EB) FUNCTION DDC_TROE(P_RATIO, F, DPRDC, TMP, RNI)
REAL(EB), INTENT(IN) :: P_RATIO, F, DPRDC, TMP
INTEGER, INTENT(IN) :: RNI
REAL(EB) :: LOGPR, LOGTEN, LOGFCENT, C, N, DLOGPRDC, DPARENTDC
TYPE(REACTION_TYPE), POINTER :: RN
REAL(EB), PARAMETER :: D=0.14_EB

RN => REACTION(RNI)
LOGPR = LOG10(MAX(PR, TWO_EPSILON_EB))
LOGPR = LOG10(MAX(P_RATIO, TWO_EPSILON_EB))
LOGTEN = LOG(10.0)
IF (RN%T2_TROE <-1.E20_EB) THEN
LOGFCENT = LOG10(MAX((1 - RN%A_TROE)*EXP(-TMP*RN%RT3_TROE) + &
Expand All @@ -722,7 +722,7 @@ REAL(EB) FUNCTION DDC_TROE(PR, F, DPRDC, TMP, RNI)
RN%A_TROE*EXP(-TMP*RN%RT1_TROE) + EXP(-RN%T2_TROE/TMP),TWO_EPSILON_EB))
ENDIF

DLOGPRDC = DPRDC/PR/LOGTEN;
DLOGPRDC = DPRDC/P_RATIO/LOGTEN;
C = -0.4_EB - 0.67_EB*LOGFCENT
N = 0.75_EB - 1.27_EB*LOGFCENT

Expand All @@ -736,20 +736,20 @@ END FUNCTION DDC_TROE


!> \brief Calculate derivative of TROE function w.r.t temperature
!> \param PR is the pressure ratio.
!> \param P_RATIO is the pressure ratio.
!> \param F is the falloff function value.
!> \param DPRDT is the derivative of TROE function w.r.t temperature.
!> \param TMP is the current temperature.
!> \param RNI is the reaction index
REAL(EB) FUNCTION DDTMP_TROE(PR, F, DPRDT, TMP, RNI)
REAL(EB), INTENT(IN) :: PR, F, DPRDT, TMP
REAL(EB) FUNCTION DDTMP_TROE(P_RATIO, F, DPRDT, TMP, RNI)
REAL(EB), INTENT(IN) :: P_RATIO, F, DPRDT, TMP
INTEGER, INTENT(IN) :: RNI
REAL(EB) :: FCENT, LOGPR, LOGTEN, LOGFCENT, DFCENTDT, C, N, DCDT, DNDT, DPARENTDT, DLOGFCENTDT, DLOGPRDT
TYPE(REACTION_TYPE), POINTER :: RN
REAL(EB), PARAMETER :: D=0.14_EB

RN => REACTION(RNI)
LOGPR = LOG10(MAX(PR, TWO_EPSILON_EB));
LOGPR = LOG10(MAX(P_RATIO, TWO_EPSILON_EB));
LOGTEN = LOG(10.0);
IF (RN%T2_TROE <-1.E20_EB) THEN
FCENT = MAX((1 - RN%A_TROE)*EXP(-TMP*RN%RT3_TROE) + &
Expand All @@ -768,7 +768,7 @@ REAL(EB) FUNCTION DDTMP_TROE(PR, F, DPRDT, TMP, RNI)
N = 0.75_EB - 1.27_EB*LOGFCENT
DCDT = -0.67*DLOGFCENTDT
DNDT = -1.27*DLOGFCENTDT
DLOGPRDT = DPRDT/PR/LOGTEN
DLOGPRDT = DPRDT/P_RATIO/LOGTEN

DPARENTDT = 2.0*(LOGPR + C)/((N - D*(LOGPR + C))**2)* &
((DLOGPRDT + DCDT) - (LOGPR + C)*(DNDT - D*(DLOGPRDT + DCDT))/(N - D*(LOGPR + C)))
Expand Down
Loading