Correcting Rate Output for Non-integer Stoichiometry (#521)

* Corrected rate output files to correctly report the effect of non-integer stoichiometry (as opposed to just counting the number of occurences of species in a reaction).

* Fixed formatting to pass style test

* Added ccoeff to findReactionsWithProductOrReactant to fix failing unit test

* Replaced rate output for stoichiometry test to output correct rates

* Comment explaining why stoich is present in the reaction_frequency_pair type

src/atchem2.f90

@@ -225,7 +225,7 @@ end subroutine FCVFUN
   allocate (detailedRatesSpecies(size( detailedRatesSpeciesName )))
   allocate (reacDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( clhs, 2 )))
   allocate (prodDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( crhs, 2 )))
-  invalid_reaction_frequency_pair = reaction_frequency_pair(-1_NPI, 0_NPI)
+  invalid_reaction_frequency_pair = reaction_frequency_pair(-1_NPI, 0_NPI, 0.0_DP)
   reacDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair
   prodDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair
   allocate (reacDetailedRatesSpeciesLengths(size( detailedRatesSpeciesName )))
@@ -243,8 +243,10 @@ end subroutine FCVFUN
   !
   ! Fill the remaining elements of each row of reac/prodDetailedRatesSpecies with the
   ! numbers of the reactions in which that species appears as a reactant/product respectively
-  call findReactionsWithProductOrReactant( detailedRatesSpecies, clhs, reacDetailedRatesSpecies, reacDetailedRatesSpeciesLengths )
-  call findReactionsWithProductOrReactant( detailedRatesSpecies, crhs, prodDetailedRatesSpecies, prodDetailedRatesSpeciesLengths )
+  call findReactionsWithProductOrReactant( detailedRatesSpecies, clhs, clcoeff, reacDetailedRatesSpecies, &
+                                           reacDetailedRatesSpeciesLengths )
+  call findReactionsWithProductOrReactant( detailedRatesSpecies, crhs, crcoeff, prodDetailedRatesSpecies, &
+                                           prodDetailedRatesSpeciesLengths )
   write (*, '(A, I0)') ' Species requiring detailed rate output (number of species found): ', size( detailedRatesSpeciesName )
   write (*,*)
 

src/configFunctions.f90

@@ -83,14 +83,16 @@ end function getIndexWithinList
   ! At the end of this function, each row of r is associated to one of
   ! the species in rSpecies. Each element of the row corresponds to a
   ! reaction in which this species appears as a product or reactant.
-  ! Each element consists of %reaction holding the reaction number, and
+  ! Each element consists of %reaction holding the reaction number,
   ! %frequency holding the number of occurences of that species in that
-  ! reaction. chs contains the product or reactant information, and
+  ! reaction, and %stoich holding the sum of the stoichiometric coefficients 
+  ! for that species in that reaction. chs contains the product or reactant information, and
   ! arrayLen is used to keep track of how many are present in each row.
-  subroutine findReactionsWithProductOrReactant( rSpecies, chs, r, arrayLen )
+  subroutine findReactionsWithProductOrReactant( rSpecies, chs, ccoeff, r, arrayLen )
     use types_mod
 
-    integer(kind=NPI), intent(in) :: rSpecies(:), chs(:,:)
+    integer(kind=NPI), intent(in) :: rSpecies(:), chs(:,:) 
+    real(kind=DP), intent(in) :: ccoeff(:)
     type(reaction_frequency_pair), intent(inout) :: r(:,:)
     integer(kind=NPI), intent(out) :: arrayLen(:)
     integer(kind=NPI) :: rCounter, i, j
@@ -125,6 +127,7 @@ subroutine findReactionsWithProductOrReactant( rSpecies, chs, r, arrayLen )
           end if
           r(i, rCounter)%reaction = chs(1, j)
           r(i, rCounter)%frequency = r(i, rCounter)%frequency + 1_NPI
+          r(i, rCounter)%stoich = r(i, rCounter)%stoich + ccoeff(j)
         end if
       end do
       arrayLen(i) = rCounter

src/dataStructures.f90

@@ -43,6 +43,9 @@ module types_mod
   type reaction_frequency_pair
     integer(kind=NPI) :: reaction
     integer(kind=NPI) :: frequency
+    ! stoich was added in Jan 2025 by Alfred Mayhew to allow proper output of rates for reactions
+    ! involving non-integer stoichiometry (instead of using frequency)
+    real(kind=DP) :: stoich
   end type reaction_frequency_pair
 
   interface operator (==)

src/outputFunctions.f90

@@ -183,10 +183,12 @@ end subroutine outputPhotolysisRates
   ! reaction, a string representing that reaction.
   pure function getReaction( speciesNames, reactionNumber ) result ( reaction )
     use types_mod
-    use reaction_structure_mod, only : clhs, crhs
+    use reaction_structure_mod, only : clhs, crhs, clcoeff, crcoeff
     use storage_mod, only : maxSpecLength, maxReactionStringLength
 
     character(len=maxSpecLength) :: reactants(10), products(10)
+    character(len=4) :: reactCoeffs(10), prodCoeffs(10)
+    character(len=4) :: tmpCoeffStr
     character(len=maxSpecLength), intent(in) :: speciesNames(*)
     integer(kind=NPI) :: i, numReactants, numProducts
     integer(kind=NPI), intent(in) :: reactionNumber
@@ -200,12 +202,21 @@ pure function getReaction( speciesNames, reactionNumber ) result ( reaction )
       if ( clhs(1, i) == reactionNumber ) then
         numReactants = numReactants + 1
         reactants(numReactants) = speciesNames(clhs(2, i))
+
+        !Store the reaction coefficient (stoichiometry) string for this reactant. If
+        !the coefficient is 1, then just use an empty string
+        if ( clcoeff(i) == 1.0 ) then
+          reactCoeffs(numReactants) = ''
+        else
+          write(tmpCoeffStr, '(F4.2)') clcoeff(i)
+          reactCoeffs(numReactants) = tmpCoeffStr
+        end if
       end if
     end do
 
     reactantStr = ' '
     do i = 1, numReactants
-      reactantStr = trim( adjustl( trim( reactantStr ) // trim( reactants(i) ) ) )
+      reactantStr = trim( adjustl( trim( reactantStr ) // trim( reactCoeffs(i) ) // trim( reactants(i) ) ) )
       if ( i < numReactants ) then
         reactantStr = trim( reactantStr ) // '+'
       end if
@@ -222,12 +233,21 @@ pure function getReaction( speciesNames, reactionNumber ) result ( reaction )
       if ( crhs(1, i) == reactionNumber ) then
         numProducts = numProducts + 1
         products(numProducts) = speciesNames(crhs(2, i))
+
+        !Store the reaction coefficient (stoichiometry) string for this product. If
+        !the coefficient is 1, then just use an empty string
+        if ( crcoeff(i) == 1.0 ) then
+          prodCoeffs(numProducts) = ''
+        else
+          write(tmpCoeffStr, '(F4.2)') crcoeff(i)
+          prodCoeffs(numProducts) = tmpCoeffStr
+        end if
       end if
     end do
 
     productStr = ' '
     do i = 1, numProducts
-      productStr = trim( adjustl( trim( productStr ) // trim( products(i) ) ) )
+      productStr = trim( adjustl( trim( productStr ) // trim( prodCoeffs(i) ) //  trim( products(i) ) ) )
       if ( i < numProducts ) then
         productStr = trim( productStr ) // '+'
       end if
@@ -295,7 +315,7 @@ subroutine outputRates( rSpecies, r, arrayLen, t, p, flag )
           write (output_file_number, '(ES15.6E3, I14, A52, I15, ES15.6E3, A, A)') t, rSpecies(i), &
                                                                                   trim( speciesNames(rSpecies(i)) ), &
                                                                                   r(i, j)%reaction, &
-                                                                                  r(i, j)%frequency * p(r(i, j)%reaction), &
+                                                                                  r(i, j)%stoich * p(r(i, j)%reaction), &
                                                                                   '  ', trim( reaction )
         end if
       end do