FDS Source: check atom balance code

Author: cxp484

Source/chem.f90

@@ -845,6 +845,8 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
 REAL(EB) :: H_G
 TYPE(USERDATA), TARGET :: USER_DATA
 LOGICAL :: ONLY_FIRST_STEP=.TRUE. ! Needed in CV_ONE_STEP
+LOGICAL :: RETURN_ORIG = .FALSE.
+REAL(EB) :: ZZ_CUR(N_TRACKED_SPECIES),EQUIV_CUR
 
 
 !======= INTERNALS ============
@@ -1012,8 +1014,9 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
    IF (IERR_C .NE. CV_SUCCESS) THEN
       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
+      ENDIF
+      RETURN_ORIG = .TRUE.
+      !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."
@@ -1024,17 +1027,24 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
 
          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)
+         CALL MOLAR_CONC_TO_MASS_FRAC(CVEC_C(1:N_TRACKED_SPECIES), ZZ_CUR(1:N_TRACKED_SPECIES))
+         CALL CALC_EQUIV_RATIO(ZZ_CUR(1:N_TRACKED_SPECIES), EQUIV_CUR)
          DO NS = 1, N_TRACKED_SPECIES
-            WRITE(LU_ERR,*)" ID, Y=",SPECIES_MIXTURE(NS)%ID, ZZ(NS)
+            WRITE(LU_ERR,*)" ID, Y, Y_CUR=",SPECIES_MIXTURE(NS)%ID, ZZ(NS),ZZ_CUR(NS)
          ENDDO 
-         WRITE(LU_ERR,*)" EQUIVALENCE RATIO, TMP=", EQUIV,TMP_IN
+         WRITE(LU_ERR,*)" EQUIVALENCE RATIO, TMP, EQUIV_CUR=", EQUIV,TMP_IN,EQUIV_CUR
          CALL CVODESTATS(CVODE_MEM) ! DIAGNOSTICS OUTPUT
-      ENDIF   
+      !ENDIF   
    ENDIF   
 ENDIF
 
-CC = CVEC_C(1:N_TRACKED_SPECIES) !DISCARD THE TEMPERATURE.
-TMP_OUT = CVEC_C(N_TRACKED_SPECIES+1)
+IF (.NOT. RETURN_ORIG) THEN
+   CC = CVEC_C(1:N_TRACKED_SPECIES) !DISCARD THE TEMPERATURE.
+   TMP_OUT = CVEC_C(N_TRACKED_SPECIES+1)
+ELSE
+   ! CC no change. To avoid error in elemental mass balance. 
+   TMP_OUT = TMP_IN
+ENDIF   
 IERR_C = FCVODEGETLASTSTEP(CVODE_MEM, CHEM_TIME_C(1))
 CHEM_TIME=CHEM_TIME_C(1)
 

Source/fire.f90

@@ -971,6 +971,7 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
             ZETA = 1._EB + Q_ALLOWED*DT/(RHO_IN * (SUM(SPECIES_MIXTURE%H_F*(ZZ_MIXED-ZZ_0))) + TWO_EPSILON_EB)
             ZETA = MAX(0._EB,MIN(1.0_EB,ZETA))
             Q_CAPPED_CELLS = Q_CAPPED_CELLS + 1
+            WRITE(LU_ERR,*)"Capped cell."
          ENDIF 
       ENDIF   
    ENDIF   
@@ -1013,7 +1014,7 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
 
 
 ! Compute heat release rate
-
+CALL PERFORM_ELEMENTAL_BALANCE(ZZ_GET,ZZ_0)
 Q_OUT = -RHO_IN*SUM(SPECIES_MIXTURE%H_F*(ZZ_GET-ZZ_0))/DT ! FDS Tech Guide (5.47)
 
 ! Extinction model
@@ -1064,6 +1065,113 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
 
 END SUBROUTINE COMBUSTION_MODEL
 
+SUBROUTINE PERFORM_ELEMENTAL_BALANCE(ZZ_GET, ZZ_0)
+   USE PROPERTY_DATA
+   REAL(EB), INTENT(IN) :: ZZ_GET(N_TRACKED_SPECIES),ZZ_0(N_TRACKED_SPECIES)
+   
+   INTEGER :: NE, NS
+   INTEGER :: N_ELEMENTS
+   REAL(EB), ALLOCATABLE :: ZZ_ELEM_GET(:), ZZ_ELEM_0(:)
+   REAL(EB) :: DIFF
+   LOGICAL :: IMBALANCE_FOUND
+
+   !-------------------------------
+   ! Setup
+   !-------------------------------
+   N_ELEMENTS = 118
+
+   ALLOCATE(ZZ_ELEM_GET(N_ELEMENTS))
+   ALLOCATE(ZZ_ELEM_0(N_ELEMENTS))
+
+   ZZ_ELEM_GET = 0.0_EB
+   ZZ_ELEM_0   = 0.0_EB
+   IMBALANCE_FOUND = .FALSE.
+
+   !-------------------------------
+   ! Compute elemental mass fractions
+   !-------------------------------
+   DO NS = 1, N_TRACKED_SPECIES
+      DO NE = 1, N_ELEMENTS
+
+         ZZ_ELEM_GET(NE) = ZZ_ELEM_GET(NE) + &
+              ZZ_GET(NS) * SPECIES_MIXTURE(NS)%ATOMS(NE) * ELEMENT(NE)%MASS / &
+              SPECIES_MIXTURE(NS)%MW
+
+         ZZ_ELEM_0(NE) = ZZ_ELEM_0(NE) + &
+              ZZ_0(NS) * SPECIES_MIXTURE(NS)%ATOMS(NE) * ELEMENT(NE)%MASS / &
+              SPECIES_MIXTURE(NS)%MW
+
+      END DO
+   END DO
+
+   !-------------------------------
+   ! Compare elemental balances
+   !-------------------------------
+   DO NE = 1, N_ELEMENTS
+
+      IF (ZZ_ELEM_GET(NE) > 0.0_EB .OR. ZZ_ELEM_0(NE) > 0.0_EB) THEN
+
+         DIFF = ABS(ZZ_ELEM_GET(NE) - ZZ_ELEM_0(NE))
+
+         IF (DIFF > 1.0E-4_EB) THEN
+            IMBALANCE_FOUND = .TRUE.
+            !WRITE(*,'(A,I3,3E14.6)') 'Imbalance in element ', NE, &
+            !     ZZ_ELEM_GET(NE), ZZ_ELEM_0(NE), DIFF
+         END IF
+
+      END IF
+
+   END DO
+
+   !-------------------------------
+   ! Optional summary message
+   !-------------------------------
+   IF (IMBALANCE_FOUND) THEN
+
+      WRITE(LU_ERR,*) ' '
+      WRITE(LU_ERR,*) '*** Elemental conservation violated ***'
+      WRITE(LU_ERR,*) ' '
+
+      !---------------------------------------
+      ! Elemental mass fractions
+      !---------------------------------------
+      WRITE(LU_ERR,'(A)') 'Elemental mass fractions (GET vs 0):'
+      WRITE(LU_ERR,'(A)') '-------------------------------------'
+
+      DO NE = 1, N_ELEMENTS
+         IF (ZZ_ELEM_GET(NE) > 0.0_EB .OR. ZZ_ELEM_0(NE) > 0.0_EB) THEN
+            WRITE(LU_ERR,'(A4,3E16.8)') TRIM(ELEMENT(NE)%ABBREVIATION), &
+                 ZZ_ELEM_GET(NE), ZZ_ELEM_0(NE), &
+                 ZZ_ELEM_GET(NE) - ZZ_ELEM_0(NE)
+         END IF
+      END DO
+
+      !---------------------------------------
+      ! Species comparison
+      !---------------------------------------
+      WRITE(LU_ERR,*) ' '
+      WRITE(LU_ERR,'(A)') 'Species with differences (ZZ_GET vs ZZ_0):'
+      WRITE(LU_ERR,'(A)') '-------------------------------------------'
+
+      DO NS = 1, N_TRACKED_SPECIES
+         IF (ABS(ZZ_GET(NS) - ZZ_0(NS)) > 1.0E-8_EB) THEN
+            WRITE(LU_ERR,'(A20,3E16.8)') TRIM(SPECIES_MIXTURE(NS)%ID), &
+                 ZZ_GET(NS), ZZ_0(NS), ZZ_GET(NS) - ZZ_0(NS)
+         END IF
+      END DO
+
+      !---------------------------------------
+      ! Totals
+      !---------------------------------------
+      WRITE(LU_ERR,*) ' '
+      WRITE(LU_ERR,'(A,2E16.8)') 'Sum ZZ_GET, ZZ_0 = ', SUM(ZZ_GET), SUM(ZZ_0)
+
+   END IF
+
+   DEALLOCATE(ZZ_ELEM_GET, ZZ_ELEM_0)
+
+END SUBROUTINE PERFORM_ELEMENTAL_BALANCE
+
 !> \brief call cvode_interface after converting mass fraction to molar concentration.
 !> \param ZZ species mass fraction array
 !> \param TMP_IN is the temperature