2
$\begingroup$

Consider the following: I would like to solve a differential equation in Fortran (with a Runge-Kutta routine), e.g. x'=-a*x, where a is a constant. I then want to solve this for different a. I now would like to write the code such that it can be applied to other problems. But in general I do not know what additional input my differential equation may need (and what size and type they have). So my question is if it is possible in Fortran to pass all additional parameters to my Runge-Kutta method (at best with a single property) so that it passes it on to my differential equation so that it can be used there? With a known set of parameters this is easy but what about an unknown set?

My approach so far is this:

PROGRAM main
  USE RKModule

  IMPLICIT NONE

  INTEGER :: l

  REAL*8, PARAMETER :: num_t = 1e3
  REAL*8, PARAMETER :: min_t = 0.0d0
  REAL*8, PARAMETER :: max_t = 10.0d0

  REAL*8, DIMENSION(1), PARAMETER :: a = 1.8d0

  REAL*8, DIMENSION(num_t) :: dGrid_t, grid_t
  COMPLEX*16, DIMENSION(num_t) :: sol
  COMPLEX*16, DIMENSION(1) :: x_in, x_out

  DO l = 1, num_t
    dGrid_t(l) = (max_t - min_t) / DBLE(num_t)
    grid_t(l) = min_t + DBLE(l-0.5d0)*dGrid_t(l)
  END DO

  x_in(1) = dcmplx(1.0d0, 0.0d0)
  x_out(1) = dcmplx(0.0d0, 0.0d0)

  DO l = 1, num_t
    sol(l) = x_in(1)
    CALL rk45(x_out, DGL, x_in, grid_t(l), dGrid_t(l), a)
    x_in = x_out
  END DO

  OPEN(UNIT = 1000, FILE = 'sol.dat', ACTION = 'WRITE')
  WRITE(1000, FMT = '(e14.7e3, 1x, e14.7e3)', ADVANCE = "NO") sol
  CLOSE(1000)

  CONTAINS

  SUBROUTINE DGL(t,f,df,vararg_in)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8, DIMENSION(1) :: b

    b = 1.0d0
    IF (PRESENT(vararg_in)) THEN
      b = vararg_in
    END IF

    df = -b*f

  END SUBROUTINE DGL

END PROGRAM main

And the Runge Kutta Routine:

MODULE RKModule
  IMPLICIT NONE

  CONTAINS

  SUBROUTINE rk45(x_out, RHS, x_in, t, dt, vararg_in)
    INTERFACE
      SUBROUTINE RHS(t,f,df,a)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: a
      END SUBROUTINE
    END INTERFACE
    REAL*8, INTENT(IN) :: t, dt 
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: x_in
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: x_out
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8 :: tHlp
    COMPLEX*16, DIMENSION(SIZE(x_in)) :: k1 ,k2, k3, k4, dx, x_hlp
    REAL*8, DIMENSION(:), ALLOCATABLE :: vararg_in_copy

    IF (PRESENT(vararg_in)) THEN
      ALLOCATE(vararg_in_copy(SIZE(vararg_in)))
      vararg_in_copy = vararg_in

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

      DEALLOCATE(vararg_in_copy)
    ELSE

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

    END IF

  END SUBROUTINE rk45

END MODULE RKModule
$\endgroup$
  • $\begingroup$ What's not working about your code above? It seems like the best approach is to use an assumed shape array as you have done. $\endgroup$ – Doug Lipinski May 23 '14 at 15:51
  • $\begingroup$ I agree with @doug-lipinski, I think you have a workable solution here although I am uncertain if you need to really make the copy here. If you want to go beyond what you have done here (a more complex argument for instance) you can pass in arbitrary data through functions without them needing to know what the data looks like by leveraging the ISO C bindings. If you are interested in this solution I can write a more in depth answer. $\endgroup$ – Kyle Mandli May 23 '14 at 20:35
  • $\begingroup$ @DougLipinski The code is working. But I like to be able to pass a generic list as vararg_in. That means I would like to pass a REAL*8, a COMPLEX*16(:,:) and other stuff all bundled in vararg_in. $\endgroup$ – DaPhil May 24 '14 at 6:19
  • $\begingroup$ @KyleMandli Yes I am interested in such a solution and would be grateful if you can supply more details. $\endgroup$ – DaPhil May 24 '14 at 6:21
  • $\begingroup$ Oh. Got it. As long as the number of types you need to support is small, it's simplest to just create multiple assumed shape arrays, one for each type. If desired, you could call the routine using keyword arguments too. Otherwise, Fortran is not the right solution if you have to have this capability. You'll have to use hacks or cross-language solutions like ISO_C_BINDING or just move to another language like C/C++ that supports variadic arguments if you have to do it this way. $\endgroup$ – Doug Lipinski May 24 '14 at 13:06
2
$\begingroup$

This question better suits to StackOverflow.

You can use class(*) or type(c_ptr). Both come from Fortran 2003. The latter is widely supported, the former only since gfortran 4.8 and other recent versions of other compilers. Note type(c_ptr) is very useful even inside Fortran, without calling any C.

The first approach:

subroutine RK4(normal arguments, vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  !if you need to preserve a copy (why?)
  allocate(vararg_copy, source=varrarg_in)

  call RHS(...,vararg_copy) !or directly vararg_in
end subroutine


subroutine DGL(...,vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  select type (vararg_in)
    select type (expected_type)
     !do something
    class default
     !some exception
  end select
end subroutine

You can pass anything to class(*), possibly a derived type, if you need variables of differing types.

The other approach:

subroutine RK4(normal arguments, vararg_in)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr

  !Doing a copy of the pointer doesn't make much sense

  call RHS(...,vararg_in_ptr)
end subroutine


subroutine DGL(...,vararg_in_ptr)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr
  type(expected_type), pointer :: vararg_in(:) !for scalar delete (:)

  !for an array
  call c_f_pointer(vararg_in_ptr, vararg_in, expected_shape)
  !for a scalar just
  call c_f_pointer(vararg_in_ptr, vararg_in)

  !do something
end subroutine
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.