Specifying Custom FORTRAN Functions in the Mechanism

AlfredMayhew commented 1 year ago

Hi, There's been a few times where I have felt it would be useful to be able to represent rates in the mechanism as a function with differing variables for different reactions. Previously, I have done the function bit beforehand and just included very long rate expressions in my mechanism files (with all of the calculation steps expanded into one line).

I have made some changes that I believe will allow a user to define custom functions. If this is something that you think would be valuable to add to AtChem2 then I'm happy to make any further changes and write a little in the documentation to explain this. Hopefully this wouldn't change the user experience for the majority of users (the customRateFuncs.f90 file can be left as its default of an empty module if it isn't needed), but adds some quite powerful flexibility for people wanting to use more complex rate expressions. I'll try and outline the changes below. Please let me know if there's any way I can improve the implementation.

Created a new file model/configuration/customRateFuncs.f90. This would be where users can specify their functions withing the custom_functions_mod module.
Added a section to build/mech_converter.py that looks for "function function_name (" and then adds function_name to the list of reserved names. I think that this bit could be improved. Do you know another way of parsing customRateFuncs.f90 to look for function names? I'm also not sure if I should be including subroutines here too? I don't really know FORTRAN very well, so I'm not 100% sure of the difference between functions and subroutines.
Added the customRateFuncs.f90 to the makefile skeleton so the build stage works properly.

I have tested this with a basic function that multiplies numbers together, so that my customRateFuncs.f90 file looks like this:

module custom_functions_mod
  implicit none

contains

  ! -----------------------------------------------------------------
  ! Multiply two numbers
  pure function multiplyNumbers( inNum1, inNum2 ) result ( outNum )
    real, intent(in) :: inNum2, inNum1
    real :: outNum

    outNum = inNum1 * inNum2

    return
  end function multiplyNumbers

end module custom_functions_mod

I can then include the following reaction in a mechanism file: % multiplyNumbers(2.0, 0.5) : = A ;. The model runs fine and A seems to be produced at the correct rate! I've also tested inputting TEMP into this function, which seems to work (provided I change the type of inNum2). Being able to input the environment variables into these functions adds even more flexibility.

Again, let me know if you can think of any ways to improve this and if you think it would be suitable for incorporation into the model. Thanks, Alfie

rs028 commented 1 year ago

@AlfredMayhew Thank you for this. It looks like a really good addition to the model. I need some time to understand the changes a think how to better integrate them, but it's definitely something useful!

AlfredMayhew commented 1 year ago

Of course, no worries. Let me know if there's anything I can do to help.

rs028 commented 9 months ago

@AlfredMayhew in order to review this PR could you please first update to the lastet version in the master branch? Thanks

AlfredMayhew commented 9 months ago

@AlfredMayhew in order to review this PR could you please first update to the lastet version in the master branch? Thanks

I think I've done this now? Is it just the master branch that needs updating? I can't seem to sync the custom_functions branch easily because of conflicts. If you need more then just let me know.

rs028 commented 9 months ago

I can't see if you done it. The usual procedure for me is:

git checkout master
git pull upstream master
git checkout custom_functions
git rebase master
git push --force

When you rebase it will ask you to resolve the conflicts (if any). This will bring the branch even with the current master branch and apply your changes on top.

AlfredMayhew commented 9 months ago

Thanks, I think I've done that now. I was a bit confused about resolving conflicts but I think I've done it now. I've run a test and it still seems to be working as it was but now GitHub is showing no conflicts. Let me know if I've missed anything!

rs028 commented 9 months ago

I am guessing that the tests fail because the makefile is looking for customRateFuncs.f90 and it doesn't exist. The easiest is probably to copy the skeleton file into the configuration directory of each test, then hopefully it should pass the test suite.

rs028 commented 9 months ago

Two more things.

The example you provide here could be converted into a proper model test
can you add some explanatory text to the manual?

codecov[bot] commented 9 months ago

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Comparison is base (63cf98e) 52.05% compared to head (d4042c3) 52.05%.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #507 +/- ## ======================================= 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

A better way to add a test (instead of using the function to multiply two numbers as I suggested above) is to use model_tests/spec_model_1 as base case. You change the mechanism so that one of the complex rate coefficients (e.g. KMT15) is implemented in the customRateFuncs.f90 file. The results of the test should be identical to the results of spec_model_1.

AlfredMayhew commented 9 months ago

I have added a test replacing KMT15. As with my latest post to #514, this test currently fails because the produced output is slightly different from the output of spec_model_1. Here, the biggest difference is for C2H4 of 1000 molecules cm-3, which is of course very small. I'm also reassured that this difference doesn't compound. The timestep after a difference between the two models shows no difference between the two models.

If you're happy then I can update the test to use my model output. I'll also add some text to the manual at some point too.

rs028 commented 9 months ago

This looks good to me. I suggest we first merge PR #514 (you should be able to, let me know if you aren't), then rebase this branch. If possible you can add some text to the manual - or just send it to me, I am revising the manual anyways so I can incorporate it into next version.

AlfredMayhew commented 9 months ago

I have updated this branch to work with the new stoichiometry branch. Let me know if there are any other changes you think would be beneficial. In terms of documentation, I will send you some modifications to the manual separately that I think would be helpful.

rs028 commented 9 months ago

I have added a couple more comments. This needs to be rebased again after PR #515 has been merged.

AlfredMayhew commented 9 months ago

Ok, that should be those changes made, including the rebase. I can see a grayed out merge pull request button, so I should be able to merge this once you approve it. Let me know if there's anything else I need to change!

AtChem / AtChem2

Specifying Custom FORTRAN Functions in the Mechanism #507

Codecov Report