Pass Module Name as Function Input in Fortran

Question

I am interested in writing a function that takes as one of its inputs a module name to use. For instance, I've written a Runge Kutta 4th order integrator that is used to solve a system of ODEs. I would prefer to write the actual Runge Kutta integration function so that there is an input to the function that defines the system of ODEs.

Here's an example module containing a function defining a system of ODEs

MODULE DIFF_EQ
! Description: This module contains the function that defines the
!              system of ODEs to be solved.

CONTAINS

    FUNCTION YDOT(t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    REAL*8, INTENT(IN)               :: t
    REAL*8, DIMENSION(2), INTENT(IN) :: y
    REAL*8, DIMENSION(3), INTENT(IN) :: PP
    REAL*8, DIMENSION(2)             :: yd

    yd = [y(2), &
          -PP(1)*SIN(y(1)) + SIN(PP(2)*t) + PP(3)]

    END FUNCTION YDOT

END MODULE DIFF_EQ

And here is the Runge Kutta integration module

MODULE RUNGE_KUTTA
! Description: This module contains the function that implements the Runge Kutta 4th 
!              order integration scheme.

CONTAINS

    FUNCTION RK4(Neq, tspan, y0, dt, DEQ, PP) RESULT(t_y)
    ! Description: This function implements the Runge Kutta 4th order integration scheme.
    !
    ! Inputs: Neq   - number of equations in system of ODEs (integer, scalar)
    !         tspan - [t0, tF] where t0 is start time and tF is end time (real, array - 2)
    !         y0    - [y1, y2, ..., yNeq] @ t0 (real, array - Neq)
    !         dt    - time step size such that t1 = t0 + dt (real, scalar)
    !         PP    - passed parameters/constants (real, array - variable)
    !
    ! Outputs: t_y  - time and solution (real, matrix - n x Neq + 1)

    USE DIFF_EQ

    IMPLICIT NONE

    INTEGER, INTENT(IN)                  :: Neq
    REAL*8, DIMENSION(2), INTENT(IN)     :: tspan
    REAL*8, DIMENSION(Neq), INTENT(IN)   :: y0
    REAL*8                               :: dt
    REAL*8, EXTERNAL                     :: DEQ
    REAL*8, DIMENSION(:), INTENT(IN)     :: PP
    INTEGER                              :: n, i
    REAL*8, DIMENSION(:,:), ALLOCATABLE  :: t_y
    REAL*8, DIMENSION(Neq)               :: k1, k2, k3, k4

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)
    ALLOCATE(t_y(n, Neq+1))

    IF (MOD(tspan(2) - tspan(1), dt) .LT. 0.000000000000001D0) THEN
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n)]
    ELSE
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n-1), tspan(2)]
    ENDIF

    t_y(1, 2:Neq+1) = y0

    PRINT *, t_y(1, 1:Neq+1)

    DO i = 2, n
        dt = t_y(i, 1) - t_y(i-1, 1)
        k1 = DEQ(t_y(i-1, 1),            t_y(i-1, 2:Neq+1),               PP)
        k2 = DEQ(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k1, PP)
        k3 = DEQ(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k2, PP)
        k4 = DEQ(t_y(i-1, 1) + dt,       t_y(i-1, 2:Neq+1) + dt*k3,       PP)
        t_y(i, 2:Neq+1) = t_y(i-1, 2:Neq+1) + 0.16666666666666667D0*dt*(k1 + 2.0D0*k2 + 2.0D0*k3 + k4)

        PRINT *, t_y(i, 1:Neq+1)
    ENDDO

    END FUNCTION RK4

END MODULE RUNGE_KUTTA

And the main program that takes the inputs and calls the functions

PROGRAM RK4_TEST

    USE DIFF_EQ
    USE RUNGE_KUTTA

    IMPLICIT NONE

    INTEGER                             :: Neq
    REAL*8, DIMENSION(2)                :: tspan
    REAL*8, DIMENSION(2)                :: y0
    REAL*8                              :: dt
    REAL*8, DIMENSION(3)                :: PP
    INTEGER                             :: n
    REAL*8, DIMENSION(:,:), ALLOCATABLE :: t_y

    Neq   = 2
    tspan = [0.0D0, 1.0D0]
    y0    = [1.0D0, 0.0D0]
    dt    = 0.1D0
    PP    = [1.0D0, 5.0D0, 0.0D0]

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)

    ALLOCATE(t_y(n, Neq+1))

    t_y = RK4(Neq, tspan, y0, dt, YDOT, PP)

END PROGRAM RK4_TEST

Using gfortran, this is compiled by

gfortran DIFF_EQ.f90 RUNGE_KUTTA.f90 RK4_TEST.f90 -o MAIN.exe

It compiles without error but I receive a Segmentation fault (core dumped) message when I try to run it.

Thanks for any help provided.

score 4 · Answer 1 · edited May 23 '17 at 12:28

You can pass a function to your solver. That should be sufficient for making your solver general, able to handle various functions. Make your solver take a function as an argument, then call it with actual argument of the particular function that you wish to solve. Or you can call it with a function pointer, with each call having the pointer pointing to the function that you wish to solve. In your example, you don't pass modules similar to DIFF_EQ, you pass functions similar to YDOT. Here is an example with function pointers: Function pointer arrays in Fortran

EDIT: OK, here is a code sample. I have shown how a second function can be passed to the solver. By using an interface statement to declare the function in the solver, the module doesn't have to be used in the solver. The functions could easily be in different modules and the source code of the solver would not need to be modified, just the source code of the calling program.

MODULE DIFF_EQ

   use ISO_FORTRAN_ENV

! Description: This module contains the function that defines the
!              system of ODEs to be solved.

CONTAINS

    FUNCTION YDOT(t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    real (real64), INTENT(IN) :: t
    real (real64), INTENT(IN) :: y(2)
    real (real64), INTENT(IN) :: PP(3)
    real (real64)             :: yd(2)

    yd = [y(2), -PP(1)*SIN(y(1)) + SIN(PP(2)*t) + PP(3)]

    END FUNCTION YDOT

    FUNCTION YDOT2 (t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    real (real64), INTENT(IN) :: t
    real (real64), INTENT(IN) :: y(2)
    real (real64), INTENT(IN) :: PP(3)
    real (real64)             :: yd(2)

    yd = [y(2), -PP(1)*cos(y(1)) + cos(PP(2)*t) + PP(3)]

    END FUNCTION YDOT2

END MODULE DIFF_EQ

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

MODULE RUNGE_KUTTA

   use ISO_FORTRAN_ENV
! Description: This module contains the function that implements the Runge Kutta 4th
!              order integration scheme.

CONTAINS

    FUNCTION RK4(DEQ_func, Neq, tspan, y0, dt, PP) RESULT(t_y)
    ! Description: This function implements the Runge Kutta 4th order integration scheme.
    !
    ! Inputs: Neq   - number of equations in system of ODEs (integer, scalar)
    !         tspan - [t0, tF] where t0 is start time and tF is end time (real, array - 2)
    !         y0    - [y1, y2, ..., yNeq] @ t0 (real, array - Neq)
    !         dt    - time step size such that t1 = t0 + dt (real, scalar)
    !         PP    - passed parameters/constants (real, array - variable)
    !
    ! Outputs: t_y  - time and solution (real, matrix - n x Neq + 1)


    IMPLICIT NONE

   interface

      function  DEQ_func ( t, y, pp )  result (yd)

         use ISO_FORTRAN_ENV
         real (real64), INTENT(IN) :: t
         real (real64), dimension (2), INTENT(IN) :: y
         real (real64), dimension (3), INTENT(IN) :: PP
         real (real64), dimension (2)             :: yd

      end function DEQ_func

   end interface

    INTEGER, INTENT(IN)                  :: Neq
    real (real64), DIMENSION(2), INTENT(IN)     :: tspan
    real (real64), INTENT(IN)                   :: y0(Neq)
    real (real64)                               :: dt
    real (real64), INTENT(IN)                   :: PP(:)
    INTEGER                              :: n, i
    real (real64), DIMENSION(:,:), ALLOCATABLE  :: t_y
    real (real64)                               :: k1(Neq), k2(Neq), k3(Neq), k4(Neq)

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)
    ALLOCATE(t_y(n, Neq+1))

    IF (MOD(tspan(2) - tspan(1), dt) .LT. 0.000000000000001D0) THEN
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n)]
    ELSE
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n-1), tspan(2)]
    ENDIF

    t_y(1, 2:Neq+1) = y0

    PRINT *, t_y(1, 1:Neq+1)
    DO i = 2, n
        dt = t_y(i, 1) - t_y(i-1, 1)
        k1 = DEQ_func(t_y(i-1, 1),              t_y(i-1, 2:Neq+1),                 PP)
        k2 = DEQ_func(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k1, PP)
        k3 = DEQ_func(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k2, PP)
        k4 = DEQ_func(t_y(i-1, 1) + dt,         t_y(i-1, 2:Neq+1) + dt*k3,         PP)
        t_y(i, 2:Neq+1) = t_y(i-1, 2:Neq+1) + 0.16666666666666667D0*dt*(k1 + 2.0D0*k2 + 2.0D0*k3 + k4)

        PRINT *, t_y(i, 1:Neq+1)
    ENDDO

    END FUNCTION RK4

END MODULE RUNGE_KUTTA

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

PROGRAM RK4_TEST

    USE RUNGE_KUTTA
    use DIFF_EQ

    IMPLICIT NONE

    INTEGER                             :: Neq
    real (real64), DIMENSION(2)                :: tspan
    real (real64), DIMENSION(2)                :: y0
    real (real64)                              :: dt
    real (real64), DIMENSION(3)                :: PP
    INTEGER                             :: n
    real (real64), DIMENSION(:,:), ALLOCATABLE :: t_y

    Neq   = 2
    tspan = [0.0D0, 1.0D0]
    y0    = [1.0D0, 0.0D0]
    dt     = 0.1D0
    PP    = [1.0D0, 5.0D0, 0.0D0]

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)

    ALLOCATE(t_y(n, Neq+1))

    write (*, '( // "calling to solve ydot:")' )
    t_y = RK4( ydot, Neq, tspan, y0, dt, PP)

    write (*, '( // "calling to solve ydot2:")' )
    t_y = RK4( ydot2, Neq, tspan, y0, dt, PP)

END PROGRAM RK4_TEST

So would you suggest that I contain all of my functions defining systems of ODEs such as `YDOT` within the `DIFF_EQ` module and then pass the function of ODEs I want to solve at that time? For instance, `DIFF_EQ` would contain `YDOT_1`, `YDOT_2`, ...., `YDOT_N`. That way, my Runge Kutta solver would always `USE DIFF_EQ` but the function I pass to the solver would be different. — Stephen, Jul 04 '13 at 06:24
Yes. One or more modules containing the functions. Pass whichever function you wish to solve to the Runge Kutta solver, either directly or via a function pointer. The Runge Kutta solver can either `use` the module, or it can have an interface statement declaring the function is solves. The second way it is more general, if later you wish to solve functions in a different module. — M. S. B., Jul 04 '13 at 13:35
thanks so much for the code sample. I had read about the interface block but didn't quite understand how it fit in for my needs. Your code sample should clear a lot of things up for me. — Stephen, Jul 06 '13 at 06:47

Pass Module Name as Function Input in Fortran

1 Answers1