Convert to using function parser rather than hardcoding and recompiling

[spco]

• [ ] Separate out functions from mechanism-rates.f90 and tools/mechanism-rate-coefficients.f90 that are not generated from the downloaded MCM file
• [ ] In converter, change the Generic Rate Coefficients and Complex Reactions sections to replace the variable names with elements of q. e.g. KRO2NO becomes q(1) and so on. Replace references to these from the Reaction Definitions section.
• [ ] Add function parser into the repo, and add to Makefile
• [ ] tools/facsimileToFortranMechanismConverter.py needs to be superceded by a version that converts to the correct syntax for the function parser to understand - this should create 2 files - one for the generic/complex sections, which will contain one line per element of q, and another to hold the rows for p. Keep the useful information in comments.
• [ ] add function parser initialisation stage
• [ ] change internals of solver to use function parser rather than the hardcoded functions
• [ ] remove hardcoded functions, and corresponding Fortran files

[spco]

This is to be done only after the pre-function parser milestone is passed.

[spco]

Fparser can't natively handle an expression like KR17 = K170/K17I because we need to declare KR17 as a variable at compile time. This affects our ability to compute the Generic Rate Coefficients and Complex Reactions section of gen/mechanism-rate-coefficients.f90 with the parser - we'd need to know the names of these at compile time, which isn't guaranteed.

I think the way around this is to do a little more in the Python stage before running atchem. Suppose we have the following gen/mechanism-rate-coefficients.f90 file:

KCH3O2 = 1.03e-13_DP*EXP(365.0_DP/TEMP)
K298CH3O2 = 3.5e-13_DP
K14ISOM1 = 3.00e7_DP*EXP(-5300_DP/TEMP)
!*;
!* Complex reactions ;
!*;
KD0 = 1.10e-05_DP*M*EXP(-10100_DP/TEMP)
KDI = 1.90e17_DP*EXP(-14100_DP/TEMP)
KRD = KD0/KDI
!****************************************************** ;
!*;
!* Reaction definitions. ;
!*;
p(1) = 5.6D-34*N2*(TEMP/300)**(-2.6)*O2   !% 5.6D-34*N2*(TEMP/300)@-2.6*O2 : O = O3 ;
p(2) = KMT01   !% KMT01 : O + NO = NO2 ;
p(3) = KMT02   !% KMT02 : O + NO2 = NO3 ;

In the python stage, we need to do more processing so that instead we get the following:

q(1) = 1.03e-13_DP*EXP(365.0_DP/TEMP)
q(2) = 3.5e-13_DP
q(3) = 3.00e7_DP*EXP(-5300_DP/TEMP)
!*;
!* Complex reactions ;
!*;
q(4)KD0 = 1.10e-05_DP*M*EXP(-10100_DP/TEMP)
q(5)KDI = 1.90e17_DP*EXP(-14100_DP/TEMP)
q(6)KRD = q(4)/q(5)
q(7)KMT01 = 1.0e-05_DP
q(8)KMT02 = q(4)/q(7)
!****************************************************** ;
!*;
!* Reaction definitions. ;
!*;
p(1) = 5.6D-34*N2*(TEMP/300)**(-2.6)*O2   !% 5.6D-34*N2*(TEMP/300)@-2.6*O2 : O = O3 ;
p(2) = q(7)   !% KMT01 : O + NO = NO2 ;
p(3) = q(8)   !% KMT02 : O + NO2 = NO3 ;

I think we just need to declare real(kind=DP) :: N2, O2, M, RH, H2O, DEC, BLH, DILUTE, JFAC, ROOFOPEN and then read and allocate the correct number of elements for q, (as we do for pcurrently).

HOWEVER... It also turns out that you can't use the function parser to read from array elements - only specified variables. I think the way to get around that is to hack the parser so that it takes as input the array q that we want it to be able to read, and uses that when it hits p(1) in the same way that intrinsics are called. This may take a bit of time, but it should work...

I'll be working on a mock-up small prototype to test the feasibility of that all.

[rs028]

As far as the first point is concerned, I think it would be good to do it anyways. It seems a cleaner way to do things (and maybe more robust).

I am a bit uncomfortable with the second point. Modifying the parser may cause some difficulty in keeping everything up to date (not to mention installation for unexperienced users). But I guess it depends on how complicated that would be.

Is there only one parser that can be used?

[spco]

There are 3 different versions of function parser: https://sourceforge.net/projects/fparser/, https://github.com/jacobwilliams/fortran_function_parser, and https://github.com/jacopo-chevallard/FortranParser. I'm looking into the latter mostly, as it is in a more modern style than the first, and slightly more recently developed than the second. They are all very similar, as the latter 2 are based on the first by Roland Schmehl, and this is itself based on the C++ parser at http://warp.povusers.org/FunctionParser/ .

As none of them are being actively developed that I can see, my intention would be to just bundle a modified distribution in with ours. This is possible under the BSD licensing, so long as we keep the necessary headers. So we shouldn't fall behind their development (as there isn't any), and we can modify it to our specific usecase. It just becomes a couple more files in our src/ directory, so it's nothing more for the user to do or install.

[rs028]

Okay then. But I think we need to decide on the priorities, depending on your time, with a view to close PR #303 and have a stable release.

[spco]

Once #347 is merged, we are halfway there already. After that, we need to:

• [ ] Rather than generate gen/mechanism-rate-coefficients.f90, create 2 files (perhaps mechanism.gc for the generic/complex reactions and mechanism.rates for the reaction rates, which can live in model/configuration along with the rest. These should not contain the q(i) = or p(i) = parts of each line, but be just the string of the RHS.
• [ ] Add FunctionParser to the repo (modified to handle q and J correctly - I have a prototype at https://github.com/spco/AtChem2/tree/function_parser which implements the q part, and J will be handled similarly when I get back to it)
• [ ] Add FunctionParser files to the build process
• [ ] Add functions to read each of mechanism.rates and mechanism.gc in as strings and pass them to the function parser initialisation stage.
• [ ] Rather than include gen/mechanism-rates-coefficients.f90, instead call the function parser to execute each of the functions in turn (essentially this does the same as the contents of gen/mechanism-rates-coefficients.f90).
• [ ] Sit back and enjoy the lack of compilation between runs and sped-up tests OR iron out all the bits I've missed, and then cry about how much the performance is affected 😆

[rs028]

Fantastic!

One minor point. mechanism.gc is still about rate coefficients, just different ones than in mechanism.rates. Maybe the naming structure should be changed

old name	new name	alternative new name
mechanism.species	mechanism.species	species.mecha
mechanism.reac	mechanism.reactants	reactants.mecha
mechanism.prod	mechanism.products	products.mecha
mechanism.gc	mechanism.generic-ratecoeff	ratecoeff-generic.mecha
mechanism.rates	mechanism.ratecoeff	ratecoeff.mecha

or something like that (these are just suggestions).

[spco]

This has been superceded by the use of a shared library in #370.

[spco]

As #370 has been merged, I'm going to close this. I'll leave the branch on my fork at https://github.com/spco/AtChem2/tree/function_parser around for a while as it may be useful for future reference.