Merge pull request #15538 from cxp484/FireX

FireX: Merge with firemodels/master

Author: cxp484

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -9140,7 +9140,8 @@ \section{Flux Limiters}
 Central differencing                & \ct{'CENTRAL'}  \\
 Godunov                             & \ct{'GODUNOV'}  \\
 Superbee (VLES, SVLES default)      & \ct{'SUPERBEE'}  \\
-CHARM (DNS, LES default)            & \ct{'CHARM'}  \\
+CHARM (DNS, LES default)            & \ct{'CHARM'}  \\ 
+\hline
 \end{tabular}
 \end{table}
 
@@ -9171,9 +9172,9 @@ \subsection{Density}
 
 The density of the gas has a natural lower bound of zero, but if the density in a cell decreases to nearly zero, the temperature would then increase to an extremely high value due to the equation of state. Thus, by default, the density is kept within the following range:
 \be
-   \min \left( 0.1 \, \rho_\infty , \frac{2 \, p_\infty \, W_{\min}}{\R \, T_{\max}} \right) < \rho <  \frac{2 \, p_\infty \, W_{\max}}{\R \, T_{\min}}
+   0.01 \; \mbox{kg/m$^3$} < \rho < 100 \; \mbox{kg/m$^3$}
 \ee
-where $W_{\min}$ and $W_{\max}$ are the minimum and maximum values of the molecular weight of the tracked gas species in units of g/mol, and $\R$ is the universal gas constant, 8314.5~J/(kmol$\cdot$K). $T_{\min}$ and $T_{\max}$ are described above. To override the limits of density, specify \ct{MINIMUM_DENSITY} and/or \ct{MAXIMUM_DENSITY} on the \ct{CLIP} line in units of kg/m$^3$.
+To override the limits of density, specify \ct{MINIMUM_DENSITY} and/or \ct{MAXIMUM_DENSITY} on the \ct{CLIP} line in units of kg/m$^3$.
 
 Clipping of density and mass fractions violates mass conservation, so it is preferable to avoid clipping if possible.  As discussed in Sec.~\ref{info:Stability_Constraints}, the time step is set to adhere to the CFL constraints of the flow field.  The proper \ct{DT} combined with flux limiters generally avoids the need for clipping.  Beyond this, FDS then employs a mass redistribution scheme, as discussed in FDS Technical Guide~\cite{FDS_Math_Guide}.  If this fails, there is yet one more attempt to avoid clipping---the time step is decreased by 10~\% ($\dt_{\mathrm{new}} = 0.9 \,\dt$) and the scalar transport equations are reiterated.  This process is carried out a maximum of \ct{CLIP_DT_RESTRICTIONS_MAX} times; the default is 5.  In some very extreme circumstances, this loop can drive the time step into numerical instability range ($\dt/\dt_{\mathrm{init}} < \mbox{\ct{LIMITING_DT_RATIO}}$).  You can control the max number of time step restrictions by setting the parameter \ct{CLIP_DT_RESTRICTIONS_MAX} on the \ct{CLIP} line (set to 0 to bypass the algorithm altogether). The number of restrictions (if any) is noted in the \ct{CHID.out} file for a given time step.
 

Source/chem.f90

@@ -807,6 +807,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
                         TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION)
 USE PHYSICAL_FUNCTIONS, ONLY : MOLAR_CONC_TO_MASS_FRAC, CALC_EQUIV_RATIO, GET_ENTHALPY, GET_MOLECULAR_WEIGHT
 USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER
+USE CHEMCONS, ONLY: CVODE_WARNING_CELLS,CVODE_ERR_CODE_MIN,CVODE_ERR_CODE_MAX
 USE GLOBAL_CONSTANTS
 USE FCVODE_MOD                 ! FORTRAN INTERFACE TO CVODE
 USE FSUNDIALS_CONTEXT_MOD      ! FORTRAN INTERFACE TO SUNCONTEXT
@@ -1015,15 +1016,18 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
    ENDIF
 
    IF (IERR_C .NE. CV_SUCCESS) THEN
-      IF (IERR_C == CV_TOO_MUCH_WORK) THEN
-         WRITE(LU_ERR,'(A, 2E18.8, I8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT, MAXTRY=", CUR_CFD_TIME, &
-                     (TEND-TCUR), MAXTRY, ". If the warning persists, reduce the timestep."
-      ELSE
-         WRITE(LU_ERR,'(A, I4, A, 2E18.8, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, &
-                        " and CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), ". If the warning persists, reduce the timestep."
+      IF(IERR_C>=CVODE_ERR_CODE_MIN .AND. IERR_C<=CVODE_ERR_CODE_MAX) THEN
+         CVODE_WARNING_CELLS(IERR_C) = CVODE_WARNING_CELLS(IERR_C) + 1
       ENDIF   
-
       IF (DEBUG) THEN
+         IF (IERR_C == CV_TOO_MUCH_WORK) THEN
+            WRITE(LU_ERR,'(A, 2E18.8, I8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT, MAXTRY=", CUR_CFD_TIME, &
+                        (TEND-TCUR), MAXTRY, ". If the warning persists, reduce the timestep."
+         ELSE
+            WRITE(LU_ERR,'(A, I4, A, 2E18.8, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, &
+                           " and CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), ". If the warning persists, reduce the timestep."
+         ENDIF
+
          CALL MOLAR_CONC_TO_MASS_FRAC(CC(1:N_TRACKED_SPECIES), ZZ(1:N_TRACKED_SPECIES))
          CALL CALC_EQUIV_RATIO(ZZ(1:N_TRACKED_SPECIES), EQUIV)
          DO NS = 1, N_TRACKED_SPECIES

Source/cons.f90

@@ -968,6 +968,10 @@ MODULE CHEMCONS
 INTEGER  :: MAX_CVODE_SUBSTEPS=100000
 INTEGER  :: CVODE_MAX_TRY=4
 INTEGER  :: CVODE_ORDER=0
+INTEGER  :: CVODE_ERR_CODE_MIN=-100
+INTEGER  :: CVODE_ERR_CODE_MAX=100
+INTEGER  :: CVODE_WARNING_CELLS(-100:100)! Index of the array is error code and value is Cell count
+CHARACTER(LEN=100) :: CVODE_WARN_MESSAGES(-100:100)
 
 ! FOR WRITING CVODE SUBSTEPS
 LOGICAL  :: WRITE_CVODE_SUBSTEPS = .FALSE.
@@ -988,4 +992,12 @@ MODULE CHEMCONS
 INTEGER :: N_IGNITION_ZONES = 0
 TYPE(IGNITION_ZONE_TYPE), DIMENSION(MAX_IGNITION_ZONES) :: IGNITION_ZONES !< Coordinates of ignition zones
 
+CONTAINS
+   SUBROUTINE INIT_CVODE_WARN_MESSAGES()
+      CVODE_WARN_MESSAGES = 'CVODE didn''t finish ODE solution with this code.'
+      CVODE_WARN_MESSAGES(-1) = 'CVODE took all internal substeps.'
+      CVODE_WARN_MESSAGES(-3) = 'Minimum step size was reached.'
+      CVODE_WARN_MESSAGES(-4) = 'Convergence test failure.'
+   END SUBROUTINE INIT_CVODE_WARN_MESSAGES
+
 END MODULE CHEMCONS

Source/fire.f90

@@ -41,6 +41,7 @@ SUBROUTINE COMBUSTION_LOAD_BALANCED(T,DT)
 
 USE SOOT_ROUTINES, ONLY: SOOT_SURFACE_OXIDATION
 USE COMP_FUNCTIONS, ONLY: CURRENT_TIME
+USE CHEMCONS, ONLY: CVODE_WARNING_CELLS,INIT_CVODE_WARN_MESSAGES
 REAL(EB), INTENT(IN) :: T,DT
 INTEGER :: NM,ICC,JCC
 REAL(EB) :: TNOW
@@ -55,12 +56,14 @@ SUBROUTINE COMBUSTION_LOAD_BALANCED(T,DT)
    COMBUSTION_INIT = .TRUE.
    NVAR_TO_SEND = N_TRACKED_SPECIES +7 ! NS, TEMP, RHO, PRES, MU, DELTA, VOL, IGN_ZN
    NVAR_TO_RECEIVE = N_TRACKED_SPECIES +4 ! NS, Q_OUT, MIX_TIME_OUT, CHI_R_OUT, CHEM_SUBIT_TMP_OUT
+   CALL INIT_CVODE_WARN_MESSAGES()
 ENDIF
 
 DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
    CALL POINT_TO_MESH(NM)
    Q     = 0._EB
    CHI_R = 0._EB
+   CVODE_WARNING_CELLS = 0
    IF (CC_IBM) THEN
       DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
          DO JCC=1,CUT_CELL(ICC)%NCELL
@@ -109,6 +112,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT)
 REAL(EB) :: MIX_TIME_OUT, Q_OUT, CHI_R_OUT, MYTEMP, MYRHO, MYMU, DELTA, VOL, TNOW2
 INTEGER :: IGN_ZN
 
+
 Q_EXISTS =  .FALSE.
 
 !------
@@ -329,6 +333,13 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT)
    ENDDO
 ENDIF   
 
+! Output accumulated CVODE warnings (instead of cell by cell) in verbose mode
+#ifdef WITH_SUNDIALS
+IF (COMBUSTION_ODE_SOLVER == CVODE_SOLVER .AND. VERBOSE) THEN
+   CALL DUMP_CVODE_WARNING_SUMMARY(NCHEM_ACTIVE_CELLS_AND_CC)
+ENDIF
+#endif
+
 END SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED
 
 
@@ -1271,6 +1282,60 @@ SUBROUTINE CALC_AFT_REAC_AND_PROD(ZZ,ZZ_REAC,ZZ_PROD)
 
 END SUBROUTINE CALC_AFT_REAC_AND_PROD
 
+!> \brief Dump CVODE warning summary
+SUBROUTINE DUMP_CVODE_WARNING_SUMMARY(NCHEM_ACTIVE_CELLS_AND_CC)
+   USE CHEMCONS, ONLY: CVODE_WARNING_CELLS,CVODE_ERR_CODE_MIN,CVODE_ERR_CODE_MAX,CVODE_WARN_MESSAGES
+
+   INTEGER, INTENT(IN) :: NCHEM_ACTIVE_CELLS_AND_CC
+   INTEGER :: NCHEM_ACTIVE_CELLS_AND_CC_GLOBAL,IERR, CVODE_ERR_CODE, TOT_WARN_CELLS
+   REAL(EB) :: PCNT
+   CHARACTER(LEN=16) :: CELL_STR, PCNT_STR
+   CHARACTER(LEN=40) :: CELL_INFO
+
+   IF (N_MPI_PROCESSES > 1) THEN
+      CALL MPI_ALLREDUCE(MPI_IN_PLACE, CVODE_WARNING_CELLS, SIZE(CVODE_WARNING_CELLS), &
+                         MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, IERR)
+
+      CALL MPI_ALLREDUCE(NCHEM_ACTIVE_CELLS_AND_CC, NCHEM_ACTIVE_CELLS_AND_CC_GLOBAL, &
+                         1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, IERR)
+   ELSE
+      NCHEM_ACTIVE_CELLS_AND_CC_GLOBAL = NCHEM_ACTIVE_CELLS_AND_CC
+   END IF
+   TOT_WARN_CELLS = SUM(CVODE_WARNING_CELLS)
+
+   IF (TOT_WARN_CELLS > 0 .AND. NCHEM_ACTIVE_CELLS_AND_CC_GLOBAL > 0 .AND. MY_RANK == 0) THEN
+
+      WRITE(LU_ERR, '(A)') '--------------- CVODE Warning Summary ----------------'
+      WRITE(LU_ERR, '(A)') ' Code     Cells (% Chem)             Message'
+      
+      DO CVODE_ERR_CODE = CVODE_ERR_CODE_MIN, CVODE_ERR_CODE_MAX
+         IF (CVODE_WARNING_CELLS(CVODE_ERR_CODE) > 0) THEN
+
+            PCNT = REAL(CVODE_WARNING_CELLS(CVODE_ERR_CODE)) / &
+                   REAL(NCHEM_ACTIVE_CELLS_AND_CC_GLOBAL) * 100.0_EB
+
+            ! Convert to strings and left-align
+            WRITE(CELL_STR, '(I0)') CVODE_WARNING_CELLS(CVODE_ERR_CODE)
+            WRITE(PCNT_STR, '(F6.2)') PCNT
+            CELL_STR = ADJUSTL(CELL_STR)
+            PCNT_STR = ADJUSTL(PCNT_STR)
+
+            ! Combine count + percentage into one left-aligned string
+            WRITE(CELL_INFO, '(A, " (", A, " %)")') TRIM(CELL_STR), TRIM(PCNT_STR)
+            CELL_INFO = ADJUSTL(CELL_INFO)
+
+            ! Print one line of summary
+            WRITE(LU_ERR, '(1X, I4, 5X, A, 5X, A)') CVODE_ERR_CODE, &
+                 TRIM(CELL_INFO), TRIM(CVODE_WARN_MESSAGES(CVODE_ERR_CODE))
+         END IF
+      END DO
+
+      WRITE(LU_ERR, '(A)') '------------------------------------------------------'
+
+   END IF
+
+END SUBROUTINE DUMP_CVODE_WARNING_SUMMARY
+
 #endif
 
 SUBROUTINE CHECK_AUTO_IGNITION(EXTINCT,TMP_IN,AIT,IIC,JJC,KKC,REAC_INDEX)