Dilute fix

[bsn502]

New files: ./dilute.py
Modified files: ./build/build_atchem2.sh

The new code uses the "DILUTE" number in environmentVariables.config to amend the mechanism provided. Tests have been successful and the code is very useful, especially when using a larger mechanism file.

Current issues: the model needs to be re-compiled if changes are made to the dilution rate; the changes currently only work if model folder is called "model".

An example using the ethane subset from the mcm and an initial concentration of ethane shows how its concentration changes under three different DILUTE rates.

[rs028]

This partially addresses the discussion on DILUTE in #355. It may eliminate the need to set this variable to NOTUSED and we can force the user to set its value to either zero or a dilution rate. If all species are diluted to zero this effectively means there is no dilution. However, this approach doubles the size of the chemical mechanism so the question is: can it slow down the execution? @spco any thoughts?

[spco]

I think I understand how this works. It makes sense - whatever we do, at some point if DILUTE is set to non-zero then we're going to have to iterate over all species and apply the factor. So this mechanism seems reasonable.

A couple of comments:

(1) This seems like it will only work for a fixed value of DILUTE, but not for a constrained value. Is that an issue? I don't see a simple way to extend this mechanism to a constrained value.

(2) Currently this only checks for NOTUSED. Could we also check for 0? If DILUTE is 0 (also handle 0.0 and 0.0e0 etc by converting the string to real), then we don't want to bother adding all these extra rows, in just the same way as NOTUSED. That will mean there's no overhead in the common can where no dilution is denoted by 0. Thinking about it, we could go further and outlaw NOTUSED as a valid value - insist that the input is 0 instead. This would require modification to the Fortran parsing of the environment variables, but simplify the possibilities for users.

[spco]

Actually, the more I think about this, the more I think this could be changed to avoid the double call to mech_converter.py, and make it cleaner overall. If the contents of dilute.py can be added directly to mech_converter.py, that makes everything a lot neater.

• Loop over speciesList to add the dilution reactions directly on about line 405 of mech_converter.py.
• You would need to add another input argument to mech_converter.py to point to the environmentVariables.config file to get the value of DILUTE. Put it second (behind the argument for the .fac file) as the others are optional.
• Edit build_atchem2.sh to fit that definition.

This removes the need to call mech_converter.py twice (which makes the flow easier to understand) and also removes the need to make a copy of the .fac file (and the resulting cleanup). Does that make sense?

[rs028]

(1) This seems like it will only work for a fixed value of DILUTE, but not for a constrained value. Is that an issue? I don't see a simple way to extend this mechanism to a constrained value.

At the moment DILUTE can be a number, NOTUSED or constrained. To be honest, I can't see a situation where DILUTE needs to be constrained, so I am not too worried about this now. That being said it doesn't hurt to have this possibility. The question is whether this approach (even if re-implemented as you suggest in #408 (comment)) can work or not when DILUTE is constrained.

(2) Currently this only checks for NOTUSED. Could we also check for 0? If DILUTE is 0 (also handle 0.0 and 0.0e0 etc by converting the string to real), then we don't want to bother adding all these extra rows, in just the same way as NOTUSED. That will mean there's no overhead in the common can where no dilution is denoted by 0. Thinking about it, we could go further and outlaw NOTUSED as a valid value - insist that the input is 0 instead. This would require modification to the Fortran parsing of the environment variables, but simplify the possibilities for users.

I agree. This is sensible. How hard is it to remove the NOTUSED option? Then DILUTE can only be one of these values:

1. 0, in which case dilute.py is not called and nothing is added to the mechanism
2. a real in which case dilute.py is called and the dilution equations are added to the mechanism
3. constrained in which case dilute.py is called and the dilution equations are added to the mechanism but the dilution rate is read from a file in the environmental constraints directory (would it work? if not it's not essential, see above comment)

Does it make sense?

[spco]

How hard is it to remove the NOTUSED option?

Good point. I don't think it would be too hard to remove the NOTUSED option - it's mostly just the handling at

AtChem2/src/inputFunctions.f90

Line 869 in a74ec61

case ('TEMP', 'RH', 'H2O', 'PRESS', 'BLHEIGHT', 'DILUTE', 'DEC')

that would need changing.

But, that actually leads us on to thinking about what DILUTE now is. Can it be references like a species in other chemical reactions? I'm guessing not. In which case, it's no longer handled in the same way as the other environment variables (EVs), and should in some ways be removed from a lot of the EV functions.

Currently, the whole process ends up with a value for dilute which can be used in the mechanism, but if it's just going to be used as a universal factor on all species as implemented in this PR, then we're no longer interested in re-evaluating that in the same way as the other EVs. So we might end up ripping a fair amount out (but just the lines that handle dilute).

To extend this to handle a constrained dilute, I think we can do something similar to what's done for the photolysis rates J<1> etc. Rather than printing

0.1 ; O2 = LOSS
0.1 ; H20 = LOSS
0.1 ; CO2 = LOSS

etc, where DILUTE = 0.1 in environmentVariables.config, instead we could print

dilute ; O2 = LOSS
dilute ; H20 = LOSS
dilute ; CO2 = LOSS

and then evaluate dilute as we currently do at each timestep. Does that make sense?

[spco]

Oh, and in fact, rather then the RHS saying LOSS or DILUTE, shouldn't it just be a blank, i.e.

dilute ; 02 = 

[codecov]

Codecov Report

Merging #408 into master will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #408   +/-   ##
=======================================
  Coverage   88.84%   88.84%           
=======================================
  Files          17       17           
  Lines        2259     2259           
=======================================
  Hits         2007     2007           
  Misses        252      252

Flag	Coverage Δ
#build	62.55% <ø> (ø)	⬆️
#tests	88% <ø> (ø)	⬆️
#unittests	35.81% <ø> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a74ec61...35c7274. Read the comment docs.

[spco]

Ok, please see #412 for what I believe is a working version. It may be helpful to manually grab that commit's contents and use it to replace the contents here (that way we easily keep all this conversation with the eventually merged PR).

I have only briefly tested this, because I don't have any real testcases where DILUTE is not set to NOTUSED (it would be very useful to such a testcase to the behaviour testsuite - please feel free to add one!). So please let me know if this doesn't work for any reason. The only slight wrinkle I see is that it doesn't update the line near the bottom of mechanism.f90:

!* End of Subset.  No. of Species = 30, No. of Reactions = 71 ;

to reflect the new reactions created. But that line is not used anywhere, it's just cruft from the original file, so I don't see an issue. (In the case I used, the count of species is out by one anyway! - not sure why).

[rs028]

Hi @spco, for the purpose of updating the manual, is this correct?

• adds DILUTE process for all species in the chemical mechanism, if DILUTE in environmentVariables.config is NOT set to NOTUSED, otherwise leaves the mechanism as is.

• If DILUTE is set to zero, the chemical mechanism is modified (i.e., 0 is not the same as NOTUSED).

• in principle, DILUTE can be constrained.

[spco]

1 and 2 are correct, but DILUTE cannot be constrained as of now, because the current implementation hard-codes the value from environmentVariables.config into mechanism.f90. You've highlighted that I could tweak things easily to allow the use of constrained values though - I'll make a change to #412 to reflect that.

[spco]

Hi @bsn502 - you'll see that the full test fails on both the Linux and Mac builds. This is discussed at #412 (comment) . @rs028 how would you like to proceed with this? Until it's fixed, Travis will always fail. Options I see are to
(a) tweak full to give a problem that can actually be solved rather than stalling;
(b) remove full from the testsuite by deleting its folder;
(c) half-remove full by leaving the folder there, but tweaking the test runner script to not run full;
(d) something else.

(a) is possibly a lot of work, (b) means we lose the big test, (c) means full is there for future use if we wanted to reinstate it with a bit of work. I guess both (b) and (c) might need some changes to the manual one way or the other.

Personally, maybe (b) is the right option for now - we can always manually reinstate full from earlier commits if need be, if it will add value - what do you all think?

[spco]

Oh, and thanks for all your continuing work @bsn502 !

[rs028]

@spco yes I think that is sensible. We want to revamp the testsuite anyways at some point.

@bsn502 Can you delete the full/ directory in travis/tests/ please?

I guess we don't need dilute.py anymore, as well. And build/build_atchem2.sh should revert to its original version?

outdated

[rs028]

Looks okay to me.
@spco One question: how do we deal with acknowledgements? I think Beth name should be on this file. The copyright notice is probably not the right place. Previously I have added a small comment (see line 15 in tools/plot/plot-atchem2_v3.py)

[spco]

Sure, I think either solution would be fine - perhaps the acknowledgement line is the right way to go?

[spco]

Probably best to squash and merge given the convoluted nature of the commit history in this PR :)

[rs028]

@bsn502 can you add an acknowledgement line to the header similar to line 15 in tools/plot/plot-atchem2_v3.py?

[rs028]

Is this ready to be merged? @spco does it need a rebase?

[spco]

Looks like it to me. I would squash and merge using the GitHub GUI, given the slightly circuitous commit history to what is essentially changes to one file and the removal of a directory.