2017-07-14 49 views
0

我是OpenMP的新手。我想使用並行do循環來解決一系列參數值的僵硬ODE系統。我在下面給出的Fortran中使用以下代碼。但是,我不知道是否允許在並行循環內部調用僵硬的求解器(作爲子例程)?另外,我想在返回主程序之前將時間序列數據寫入子程序中帶有諸如「r_value_s__value.txt」等文件名的文件中。任何人都可以幫忙以下是代碼和錯誤。我用gfortran與標誌-fopenmp進行編譯。OpenMP參數並行掃描

PROGRAM OPENMP_PARALLEL_STIFF 

      USE omp_lib 
      IMPLICIT NONE 

      INTEGER :: I, J 
      INTEGER, PARAMETER :: RTOT=10, STOT=15 
      INTEGER :: TID 
      INTEGER, PARAMETER :: NUM_THREADS=8 
      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,TID) 

      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 for the 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=3011+II+(JJ-1)*RTOT 
      IND=TID+3011+II+(JJ-1)*RTOT 
      WRITE (*,12)SERIAL_NUMBER,TID 
    12 FORMAT ("SL. NO. ",I5," THREAD NO.",I3) 

      r=(II-1)*dr 
      s=s0+JJ*ds 

      EPS = 1.0D-9 

!   Open the output file: 

      WRITE (filename,93)r,s 
    93 FORMAT ("r_",f6.4,"_s_",f4.2,".txt") 
      OPEN (UNIT=IND,FILE=filename,STATUS='UNKNOWN',ACTION='WRITE') 

!  Parameters for the stiff ODE system 

      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 

!   S and r dependent parameters 

      e1s = e1/s; e2s = e2/(s**2); e3s = e3/s; e4s = e4/s 
      dels = del*s; rs = r*s 
      e1v = e1/v;  e2v = e2/(v**2); e3v = e3/v; e4v = e4/v 
      delv = del*v;  rv = r*v 

!   SET INITIAL PARAMETERS for INTEGRATION ROUTINES 

      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: USING MODULE USEFUL_PARAMETERS_N_FUNC 

      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) 

!   Integration: 

      DO I=1,NTOUT 

      TOUT = (I-1)*DELTAT 

      CALL DVODE_F90(F_FUNC,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS) 

!   Stop the integration in case of an error 

      IF (ISTATE<0) THEN 
      WRITE (*,*)"ISTATE ", ISTATE 
      STOP 
      END IF 

!   WRITE DATA TO FILE 

      WRITE (IND,*) TOUT,T, Y(NEQ-2) 

     END DO 


     CLOSE(UNIT=IND) 

     RETURN 
    END SUBROUTINE STIFF_DRIVER 

At line ** of file openmp_parallel_stiff.f90 (unit = 3013) Fortran runtime error: File already opened in another unit

+1

我懷疑在'stiff_driver'中,對於兩個不同的'II'和'JJ'對,'r'和's'的值是相同的(對4位和2位小數)。你可以檢查嗎? – francescalus

+0

由於您的代碼當前有效,因此一個線程不會同時調用'STIFF_DRIVER',因此您可以使用更簡單的方案來導出單元編號,例如, '3010 + tid'。使用Fortran 2008編譯器,也可以執行OPEN(NEWUNIT = ind,...).... CLOSE(UNIT = ind)'。調用'OPEN(NEWUNIT = ind,...)'將文件連接到一個未使用的單元號,並在'ind'中返回後者。如果Fortran運行時不是線程安全的,則應該將「OPEN」和「CLOSE」語句包含在關鍵部分中。 –

回答

0

的問題是,您選擇的格式:f6.4r將溢出的r>=10。然後,對於所有線程上的r>=10的所有值,輸出將爲六個星號******(取決於編譯器)。 s也是如此。

我會建議限制/檢查這些值的範圍或擴展格式以兌現更多數字。


如@francescalus所提到的,另一種可能性是命中的IIJJ其中rs是相同的組合。

只是爲了好玩 - 讓我們做數學題:

r=(II-1)*dr 
s=s0+JJ*ds 

r=s如下

(II-1)*dr = s0+JJ*ds 

II = 1 + s0/dr + JJ*ds/dr 

使用常量s0=0.450D0, dr=1.0D-4, ds=1.0D-2產量

II = 4501 + JJ*10 

因此,每當這個組合對於兩個(或更多)線程都是真實的,你就會遇到觀察到的問題。

這種情況下的簡單解決方案:將線程號添加到文件名。

+0

謝謝@francescalus和亞歷山大。 FORMAT是你正確指出的問題。 –