[FEATURE REQUEST] Add toggle for calling UPDATE_SUN, UPDATE_RCONST within integrator step

[laestrada]

Original issue thread found here. Thread content copied below: yantosca commented on Jun 7 RolfSander wrote:

"Indeed, calling UPDATE_RCONST several times (i.e., once for each KPP substep) is often not necessary. It's good to have the option to switch that off. Maybe this could be implemented via a currently not-used element in the ICNTRL array, e.g.:"

IF (ICNTRL(15)==1) THEN
    CALL Update_SUN()
    CALL Update_RCONST()
ENDIF

yantosca commented on Jun 7 RolfSander, thanks for this. I had removed those calls from the integrator files but this is an easy way to be able to toggle these calls. I'll add this into the next KPP update.

RolfSander commented on Jun 8 After taking a closer look at these subroutines, I think that Update_SUN is probably never necessary. I guess it would be sufficient to use the ICNTRL(15) switch only for Update_RCONST:

IF (ICNTRL(15)==1) THEN
    CALL Update_RCONST()
ENDIF

yantosca commented on Jun 10 This is now done in commit e7540d3. All of the F90 integrator modules have this IF statement that can toggle UPDATE_RCONST if the user wishes. This will go into the next KPP version. I will close out this issue now.

yantosca commented on Jun 11 My updates actually broke KPP so I went back to the GC_updates branch and started a new dev branch. In commit 342a4ac I added the call to Update_RCONST if ICNTRL(15)==1 only to the F90 Rosenbrock integrator file (int/rosenbrock.f90) since that is what GEOS-Chem users. As time allows, I'll make the similar edits into the other F90 integrator files.

It was more complex than I thought, I had to pass a logical flag with the value of (ICNTRL(15)==1) that down a few levels from the ros_Integrator routine. Straightforward but I didn't have time to do that for the rest of the code. I'll reopen this issue to denote that this still has to be done for the other integrator files.

RolfSander commented on Jun 13 Hmm, I didn't realize that this would be so complex. An alternative for a quick implementation in the other solvers might be to declare

LOGICAL :: Do_Update_Rates

as a module-wide variable, i.e., before the CONTAINS statement of the module

[yantosca]

Thinking out loud... wondering if it would be better to use C-preprocessor switches to shunt this behavior as such:

#ifndef NO_UPDATE_SUN_IN_INTEGRATOR
   CALL UPDATE_SUN()
#endif
#ifndef NO_UPDATE_RCONST_IN_INTEGRATOR
   CALL Update_RCONST()
#endif

My thinking is that using Fortran IF statements would incur some extra CPU overhead every time the integrator is called. For models like GEOS-Chem, the integrator is called for each grid box in the chemistry domain multiplied by the number of internal KPP timesteps, so the time spent in IF statements might add up to non-negligible amount of time. Using #ifdefs would also be less intrusive to the existing code.

Thoughts? @msl3v, @jimmielin, @laestrada, @RolfSander

[jimmielin]

Hi Bob, I agree preprocessor statements may be more performance efficient. Eliminating code at the compiler level is going to be faster than IF-branching at runtime for sure, if there is read code that GEOS-Chem or another model will certainly never use.

[msl3v]

Agree with the performance efficiency perspective. It would be extremely useful, though, if we could update 'some' rates, and just not all of them. This would permit concentration-based rate calculations essential for several parameterized rate forms.

There's always the option to write vectors to *_Parameters.[c,F,F90] controlling which rates are updated. If empty vector, no calcs, toggled via a preprocessor.

[RolfSander]

Yes, it will be faster with C-preprocessor switches, but I wonder if it really matters. Maybe run some benchmark tests first, with and without the Fortran IF, and check if the gain in performance is significant.

Also, if we use C-preprocessor directives, we lose the option to control this at run time. For example, one may want to activate Update_RCONST() only during sunrise and sunset, when the changes are big.

[msl3v]

The IF toggle makes the most sense. The condition of sunrise & sunset is a good example. I'll do a performance test of this with an IF statement in a small boxmodel under development. I agree it will likely be a small impact.

[msl3v]

P.S. if controlling through ICNTRL(), what abpout the following? ICNTRL(X) = 1 ! Update_RCONST() ICNTRL(X) = 2 ! Update_PHOTO() ICNTRL(X) = 3 ! Update both

[RolfSander]

Thanks, Mike, for performing the tests!

I like ICNTRL(X) = 1 for RCONST and 2 for PHOTO. I don't think "3 for both" is necessary because the subroutine Update_RCONST already includes Update_PHOTO.

Before choosing the array index X for ICNTRL(X), let's check that it is not used by any of the integrators.

BTW: There is also a subroutine called Update_SUN. I think it exists only for testing purposes, and no-one uses it. Maybe we can delete it.

[yantosca]

I like ICNTRL(X) = 1 for RCONST and 2 for PHOTO. I don't think "3 for both" is necessary because the subroutine Update_RCONST already includes Update_PHOTO.

We can go with this. I can now start working on this.

Before choosing the array index X for ICNTRL(X), let's check that it is not used by any of the integrators.

I think we said ICNTRL(15) is OK, but I will confirm.

BTW: There is also a subroutine called Update_SUN. I think it exists only for testing purposes, and no-one uses it. Maybe we can delete it.

I'm not sure if any of the CI tests use it. I'll double check that too.

[yantosca]

Grepping for ICNTRL in the int folder, we get

$ grep -in "icntrl(15)" *.f90
$ grep -in "icntrl(16)" *.f90
$ grep -in "icntrl(17)" *.f90
$ grep -in "icntrl(18)" *.f90
$ grep -in 'icntrl(19)" *.f90

so these are the only slots of ICNTRL that are free in all integrators. I'll use ICNTRL(15).

[yantosca]

@RolfSander @msl3v: question: most of the integrators have calls to Update_RCONST and Update_PHOTO in both the KPP FUN_CHEM and JAC_CHEM routines. I'm assuming that we will want to also test for ICNTRL(15) in JAC_CHEM as well, correct?

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE JAC_CHEM (T, V, JF)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    USE KPP_ROOT_Parameters
    USE KPP_ROOT_Global
    USE KPP_ROOT_JacobianSP
    USE KPP_ROOT_Jacobian, ONLY: Jac_SP
    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO

    IMPLICIT NONE

    KPP_REAL :: V(NVAR), T, Told
#ifdef FULL_ALGEBRA    
    KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
    INTEGER :: i, j 
#else
    KPP_REAL :: JF(LU_NONZERO)
#endif   

    !Told = TIME
    !TIME = T
    !CALL Update_SUN()
    !CALL Update_RCONST()
    !CALL Update_PHOTO()
    !TIME = Told

#ifdef FULL_ALGEBRA    
    CALL Jac_SP(V, FIX, RCONST, JV)
    DO j=1,NVAR
      DO i=1,NVAR
         JF(i,j) = 0.0d0
      END DO
    END DO
    DO i=1,LU_NONZERO
       JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
    END DO
#else
    CALL Jac_SP(V, FIX, RCONST, JF) 
#endif   

    Njac=Njac+1

  END SUBROUTINE JAC_CHEM

[RolfSander]

I'm not familiar with JACCHEM (the Rosenbrock solvers don't have it). As long as the default setting for ICNTRL(15) leaves the `Update*` calls unchanged, anything is fine with me.

[msl3v]

@yantosca, I also don't recognize JAC_CHEM(). Can you say where this is?

[yantosca]

@msl3v: I think it is in kpp_dvode.f90. But now that I look at it again, it has calls to Update_SUN and Update_RCONST commented out internally, so maybe we don't need to worry about it.

Also, I've been doing:

ICNTRL(15) = 1 ! Update_RCONST() only
ICNTRL(15) = 2 ! Update_RCONST() then call Update_PHOTO()
ICNTRL(15) = 3 ! Update_SUN()    then call Update_RCONST()

The Update_SUN() function is needed for the C-I tests to run, as it is called from the drv/general*.f90 driver routines.

[RolfSander]

Sounds good. And with ICNTRL(15) = 0, the same UPDATE_* subroutines are called as in the current version of our integrators, whatever that may be. Right?

[yantosca]

Right now I have ICNTRL(15) = 0 turning off calls to the Update_* routines in the integrators but I of course can change that. Maybe this should be:

ICNTRL(15) = -1 : Turn off all Update_* functions w/in integrators
ICNTRL(15) = 0 : Status quo, call all Update_* functions as in the current code
ICNTRL(15) = 1: Update_RCONST only
ICNTRL(15) = 2: Update RCONST then Update_PHOTO
ICNTRL(15) = 3: Update_SUN then Update_RCONST

I think we can do that with a couple of IF statements:

IF ( ICNTRL(15) == 3 ) CALL Update_SUN()
IF ( ICNTRL(15) > -1 ) CALL Update_RCONST()
IF ( ICNTRL(15) == 2 ) CALL Update_PHOTO()

so Update_RCONST is called for any positive value for ICNTRL(15), and the other update functions are only called for specific values of ICNTRL(15).

I think the IF statements (without ELSEs) should be pretty efficient. It's the ELSE block where you get bogged down as the CPU has to try to predict which branch to take.

Also I will run diff tests w/r/t the prior dev branch to make sure the C-I tests give the same results.

[yantosca]

Oops that will work for the Update_RCONST but not the others. Let me think about this a bit.

[RolfSander]

Let me think about this a bit.

Yes, a "bit" could be the solution :-)

If we think of ICNTRL(15) as a three-bit number, we can assign each bit to one of the UPDATE_* calls, e.g.:

Update_RCONST() = 1 = OOI
Update_PHOTO()  = 2 = OIO
Update_SUN()    = 4 = IOO

Then, if the user defines for example ICNTRL(15) = 5 = IOI, then the calls to Update_SUN and Update_RCONST are activated.

In Fortran90, I think this could be coded as:

LOGICAL :: L_Update_RCONST, L_Update_PHOTO, L_Update_SUN

L_Update_RCONST = IAND(ICNTRL(15), 1))
L_Update_PHOTO  = IAND(ICNTRL(15), 2))
L_Update_SUN    = IAND(ICNTRL(15), 4))

IF (L_Update_RCONST) CALL Update_RCONST
IF (L_Update_PHOTO ) CALL Update_PHOTO
IF (L_Update_SUN   ) CALL Update_SUN

For ICNTRL(15)=0 and ICNTRL(15)=-1, special treatment would be necessary:

IF (ICNTRL(15)=0) THEN
  L_Update_RCONST = ...  !  default value
  L_Update_PHOTO  = ...  !  default value
  L_Update_SUN    = ...  !  default value
ENDIF

IF (ICNTRL(15)<0) THEN
  L_Update_RCONST = .FALSE. 
  L_Update_PHOTO  = .FALSE. 
  L_Update_SUN    = .FALSE. 
ENDIF

[yantosca]

I was just looking into this. The only thing is I don't want to add a lot of IF's into the integrators as that will slow things down. Let me percolate on this...

[RolfSander]

Indeed, lots of IFs in the main loop should be avoided. Most of the commands, however, can be put into the initialization.

Then we only have one IF (L_Update_*) for each Update_* inside the time loop.

In addition, accessing a logical for IF(L_Update_*) is probably faster than selecting an element out of an array for IF(ICNTRL(15)).

[RolfSander]

I think I was a bit imprecise in my previous post. We do need to evaluate ICNTRL(15) in every KPP call because the user may want to change it during the model run (e.g., at sunrise or sunset).

However, during the KPP-calculated substeps, ICNTRL(15) does not change, and we could use LUPDATE*.

[yantosca]

So what I have been working on now looks like this:

SUBROUTINE INTEGRATE( TIN, TOUT, &
  ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )

   IMPLICIT NONE

   KPP_REAL, INTENT(IN) :: TIN  ! Start Time
   KPP_REAL, INTENT(IN) :: TOUT ! End Time
   ! Optional input parameters and statistics
   INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
   KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
   INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
   KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
   INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U

   KPP_REAL :: RCNTRL(20), RSTATUS(20)
   INTEGER       :: ICNTRL(20), ISTATUS(20), IERR

   INTEGER, SAVE :: Ntotal = 0

   ICNTRL(:)  = 0
   RCNTRL(:)  = 0.0_dp
   ISTATUS(:) = 0
   RSTATUS(:) = 0.0_dp

    !~~~> fine-tune the integrator:
   ICNTRL(1) = 0        ! 0 - non-autonomous, 1 - autonomous
   ICNTRL(2) = 0        ! 0 - vector tolerances, 1 - scalars

   ! If optional parameters are given, and if they are >0,
   ! then they overwrite default settings.
   IF (PRESENT(ICNTRL_U)) THEN
     WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
   END IF
   IF (PRESENT(RCNTRL_U)) THEN
     WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
   END IF

   ! Determine the settings of the L_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)
   CALL Integrator_Update_Options( ICNTRL(15) )

   ! Call the Rosenbrock solver
   CALL Rosenbrock(NVAR,VAR,TIN,TOUT,   &
         ATOL,RTOL,                &
         RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)

with:

SUBROUTINE Integrator_Update_Options( option )

!    option      -> determine whether to call Update_RCONST, Update_PHOTO,
!                   and Update_SUN from within the integrator
!        = -1 :   Do not call Update_* functions within the integrator
!        =  0 :   Status quo: Call whichever functions are normally called
!        =  1 :   Call Update_RCONST from within the integrator
!        =  2 :   Call Update_PHOTO from within the integrator
!        =  3 :   Call Update_RCONST and Update_PHOTO from within the int.
!        =  4 :   Call Update_SUN from within the integrator
!        =  5 :   Call Update_SUN and Update_RCONST from within the int.

!~~~> Input variable
  INTEGER, INTENT(IN) :: option

  ! Option -1: turn off all Update_* calls within the integrator
  IF ( option == -1 ) THEN
     Do_Update_RCONST = .FALSE.
     Do_Update_PHOTO  = .FALSE.
     Do_Update_SUN    = .FALSE.
     RETURN
  ENDIF

  ! Option 0: status quo: Call update functions if defined
  IF ( option == 0 ) THEN  
     Do_Update_RCONST = .TRUE.
     Do_Update_PHOTO  = .TRUE.
     Do_Update_SUN    = .TRUE. 
     RETURN
  ENDIF

  ! Otherwise determine from the value passed
  Do_Update_RCONST = ( IAND( option, 1 ) > 0 )
  Do_Update_PHOTO  = ( IAND( option, 2 ) > 0 )
  Do_Update_SUN    = ( IAND( option, 4 ) > 0 )

END SUBROUTINE Integrator_Update_Options

and the Do_Update_* variables are now global module variables. So we do allow the user to change them at each call to the Integrate function but not during individual calls to FunTemplate or JacTemplate. I think this is the desired behavior that you were describing.

[RolfSander]

Looks good!

• If SUBROUTINE Integrator_Update_Options is the same for all integrators, we may be able to put it into util/util.f90.

• Although you don't mention option=6 and option=7 in the comment, they are also valid possibilities.

[yantosca]

Nice! I'll put it into util/util.f90.

Option 6 and 7 are valid possibilities but I don't think any of the integrators have them.

[yantosca]

Updated comments in commit f34285ad87916c279edccf426358d23dc9f74c91. Now merged into dev.