How to pass a constant value in DVODE code, specially in user defined function “F”? 1

[Priyabrata Mohapatra]

Hello,
I'm using DVODE fortran 90 version. I'm trying to pass a constant value in variable form to the user defined function "F" [in this case it's FEX]. I'm experimenting on example problem,

dy1/dt = -.04*y1 + 1.e4*y2*y3
dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
dy3/dt = 3.e7*y2**2

on the interval from t = 0.0 to t = 4.e10, with initial conditions
y1 = 1.0, y2 = y3 = 0. The problem is stiff.

So example code given is,

PROGRAM EXAMPLE
!-----------------------------------------------------------------------
! EXAMPLE PROBLEM
!
! The following is a simple example problem, with the coding
! needed for its solution by DVODE. The problem is from chemical
! kinetics, and consists of the following three rate equations:
! dy1/dt = -.04*y1 + 1.e4*y2*y3
! dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
! dy3/dt = 3.e7*y2**2
! on the interval from t = 0.0 to t = 4.e10, with initial conditions
! y1 = 1.0, y2 = y3 = 0. The problem is stiff.
!
! The following coding solves this problem with DVODE, using MF = 21
! and printing results at t = .4, 4., ..., 4.e10. It uses
! ITOL = 2 and ATOL much smaller for y2 than y1 or y3 because
! y2 has much smaller values.
! At the end of the run, statistical quantities of interest are
! printed. (See optional output in the full description below.)
!
EXTERNAL FEX, JEX
DOUBLEPRECISION ATOL, RPAR, RTOL, RWORK, T, TOUT, Y
DIMENSION Y (3), ATOL (3), RWORK (67), IWORK (33)
NEQ = 3
Y (1) = 1.0D0
Y (2) = 0.0D0
Y (3) = 0.0D0
T = 0.0D0
TOUT = 0.4D0
ITOL = 2
RTOL = 1.D-4
ATOL (1) = 1.D-8
ATOL (2) = 1.D-14
ATOL (3) = 1.D-6
ITASK = 1
ISTATE = 1
IOPT = 0
LRW = 67
LIW = 33
MF = 21
DO 40 IOUT = 1, 12
CALL DVODE (FEX, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, &
ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF, RPAR, IPAR)
WRITE (6, 20) T, Y (1), Y (2), Y (3)
20 FORMAT (' At t =',D12.4,' y =',3D14.6)
IF (ISTATE.LT.0) GOTO 80
40 TOUT = TOUT * 10.
WRITE (6, 60) IWORK (11), IWORK (12), IWORK (13), IWORK (19), &
IWORK (20), IWORK (21), IWORK (22)
60 FORMAT(/' No. steps =',I4,' No. f-s =',I4, &
& ' No. J-s =',I4,' No. LU-s =',I4/ &
& ' No. nonlinear iterations =',I4/ &
& ' No. nonlinear convergence failures =',I4/ &
& ' No. error test failures =',I4/)
STOP
80 WRITE (6, 90) ISTATE
90 FORMAT(///' Error halt: ISTATE =',I3)
STOP
END PROGRAM EXAMPLE

SUBROUTINE FEX (NEQ, T, Y, YDOT, RPAR, IPAR)
DOUBLEPRECISION RPAR, T, Y, YDOT
DIMENSION Y (NEQ), YDOT (NEQ)
YDOT (1) = - .04D0 * Y (1) + 1.D4 * Y (2) * Y (3)
YDOT (3) = 3.D7 * Y (2) * Y (2)
YDOT (2) = - YDOT (1) - YDOT (3)
RETURN
END SUBROUTINE FEX

SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD, RPAR, IPAR)
DOUBLEPRECISION PD, RPAR, T, Y
DIMENSION Y (NEQ), PD (NRPD, NEQ)
PD (1, 1) = - .04D0
PD (1, 2) = 1.D4 * Y (3)
PD (1, 3) = 1.D4 * Y (2)
PD (2, 1) = .04D0
PD (2, 3) = - PD (1, 3)
PD (3, 2) = 6.D7 * Y (2)
PD (2, 2) = - PD (1, 2) - PD (3, 2)
RETURN
END SUBROUTINE JEX

But now, I want to put a variable "C" in the equation,
dy1/dt = C*y1 + 1.e4*y2*y3, where C=-0.04

The reason, I'm trying to do is, I'm trying to solve a chemical kinetic equation, where I've to update the coefficient of variable (y1, y2..) at every time step, for a particular node. That means I need to calculate the coefficient in main program then pass to the user defined function "F", [in this case it's FEX], for my problem."


Now I tried to pass a constant value to the program, but I failed to introduce a variable. Please help, how to introduce a variable in function.

I posted the link for DVODE f90 version.

 

[mikrom]

To pass the constant C from the main program to the user defined subroutine, I would use COMMON BLOCK.

 

[mikrom]

Here is the example, how to pass the variable C from the main program to the user defined function FOO() using COMMON BLOCK

[COLOR=#800080]program[/color] common_example
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]external[/b][/color] FOO
  [COLOR=#2e8b57][b]integer[/b][/color] :: j, c
[COLOR=#2e8b57][b]  real[/b][/color] :: x, y
  [COLOR=#2e8b57][b]common[/b][/color] [COLOR=#a52a2a][b]/[/b][/color]globals[COLOR=#a52a2a][b]/[/b][/color]c
  
  [COLOR=#a52a2a][b]do[/b][/color] j[COLOR=#a52a2a][b]=[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]5[/color]
    [COLOR=#0000ff]! changing the constant used in F in the main program loop[/color]
    c [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]10[/color] [COLOR=#a52a2a][b]*[/b][/color] j
    
    x [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color] [COLOR=#a52a2a][b]*[/b][/color] j [COLOR=#0000ff]! input  value of the subroutine SOLVE[/color]
    y [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]0[/color]       [COLOR=#0000ff]! result value of the subroutine SOLVE[/color]
    [COLOR=#0000ff]! call solver routine[/color]
    [COLOR=#008b8b]call[/color] SOLVE(FOO, x, y)

    [COLOR=#0000ff]! print iteration result[/color]
    [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#ff00ff]"('# j = ', I2, ': x = ', f5.2, '  ==> y = ', f6.2)"[/color]) j, x, y  
    [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color])
  [COLOR=#a52a2a][b]end do[/b][/color]

[COLOR=#800080]end program[/color]

[COLOR=#0000ff]! solver subroutine[/color]
[COLOR=#800080]subroutine[/color] SOLVE(F, x_in, y_out)
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]external[/b][/color] FOO
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x_in
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]out[/b][/color]) :: y_out

  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) [COLOR=#ff00ff]'* calling SOLVE('[/color], x_in, [COLOR=#ff00ff]','[/color], y_out, [COLOR=#ff00ff]')'[/color]
   
  [COLOR=#0000ff]! call user function [/color]
  [COLOR=#008b8b]call[/color] F(x_in, y_out)
  [COLOR=#0000ff]! compute solution[/color]
  y_out [COLOR=#a52a2a][b]=[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#a52a2a][b]*[/b][/color] y_out

  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) [COLOR=#ff00ff]'* computed solution: '[/color], y_out 
[COLOR=#800080]end subroutine[/color]

[COLOR=#0000ff]! subroutine to provide the user defined function f(x)[/color]
[COLOR=#800080]subroutine[/color] FOO(X_input, Y_output)
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: X_input
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]out[/b][/color]) :: Y_output   
  [COLOR=#2e8b57][b]integer[/b][/color] :: c
  [COLOR=#2e8b57][b]common[/b][/color] [COLOR=#a52a2a][b]/[/b][/color]globals[COLOR=#a52a2a][b]/[/b][/color]c
 
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) [COLOR=#ff00ff]'  ** user subroutine FOO() begin :'[/color]
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) [COLOR=#ff00ff]'  **   c = '[/color], c

  [COLOR=#0000ff]! functional value   [/color]
  y_output [COLOR=#a52a2a][b]=[/b][/color] c [COLOR=#a52a2a][b]*[/b][/color] x_input
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#ff00ff]"('   **   f(', f5.2, ') = ', f6.2)"[/color]) x_input, y_output    
  [COLOR=#a52a2a][b]write[/b][/color]([COLOR=#a52a2a][b]*[/b][/color],[COLOR=#a52a2a][b]*[/b][/color]) [COLOR=#ff00ff]'  ** user subroutine FOO() end.'[/color]
[COLOR=#800080]end subroutine[/color]

Output:

mikrom@mikrom-Lenovo-S500 ~/Work/fortran $ gfortran common_example.f95 -o common_example
mikrom@mikrom-Lenovo-S500 ~/Work/fortran $ ./common_example
 * calling SOLVE(   1.00000000     ,   0.00000000     )
   ** user subroutine FOO() begin :
   **   c =           10
   **   f( 1.00) =  10.00
   ** user subroutine FOO() end.
 * computed solution:    20.0000000    
# j =  1: x =  1.00  ==> y =  20.00

 * calling SOLVE(   2.00000000     ,   0.00000000     )
   ** user subroutine FOO() begin :
   **   c =           20
   **   f( 2.00) =  40.00
   ** user subroutine FOO() end.
 * computed solution:    80.0000000    
# j =  2: x =  2.00  ==> y =  80.00

 * calling SOLVE(   3.00000000     ,   0.00000000     )
   ** user subroutine FOO() begin :
   **   c =           30
   **   f( 3.00) =  90.00
   ** user subroutine FOO() end.
 * computed solution:    180.000000    
# j =  3: x =  3.00  ==> y = 180.00

 * calling SOLVE(   4.00000000     ,   0.00000000     )
   ** user subroutine FOO() begin :
   **   c =           40
   **   f( 4.00) = 160.00
   ** user subroutine FOO() end.
 * computed solution:    320.000000    
# j =  4: x =  4.00  ==> y = 320.00

 * calling SOLVE(   5.00000000     ,   0.00000000     )
   ** user subroutine FOO() begin :
   **   c =           50
   **   f( 5.00) = 250.00
   ** user subroutine FOO() end.
 * computed solution:    500.000000    
# j =  5: x =  5.00  ==> y = 500.00

 

[salgerman]

hhhmmm...if I understand correctly, I am thinking you have several choices.

First, if you have the source for DVODE (which I think you do), you can make changes to accommodate the additional constant in the call to DVODE; after all, it seem you have already made some changes to the expected list of arguments to FEX: RPAR and IPAR do not seem to show in the first google search for DVODE.

Second, you can simply place your FEX in a module, declare the variable C in it and use the module in your program and set the value of C there.

Third, I don't know enough to know whether this affects other things, but you could "increase" the size of your system from Y(1..3) to Y(1..4) and use the fourth "equation" as some kind of dummy equation to keep passing a constant value and use it in the spot for C.

 

[mikrom]

...
it seem you have already made some changes to the expected list of arguments to FEX: RPAR and IPAR do not seem to show in the first google search for DVODE
...

Hi salgerman,
If you look in the source which is available here:

http://www.netlib.org/ode/vode.f

you would see that OP tried only example program which is provided the source introductory comments.

 

[salgerman]

I see.

Also, somehow, I ended up posting without seeing the previous post.

I do know that I on purpose mentioned the use of a module instead of a common block...we need to stop using common blocks. Besides, the OP said the program is in F90.

 

[Priyabrata Mohapatra]

Hi mikrom,
It's a great help. Thanks. It's showing good result now. And I'm happy that I could able to pass the variable now. Will keep posting if I get any further issue.

Once again, thanks a lot.

 

[mikrom]

Hi Priyabrata Mohapatra, you are welcome.