Integrator, Rosenbrock 3

[Giorgio4766]

Ask a question about KPP:

Hi everyone, I am setting up an experiment where I need to be sure that I am using the Rosenbrock 3 integrator in order to make results comparable to another dataset. In the .def file I specified the following command,

INTEGRATOR rosenbrock

but that should refer to "rodas3" which I am not sure refers to Rosenbrock 3. Can you please suggest me what #INTEGRATOR command is needed to use the Ros3 integrator?

Thank you!

[yantosca]

Thanks for writing @Giorgio4766. Options specific to each integrator are included in the comments in the Fortran integrator modules (e.g. ROOT_Integrator.F90). You can use ICNTRL and RCNTRL to pass options. For the Rosenbrock integrator, ICNTRL(3) is the setting for the selection of numerical method:

!
!    Note: For input parameters equal to zero the default values of the
!       corresponding variables are used.
!
!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
!
!    ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
!              = 1: AbsTol, RelTol are scalars
!
!    ICNTRL(3)  -> selection of a particular Rosenbrock method
!        = 0 :    Rodas3 (default)
!        = 1 :    Ros2
!        = 2 :    Ros3
!        = 3 :    Ros4
!        = 4 :    Rodas3
!        = 5 :    Rodas4
!
!    ICNTRL(4)  -> maximum number of integration steps
!        For ICNTRL(4)=0) the default value of 200000 is used
!
!    ICNTRL(12)  -> use auto-reduce solver? set threshold in RCNTRL(12)
!    ICNTRL(13)  -> ... append slow species when auto-reducing?
!    ICNTRL(14) -> choose a target species instead for determining threshold?
!                  if yes, specify idx. then RCNTRL(12) is obsolete.
!
!    ICNTRL(15) -> Toggles calling of Update_* functions w/in the integrator
!        = -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 the int.
!
!    ICNTRL(16) -> 
!        = 0 : allow negative concentrations (default)
!        = 1 : set negative concentrations to zero
!
!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
!          It is strongly recommended to keep Hmin = ZERO
!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
!    RCNTRL(3)  -> Hstart, starting value for the integration step size
!
!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
!                          (default=0.1)
!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
!         than the predicted value  (default=0.9)
!
!    RCNTRL(12) -> threshold for auto-reduction (req. ICNTRL(12)) (default=100)
!    RCNTRL(14) -> AR threshold ratio (default=0.01)

(note: rosenbrock_autoreduce shown above... there are some extra settings in ICNTRL(12:14) and RCNTRL(12,14) but the rest is the same.

[Giorgio4766]

Dear Yantosca, thanks a lot for your quick reply. I see that Integrator.f90 file contains, at the beginning of the file, ICNTRL is initialized at zero: ICNTRL(:) = 0

meaning that also ICNTRL(3)=0 so Rodas3 is also set as integrator. To change the integrator to Ros3 I would do this

!~~~>   Initialize the particular Rosenbrock method selected

     ICNTRL(3)=2
  SELECT CASE (ICNTRL(3))
     CASE (1)
       CALL Ros2
     CASE (2)
       CALL Ros3
     CASE (3)
       CALL Ros4
     CASE (0,4)
       CALL Rodas3
     CASE (5)
       CALL Rodas4
     CASE DEFAULT
       PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) 
       CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
       RETURN
   END SELECT

fixing ICNTRL(3)=2 directly in the Integrator.f90 file. Would that be correct?

Thanks,

Giorgio.

[yantosca]

Thanks for replying @Giorgio4766. So what happens is that in routine INTEGRATE the ICNTRL and RCNTRL values are passed in as ICNTRL_U and RCNTRL_U:

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

and then further down in the file nonzero values of ICNTRL_U and RCNTRL_U are copied to ICNTRL and RCNTRL:

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

   !~~~> fine-tune the integrator:
   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
   IF ( PRESENT( RCNTRL_U ) ) THEN
      WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U
   ENDIF

so whatever you specify in ICNTRL and RCNTRL from the calling routine to the routine INTEGRATE will be preserved. There should be no need to hardwire the value in the *_Integrator.F90 file.

You can also use a debug print statement in *_Integrator.F90 file to verify that the ICNTRL setting is what it should be.

[yantosca]

@Giorgio4766: So you would want to call Integrate like this:

INTEGER :: ICNTRL(20)
...
ICNTRL = 0
ICNTRL(3) = 2
CALL Integrate( Tin, Tout, ICNTRL_U=ICNTRL, ... )

[Giorgio4766]

Dear Bob, thank you a lot for your timely replies. I will then act on the general.f90 file where the SUBROUTINE INTEGRATE is called and add the crucial command ICNTRL(3) = 2 . I will let you know how it goes then.

I am glad to read that you are based in Cambridge (MA). I was there last summer and have great memories from that trip.

Danke,

Cheers,

Giorgio.

[Giorgio4766]

Dear Bob, I am afraid I was not still able to configure the solver to Rosenbrock3, (Ros3). At least I am not sure it went trough.

I exectued the command ./../bin/kpp dynho.kpp (dynho is the prefix of the experiment). I then checked that CALL Integrate( Tin, Tout, ICNTRL_U .. ) is present in general.f90 and in ./dynho_Main.f90 . I operated only on dynho_Main.f90 to include

INTEGER :: ICNTRL(20) ... ICNTRL = 0 ICNTRL(3) = 2 CALL Integrate( Tin, Tout, ICNTRL_U=ICNTRL, ... )

as you previously suggested. I then compiled with no errors. Unfortunately I don't know how to check whether or not ICNTRL, and ICNTRL_U has been made equal to ICNTRL(3)=2.

I hope you can help with that!

Thank you

[yantosca]

@Giorgio: you can put a print statement to check the values of ICNTRL in dyno_Integrator.F90. Then if you see its' OK you can delete it.

[Giorgio4766]

Dear Bob, done and dusted! I changed the third value of ICNTRL_U, which refers to the selection of the Rosenbrock solver, in _main.f90 . I then printed both values of ICNTRL and ICNTRL_U in the _Integrator.f90 and both show the wanted value (=2).

Thanks a lot for your help!