AtChem / AtChem2

Atmospheric chemistry box-model for the MCM
MIT License
59 stars 23 forks source link

Adding Stoichiometric Coefficients #514

Closed AlfredMayhew closed 9 months ago

AlfredMayhew commented 9 months ago

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:

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[bot] commented 9 months ago

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 :warning:
Additional details and impacted files ```diff @@ 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 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

rs028 commented 9 months ago

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 commented 9 months ago

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.