*********************************************************************** * Benchmark Program using SHUS(Simple High Resolution Upwind Scheme) * * Eiji SHIMA ALL RIGHT RESERVED, 2003 * *********************************************************************** PARAMETER (NX=2001,NY=2001,IMAX=500,JMAX=500,NUM=10) IMPLICIT REAL*8(A-H,O-Z) REAL*8 RO(0:NX,0:NY),ROU(0:NX,0:NY),ROV(0:NX,0:NY),E(0:NX,0:NY) REAL*8 XNI(0:NX,0:NY),YNI(0:NX,0:NY),DLI(0:NX,0:NY) REAL*8 XNJ(0:NX,0:NY),YNJ(0:NX,0:NY),DLJ(0:NX,0:NY) REAL*8 FRO(0:NX,0:NY),FROU(0:NX,0:NY),FROV(0:NX,0:NY), & FE(0:NX,0:NY) REAL*4 time0(2),TIME1(2),DTIME DO 100 J=0,NY DO 100 I=0,NX RO(I,J)=1.0 ROU(I,J)=0.1 ROV(I,J)=0.0 E(I,J) =1.0 XNI(I,J)=1.0 YNI(I,J)=0.0 DLI(I,J)=1.0 XNJ(I,J)=0.0 YNJ(I,J)=1.0 DLJ(I,J)=1.0 100 CONTINUE GAMM =1.4 GAMM4=1.0/(GAMM-1.0) GAMM1=GAMM-1.0 ONE=1.0 DT =0.001 ****** TITLE ****** WRITE(*,*) WRITE(*,'(''##### BENCHMARK OF COMPRESSIBLE CFD ######'')') WRITE(*,'(''##### SHUS VER.2 #########################'')') WRITE(*,'(''ARRAY SIZE (0:NX,0:NY)=(0:'',I4, & '',0:'',I4,'')'')')NX,NY WRITE(*,'('' (IMAX,JMAX)=('',I4, & '',''I4,'')'')')IMAX,JMAX WRITE(*,'(''ITERATION='',I4)')NUM DO 200 L=1,3 * ***** START TIMING ****** * **** FOR UNIX WS ****** EE0=DTIME(TIME0) * **** FOR FUJITSU FRT COMPILER ***** CF CALL TIME(IEE0) CF CALL CLOCKM(IE0) DO 10 N=1,NUM DO 20 J=1,JMAX DO 20 I=1,IMAX * ******* EULER FLUX BY SHUS ************* * ******* NO. OF FLOATING POINT OPERATION IN THIS LOOP = 103 * ******* SET VARIABLES ****** RO0 = RO(I,J) ROU0 = ROU(I,J) ROV0 = ROV(I,J) E0 = E(I,J) RO6 = RO(I,J+1) ROU6 = ROU(I,J+1) ROV6 = ROV(I,J+1) E6 = E(I,J+1) XN = XNJ(I,J+1) YN = YNJ(I,J+1) DS = DLJ(I,J+1) * ******* AT RIGHT ****** RO0I = 1.0/RO0 U0 = ROU0*RO0I V0 = ROV0*RO0I P0 = GAMM1*(E0-0.5*RO0*(U0**2+V0**2)) VN0 = XN*U0+YN*V0 H0 = (E0+P0)*RO0I FRO0 = RO0*VN0 * ******* AT LEFT ****** RO6I = 1.0/RO6 U6 = ROU6*RO6I V6 = ROV6*RO6I P6 = GAMM1*(E6-0.5*RO6*(U6**2+V6**2)) VN6 = XN*U6+YN*V6 H6 = (E6+P6)*RO6I FRO6 = RO6*VN6 * ***** SET AVERAGE VARIABLES ROC = (RO0+RO6) VN = 0.5*(VN0+VN6) C = SQRT(GAMM*(P0+P6)/ROC) ROC = ROC*0.5 CI = 1.0/C AMC = VN*CI * ****** DIFFERENCE OF PRIMITIVE QUANTITIES ***** DU = VN6-VN0 DRO = RO6-RO0 DP = P6 -P0 * ****** FLUXES ****** ABS1 = ABS(AMC+1.0) ABS2 = ABS(AMC-1.0) FRO1 = ABS(VN)*DRO+0.5*(ABS1-ABS2)*ROC*DU & +0.5*(ABS1+ABS2-2.0*ABS(AMC))*DP*CI FROC = 0.5*(FRO0+FRO6-FRO1) FF0 = 0.5*(FROC+ABS(FROC)) FF6 = FROC-FF0 FROUC = FF0*U0+FF6*U6 FROVC = FF0*V0+FF6*V6 FEC = FF0*H0+FF6*H6 * ******* PRESSURE ****** AM0 = VN0*CI AM01 = MAX(-ONE,MIN(ONE,AM0)) BT0 = (-AM01+2.0)*(AM01+1.0)**2 AM6 = VN6*CI AM61 = MAX(-ONE,MIN(ONE,AM6)) BT6 =-(-AM61-2.0)*(AM61-1.0)**2 PC = 0.25*(P0*BT0+P6*BT6) * ****** TOTAL FLUX ****** FRO(I,J+1) = DS*FROC FROU(I,J+1) = DS*(FROUC+XN*PC) FROV(I,J+1) = DS*(FROVC+YN*PC) FE(I,J+1) = DS*FEC 20 CONTINUE DO 50 J=1,JMAX DO 50 I=1,IMAX * ******* NO. OF FLOATING POINT OPERATION IN THIS LOOP = 12 FRO1 =FRO(I,J+1)*DT FROU1=FROU(I,J+1)*DT FROV1=FROV(I,J+1)*DT FE1 =FE(I,J+1)*DT RO(I,J) =RO(I,J) -FRO1 ROU(I,J) =ROU(I,J) -FROU1 ROV(I,J) =ROV(I,J) -FROV1 E(I,J) =E(I,J) -FE1 RO(I,J+1) =RO(I,J+1) +FRO1 ROU(I,J+1)=ROU(I,J+1)+FROU1 ROV(I,J+1)=ROV(I,J+1)+FROV1 E(I,J+1) =E(I,J+1) +FE1 50 CONTINUE 10 CONTINUE ****** STOP TIMING ****** * **** FOR UNIX WS ****** EE1=DTIME(TIME1) * **** FOR FUJITSU FRT COMPILER ***** CF CALL TIME(IEE1) CF CALL CLOCKM(IE1) CF EE1 = 1.0E-3*FLOAT(IEE1-IEE0) CF TIME1(1) = 1.0E-3*FLOAT(IE1-IE0) GOSA=0 DO 30 J=1,JMAX DO 30 I=1,IMAX GOSA=GOSA+FRO(I,J)+FROU(I,J)+FROV(I,J)+FE(I,J) 30 CONTINUE WRITE(*,'(''##### TRIAL '',I4,'' #####'')') L NFLOP=IMAX*JMAX*NUM*(103+12) IF(EE1.NE.0.0) XFLOP=NFLOP/EE1*1.0E-6 IF(TIME1(1).NE.0.0) XFLOP2=NFLOP/TIME1(1)*1.0E-6 WRITE(*,6001) TIME1(1), EE1, GOSA WRITE(*,6004) XFLOP2, XFLOP 6001 FORMAT(' CPUTIME=',F18.15,3X,'ELAPSED=',F18.15,3X,'GOSA=',E13.4) 6004 FORMAT(' MFLOPS(CPUTIME)=',F10.5,7X,'MFLOPS(ELAPSTIME)=',F10.5) 200 CONTINUE END