AtChem / AtChem2

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

Dilute fix #408

Closed bsn502 closed 4 years ago

bsn502 commented 5 years ago

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.

dilute_test_v2

rs028 commented 5 years ago

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 commented 5 years ago

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.

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 commented 5 years ago

(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 https://github.com/AtChem/AtChem2/pull/408#issuecomment-546301638) 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 commented 5 years ago

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 https://github.com/AtChem/AtChem2/blob/a74ec610bfa3ecb5d32fda1e484ded345989cb4d/src/inputFunctions.f90#L869 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 commented 5 years ago

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

dilute ; 02 = 
codecov[bot] commented 5 years ago

Codecov Report

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

Impacted file tree graph

@@           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% <ø> (ø) :arrow_up:
#tests 88% <ø> (ø) :arrow_up:
#unittests 35.81% <ø> (ø) :arrow_up:

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 commented 4 years ago

Ok, please see https://github.com/AtChem/AtChem2/pull/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 commented 4 years ago

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

spco commented 4 years ago

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 commented 4 years ago

Hi @bsn502 - you'll see that the full test fails on both the Linux and Mac builds. This is discussed at https://github.com/AtChem/AtChem2/pull/412#issuecomment-571563440 . @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 commented 4 years ago

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

rs028 commented 4 years ago

@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?

spco commented 4 years ago

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

spco commented 4 years ago

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

rs028 commented 4 years ago

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

rs028 commented 4 years ago

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

spco commented 4 years ago

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.