KineticPreProcessor / KPP

The KPP kinetic preprocessor is a software tool that assists the computer simulation of chemical kinetic systems
GNU General Public License v3.0
21 stars 11 forks source link

Toggle calling of Update_RCONST, Update_PHOTO, Update_SUN within integrators -- Closes #9 #24

Closed yantosca closed 2 years ago

yantosca commented 2 years ago

Overview

This PR implements the feature request in #9. We now can use the setting of KPP runtime option ICNTRL(15) to determine whether or not rate constants can be updated within the integrator itself.

    ! Determine the settings of the Do_Update_* flags, which determine
    ! whether or not we need to call Update_* routines in the integrator
    ! (or not, if we are calling them from a higher-level)
    ! ICNTRL(15) = -1 ! Do not call Update_* functions within the integrator
    !            =  0 ! Status quo
    !            =  1 ! Call Update_RCONST from within the integrator
    !            =  2 ! Call Update_PHOTO from within the integrator
    !            =  3 ! Call Update_RCONST and Update_PHOTO from w/in the int.
    !            =  4 ! Call Update_SUN from within the integrator
    !            =  5 ! Call Update_SUN and Update_RCONST from within the int.   
    !            =  6 ! Not implemented```

NOTE: For global models such as GEOS-Chem, you will want to use ICNTRL(15) = -1, as rates are computed outside of the integrator in order to avoid computational overhead. But for other implementations using KPP (such as box models), it may be advantageous to be able to update rates at each call to the integrator.

Validation

I ran the continuous integration tests manually at these 2 commits:

  1. Dev (with modifications) @ commit a9d2794e7a7611092edc0615a68dca8ebe4b6d29 : Dev CI test output log

  2. Ref (parent commit) @ commit 109ed9213860bc4357b45170e4eb92766681b407 : Ref CI test output log

These files are identical, thus proving that the default settings (ICNTRL(15) = 0) yield the same results.

yantosca commented 2 years ago

Comments updated! Will merge on top of dev shortly.

@RolfSander, feel free to pull in the patch-memory branch to dev as well.

RolfSander commented 2 years ago

@yantosca, I'm afraid I have to re-open this discussion because I realized that we need to distinguish between ICNTRL and ICNTRL_U.

ICNTRL defines the integrator-specific default which needs to be added for all compilers. This can be done in the "fine-tune the integrator" section. To keep the current behaviour unchanged, I think the correct values are:

ICNTRL(15) = 0 for kpp_radau5.f90 and kpp_seulex.f90

ICNTRL(15) = 7 for kpp_dvode.f90, kpp_lsode.f90, runge_kutta_adj.f90, runge_kutta.f90 and runge_kutta_tlm.f90

ICNTRL(15) = 5 for all other integrators

In contrast to ICNTRL, the array ICNTRL_U contains the user-supplied values. If the user provides a zero, the default is kept. This is done with the following line:

WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)

Here we have the next prolem: When the user supplies ICNTRL_U(15)=-1, the value in ICNTRL is not changed at all. So we should change this line to:

WHERE(ICNTRL_U(:) /= 0) ICNTRL(:) = ICNTRL_U(:)

Once these corrections are made, I think we don't need the "IF (option==0)" block in util.f90 anymore.

Sorry I didn't see all this earlier...

yantosca commented 2 years ago

@yantosca, I'm afraid I have to re-open this discussion because I realized that we need to distinguish between ICNTRL and ICNTRL_U.

ICNTRL defines the integrator-specific default which needs to be added for all compilers. This can be done in the "fine-tune the integrator" section. To keep the current behaviour unchanged, I think the correct values are:

ICNTRL(15) = 0 for kpp_radau5.f90 and kpp_seulex.f90

ICNTRL(15) = 7 for kpp_dvode.f90, kpp_lsode.f90, runge_kutta_adj.f90, runge_kutta.f90 and runge_kutta_tlm.f90

ICNTRL(15) = 5 for all other integrators

In contrast to ICNTRL, the array ICNTRL_U contains the user-supplied values. If the user provides a zero, the default is kept. This is done with the following line:

WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)

Here we have the next prolem: When the user supplies ICNTRL_U(15)=-1, the value in ICNTRL is not changed at all. So we should change this line to:

WHERE(ICNTRL_U(:) /= 0) ICNTRL(:) = ICNTRL_U(:)

Once these corrections are made, I think we don't need the "IF (option==0)" block in util.f90 anymore.

Sorry I didn't see all this earlier...

No worries @RolfSander I'll make these fixes.

yantosca commented 2 years ago

Also @RolfSander, I noticed that some integrators had the calls to the Update* commands commented out. Should we uncomment those all? I didn't want to break any functionality.

yantosca commented 2 years ago

This should now be fixed with commit f3a92ba292a200dafacc9b3c8b17f79e46ab9499

RolfSander commented 2 years ago

Yes, I think we should uncomment all calls to the Update_* commands. This should not break any functionality as long as ICNTRL(15) is defined correctly.

I've done that. I noticed that a few of the integrators had

    !~~~> if optional parameters are given, and if they are >0,
    !     then use them to overwrite default settings
    IF (PRESENT(ICNTRL_U)) ICNTRL(1:20) = ICNTRL_U(1:20)
    IF (PRESENT(RCNTRL_U)) RCNTRL(1:20) = RCNTRL_U(1:20)

instead of the WHERE statements. That will just overwrite everything regardless of whether it's /= 0, so it's not quite an apples-to-apples comparison w/ the new code. Just FYI.