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

laestrada · 2021-10-13T14:27:34Z

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

The text was updated successfully, but these errors were encountered:

yantosca · 2021-11-01T19:22:09Z

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 · 2021-11-01T19:24:57Z

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 · 2021-11-01T19:47:01Z

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 · 2021-11-02T08:25:48Z

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 · 2021-11-11T13:51:40Z

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 · 2021-11-11T14:14:37Z

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 · 2021-11-11T15:34:07Z

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 · 2022-04-14T14:33:28Z

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 · 2022-04-14T14:41:24Z

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 · 2022-04-14T15:18:19Z

@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 · 2022-04-14T15:29:02Z

I'm not familiar with JAC_CHEM (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 · 2022-04-14T21:48:33Z

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

yantosca · 2022-04-14T22:10:39Z

@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 · 2022-04-15T08:55:00Z

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 · 2022-04-15T13:28:43Z

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 · 2022-04-15T13:31:49Z

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

RolfSander · 2022-04-15T14:00:31Z

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 · 2022-04-15T14:19:09Z

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 · 2022-04-15T14:26:49Z

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 · 2022-04-15T14:47:01Z

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 L_UPDATE_*.

yantosca · 2022-04-15T15:22:44Z

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 · 2022-04-15T16:22:35Z

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 · 2022-04-15T16:37:14Z

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 · 2022-04-19T13:54:49Z

Updated comments in commit f34285a. Now merged into dev.

We had passed Do_Update_Rates to down to FunTemplate and JacTemplate as a crude way of trying to toggle Update_RCONST on/off. This is now superseded by the updates in KineticPreProcessor/KPP PR #9. We have now deleted this obsolete code from int/rosenbrock.f90. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

laestrada added the feature New feature or request label Oct 13, 2021

yantosca added the integrators Related to numerical integrators label Oct 13, 2021

yantosca mentioned this issue Apr 15, 2022

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

Merged

yantosca linked a pull request Apr 15, 2022 that will close this issue

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

Merged

yantosca closed this as completed Apr 19, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

laestrada commented Oct 13, 2021 •

edited

Loading

yantosca commented Nov 1, 2021

jimmielin commented Nov 1, 2021

msl3v commented Nov 1, 2021 •

edited

Loading

RolfSander commented Nov 2, 2021

msl3v commented Nov 11, 2021

msl3v commented Nov 11, 2021

RolfSander commented Nov 11, 2021

yantosca commented Apr 14, 2022

yantosca commented Apr 14, 2022 •

edited

Loading

yantosca commented Apr 14, 2022

RolfSander commented Apr 14, 2022

msl3v commented Apr 14, 2022

yantosca commented Apr 14, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022 •

edited

Loading

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

yantosca commented Apr 19, 2022

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

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

Comments

laestrada commented Oct 13, 2021 • edited Loading

yantosca commented Nov 1, 2021

jimmielin commented Nov 1, 2021

msl3v commented Nov 1, 2021 • edited Loading

RolfSander commented Nov 2, 2021

msl3v commented Nov 11, 2021

msl3v commented Nov 11, 2021

RolfSander commented Nov 11, 2021

yantosca commented Apr 14, 2022

yantosca commented Apr 14, 2022 • edited Loading

yantosca commented Apr 14, 2022

RolfSander commented Apr 14, 2022

msl3v commented Apr 14, 2022

yantosca commented Apr 14, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022 • edited Loading

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

RolfSander commented Apr 15, 2022

yantosca commented Apr 15, 2022

yantosca commented Apr 19, 2022

laestrada commented Oct 13, 2021 •

edited

Loading

msl3v commented Nov 1, 2021 •

edited

Loading

yantosca commented Apr 14, 2022 •

edited

Loading

yantosca commented Apr 15, 2022 •

edited

Loading