Adding Stoichiometric Coefficients

[AlfredMayhew]

I have looked into adding stoichiometric coefficients following discussions in #513. It turns out there was already a coefficient assigned in src/inputFunctions.f90 for each species in each reaction, it was just being set to 1 every time. This coefficient of 1 was then already included in calculations in the resid subroutine in src/solverFunctions.f90. This made it fairly easy to modify the code to have variable coefficients!

I have made the following adjustments:

• Check for stoichiometric coefficients when reading the species (reactants and products) in each reaction in build/mech_converter.py. If no coefficient is given then assign a coefficient of 1.
• Carry these coefficients through as a new column in the files mechanism.reac and mechanism.prod.
• Read the coefficients in src/inputFunctions.f90 instead of assigning a value of 1 for all species.

I've run some models to test this and it seems to be behaving well. I've compared against the previous version with species duplicated and the results are identical. I've also tested some edge-cases like setting the coefficient to 0 and non-integer reactant coefficients, which also work as intended.

Let me know if there's any more changes I need to make, otherwise this resolves #513.

[codecov]

Codecov Report

Attention: 2 lines in your changes are missing coverage. Please review.

Comparison is base (d3ce21b) 52.05% compared to head (0cab63b) 52.05%.

Files	Patch %	Lines
src/inputFunctions.f90	66.66%	2 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #514   +/-   ##
=======================================
  Coverage   52.05%   52.05%           
=======================================
  Files          17       17           
  Lines        2096     2096           
  Branches      166      166           
=======================================
  Hits         1091     1091           
  Misses        933      933           
  Partials       72       72           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[rs028]

I haven't tested myself but it looks okay to me.

A suggestion for a test. Using model_tests/spec_model_1 as base case, you can change couple of reactions to use stoichiometric coefficients. for example:

HO2 + HO2 = H2O2
HCHO = CO + HO2 + HO2

The results of the test should be identical to the results of spec_model_1.

[AlfredMayhew]

I added a test of the stoichiometry here, but it is failing. To make the test, I copied the spec_model_1 directory and manually edited the mechanism file along with mechanism.prod.cmp and mechanism.react.cmp to replace any duplicated species with a single species with a stoichiometric coefficient.

Looking at the files from the test, there only seems to be some small differences in the reactionRates files. When comparing the species concentrations, there are no large differences between the two models (the largest difference is 10 molecules cm-3 of HO2, and the largest % difference is 2E-4%), but maybe this falls outside of the tolerance of numdiff?

I don't think it's unreasonable to accept small differences between using the split mechanism and the stoichiometric coefficient mechanism since the problem being passed to the solver is slightly different, but I wanted to check that it would be OK to adjust this test to use the output from my stoichiometric model rather than comparing to the exact output from spec_model_1?

Note that I have also made changes to the .kpp -> .fac conversion ( c34cdaa ) to allow my kpp test to complete as my filepath contains points ('.'), so was causing errors in trying to create the new kpp file.

[rs028]

I am satisfied that the differences between spec_model_1 and spec_model_stoich are negligible. This can be merged, thank you!