Open CRquantum opened 2 years ago
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
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.
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:
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!
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.
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.
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 Androidhttps://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 GitHubhttps://github.com/princemahajan/FLINT/issues/11#issuecomment-1007815085, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AMK5YH3P6JAELMGKENDQZJDUU5YPNANCNFSM5LOEOC4A. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://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: @.***>
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
Hi bro,
A quick question, a piece of my code is like below,
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.
Thanks a lot! Best regards, Rong