Merge pull request #14804 from cxp484/master

FDS Source: Handle reaction exponent <1.

Author: cxp484

Source/chem.f90

@@ -101,7 +101,7 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
           K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PR
 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
+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)
@@ -144,7 +144,14 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
    ! MULTIPLY WITH MOLAR_CONCENTRATION ^ STOICHIOMETRIC_COEFF
    DO NS=1,RN%N_SPEC
       IF (CVEC(YP2ZZ(RN%N_S_INDEX(NS))) < MIN_SPEC(YP2ZZ(RN%N_S_INDEX(NS)))) CYCLE REACTION_LOOP
-      R_F = R_F*(CVEC(YP2ZZ(RN%N_S_INDEX(NS))))**RN%N_S(NS) 
+      EXPONENT = RN%N_S(NS)
+      CONC = CVEC(YP2ZZ(RN%N_S_INDEX(NS)))
+      IF (EXPONENT < 1._EB) THEN
+         R_F = R_F * CONC**(EXPONENT - 1._EB)
+         R_F = R_F * CONC
+      ELSE
+          R_F = R_F * CONC**EXPONENT
+      END IF
    ENDDO
 
    ! CALCULATE B_I BASED ON TYPE OF REACTION
@@ -207,7 +214,6 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC, TN, USER_DATA)
 ! PRESSURE DERIVATIVE (CONSTANT PRESSURE ASSUMPTION)
 FVEC(N_TRACKED_SPECIES+2) = 0._EB
 
-
 END SUBROUTINE DERIVATIVE
 
 
@@ -336,11 +342,11 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,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
-REAL(EB) :: ZZ(N_TRACKED_SPECIES), CP_I(N_TRACKED_SPECIES), HS_I(N_TRACKED_SPECIES)    
+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
 REAL(EB) :: VOL_CHANGE_TERM1, VOL_CHANGE_TERM2, VOL_CHANGE_TERM3, VOL_CHANGE_TERM, SUM_DOMEGA_DOT_BY_DT, MW0, MW, &
-            SUM_CP_ZZ0,SUM_H_ZZ0
+            SUM_CP_ZZ0,SUM_H_ZZ0, EXPONENT, CONC, CONC_EXP
 INTEGER :: I,NS, NS1, NS2, ITMP
 TYPE(REACTION_TYPE), POINTER :: RN
 
@@ -391,10 +397,8 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)
       DKCDTBYKC = ((RN%DELTA_G(MIN(I_MAX_TEMP,NINT(TMP)))*TMPI + RN%DELTA_S(MIN(I_MAX_TEMP,NINT(TMP))))+RN%C0_EXP)*TMPI 
    ENDIF
 
-   ! MULTIPLY WITH MOLAR_CONCENTRATION ^ STOICHIOMETRIC_COEFF
    DO NS=1,RN%N_SPEC
       IF (CVEC(YP2ZZ(RN%N_S_INDEX(NS))) < MIN_SPEC(YP2ZZ(RN%N_S_INDEX(NS)))) CYCLE REACTION_LOOP
-      R_F = R_F*(CVEC(YP2ZZ(RN%N_S_INDEX(NS))))**RN%N_S(NS)
    ENDDO
 
 
@@ -431,13 +435,41 @@ SUBROUTINE JACOBIAN(CVEC,FVEC,JMAT,TN,USER_DATA)
 
    !Contribution of qi
    DO NS1 = 1, RN%N_SPEC
+      ! MULTIPLY WITH MOLAR_CONCENTRATION ^ STOICHIOMETRIC_COEFF
+      CONC_EXP = 1.0_EB
+      DO NS2=1,RN%N_SPEC
+         EXPONENT = RN%N_S(NS2)
+         CONC = CVEC(YP2ZZ(RN%N_S_INDEX(NS2)))
+         IF (NS2 == NS1) THEN
+            CONC_EXP = CONC_EXP * MERGE((CONC+TWO_EPSILON_EB)**(EXPONENT - 1._EB), CONC**(EXPONENT - 1._EB), EXPONENT < 1._EB)
+         ELSE
+            IF (EXPONENT < 1._EB) THEN
+               CONC_EXP = CONC_EXP * CONC * (CONC+TWO_EPSILON_EB)**(EXPONENT - 1._EB)
+            ELSE
+               CONC_EXP = CONC_EXP * CONC**EXPONENT
+            END IF
+         ENDIF
+      ENDDO
+
       DO NS=1,RN%N_SMIX_FR
-         DCVEC1 = R_F*RN%NU_NN(RN%NU_INDEX(NS))*RN%N_S(NS1)/CVEC(YP2ZZ(RN%N_S_INDEX(NS1)))
+         DCVEC1 = R_F*CONC_EXP*RN%NU_NN(RN%NU_INDEX(NS))*RN%N_S(NS1)
          JMAT((YP2ZZ(RN%N_S_INDEX(NS1))),RN%NU_INDEX(NS)) = &
          JMAT((YP2ZZ(RN%N_S_INDEX(NS1))),RN%NU_INDEX(NS))+ DCVEC1 
       ENDDO   
    ENDDO 
 
+   ! CALCULATE THE REACTION RATE
+   DO NS=1,RN%N_SPEC
+      EXPONENT = RN%N_S(NS)
+      CONC = CVEC(YP2ZZ(RN%N_S_INDEX(NS)))
+      IF (EXPONENT < 1._EB) THEN
+         R_F = R_F * CONC**(EXPONENT - 1._EB)
+         R_F = R_F * CONC
+      ELSE
+          R_F = R_F * CONC**EXPONENT
+      END IF
+   ENDDO
+
    ! Add contribution of C_I
    IF (RN%THIRD_BODY) THEN
       DO NS = 1,N_TRACKED_SPECIES
@@ -792,6 +824,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 REAL(C_DOUBLE)                 :: CVEC_C(N_TRACKED_SPECIES+2)    ! N_SP + 2 (FOR TEMPERATURE AND PRESSURE)
 REAL(C_DOUBLE)                 :: ATOLVEC_C(N_TRACKED_SPECIES+2) ! N_SP + 2 
 INTEGER(C_LONG)                :: MAXSTEPS_C     ! MAXIMUM NUMBER OF INTERNAL STEPS
+INTEGER(C_INT)                 :: MAXORD_C       ! Maximum number of order.
 INTEGER(C_INT64_T)             :: NEQ
 REAL(C_DOUBLE)                 :: CHEM_TIME_C(1)    ! OUTPUT CHEMICAL TIME
 
@@ -890,6 +923,15 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
    STOP 1
 END IF
 
+! SET MAX ORDER
+MAXORD_C = 5
+IF (IS_EXPONENT_LT_1) MAXORD_C = 1
+IERR_C = FCVODESETMAXORD(CVODE_MEM, MAXORD_C)
+IF (IERR_C /= 0) THEN
+   WRITE(LU_ERR,*) 'ERROR IN FCVODESETMAXORD, IERR = ', IERR_C, '; HALTING'
+   STOP 1
+END IF
+
 ! SET ERROR HANDLER
 IERR_C = FCVODESETERRHANDLERFN(CVODE_MEM, C_FUNLOC(FDS_CVODE_ERR_HANDLER), C_NULL_PTR)
 IF (IERR_C /= 0) THEN
@@ -920,7 +962,8 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
    ONLY_FIRST_STEP = .TRUE.
    IF (WRITE_SUBSTEPS) THEN ! This WRITE_SUBSTEPS is only true for few verification cases.
       ONLY_FIRST_STEP = .FALSE.
-      ALLOCATE(CVODE_SUBSTEP_DATA((CVODE_MAX_TRY+1)*MAX_CVODE_SUBSTEPS, N_TRACKED_SPECIES+4))
+      IF (.NOT. ALLOCATED(CVODE_SUBSTEP_DATA)) &
+          ALLOCATE(CVODE_SUBSTEP_DATA((CVODE_MAX_TRY+1)*MAX_CVODE_SUBSTEPS, N_TRACKED_SPECIES+4))
    ENDIF
 
    SUBSTEP_COUNT = 0
@@ -941,11 +984,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
       ENDIF  
    END DO
 
-   IF (WRITE_SUBSTEPS) THEN
-      TOTAL_SUBSTEPS_TAKEN = SUBSTEP_COUNT
-      STOP_STATUS=CVODE_SUBSTEP_STOP
-      RETURN
-   ENDIF   
+
 ENDIF
 
 IF (IERR_C /= 0) THEN
@@ -999,6 +1038,12 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 CALL FN_VDESTROY(SUNATOL)
 IERR_C = FSUNCONTEXT_FREE(SUNCTX)
 
+IF (WRITE_SUBSTEPS) THEN
+   TOTAL_SUBSTEPS_TAKEN = SUBSTEP_COUNT
+   STOP_STATUS=CVODE_SUBSTEP_STOP
+   RETURN
+ENDIF
+
 
 END SUBROUTINE CVODE_SERIAL
 

Source/cons.f90

@@ -894,6 +894,7 @@ MODULE CHEMCONS
 INTEGER  :: MAX_CVODE_SUBSTEPS=100000
 REAL(EB) :: MAX_CHEM_TIME=1.E-6_EB
 INTEGER  :: CVODE_MAX_TRY=4
+LOGICAL  :: IS_EXPONENT_LT_1 = .FALSE.
 
 ! FOR WRITING CVODE SUBSTEPS
 LOGICAL  :: WRITE_CVODE_SUBSTEPS = .FALSE.

Source/read.f90

@@ -5061,7 +5061,7 @@ END SUBROUTINE READ_REAC
 
 
 SUBROUTINE PROC_REAC_1
-USE CHEMCONS, ONLY: YP2ZZ
+USE CHEMCONS, ONLY: YP2ZZ, IS_EXPONENT_LT_1
 USE PROPERTY_DATA, ONLY : PARSE_EQUATION, SHUTDOWN_ATOM
 REAL(EB) :: MASS_PRODUCT,MASS_REACTANT,REACTION_BALANCE(118),NU_Y(N_SPECIES)
 INTEGER :: NS,NS2, NR
@@ -5284,6 +5284,15 @@ SUBROUTINE PROC_REAC_1
       ENDIF
    ENDDO
 
+   IF(TRIM(ODE_SOLVER)=='CVODE' .AND. .NOT. IS_EXPONENT_LT_1) THEN
+      DO NS=1,RN%N_SPEC
+         IF (RN%N_S(NS) .LT. 1) THEN
+            IS_EXPONENT_LT_1 = .TRUE.
+            WRITE(LU_ERR,*)"Info: CVODE solver order set to 1 because one of the reaction involves exponent < 1.0. "
+            EXIT
+         ENDIF
+      ENDDO
+   ENDIF
 
    ! Normalize the stoichiometric coefficients by that of the fuel.
    RN%NU_FUEL_0 = -RN%NU(RN%FUEL_SMIX_INDEX)