erkvar%Integrate will change the initial time and final time?

[CRquantum]

Hi bro,

A quick question, a piece of my code is like below,

    write(6,*) 'before =', x0_test,xf_test
    call erkvar%Init(me, MAX_STEPS, Method=ERK_DOP853, ATol=[atol], RTol=[rtol],&
        InterpOn=.false., EventsOn=.FALSE. )
    if (erkvar%status == FLINT_SUCCESS) then   
        stiffstatus = stifftestval
        stepsz = stepsz0       
        call erkvar%Integrate(x0_test, y0_test, xf_test, yf_test, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
            IntStepsOn=.false.,Xint = Xint, Yint = Yint, &
            EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus)        
    endif  
    write(6,*) 'after =', x0_test,xf_test

You know the x0_test and xf_test are the initial time and final time for the integral.
I expect that before and after call erkvar%Integrate, x0_test and xf_test will not change.

However, I notice that sometimes x0_test and xf_test will be changed after call erkvar%Integrate .

I checked erkvar%status, and, errstr.
I found erkvar%status = 5 , and
errstr is "Maximum integration steps reached"

It looks like, in this case FLINT did not complete the intergration, so x0_test and xf_test are somehow changed.

In such case, increasing max_steps does not really help.

1. I wonder, so in such case, do you have some suggestions? Like how to complete the integral correctly? Do I use a different integral method instead of the default ERK_DOP853?
2. If FLINT cannot complete the integral correctly, does that basically mean that the ODE does not have a numerical solution?

Thanks a lot!
Best regards,
Rong

[CRquantum]

Oh, one more question,
if I want to let the program stop when ODEs are stiff, does the following script looks correct?
Thanks a lot in advance!

    call erkvar%Init(me, MAX_STEPS, Method=ERK_DOP853, ATol=[atol], RTol=[rtol],&
        InterpOn=.false., EventsOn=.FALSE. )
    if (erkvar%status == FLINT_SUCCESS) then   
        stiffstatus = stifftestval
        stepsz = stepsz0       
        call erkvar%Integrate(x0_test, y0_test, xf_test, yf_test, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
            IntStepsOn=.false.,Xint = Xint, Yint = Yint, &
            EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus)        
    endif  
    if (stiffstatus == -1) then 
        write(6, *) 'problem is stiff'
        stop
    endif

[princemahajan]

If you are seeing this error message "Maximum integration steps reached" then it means FLINT is returning the last integrated solution in yf_test and the corresponding time or independent variable in xf_test.

If you increase the max. steps allowed and still seeing the same message then it is just probably because it needs more steps to reach the final time. You can check this by comparing the returned xf_test value before and after increasing the max steps to see if the integration was advanced further at all or not.

[princemahajan]

for setting the stiffness test option, see the FLINT documentation here what value to pass to FLINT and what it will return. Here is excerpt from that:

---

[in,out]	stifftest	FLINT will use Hairer's DOP853 alogorithm to test for stiffness of the differential equations if StffTest == 1 or StffTest == 2. Test is disabled by default or if Stifftest == 0.StiffTest == 1: If DiffEq are detected to be stiff then -1 is returned in StiffTest and the integration halts at the current step with appropriate status code.StiffTest == 2: If DiffEq are detected to be stiff then -1 is returned in StiffTest but the integration continues.

[CRquantum]

Thank you my man!
Yeah, if I increase the max_steps to a very big number, I do see the xf_test increased.
eg, I want the range to be [24,48], with max_steps=10^5, xf_test stopped at about 3.5. With max_steps=10^8, xf_test increased to 35. However, in this case, it may show
erkvar%status = 9, Step size reduced to very small value
So the integration again may fail. because increasing max_steps does not help anymore.

In this case, do you have any suggestions?

LIke perhaps this indicate it is a stiff problem or something?

Or, do I decrease the integral range?

Thanks much in advance!

Happy new year by the way! LOL!

[princemahajan]

Happy new year, pal!
Yes your diffeq does seem like a stiff system. Use the stiffest option with value 1 and this will halt the integration as soon as FLINT detects the diffeq to be stiff from that point on. If that is the case, you just have to try solvers for the stiff system. FLINT as of now has only nonstiff solvers.

[CRquantum]

Thank you bro!

May I ask, if it is a stiff problem, does stiffstatus have to return exactly -1?

I mean, like,

stiffstatus = 1   
call erkvar%Integrate(x0_test, y0_test, xf_test, yf_test, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
            IntStepsOn=.false.,Xint = Xint, Yint = Yint, &
            EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus)        

You see, I set stiffstatus to 1 and send it to StiffTest in the call erkvar%Integrate(...)
So I expect that the integral should halt if stiff problem is detected.

Now, after call erkvar%Integrate(...), stiffstatus is still 1 actually. However the error code shows "Maximum integration steps reached" or whatever.
In this case, does that mean the problem is actually non-stiff?

I mean does stiffstatus have to return exactly -1 so that it is a stiff system?

But anyway, I can use DVODE for stiff problem.

[princemahajan]

Can you post your differential equations here? It may be possible that your FLINT stiff test (which is basically same as in Hairers codes) is not able to detect stiffness for your equations if it is not returning -1 in stifftest. If you use a stiff solver and it has no issues then it will be clear that yours is a stiff system and FLINT stiff test is not working as it should for these equations. Get Outlook for Android<https://aka.ms/ghei36>

…

________________________________ From: Chen Rong ***@***.***> Sent: Friday, January 7, 2022 5:21:26 PM To: princemahajan/FLINT ***@***.***> Cc: Bharat Mahajan ***@***.***>; Comment ***@***.***> Subject: Re: [princemahajan/FLINT] erkvar%Integrate will change the initial time and final time? (Issue #11) Thank you bro! May I ask, if it is a stiff problem, does stiffstatus have to return exactly -1? I mean, like, stiffstatus = 1 call erkvar%Integrate(x0_test, y0_test, xf_test, yf_test, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, & IntStepsOn=.false.,Xint = Xint, Yint = Yint, & EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus) So set stiffstatus to 1 and send it to StiffTest in the call erkvar%Integrate(...) So I expect that the integral should halt if stiff problem is detected. Now, after call erkvar%Integrate(...), stiffstatus is still 1 actually. However the error code shows "Maximum integration steps reached" or whatever. In this case, does that mean the problem is actually non-stiff? I mean does stiffstatus have to return exactly -1 so that it is a stiff system? But anyway, I can use DVODE for stiff problem. — Reply to this email directly, view it on GitHub<#11 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AMK5YH3P6JAELMGKENDQZJDUU5YPNANCNFSM5LOEOC4A>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you commented.Message ID: ***@***.***>

[CRquantum]

Thank you very much.
Sure, if you have time you could download the whole VS project illustration code here,
https://gitlab.com/CRquantum/flint_check
It is VS2017 + Intel OneAPI.

Or, you could simply save the text below as a f90 file as a drive file (combine it with all yours files in the src folder in FLINT), just like your test_CR3BP.f90.

module MyDiffEq

    use FLINT
    use, intrinsic :: IEEE_ARITHMETIC
    implicit none 
    
    ! user diff eq system
    type, extends(DiffEqSys) :: PM_ODE_Sys
        !Prim
        real(wp) :: Ka
        real(wp) :: Vmax0
        real(wp) :: Km
        real(wp) :: Vc0
        real(wp) :: FA1
        real(wp) :: KCP
        real(wp) :: KPC
        !Cov
        real(wp) :: wt
        real(wp) :: age
        !Sec
        real(wp) :: VM  
        real(wp) :: V 
    contains
        procedure :: F => YP
    end type PM_ODE_Sys    
contains      
    function YP(me, X, Y, Params)
    implicit none
    intrinsic :: size  
    class(PM_ODE_Sys), intent(inout) :: me !< Differential Equation object
    real(WP), intent(in) :: X
    real(WP), intent(in), dimension(me%n) :: Y 
    real(WP), intent(in), dimension(:), optional :: Params
    real(WP), dimension(size(Y)) :: YP
    !Diffeqs
    !write (6,*) 'ODEs called'
    YP(1) = -me%Ka*Y(1)
    YP(2) = me%Ka*Y(1) - me%VM/(me%Km*me%V+Y(2))*Y(2) - me%KCP*Y(2) + me%KPC*Y(3)
    YP(3) = me%KCP*Y(2) - me%KPC*Y(3)    
    !write (6,*) 'ODEs = ', YP
    return
    end function YP
          
end module MyDiffEq

module model
    use FLINT
    use MyDiffEq   
    implicit none
    ! Turn on the events
    logical, private, parameter :: EventsEnabled = .FALSE.
    logical, private, dimension(3), parameter :: EvMask = [.FALSE.,.FALSE.,.FALSE.] 
    ! Event step size
    real(WP), private, parameter :: evstepsz = 0.00001_WP
    ! event tolerance
    real(WP), private, dimension(3) :: evtol = [0.001_WP,1.0e-9_WP,1.0e-9_WP]
    ! scalar tolerance
    real(kind=wp), private, parameter :: rtol   = 1.0e-4_WP  ! originally it is 10^-11.
    real(kind=wp), private, parameter :: atol   = rtol*1.0e-3_WP    
    ! max no. of steps
    integer, private, parameter :: MAX_STEPS = 10**4 ! decreae this can increase speed.
    logical, private, parameter :: CONST_STEPSZ = .FALSE. 

    contains
    
    subroutine my_ode_init(me)
    class(PM_ODE_Sys) :: me
    
    me%n = 3 ! # of ODEs
    me%m = 0 ! # of events, 0, so turn it off.
    
    me%Ka = 6.51750402383093
    me%Vmax0 = 23.1441646658992
    me%Km = -4.83193907550415
    me%Vc0 = 3.48392293236760
    me%FA1 = 0.543133339226148
    me%KCP = 5.80211181194353 
    me%KPC = 8.18019402422159
    
    me%wt = 59.05249
    me%age = 17.57132
    
    me%VM = me%Vmax0 * me%wt**0.75
    me%V = me%Vc0 * me%wt   
    return
    end subroutine my_ode_init
    
    subroutine test()
    real(kind=wp) :: x0_test, xf_test
    real(kind=wp), dimension(3) :: y0_test, yf_test ! here neq=3    
    real(kind=wp) :: stepsz0, stepsz
    integer :: stifftestval, stiffstatus
    real(kind=wp), dimension(:), allocatable, target :: Xint, Xarr
    real(kind=wp), dimension(:,:), allocatable, target :: Yint, Yarr
    real(kind=wp), allocatable, dimension(:,:) :: EventStates
    character(len=:), allocatable :: errstr
    integer :: fcalls, naccpt, nrejct
    type(ERK_class) erkvar
    type(PM_ODE_Sys) :: My_obj    
    
    stepsz0 = 0.0E-3_wp    ! let FLINT compute the initial step-size    
    stifftestval = 1  ! check for stiffness and stop integration if stiff 
    
    stiffstatus = stifftestval
    stepsz = stepsz0   
    
    x0_test = 2.0_wp
    xf_test = 24.0_wp
    y0_test = [0.0_wp, 155.921971503913_wp, 101.637238186467_wp]
         
    write(6,*) 'integral range [x0, xf] =', x0_test, xf_test
    write(6,*) 'inital condition y0 =', y0_test
    
    call my_ode_init(My_obj) ! initalize the ODE system My_obj
    
    write(6,*) My_obj%n, My_obj%m, My_obj%Ka, My_obj%Vmax0, My_obj%Km, My_obj%Vc0, My_obj%FA1, My_obj%KCP, My_obj%KPC &
        , My_obj%wt, My_obj%age, My_obj%VM, My_obj%V 
    write(6,*) My_obj%F(0.0_WP,y0_test)
    
    call erkvar%Init(My_obj, MAX_STEPS, Method=ERK_DOP853, ATol=[atol], RTol=[rtol],&
        InterpOn=.false., EventsOn=.FALSE. ) 
    
    if (erkvar%status == FLINT_SUCCESS) then  
        write(6,*) 'FLINT Init success'      
        call erkvar%Integrate(x0_test, y0_test, xf_test, yf_test, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
            IntStepsOn=.false.,Xint = Xint, Yint = Yint, &
            EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus)   
        if (erkvar%status == FLINT_SUCCESS) then   
            write(6,*) 'FLINT success'           
            call erkvar%Info(stiffstatus, errstr, nAccept=naccpt, nReject=nrejct, nFCalls=fcalls)        
            print *, fcalls, naccpt, nrejct        
        else
            write(6,*) 'FLINT fail with status: ', erkvar%status
            write(6,*) 'integral range [x0, xf] =', x0_test, xf_test
            write(6,*) 'y0 =', y0_test
            write(6,*) 'yf =', yf_test
            error stop
        endif         
          
    else
        write(6,*) 'init fail'
    endif    

    return
    end subroutine test

end module model
    

    program TestFLINT_CR
    use model   
    implicit none  
    call test()
    stop ('Program ends normally!')
    end program TestFLINT_CR

It is just a 3 ODE system,

#Diffeq
XP(1) = -Ka*X(1)
XP(2) = Ka*X(1)-VM/(Km*V+X(2))*X(2)-KCP*X(2)+KPC*X(3)
XP(3) = KCP*X(2) - KPC*X(3)

where all the parameters are defined below. They are initialized in subroutine my_ode_init(me).

#Primary parameters
Ka
Vmax0
Km
Vc0
FA1 ! this is not used.
KCP
KPC

#Covariates
wt
age

#Secondary parameters
VM = Vmax0 * wt**0.75
V = Vc0 * wt

It looks like these 3 ODEs will show status 5 error, looks like a stiff problem.
However DVODE can handle this problem without error and can do interpolation too, see
https://fortran-lang.discourse.group/t/has-anyone-used-dvode-and-is-it-good/2281/23?u=crquantum

Wish FLINT could handle stiff problem better :)

Best regards,
Rong