Bug: Forward Euler integrator does not return ending timestep via RSTATUS · Issue #135 · KineticPreProcessor/KPP

While helping @cschooling debug a problem with the GEOS-Chem carbon gases mechanism in box model mode, I discovered that the Forward Euler integrator does not return the ending timestep to the calling program.

In the Integrate routine in the int/feuler.f90 module, the ending timestep (TEND) should be saved as RSTATUS(1) but is not:

KPP/int/feuler.f90

Lines 36 to 117 in b717a0d

    
           SUBROUTINE Integrate( TIN,       TOUT,      ICNTRL_U , RCNTRL_U,           & 
        
                                 ISTATUS_U, RSTATUS_U, IERR_U                        ) 
        
             ! ICNTRL(16) 
        
             ! 0 -> do nothing. 
        
             ! 1 -> set negative values to zero 
        
             ! 2 -> return with error code 
        
             ! 3 -> stop at negative 
        
             ! 
        
             ! ICNTRL(17) = verbose error output 
        
             !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
        
             ! User interface routine to the KPP forward Euler integrator 
        
             !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
        
             USE KPP_ROOT_Util, ONLY : Integrator_Update_Options 
        
             !~~~> Inputs 
        
             KPP_REAL, INTENT(IN)            :: TIN           ! Start Time 
        
             KPP_REAL, INTENT(IN)            :: TOUT          ! End Time 
        
             INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)  ! Input options 
        
             KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)  ! Input options 
        
             !~~~> Outputs 
        
             INTEGER,  INTENT(IN),  OPTIONAL :: ISTATUS_U(20) ! Returned status values 
        
             KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) ! Returned status values 
        
             INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U        ! Error code 
        
             ! Local variables 
        
             INTEGER                         :: ICNTRL(20) 
        
             KPP_REAL                        :: RSTATUS(20) 
        
             INTEGER                         :: IERR 
        
             !~~~> Zero input and output arrays for safety's sake 
        
             ICNTRL     = 0 
        
             RSTATUS    = 0.0_dp 
        
             !~~~> fine-tune the integrator: 
        
             ICNTRL(1 ) = 0       ! Verbose error output 
        
             ICNTRL(2 ) = 0       ! Stop upon negative results 
        
             ICNTRL(15) = 5       ! Call Update_SUN and Update_RCONST from w/in the int. 
        
            !~~~> if optional parameters are given, and if they are /= 0, 
        
            !     then use them to overwrite default settings 
        
            IF ( PRESENT( ICNTRL_U ) ) THEN 
        
               WHERE( ICNTRL_U /= 0 ) ICNTRL = ICNTRL_U 
        
            ENDIF 
        
            !~~~> 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 ! Call Update_SUN and Update_PHOTO from within the int. 
        
            !            =  7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int. 
        
            CALL Integrator_Update_Options( ICNTRL(15),          & 
        
                                            Do_Update_RCONST,    & 
        
                                            Do_Update_PHOTO,     & 
        
                                            Do_Update_Sun       ) 
        
            !~~~> In order to remove the prior EQUIVALENCE statements (which 
        
            !~~~> are not thread-safe), we now have declared VAR and FIX as 
        
            !~~~> threadprivate pointer variables that can point to C. 
        
            VAR => C(1:NVAR) 
        
            FIX => C(NVAR+1:NSPEC) 
        
            !~~~> Call the integrator 
        
            CALL ForwardEuler( NVAR, C(1:NVAR), TIN, TOUT, ICNTRL, IERR ) 
        
            !~~~> Free pointers 
        
            VAR => NULL() 
        
            FIX => NULL() 
        
            !~~~> Return error status (NOTE: ISTATUS_U does nothing) 
        
            IF ( PRESENT( IERR_U    ) ) IERR_U    = IERR 
        
            IF ( PRESENT( RSTATUS_U ) ) RSTATUS_U = RSTATUS 
        
           END SUBROUTINE Integrate

The problem is that the general box-model driver relies on the value of RSTATUS(1) as returned by the integrator in order to get the ending timestep of the integration. This value is then compared against TEND to see if it is time to exit the loop, as shown here:

KPP/drv/general.f90

Lines 35 to 61 in b717a0d

    
           !~~~> Time loop 
        
                 T = TSTART 
        
           kron: DO WHILE (T < TEND) 
        
                   TIME = T 
        
           !~~~> Write values of monitored species at each iteration 
        
                   CALL GetMass( C, DVAL ) 
        
                   WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,       & 
        
                              ( TRIM(SPC_NAMES(MONITOR(i))),           & 
        
                                C(MONITOR(i))/CFACTOR, i=1,NMONITOR ) 
        
                   CALL SaveData() 
        
           !~~~> Update sunlight intensity and reaction rates 
        
           !~~~> before calling the integrator. 
        
                   CALL Update_SUN() 
        
                   CALL Update_RCONST() 
        
           !~~~> Call the integrator 
        
                   CALL INTEGRATE( TIN       = T,        & 
        
                                   TOUT      = T+DT,     & 
        
                                   ICNTRL_U  = ICNTRL,   & 
        
                                   RCNTRL_U  = RCNTRL,   & 
        
                                   RSTATUS_U = RSTATE   ) 
        
                   T = RSTATE(1) 
        
                 END DO kron

But because the RSTATUS(1) is not returned, then the loop will execute infinitely.

I will make a PR to fix this and will also add a C-I test for the Forward Euler integrator in a subsequent PR.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bug: Forward Euler integrator does not return ending timestep via RSTATUS #135

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	SUBROUTINE Integrate( TIN, TOUT, ICNTRL_U , RCNTRL_U, &
	ISTATUS_U, RSTATUS_U, IERR_U )
	! ICNTRL(16)
	! 0 -> do nothing.
	! 1 -> set negative values to zero
	! 2 -> return with error code
	! 3 -> stop at negative
	!
	! ICNTRL(17) = verbose error output

	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	! User interface routine to the KPP forward Euler integrator
	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	USE KPP_ROOT_Util, ONLY : Integrator_Update_Options

	!~~~> Inputs
	KPP_REAL, INTENT(IN) :: TIN ! Start Time
	KPP_REAL, INTENT(IN) :: TOUT ! End Time
	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) ! Input options
	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) ! Input options

	!~~~> Outputs
	INTEGER, INTENT(IN), OPTIONAL :: ISTATUS_U(20) ! Returned status values
	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) ! Returned status values
	INTEGER, INTENT(OUT), OPTIONAL :: IERR_U ! Error code

	! Local variables
	INTEGER :: ICNTRL(20)
	KPP_REAL :: RSTATUS(20)
	INTEGER :: IERR

	!~~~> Zero input and output arrays for safety's sake
	ICNTRL = 0
	RSTATUS = 0.0_dp

	!~~~> fine-tune the integrator:
	ICNTRL(1 ) = 0 ! Verbose error output
	ICNTRL(2 ) = 0 ! Stop upon negative results
	ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int.

	!~~~> if optional parameters are given, and if they are /= 0,
	! then use them to overwrite default settings
	IF ( PRESENT( ICNTRL_U ) ) THEN
	WHERE( ICNTRL_U /= 0 ) ICNTRL = ICNTRL_U
	ENDIF

	!~~~> 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 ! Call Update_SUN and Update_PHOTO from within the int.
	! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int.
	CALL Integrator_Update_Options( ICNTRL(15), &
	Do_Update_RCONST, &
	Do_Update_PHOTO, &
	Do_Update_Sun )

	!~~~> In order to remove the prior EQUIVALENCE statements (which
	!~~~> are not thread-safe), we now have declared VAR and FIX as
	!~~~> threadprivate pointer variables that can point to C.
	VAR => C(1:NVAR)
	FIX => C(NVAR+1:NSPEC)

	!~~~> Call the integrator
	CALL ForwardEuler( NVAR, C(1:NVAR), TIN, TOUT, ICNTRL, IERR )

	!~~~> Free pointers
	VAR => NULL()
	FIX => NULL()

	!~~~> Return error status (NOTE: ISTATUS_U does nothing)
	IF ( PRESENT( IERR_U ) ) IERR_U = IERR
	IF ( PRESENT( RSTATUS_U ) ) RSTATUS_U = RSTATUS

	END SUBROUTINE Integrate

	!~~~> Time loop
	T = TSTART
	kron: DO WHILE (T < TEND)

	TIME = T

	!~~~> Write values of monitored species at each iteration
	CALL GetMass( C, DVAL )
	WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T, &
	( TRIM(SPC_NAMES(MONITOR(i))), &
	C(MONITOR(i))/CFACTOR, i=1,NMONITOR )
	CALL SaveData()

	!~~~> Update sunlight intensity and reaction rates
	!~~~> before calling the integrator.
	CALL Update_SUN()
	CALL Update_RCONST()

	!~~~> Call the integrator
	CALL INTEGRATE( TIN = T, &
	TOUT = T+DT, &
	ICNTRL_U = ICNTRL, &
	RCNTRL_U = RCNTRL, &
	RSTATUS_U = RSTATE )
	T = RSTATE(1)

	END DO kron

Bug: Forward Euler integrator does not return ending timestep via RSTATUS #135

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions