I am new to OpenMP. I want to solve a stiff ODE system for a range of parameter values using parallel do loops. I use the following Pseudocode in FORTRAN given below. However, I do not know whether calling a stiff solver(as a subroutine) inside a parallel do loop is allowed or not? Also, I want to write the time series data into files with filenames such as "par_value.txt" in the subroutine before the return to the main program. Can anyone help.
PROGRAM OPENMP_PARALLEL_STIFF
!
USE omp_lib
!
IMPLICIT NONE
!
INTEGER :: I, J
INTEGER, PARAMETER :: RTOT=10, STOT=15
INTEGER :: TID
INTEGER, PARAMETER :: NUM_THREADS=3
DOUBLE PRECISION :: T_INITIAL, T_FINAL
CALL OMP_SET_NUM_THREADS(NUM_THREADS)
CALL CPU_TIME(T_INITIAL)
!
PRINT*, "TIME INITIAL ",T_INITIAL
!
!$OMP PARALLEL DO PRIVATE(I,J)
DO I=1,RTOT
DO J=1,STOT
TID=OMP_GET_THREAD_NUM()
CALL STIFF_DRIVER(TID,I,J,RTOT,STOT)
END DO
END DO
!$OMP END PARALLEL DO
!
CALL CPU_TIME(T_FINAL)
!
PRINT*, "TIME FINAL ",T_FINAL
!
PRINT*, "TIME ELAPSED ",(T_FINAL-T_INITIAL)/NUM_THREADS
END PROGRAM OPENMP_PARALLEL_STIFF
SUBROUTINE STIFF_DRIVER(TID,II,JJ,RTOT,STOT)
!
USE USEFUL_PARAMETERS_N_FUNC
USE DVODE_F90_M
! -----------------------------------------------------------------
! Type declarations:
! -----------------------------------------------------------------
IMPLICIT NONE
! -----------------------------------------------------------------
! Number of odes and number of event functions for this problem:
! -----------------------------------------------------------------
INTEGER :: SERIAL_NUMBER, TID
INTEGER :: II, JJ, RTOT, STOT, IND
INTEGER :: J, NTOUT
INTEGER :: ITASK, ISTATE, ISTATS, I
! -----------------------------------------------------------------
! parameters : declaration
! -----------------------------------------------------------------
DOUBLE PRECISION, PARAMETER :: s0=0.450D0, dr=1.0D-4, ds=1.0D-2
! -----------------------------------------------------------------
DOUBLE PRECISION, DIMENSION(NEQ) :: Y, YOUT
DOUBLE PRECISION :: ATOL, RTOL, RSTATS, T, TOUT, EPS, TFINAL, DELTAT
DIMENSION :: RSTATS(22), ISTATS(31)
DOUBLE PRECISION :: bb, cc, ba, ba1, eta
CHARACTER(len=45) :: filename
! ---------------------------------------------------------------
!
TYPE (VODE_OPTS) :: OPTIONS
!
! -----------------------------------------------------------------
SERIAL_NUMBER=11+II+(JJ-1)*RTOT
! IND=TID+11+II+(JJ-1)*RTOT
WRITE (*,12)SERIAL_NUMBER,TID
12 FORMAT ("SL. NO. ",I3," THREAD NO.",I3)
! -----------------------------------------------------------------
r=(II-1)*dr
s=s0+JJ*ds
! -----------------------------------------------------------------
EPS = 1.0D-9
! ------------------------------------------------------------------
! Open the output file and the plot files:
! -----------------------------------------------------------------
WRITE (filename,93)r,s
93 FORMAT ("r_",f6.4,"_s_",f4.2,".txt")
OPEN (UNIT=NEWUNIT(IND),FILE=filename,STATUS='UNKNOWN')
! -----------------------------------------------------------------
EPS = 1D-9
! -----------------------------------------------------------------
! Parameters for the stiff ODE system
! -----------------------------------------------------------------
noise = 0.0010D0
pi = 4.0D0*atan(1.0D0)
q0 = 0.60D0; v = 3.0D0
Va = 20.0D-4; Vs = 1.0D-1
e1 = 1.0D-1; e2 = 1.10D-5; e3 = 2.3D-3; e4=3.0D-4
del = 1.7D-4; mu = 5.9D-4
al = 1.70D-4; be = 8.9D-4; ga = 2.5D-1
! ---------------------------------------------------------------
! Surrounding
! ---------------------------------------------------------------
e1s = e1/s; e2s = e2/(s**2); e3s = e3/s; e4s = e4/s
dels = del*s; rs = r*s
! --------------------------------------------------------------
! discrete subsystems
! --------------------------------------------------------------
e1v = e1/v; e2v = e2/(v**2); e3v = e3/v; e4v = e4/v
delv = del*v; rv = r*v
! ---------------------------------------------------------------
! Assign values for q
! ---------------------------------------------------------------
CALL RANDOM_SEED
DO I = 1, n
CALL RANDOM_NUMBER(HARVEST=bb)
IF(bb <= 0.0D0)THEN
bb = 0.00010D0
END IF
CALL RANDOM_NUMBER(HARVEST=cc)
ba = -2.0D0*DLOG(bb)
ba1 = DCOS(2.0D0*pi*cc)
eta = noise*DSQRT(ba)*ba1
! ---------------------------------------------------------------
q(I) = q0 + eta
! ---------------------------------------------------------------
END DO
! -----------------------------------------------------------------
! Set the integration parameters:
! -----------------------------------------------------------------
T = 0.0D0
TFINAL = 200.0D0
DELTAT = 0.10D0
NTOUT = INT(TFINAL/DELTAT)
RTOL = EPS
ATOL = EPS
ITASK = 1
ISTATE = 1
! -----------------------------------------------------------------
! Set the initial conditions:
! -----------------------------------------------------------------
CALL Y_INITIAL(NEQ,Y)
! -----------------------------------------------------------------
! Set the VODE_F90 options:
! -----------------------------------------------------------------
OPTIONS = SET_OPTS(DENSE_J=.TRUE.,USER_SUPPLIED_JACOBIAN=.FALSE., &
RELERR=RTOL,ABSERR=ATOL,MXSTEP=100000)
! -----------------------------------------------------------------
! Perform the integration:
! -----------------------------------------------------------------
DO I=1,NTOUT
! -----------------------------------------------------------------
TOUT = (I-1)*DELTAT
! -----------------------------------------------------------------
CALL DVODE_F90(F_FUNC,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
! -----------------------------------------------------------------
! Stop the integration if an error occurred:
! -----------------------------------------------------------------
IF (ISTATE<0) THEN
WRITE (*,*)"ISTATE ", ISTATE
STOP
END IF
! -----------------------------------------------------------------
! WRITE DATA TO FILE
! -----------------------------------------------------------------
WRITE (IND,*) TOUT,T, x_mean, Y(NEQ-2)
END DO
! -----------------------------------------------------------------
CLOSE(UNIT=IND)
END SUBROUTINE STIFF_DRIVER
I also used the following function in the modules USEFUL_PARAMETERS_N_FUNC to obtain unit number for opening files in the called subroutine to write the time series data.
INTEGER FUNCTION NEWUNIT(UNIT)
INTEGER, INTENT(OUT), OPTIONAL :: UNIT
! local
INTEGER, PARAMETER :: LUN_MIN=3011, LUN_MAX=4000
LOGICAL :: open_ed
INTEGER :: l
! begin
NEWUNIT=-1
DO l=LUN_MIN,LUN_MAX
INQUIRE(UNIT=l,OPENED=open_ed)
IF(.NOT.open_ed) THEN
NEWUNIT=l
EXIT
END IF
END DO
IF (PRESENT(UNIT)) UNIT=NEWUNIT
END FUNCTION NEWUNIT
However, I get on compiling and running the code
Fortran runtime error: File already opened in another unit:
A storage allocation error occurred.
The error occurred at location I1.
In the above message, I1 = 280
LEVEL = 2 in XERRDV. Stopping.
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
