Adaptive solver / mechanism auto-reduction feature

[jimmielin]

Hi all,

I'm opening up this PR for the eventual merge of the auto-reduce feature into KPP 3.0.0, which @msl3v and I have been working on.

For some reason GitHub shows a few files changed that don't belong to this feature. It'll be necessary to merge this with a little manual work. As far as I know, the changes should be in these files only:

• int/rosenbrock_autoreduce.def
• int/rosenbrock_autoreduce.f90
• src/code.h (support LOGICAL)
• src/code_f90.c (support LOGICAL)
• src/gen.c
• src/scan.l
• src/scan.y
• src/scanner.c

To enable this, run kpp with an input file that has #AUTOREDUCE on. The #INTEGRATOR rosenbrock_autoreduce must be selected.

To use AR at runtime, ICNTRL(8) = .true. needs to be set. Two thresholding methods are supported:

1. RCTNRL(8) sets an absolute threshold in molec cm-3
2. or, ICNTRL(10) /= 0 chooses a target species ID to base AR threshold on; RCNTRL(10) is the factor that is multiplied to the target species' max(P, L) to be used as the threshold

The # for 8, 10 here may need to be changed to be standardized across integrators during the merge, noting that GEOS-Chem uses these as well and fullchem_mod.F90 needs to be changed accordingly.

[RolfSander]

It's great to see that we are now almost ready to merge the autoreduction feature into dev! However, the pull request shows a couple of conflicts that must be resolved first. I would suggest the following sequence of actions:

1) Merge both cleanup_int and docs into dev first. (I think cleanup_int and docs can then be deleted)

2) Merge dev into autoreduce. This will allow us to work on the merge conflicts inside the autoreduce branch without negatively affecting the dev branch.

I'm open for alternative ideas though. Let me know what you think...

[jimmielin]

Hi Rolf, I agree with this course of action. I think once cleanup_int and docs are merged into dev, this PR's changes will be clearer (a lot of the conflicts seem to be with the docs, probably because I rebased autoreduce over some commit in dev).

Then we can see if autoreduce can be merged in directly into dev. I don't see why it would negatively affect dev since we took care in implementing this into a separate integrator, rosenbrock_autoreduce.f90, so the sooner it's done the better before everything starts to diverge too far. Please let me know what you think!

[RolfSander]

Unless @yantosca or @msl3v have any objections, I'm happy to merge cleanup_int and docs into dev now. See:

https://github.com/KineticPreProcessor/KPP/pull/44

Let's see if there will be any remaining conflicts...

[msl3v]

No objections. (I'm only sparsely available for about a week, so I hope to get caught up with this in any case.)

On 6/21/22 15:47, Rolf Sander wrote:

Unless @yantosca https://github.com/yantosca or @msl3v https://github.com/msl3v have any objections, I'm happy to merge cleanup_int and docs into dev now. See:

44 https://github.com/KineticPreProcessor/KPP/pull/44

Let's see if there will be any remaining conflicts...

— Reply to this email directly, view it on GitHub https://github.com/KineticPreProcessor/KPP/pull/46#issuecomment-1162259453, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVE233OOHM4UTJ5HB32JCDVQIL6FANCNFSM5ZMXXFNQ. You are receiving this because you were mentioned.Message ID: @.***>

[yantosca]

Once PR #47 is approved then we can go ahead and merge cleanup_int into dev.

[RolfSander]

We have now merged both cleanup_int and docs into dev now.

@jimmielin: You can merge dev now into your autoreduce branch and check for any potential remaining conflicts there.

[jimmielin]

Thanks. I've attempted a merge of dev into autoreduce. I've cleared all the git conflicts, but I need to make sure the AR functionality is still intact and that it builds ok. GitHub still seems to show some unrelated files in doc being updated for some reason. I'll report back...

[jimmielin]

I'm getting a weird error when I attempt to build the GEOS-Chem mechanism under the GEOS-Chem source code's KPP directory.

./build_mechanism.sh fullchem

This is KPP-2.5.0.

KPP is parsing the equation file.
Fatal error : atoms: Can't read file
Program aborted
KPP failed to build 'gckpp_Rates.F90'! Aborting.

@yantosca have you seen this before? It seems to be looking for an atoms file which GEOS-Chem does not have. Maybe a new default has been introduced which makes it fail with the GEOS-Chem gckpp.kpp?

Update: Looks like the latest build in dev fails with GEOS-Chem as well, so this was not introduced by autoreduce

[RolfSander]

Sorry if we forgot to mention this: The file atoms has been renamed to atoms.kpp. Therefore, you need #INCLUDE atoms.kpp now.

[jimmielin]

Thanks Rolf! Looks like the following changes are needed for GEOS-Chem once this is merged:

• #UPPERCASEF90 isn't enabled in the 13.4.0 version I have yet
• #INCLUDE atoms in fullchem.eqn needs to be changed to #INCLUDE atoms.kpp.

• ohreact_parser.py trips up here:

  
   '        MGLY + OH --> PCO + MCO3 + CO                                                               ' /) ! index up to 120

with the "up to 120" wording. A fix needs to be made replacing

rn = int(spl[1].split('index ')[1])

with

            try:
                rn = int(spl[1].split('index ')[1])
            except ValueError:
                rn = int(spl[1].split('index up to ')[1])

I'm building this now and testing out the new generated code inside GEOS-Chem. Once I complete a full run, it looks like the merge should work. I don't know why the CI tests are failing (I haven't touched the other integrators at all and #AUTOREDUCE is off by default).

[yantosca]

Hi @jimmielin, the CI tests are failing because VAR and FIX aren't declared as POINTER in the _Global file. It looks like the code we had in gen.c to do that got clobbered by the merge.

[jimmielin]

Thanks @yantosca, fixed. Let's see if the CI tests pass. Also I had to make a fix to OHreact_parser.py - I'll make a separate PR for that over at the GEOS-Chem repo since it's a backwards compatible change.

[yantosca]

Also in GEOS-Chem 14.0.0 I've already fixed the build_mechanism.sh and the OHreact_parser.py files for KPP 2.5.0. I didn't bother prepping GEOS-Chem 13.4.* for KPP 2.5.0 since we are in the process of releasing the GEOS-Chem 14.0.0 version.

[yantosca]

@jimmielin -- if you want to run the C-I tests locally, do this:

cd $KPP_HOME/.ci-pipelines
./ci-manual-testing-script.sh

and then that should do the same thing as the C-I tests on Azure. (It won't run them in the container but that's good enough for finding bugs).

[yantosca]

and then ./ci-manual-cleanup-script.sh to remove all of the generated output from the local C-I tests.

[jimmielin]

Thanks Bob! That'll be useful, I'm still running into a segfault in the CI test. I'll debug it locally.

I've just noticed that you made the fix in GEOS-Chem 14 as well. Looks like this will be necessary only for the 13.4.0 testing version I have. I'll try to run the full model with and without auto-reduction, but we are getting close to making this mergeable.

[yantosca]

Great. What I like to do is to run the C-I tests on the branch that I'm working on, and then again on the parent branch (requires make distclean and rebuilding kpp) and then diff them to see if there are identical results. (I pipe the output of ./ci-manual-testing-script.sh to log files).

[jimmielin]

Thanks - should be all fixed now (fingers crossed). I had a D'oh moment where I had segfaults in a weird place (beginning of a function call) and then realized I needed ulimit -s unlimited on my local system.

[jimmielin]

Looking good. I noticed that these declarations are no longer present in dev:

//===========================================================================
// KPP 2.3.2_gc, Bob Yantosca (26 Mar 2021)
// Define extra arrays and scalars for gckpp_Global.F90
//
  MW = DefvElm( "MW",   real, -NSPEC,
        "Species molecular weight [g/mole]" );
  SR_MW = DefvElm( "SR_MW", real, -NSPEC,
           "Square root of species molecular weight [g/mole]" );
  SR_TEMP = DefElm( "SR_TEMP", real,
            "Square root of Temperature [K**0.5]" );
  K300_OVER_TEMP = DefElm( "K300_OVER_TEMP", real,
               "300.0 / Temperature [K]" );
  TEMP_OVER_K300 = DefElm( "TEMP_OVER_K300", real,
               "Temperature [K] / 300.0");
  NUMDEN = DefElm( "NUMDEN", real,
           "Air number density [#/cm3]");

I'll hack around them for my 13.4 code since I saw these were moved to commonIncludeVars.H in GEOS-Chem 14. We can sort out that part when the AR code is merged into GEOS-Chem. In fact I think the merge for the AR code in GEOS-Chem would just be for the changes made in fullchem_mod.F90 and input_opt_mod.F90 etc. for the options, everything else should be "vanilla" KPP now.

[yantosca]

Hi @jimmielin, yes those variables were from the version of KPP that was kind of "hardwired" for GEOS-Chem. They are now declared locally in GEOS-Chem in the commonIncludeVars.H file.

Also I noticed a couple of unused variable warnings in the C code:

gen.c:433:5: warning: unused variable ‘flxind’ [-Wunused-variable]
  433 | int flxind[MAX_EQN];
      |     ^~~~~~
gen.c:423:8: warning: unused variable ‘j’ [-Wunused-variable]
  423 | int i, j;
      |        ^

Removing the unused variables will clear the warnings.

[jimmielin]

Thanks @yantosca. I've backported some changes from GC 14 into my copy of GC 13.4 so it compiles and runs. I've also just updated rosenbrock_autoreduce, the Rosenbrock solver with auto-reduction feature, to be in sync with other changes made in rosenbrock.f90. Pending checks on the output I think this is ready to merge into dev. Hopefully with all the changes made now, it should be painless.

Edit: There seems to be some issues in testing - I'll update once this is fully mergeable

• Box model tests passed
• still testing full model... changes to VAR to C for thread safety need to be backported so the model runs as well with kpp 2.5.0 now testing output. runs ok as expected output ok. zero diff with no-AR or AR on

[yantosca]

Thanks @jimmielin. Feel free to merge into dev pending successful checks.

[jimmielin]

Thanks @yantosca, CI outputs also verified to match (diff shows some differences that are purely cosmetic). The repository requires a review prior to merge, so please feel free to go over the changed files and merge when ready.

➜  .ci-pipelines git:(autoreduce) ./ci-manual-testing-script.sh | tee autoreduce_out.txt
➜  src git:(dev) make distclean && make clean && make
➜  .ci-pipelines git:(dev) ./ci-manual-testing-script.sh | tee dev_out.txt
➜  .ci-pipelines git:(dev) diff dev_out.txt autoreduce_out.txt 
251d250
< gfortran -cpp -O -g  -c rk_Stoichiom.f90
253a253
> gfortran -cpp -O -g  -c rk_Stoichiom.f90
610d609
< gfortran -cpp -O -g  -c rktlm_Hessian.f90
612a612
> gfortran -cpp -O -g  -c rktlm_Hessian.f90
984d983
< gfortran -cpp -O -g  -c ros_Hessian.f90
985a985
> gfortran -cpp -O -g  -c ros_Hessian.f90
1704d1703
< gfortran -cpp -O -g  -c rosadj_Jacobian.f90
1706d1704
< gfortran -cpp -O -g  -c rosadj_LinearAlgebra.f90
1707a1706,1707
> gfortran -cpp -O -g  -c rosadj_LinearAlgebra.f90
> gfortran -cpp -O -g  -c rosadj_Jacobian.f90
1987,1988d1986
< gfortran -cpp -O -g  -c rostlm_Jacobian.f90
< gfortran -cpp -O -g  -c rostlm_LinearAlgebra.f90
1990a1989,1990
> gfortran -cpp -O -g  -c rostlm_Jacobian.f90
> gfortran -cpp -O -g  -c rostlm_LinearAlgebra.f90
2482d2481
< gfortran -cpp -O -g  -c sd_LinearAlgebra.f90
2483a2483
> gfortran -cpp -O -g  -c sd_LinearAlgebra.f90
2841d2840
< gfortran -cpp -O -g  -c sdadj_Hessian.f90
2843a2843
> gfortran -cpp -O -g  -c sdadj_Hessian.f90
2929d2928
< gfortran -cpp -O -g  -c small_f90_Stoichiom.f90
2931a2931
> gfortran -cpp -O -g  -c small_f90_Stoichiom.f90
3479a3480
> gfortran -cpp -O -g  -c small_strato_Hessian.f90
3482d3482
< gfortran -cpp -O -g  -c small_strato_Hessian.f90

[RolfSander]

I would like to run a few tests with my MECCA model, using the code from the autoreduce branch but with #AUTOREDUCE OFF. Hopefully I will get (almost) the same results as before. I think I can do these tests until tomorrow.

[jimmielin]

Thanks @RolfSander. Looking forward to your tests. I am actually curious what happens in MECCA if #AUTOREDUCE on and #INTEGRATOR rosenbrock_autoreduce is enabled as well, but I guess it would be a question for another day :smiley:

[RolfSander]

Yes, F77_* can be removed from rosenbrock_autoreduce.def. I also removed the entries which are default settings anyway.

[RolfSander]

I noticed that autoreduce uses additional ICNTRL (8,9,10), RCNTRL (8,10), and RSTATUS (4) elements. Can we add short descriptions of them to 04_input_for_kpp.rst and 05_output_from_kpp.rst in docs/source/using_kpp?

[jimmielin]

Thanks @RolfSander and @yantosca for the changes. I'll update the documentation with the new switch information for ICNTRL RCNTRL and RSTATUS(4) now, will push a commit shortly.

[jimmielin]

I've updated the documentation. rosenbrock_autoreduce is a long name so I had to change the table ReST text a little to make it look tidy, but I just added a few rows.

[RolfSander]

Thanks, @jimmielin, for adding the documentation!

In rosenbrock_autoreduce, the elements 8 and 10 of RCNTRL now have a different meaning than in the other integrators. Do you think it would make sense to use the elements 12 and 13 instead for the solver-specific settings of rosenbrock_autoreduce?

[jimmielin]

I agree, @RolfSander. I was worried that the switches' meaning would be inconsistent as well. I will move both ICNTRL and RCNTRL added to positions 12, 13, and 14, for consistency. I will push a commit shortly both here and over at GEOS-Chem.

[RolfSander]

@jimmielin: Thanks! Yes, 12, 13, and 14 are suitable positions.

@yantosca: How can I activate the updated *.rst files? Or are you doing this?

[jimmielin]

Hi @RolfSander, I've updated the indices. I think that the updated .rst files only propagate once this is pushed to the branch that docs are being built from, not sure if it's main or dev.

[yantosca]

@jimmielin: For now the ReadTheDocs output is being built from the docs branch but we can change that when we release.

[RolfSander]

Oh, sorry, I didn't know this. I deleted the docs branch after it was merged into dev. I guess we'll have to recreate it (or choose another branch to build the docs.

[yantosca]

For now I can point it to the main branch, so that the docs will be pegged to KPP 2.5.0 until we merge in the new features.

[jimmielin]

Verified that the changed version of KPP works correctly with GEOS-Chem as well over at https://github.com/jimmielin/geos-chem/tree/staging/autoreducekpp. I'm happy to merge this anytime.

[RolfSander]

Merging is fine for me now.

[yantosca]

Thanks @jimmielin! We should merge. Were we going to release this as an intermediate 2.6.0 version, or should we just hold off until 3.0.0?

[RolfSander]

I think we can call it 2.6.0 now. After that we can take care of the cleanup_util branch and check the documentation.

[yantosca]

Thanks. I can do the merge & do a review of the docs before pushing a release. Can fix any last-minute issues before release.

[yantosca]

I've pushed the KPP 2.6.0 release and made some minor updates to the doc (such as adding the reference to @jimmielin's paper in prep). But I had to manually create a DOI () since there was a bug in the .zenodo.json file (now fixed in cleanup_util).

On to KPP 3.0.0!

[RolfSander]

@jimmielin: I have two questions about rosenbrock_autoreduce.f90:

• Many variables are declared as REAL(kind=dp). Is this correct or should it be KPP_REAL?

• I think cWAXPY is never used. Can it e deleted?

[jimmielin]

Hi @RolfSander,

1. I think it can be replaced with KPP_REAL.
2. Yes. The version committed here has been highly optimized and is not the original implementation by @msl3v. I removed the original implementation (which used cWAXPY and had a much clearer flow in code...) to clean up the merge.

I'm happy to work on and submit a pull request to clean up the code and some cosmetic changes to remove trailing spaces etc which trip up git merge --check.

[RolfSander]

@jimmielin: Thanks!

Since we don't have the autoreduce branch anymore, I think you can put these small changes directly into dev (if @yantosca agrees).

[yantosca]

That is fine. We can push a 2.6.1 release, which would also fix the Auto-generation of the DOI.