Ideas for KPP 4.0.0

[RolfSander]

Add your ideas for KPP 4.0.0 here. Everything we don't have the time to work on now but may be a nice addition in the future...

• Replace all BLAS/LAPACK functions (WAXPY, WSCAL, etc.) in integrator files.

• Delete #WRITE_OPT and WriteOptions() in debug.c.

• Use ICNTRL in tau_leap and gillespie integrators.

• Add Python (wrapper) and Julia as new languages to KPP.

• Have the Fortran90 code create a derived type with global variables.

• Formal uncertainty quantification and propagation through an integration step.

• Delete #LUMP?

• Delete #TRANSPORT and #TRANSPORTALL

• Update obsolete fortran code, e.g., arithmetic IF statements

• Enlarge some limits so that the complete MCM can be used

• Add #DEBUG command

• Write input file checker, see Input file checker #95

• Use Y instead of C in Update_RCONST(), see Are we using Update_RCONST() correctly? #93

...

[RolfSander]

Use ICNTRL in tau_leap and gillespie integrators

[RolfSander]

add Python and Julia as new languages to KPP

[yantosca]

Have the Fortran90 code create a derived type with global variables

[msl3v]

I have two different KPP models in development. Features for a v4.0.0 list will definitely come out of this.

[yantosca]

At IGC10 I learned that @drewpendergrass (student in my group) has a Python wrapper for KPP (that works on generated F90 code). We can try to add this to KPP 3.1.0 or KPP 4.0.0. Wrapping the F90 code would probably be the quickest way to get KPP functionality in Python and/or Julia (and it'd probably be faster than executing native code).

[RolfSander]

This sounds great! Is the interface based on f2py?

https://numpy.org/devdocs/f2py

[drewpendergrass]

Yes it's based on f2py! The wrapper is quite basic at this point, but it makes it easy to do box modeling within a Python environment -- which is much nicer than messing around in fortran.

[RolfSander]

Nice :-)

I've used f2py for other (pretty small) projects, and my experience has been very good so far!

[jimmielin]

Based on the discussion above I currently have this list that we can include in the outlook paragraph of the manuscript:
(1) adding support for modern languages such as Python and Julia;
(2) refactoring of the generated code to avoid global data structures for easier parallelization;
(3) streamlining inputs and outputs of all integrators for consistency; and
(4) supporting the adaptive solver option in other integrators and programming languages.

Feel free to add more to the list so we can include some potential future directions for the next version of KPP.

[drewpendergrass]

Just noting that I am happy to bottomline the Python/Julia wrapper in KPP4.

[jimmielin]

Thanks @drewpendergrass! It would be awesome to have KPP built with a wrapper so it can interface with modern languages more easily 😄

[RolfSander]

Thanks, @jimmielin, the list looks good. We may have to add the MCM if we are not able to finish that project (#4) for KPP-3.0.0:

(5) Improve interaction and compatibility with the Master Chemical Mechanism (http://mcm.york.ac.uk)

[msl3v]

Preprocess and make multiple mechanisms available in one set of modules. e.g by defining multiple *.eqn files in one *.kpp file, and being provided a module interface that lets the user pick which mechanism to run within a program.

[RolfSander]

@msl3v: This is very similar to what we already have in the MECCA system: A python preprocessing script that extracts selected reactions from a comprehensive kpp file. And also a selection of the desired mechanism at runtime.

[msl3v]

@RolfSander Is the python method you reference similar to how MECCA operated (or still operates?) to build mechanisms? Are the goals still the same? e.g. describe and process trop or strat, or MBL regimes?

What I was describing above would actually pre-process multiple mechanism definitions within one set of modules, e.g. multiple definitions of LU_IROW, LU_ICOL, multiple back-substitution solutions, Fun(), etc. Possibly all using the same parameter definitions (species, etc).

The goal wouldn't necessarily be achieving efficiency. Rather to expand the achievable complexity in boolean environments, e.g. in cloud, at a system boundary, etc. where process splitting would otherwise undermine the result.

edit: basically it would permit consistent step-wise changes in problem size

[msl3v]

Next request/brainstorm: create an add-on infrastructure. A standard method of adding standardized libraries/packages that expand system capacity. e.g. biogeochem functions (I've got a prototype enzyme mechanism built with KPP) that would be developed independently but explicity for KPP.

edit: another example, a package for Julia, rather than building it in directly.

[RolfSander]

@msl3v: I was referring to xmecca. This used to be a tcsh/sed/awk script mixture but now it has been converted into a consistent python script xmecca.py.

[msl3v]

Yes! I remember xmecca.

…

On 7/21/22 13:38, Rolf Sander wrote: @msl3v <https://github.com/msl3v>: I was referring to xmecca. This used to be a tcsh/sed/awk script mixture but now it has been converted into a consistent python script xmecca.py. — Reply to this email directly, view it on GitHub <#43 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVE2347W3DN7TERX4G4FUDVVGDJ7ANCNFSM5WY6NKSQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[RolfSander]

We still have a couple of WAXPY, WSCAL, etc. commands in our integrator
files. I think that these BLAS/LAPACK functions can all be replaced by
Fortran 90 commands. Something to keep in mind for KPP-4.0.0...

[RolfSander]

Delete #WRITE_OPT and WriteOptions() in debug.c because it prints some wrong values. Most probably because it depends on some use* variables which are not yet defined. The correct values can be found in the logfile.

[msl3v]

Idea: Formal uncertainty quantification and propagation through an integration step.

[RolfSander]

@msl3v: I added your suggestion to the to-do list at the top of this page.

[RolfSander]

Is anyone using the #LUMP command? If not, can we delete it? Looking at the KPP source code, I can't even see if it is already implemented properly. Maybe it doesn't work at all?

[msl3v]

Not here.

[msl3v]

Idea: variable stoichiometric coeffs

[jimmielin]

Idea: Tagged names for reactions

CAM-chem's preprocessor has this. It's handy so that reaction rates, rate constants, etc. when output in a model can have nice names instead of RxnRate_EQ001, RxnRate_EQ002, ... we have now in GEOS-Chem.

Their definitions look something like

[HO2NO2_OH]            HO2NO2 + OH  -> H2O + NO2 + O2                           ; 4.5e-13, 610 
[N_NO2a]               N + NO2  -> N2O + O                                      ; 2.9e-12, 220 
[N_NO2b]               N + NO2  -> 2*NO                                         ; 1.45e-12, 220 

We could probably add this at the end of each line, similar to how comments are done in the mechanism right now in GEOS-Chem:

O3 + MO2 = CH2O + HO2 + O2 :                 GCARR_ac(2.90d-16, -1000.0d0);                                                    {name=O3_MO2; 2014/02/03; Eastham2014; SDE}
OH + OH = H2O + O :                          1.80d-12;                                                                         {name=OH_OH; 2014/02/03; Eastham2014; SDE}

We would just need to find a place for the model to retrieve this, simple mapping of the reaction id# to the "tag"/name string.

[msl3v]

This is a good idea. CAABA/MECCA does this, though with external scripts, I think. Right now, KPP ignores anything in {} brackets so reactions can be tagged this way. It may be possible to put a toggle switch that lets KPP see the reaction tags.

[yantosca]

KPP should already be able to use equation tags in < > already, right? At least we have it in the documentation example

#EQUATIONS { Small Stratospheric Mechanism }
<R1>  O2   + hv = 2O            : (2.643E-10) * SUN*SUN*SUN;
<R2>  O    + O2 = O3            : (8.018E-17);
<R3>  O3   + hv = O   + O2      : (6.120E-04) * SUN;
<R4>  O    + O3 = 2O2           : (1.576E-15);
<R5>  O3   + hv = O1D + O2      : (1.070E-03) * SUN*SUN;
<R6>  O1D  + M  = O   + M       : (7.110E-11);
<R7>  O1D  + O3 = 2O2           : (1.200E-10);
<R8>  NO   + O3 = NO2 + O2      : (6.062E-15);
<R9>  NO2  + O  = NO  + O2      : (1.069E-11);
<R10> NO2  + hv = NO  + O       : (1.289E-02) * SUN;

[msl3v]

It would be useful to use these tags in place of a generic name in *_Parameters.*, for sure.

[RolfSander]

As suggested by @yantosca, we can probably use the already existing equation tags between < and >. I'm not sure if the underscore is already allowed but it could be added easily.

[jimmielin]

Very minor and very esoteric feature, "empty" products in reactions

CAM-chem uses a similar feature for hard coding some artificial sinks into the mechanism.
e.g., SF6 + hv = sink or soa + hv = . Tried building a mechanism with empty products and it throws a syntax error (as expected) with Misplaced ';'

Probably could just use one of the DEFFIX species as product to work around this but just wanted to put it here in case there was further interest...

[RolfSander]

For products you don't care about, you can use the pre-defined dummy
species PROD, e.g.:

SF6 + hv = PROD

https://kpp.readthedocs.io/en/stable/using_kpp/04_input_for_kpp.html#equations

[jimmielin]

Thanks @RolfSander! I should've looked closely at the documentation first, my bad 😅

[RolfSander]

We should update obsolete fortran code, e.g., arithmetic IF
statements. The gfortran compiler already complains about this. In util/blas.f90, we currently have:

IF (incX .EQ. incY) IF (incX-1) 5,20,60 

[RolfSander]

Now that #TRANSPORT is deleted, I had a look at other rarely used
options. I noticed that there is a #WRITE_OPT command which isn't even
mentioned in our readthedocs manual. When activated, it calls
WriteOptions() in debug.c. It creates logging output which is nice
but redundant. We already have the same output in GenerateLog() in
gen.c.

@yantosca: Should we delete #WRITE_OPT?

[yantosca]

Hi @RolfSander, thanks for checking the old options. I think we can delete #WRITE_OPT. I'm not sure it was in the PDF manual, since that's what I took as the starting point for the RTD manual.

I just pushed PR #88 to remove LUMP. I can make another PR to remove #WRITE_OPT plus any other options that may no longer be used.

[RolfSander]

I don't think that #WRITE_OPT has ever been in any KPP manual. Yes,
let's delete it.

There are 3 similar options: WRITE_ATM, WRITE_SPC and
WRITE_MAT. They also create potentially interesting output for
debugging. Unlike #WRITE_OPT, they are not duplicated elsewhere, thus
we may want to keep them. If we keep them, we could mention them in the
RTD manual.

While I'm looking at these options, I decided to check all of them.
Here's what I've found:

• I think these can be deleted:

  • #DEFRAD
  • #SPARSEDATA (scanner.c says that it is deprecated)
  • #USE (scanner.c says that it is deprecated)
  • #USES (seems to occur only in *.el)

• #AUTOREDUCE isn't mentioned in the manual yet.

• #DECLARE is missing in *.el

• #FAMILIES appears twice in *.el

• For these, I have no idea:

  • #INITIALIZE
  • #RUN
  • #XGRID
  • #YGRID
  • #ZGRID

[yantosca]

Thanks @RolfSander. I would keep the #WRITE_{ATM,SPC_MAT} options for now. We can get rid of the others. I'll also add a note about #AUTOREDUCE documentation (based on @jimmielin's paper).

It looks like the #INITIALIZE, etc. might have been from running KPP as a model on a 3-D grid. I think we can also let those go.

[RolfSander]

OK, sounds good. Let me know if there's anything I can do...

[RolfSander]

I forgot to mention #SETRAD. Like #DEFRAD, it is also obsolete.

[RolfSander]

I suggest to enlarge some limits so that the complete MCM can be used. Increasing MAX_NO_OF_LINES should not cause any problems. Regarding MAX_EQN, we need to check if the new value is okay for all machines (except, of course, for those machines which already have to use a smaller value). You can find my first tests in the mcm branch.

[RolfSander]

I suggest to add a #DEBUG command. When set to ON, KPP should create a log file with lots of information.

[95moon]

Dear Professor Rolf Sander,

I hope this message finds you well. I am a beginner, currently a second-year Master's student in Atmospheric Science. Based on the KPP website example, I understand that running a box model requires four initial files: eqn, spc, dev, and a KPP file. However, in the MCM example you provided, I only seem to see settings for the initial model parameters in the driver_mcm document. After following the instructions in the README and running the model, I found that the mcm.exe did not output any results, which left me quite confused. Could you please provide more details on this issue? If possible, feel free to contact me via email at [email protected]. Additionally, my programming skills are limited to Python.

Thank you very much for your time and help.

Best regards,
XinJie Liu

[yantosca]

@95moon could you please open a new issue with your question? That will allow us to keep this issue related to ideas for KPP 4.0.0.