INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Jobs

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

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

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

(OP)
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.

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

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

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

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

CODE

program common_example
  implicit none
  external FOO
  integer :: j, c
  real :: x, y
  common /globals/c
  
  do j=1, 5
    ! changing the constant used in F in the main program loop
    c = 10 * j
    
    x = 1.0 * j ! input  value of the subroutine SOLVE
    y = 0       ! result value of the subroutine SOLVE
    ! call solver routine
    call SOLVE(FOO, x, y)

    ! print iteration result
    write(*,"('# j = ', I2, ': x = ', f5.2, '  ==> y = ', f6.2)") j, x, y  
    write(*,*)
  end do

end program

! solver subroutine
subroutine SOLVE(F, x_in, y_out)
  implicit none
  external FOO
  real, intent(in) :: x_in
  real, intent(out) :: y_out

  write(*,*) '* calling SOLVE(', x_in, ',', y_out, ')'
   
  ! call user function 
  call F(x_in, y_out)
  ! compute solution
  y_out = 2 * y_out

  write(*,*) '* computed solution: ', y_out 
end subroutine

! subroutine to provide the user defined function f(x)
subroutine FOO(X_input, Y_output)
  implicit none
  real, intent(in) :: X_input
  real, intent(out) :: Y_output   
  integer :: c
  common /globals/c
 
  write(*,*) '  ** user subroutine FOO() begin :'
  write(*,*) '  **   c = ', c

  ! functional value   
  y_output = c * x_input
  write(*,"('   **   f(', f5.2, ') = ', f6.2)") x_input, y_output    
  write(*,*) '  ** user subroutine FOO() end.'
end subroutine 

Output:

CODE

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 

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

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.

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

Quote (salgerman)


...
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.

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

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.

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

(OP)
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.

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

Hi Priyabrata Mohapatra, you are welcome.

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Resources

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close