Adding Stoichiometric Coefficients (#514)

* Allowed for stoichometric coefficients in the chemical mechanism * added stoichometry to mechanism.prod and .reac files in tests to allow them to complete successfully * Added a test for the stoichometry based on the spec_model_1 test
AtChem · Jan 5, 2024 · 27ede44 · 27ede44
1 parent d3ce21b
commit 27ede44
Show file tree

Hide file tree

Showing 97 changed files with 15,893 additions and 11,521 deletions.
diff --git a/build/kpp_conversion.py b/build/kpp_conversion.py
@@ -203,7 +203,7 @@ def write_fac_file(input_file):
     print('Running write_fac_file() on: ' + str(input_file))
 
     contents1, contents2, contents3, contents4 = kpp_to_facsimile(input_file)
-    output_file = input_file.split('.')[0] + '.fac'
+    output_file = input_file.rsplit('.', 1)[0] + '.fac'
 
     with open(output_file, 'w') as file_open:
         file_open.write('\n* Generic Rate Coefficients ;\n')

diff --git a/build/mech_converter.py b/build/mech_converter.py
@@ -121,6 +121,34 @@ def tokenise_and_process(input_string, vars_dict):
     # Return the reconstructed string.
     return new_rhs
 
+def separate_stoichiometry(input_species):
+    """
+    This function takes in a string of a species from the mechanism and 
+    separates the species name from any preceeding stoichiometric coefficient.
+    This assumes that no species names will begin with a number.
+
+    Args:
+        input_species(str): a string containing a species name with a possible
+                            preceeding coefficient (e.g. H2O2, 2H2O2, 2 H2O2, or 0.5H2O2)
+
+    Returns:
+        split_spec (tuple): a tuple of a float and a string. The first item 
+                            (float) is the stoichiometric coefficient and the 
+                            second (string) is the species name.
+    """
+
+    #regex to match the potential coefficient and name sections of the input
+    in_pat = re.compile(r"^ *(\d*\.?\d*) *([a-zA-Z_].*) *$")
+    pat_match = in_pat.match(input_species)
+    if pat_match:
+        if pat_match[1]: #if there is a coefficient passed 
+            return (float(pat_match[1]), pat_match[2])
+        else: #if there is no coefficient then output an assumed coefficient of 1
+            return (1.0, pat_match[2])
+    else:
+        raise Exception(f"""Reaction species does not match the correct 
+                        format: '{input_species}'. Note that species names should 
+                        not begin with numerical characters.""")
 # ------------------------------------------------------------ #
 
 def convert_to_fortran(input_file, mech_dir, mcm_vers):
@@ -177,7 +205,7 @@ def convert_to_fortran(input_file, mech_dir, mcm_vers):
 
     # Check if the chemical mechanism file is in KPP format, in which case convert it
     # to FACSIMILE format (see documentation of `kpp_conversion.py` for more info)
-    if input_filename.split('.')[1] == 'kpp':
+    if input_filename.split('.')[-1] == 'kpp':
         input_fac = kpp_conversion.write_fac_file(os.path.join(input_directory, input_filename))
     else:
         input_fac = input_filename
@@ -399,39 +427,52 @@ def convert_to_fortran(input_file, mech_dir, mcm_vers):
 
             # Ignore empty reactantsList.
             if not reactantsList.strip() == '':
-                # Compare each reactant against known species.
+                # Compare each reactant against known species and note the 
+                # stoichometric coefficients for each reactant.
                 reactantNums = []
+                reactantStoichs = []
                 for x in reactants:
+                    #split the reactant into the species name and the stoichometric coefficient
+                    x_coeff, x_name = separate_stoichiometry(x)
+                    # Add the stoichometric coefficient to reactantStoichs
+                    reactantStoichs.append(x_coeff)
                     # If the reactant is a known species then add its number to reactantNums.
-                    if x in speciesList:
-                        reactantNums.append(speciesList.index(x)+1)
+                    if x_name in speciesList:
+                        reactantNums.append(speciesList.index(x_name)+1)
                     else:
                         # Reactant x is not a known species. Add reactant to speciesList,
                         # then add this number to reactantNums to record this reaction.
-                        speciesList.append(x)
+                        speciesList.append(x_name)
                         reactantNums.append(len(speciesList))
 
                 # Write the reactants to mech_reac_list.
-                mech_reac_list.extend([str(reactionNumber) + ' ' + str(z) \
-                                       + '\n' for z in reactantNums])
+                mech_reac_list.extend([f"{reactionNumber} {z} {y}\n" for \
+                                       y,z in zip(reactantStoichs, reactantNums)])
 
             # Ignore empty productsList.
             if not productsList.strip() == '':
-                # Compare each product against known species.
+                # Compare each product against known species and note the 
+                # stoichometric coefficients for each product.
                 productNums = []
+                productStoichs = []
                 for x in products:
+                    #split the product into the species name and the stoichometric coefficient
+                    x_coeff, x_name = separate_stoichiometry(x)
+                    # Add the stoichometric coefficient to productStoichs
+                    productStoichs.append(x_coeff)
+
                     # If the reactant is a known species then add its number to reactantNums.
-                    if x in speciesList:
-                        productNums.append(speciesList.index(x)+1)
+                    if x_name in speciesList:
+                        productNums.append(speciesList.index(x_name)+1)
                     else:
                         # Product x is not a known species. Add product to speciesList,
                         # then add this number to productNums to record this reaction.
-                        speciesList.append(x)
+                        speciesList.append(x_name)
                         productNums.append(len(speciesList))
 
                 # Write the products to mech_prod_list.
-                mech_prod_list.extend([str(reactionNumber) + ' ' + str(z) \
-                                       + '\n' for z in productNums])
+                mech_prod_list.extend([f"{reactionNumber} {z} {y}\n" for \
+                                       y,z in zip(productStoichs, productNums)])
 
     # -------------------------------------------------
 

diff --git a/src/inputFunctions.f90 b/src/inputFunctions.f90
@@ -64,6 +64,7 @@ subroutine readReactions()
 
     integer(kind=NPI) :: lhs_size, rhs_size
     integer(kind=NPI) :: k, l, count
+    real(kind=DP) :: m
     integer(kind=IntErr) :: ierr
 
     call inquire_or_abort( trim( configuration_dir ) // '/mechanism.reac', 'getReactantAndProductListSizes()')
@@ -82,16 +83,16 @@ subroutine readReactions()
     ! read data for lhs of equations
     ! lhs(1, i) contains the reaction number
     ! lhs(2, i) contains the species number of the reactant
-    ! lhs(3, i) contains 1 as a constant factor (stoichiometric coefficient)
+    ! coeff(i) contains the stoichiometric coefficient
     count = 0
     read (10,*, iostat=ierr)
-    read (10,*, iostat=ierr) k, l
+    read (10,*, iostat=ierr) k, l, m
     do while ( ierr == 0 )
       count = count + 1
       clhs(1, count) = k
       clhs(2, count) = l
-      clcoeff(count) = 1.0_DP
-      read (10,*, iostat=ierr) k, l
+      clcoeff(count) = m
+      read (10,*, iostat=ierr) k, l, m
     end do
     close (10, status='keep')
 
@@ -101,17 +102,17 @@ subroutine readReactions()
     ! read data for rhs of equations
     ! rhs(1, i) contains the reaction number
     ! rhs(2, i) contains the species number of the product
-    ! coeff(i) contains 1.0 as a constant factor (stoichiometric coefficient)
+    ! coeff(i) contains the stoichiometric coefficient
     count = 0
     ierr = 0
     read (11,*, iostat=ierr)
-    read (11,*, iostat=ierr) k, l
+    read (11,*, iostat=ierr) k, l, m
     do while ( ierr == 0 )
       count = count + 1
       crhs(1, count) = k
       crhs(2, count) = l
-      crcoeff(count) = 1.0_DP
-      read (11,*, iostat=ierr) k, l
+      crcoeff(count) = m
+      read (11,*, iostat=ierr) k, l, m
     end do
     close (11, status='keep')
 

diff --git a/tests/model_tests/INFO.md b/tests/model_tests/INFO.md
@@ -23,3 +23,5 @@ starting at 6:30 am on 9/11/2008.
 - `spec_model_1` is the base case
 
 - `spec_model_kpp` is the same as the base case but the mechanism is in KPP format.
+
+- `spec_model_stoich` is the same as the base case but the mechanism has been adjusted to include stoichometric coefficients (e.g. 'NO + NO = NO2 + NO2' becomes '2 NO = 2 NO2')
diff --git a/tests/model_tests/env_model_1/configuration/mechanism.prod.cmp b/tests/model_tests/env_model_1/configuration/mechanism.prod.cmp
@@ -1,8 +1,8 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
-2 4
-2 2
-3 1
-4 1
-5 3
-6 3
+1 2 1.0
+2 4 1.0
+2 2 1.0
+3 1 1.0
+4 1 1.0
+5 3 1.0
+6 3 1.0
diff --git a/tests/model_tests/env_model_1/configuration/mechanism.reac.cmp b/tests/model_tests/env_model_1/configuration/mechanism.reac.cmp
@@ -1,9 +1,9 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
-2 3
-3 2
-4 2
-5 2
-5 4
-6 4
-6 1
+1 1 1.0
+2 3 1.0
+3 2 1.0
+4 2 1.0
+5 2 1.0
+5 4 1.0
+6 4 1.0
+6 1 1.0
diff --git a/tests/model_tests/env_model_2/configuration/mechanism.prod.cmp b/tests/model_tests/env_model_2/configuration/mechanism.prod.cmp
@@ -1,8 +1,8 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
-2 4
-2 2
-3 1
-4 1
-5 3
-6 3
+1 2 1.0
+2 4 1.0
+2 2 1.0
+3 1 1.0
+4 1 1.0
+5 3 1.0
+6 3 1.0
diff --git a/tests/model_tests/env_model_2/configuration/mechanism.reac.cmp b/tests/model_tests/env_model_2/configuration/mechanism.reac.cmp
@@ -1,9 +1,9 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
-2 3
-3 2
-4 2
-5 2
-5 4
-6 4
-6 1
+1 1 1.0
+2 3 1.0
+3 2 1.0
+4 2 1.0
+5 2 1.0
+5 4 1.0
+6 4 1.0
+6 1 1.0
diff --git a/tests/model_tests/env_model_3/configuration/mechanism.prod.cmp b/tests/model_tests/env_model_3/configuration/mechanism.prod.cmp
@@ -1,8 +1,8 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
-2 4
-2 2
-3 1
-4 1
-5 3
-6 3
+1 2 1.0
+2 4 1.0
+2 2 1.0
+3 1 1.0
+4 1 1.0
+5 3 1.0
+6 3 1.0
diff --git a/tests/model_tests/env_model_3/configuration/mechanism.reac.cmp b/tests/model_tests/env_model_3/configuration/mechanism.reac.cmp
@@ -1,9 +1,9 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
-2 3
-3 2
-4 2
-5 2
-5 4
-6 4
-6 1
+1 1 1.0
+2 3 1.0
+3 2 1.0
+4 2 1.0
+5 2 1.0
+5 4 1.0
+6 4 1.0
+6 1 1.0
diff --git a/tests/model_tests/env_model_4/configuration/mechanism.prod.cmp b/tests/model_tests/env_model_4/configuration/mechanism.prod.cmp
@@ -1,8 +1,8 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
-2 4
-2 2
-3 1
-4 1
-5 3
-6 3
+1 2 1.0
+2 4 1.0
+2 2 1.0
+3 1 1.0
+4 1 1.0
+5 3 1.0
+6 3 1.0
diff --git a/tests/model_tests/env_model_4/configuration/mechanism.reac.cmp b/tests/model_tests/env_model_4/configuration/mechanism.reac.cmp
@@ -1,9 +1,9 @@
 4 6 7 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
-2 3
-3 2
-4 2
-5 2
-5 4
-6 4
-6 1
+1 1 1.0
+2 3 1.0
+3 2 1.0
+4 2 1.0
+5 2 1.0
+5 4 1.0
+6 4 1.0
+6 1 1.0
diff --git a/tests/model_tests/firstorder/configuration/mechanism.prod.cmp b/tests/model_tests/firstorder/configuration/mechanism.prod.cmp
@@ -1,2 +1,2 @@
 2 1 0 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
+1 2 1.0
diff --git a/tests/model_tests/firstorder/configuration/mechanism.reac.cmp b/tests/model_tests/firstorder/configuration/mechanism.reac.cmp
@@ -1,2 +1,2 @@
 2 1 0 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
+1 1 1.0
diff --git a/tests/model_tests/secondorder/configuration/mechanism.prod.cmp b/tests/model_tests/secondorder/configuration/mechanism.prod.cmp
@@ -1,3 +1,3 @@
 2 1 0 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 2
-1 2
+1 2 1.0
+1 2 1.0
diff --git a/tests/model_tests/secondorder/configuration/mechanism.reac.cmp b/tests/model_tests/secondorder/configuration/mechanism.reac.cmp
@@ -1,3 +1,3 @@
 2 1 0 numberOfSpecies numberOfReactions numberOfGenericComplex
-1 1
-1 1
+1 1 1.0
+1 1 1.0
Original file line number	Diff line number	Diff line change
Expand Up		@@ -23,3 +23,5 @@ starting at 6:30 am on 9/11/2008.
		- `spec_model_1` is the base case

		- `spec_model_kpp` is the same as the base case but the mechanism is in KPP format.

		- `spec_model_stoich` is the same as the base case but the mechanism has been adjusted to include stoichometric coefficients (e.g. 'NO + NO = NO2 + NO2' becomes '2 NO = 2 NO2')