CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Moving Mesh method for one-- C dimensional Euler equations C C Local LF fluxes C Copyright(c) by Huazhong Tang CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PARAMETER(NNX=15000) IMPLICIT REAL*8 (A-H,O-Z) COMMON/COM1/UO(3,NNX),UN(3,NNX),XO(NNX),XN(NNX) COMMON/COM2/NPX,NXI COMMON/EULER/GAMA DATA NPX,NXI/503,503/ DATA XLEFT,XRIGHT,XPO,GAMA/0.0,1.0,0.3,1.4/ C C Give a parition of the domain and initial value C DX=(XRIGHT-XLEFT)/(NPX-3) DO J=1,NPX+4 XO(J)=XLEFT+(J-3)*DX XN(J)=XO(J) ENDDO DO J=1,NPX+2 XXPO=0.5*(XO(J)+XO(J+1)) IF(XXPO.LE.XPO) THEN UO(1,J)=1.0 UO(2,J)=0.0 UO(3,J)=2.5 ELSE UO(1,J)=0.125 UO(2,J)=0.0 UO(3,J)=0.25 ENDIF IF(XXPO.LT.XPO) THEN UO(1,J)=1000.0 UO(2,J)=0.0 UO(3,J)=1000.0/(GAMA-1.0) ELSE UO(1,J)=1.0 UO(2,J)=0.0 UO(3,J)=1.0/(GAMA-1.0) ENDIF ENDDO time=0.0 C C Moving mesh according to initial value C WRITE(*,*)'IF IADAP.LT.0 THEN DO NOT MOVE MESH!!' READ(*,*) IADAP CALL AMMM(time,IADAP) C C PDEs solver C CALL PDES(IADAP) C C Output the final grid distribution C OPEN(11,FILE='gridp.dat',status='unknown') DO J=1,NPX+2 WRITE(11,*) XO(J), 0.5 ENDDO CLOSE(11) OPEN(16,FILE='ene.dat',status='unknown') OPEN(15,FILE='pre.dat',status='unknown') OPEN(12,FILE='vel.dat',status='unknown') OPEN(11,FILE='rho.dat',status='unknown') DO J=3,NPX XPO=0.5*(XO(J)+XO(J+1)) PRE=(GAMA-1.0)*(UO(3,J)-0.5*UO(2,J)**2/UO(1,J) ) WRITE(11,*) XPO,UO(1,J) WRITE(12,*) XPO,UO(2,J)/UO(1,J) WRITE(15,*) XPO,PRE WRITE(16,*) XPO, PRE/(GAMA-1.)/UO(1,J) ENDDO CLOSE(11) STOP END SUBROUTINE AMMM(time,IADAP) PARAMETER(NNX=15000) IMPLICIT REAL*8 (A-H,O-Z) COMMON/COM1/UO(3,NNX),UN(3,NNX),XO(NNX),XN(NNX) COMMON/COM2/NPX,NXI COMMON/EULER/GAMA dimension wo(nnx),wn(nnx),DUX(3,NNX) & ,DUU(NNX),DUS(NNX),DUV(NNX) C DATA ALPHA,BETA,THETA/20.0,100.0,0.9/ DATA ALPHA,BETA,THETA/1000.0,1000.0,100.0/ C FMON(X,DX,ALPHA,BETA)=SQRT(1.0+ALPHA*X*X+BETA*DX*DX ) FMON(X,Y,Z,ALPHA,BETA,THETA) & =SQRT(1.0+ALPHA*X*X+BETA*Y*Y+THETA*Z*Z ) c FMON(X,DX,ALPHA,BETA,THETA)=THETA*EXP(ALPHA*ABS(X)) c & +EXP(BETA*ABS(DX))*(1.-THETA) IF(IADAP.LT.0) RETURN !DO NOT MOVE MESH NIT=20 IF(time.le.0.0) NIT=100 OMEGA=0.0 DO J=1,NXI+2 XN(J)=XO(J) DO L=1,3 UN(L,J)=UO(L,J) ENDDO ENDDO NIA=3 NIB=NPX-1 NTA=4 NTB=NXI-1 C C Solve the mesh equation C IT=0 900 IT=IT+1 C C DEFINE MONITOR FUNCTION C SUU=0.0 SUS=0.0 DO J=NIA,NIB RHOP1=UO(1,J+1) UXP1 =UO(2,J+1)/RHOP1 PREP1=(GAMA-1.)*(UO(3,J+1)-0.5*RHOP1*UXP1*UXP1) RHOM1=UO(1,J-1) UXM1 =UO(2,J-1)/RHOM1 PREM1=(GAMA-1.)*(UO(3,J-1)-0.5*RHOM1*UXM1*UXM1) DUU(J)=(UXP1-UXM1)/2.0 DUV(J)=UXP1+UXM1-2.*UO(2,J)/UO(1,J) DUS(J)=((PREP1/RHOP1**GAMA)-(PREM1/RHOM1**GAMA))/2.0 DUV(J)=(PREP1-PREM1)/2.0 SUU=MAX(SUU,DUU(J)) SUS=MAX(SUS,DUS(J)) ENDDO DO J=NIA,NIB UU=DUU(J)!/(SUU+1.0E-10) DU=DUS(J)!/(SUS+1.0E-10) DV=DUV(J) C WO(J)=FMON(UU,DU,ALPHA,BETA) WO(J)=FMON(UU,DU,DV,ALPHA,BETA,THETA) ENDDO DO 100 IJK=1,4 DO J=NIA,NIB WN(J)=WO(J) IF(J.GT.NIA.and.J.LT.NIB) THEN WN(J)=(WO(J+1)+2.*WO(J)+WO(J-1))/4. ENDIF ENDDO DO J=NIA,NIB WO(J)=WN(J) ENDDO 100 CONTINUE DO J=NTA,NTB XN(J)=(WO(J)*XO(J+1)+WO(J-1)*XN(J-1))/(WO(J)+WO(J-1)) ENDDO c DO J=NTA,NTB c XN(J)=(1.0-OMEGA)*XN(J)+OMEGA*XO(J) c ENDDO c Wang Han c MI=10 c MN=(NPX-3)/MI c DO J=NTA-1+MN,NTB,MN c XN(J)=(WO(J)*XO(J+1)+WO(J-1)*XN(J-1))/(WO(J)+WO(J-1)) c ENDDO c XN(NTA-1)=XO(NTA-1) c XN(NTB+1)=XO(NTB+1) c c DO J=NTA-1+MN,NTB+1,MN c DX=(XN(J)-XN(J-MN))/MN c DO K=J-MN,J-1 c XN(K)=XN(J-MN)+DX*(K-J+MN) c ENDDO c ENDDO SUM=0.0 DO J=NTA,NTB SUM=SUM+ABS(XN(J)-XO(J)) ENDDO SUM=SUM/(NTB-NTA) C GOTO 2000 c c non-conservative interpolation c if( mod(IT,1000).eq.0) THEN DO L=1,3 DO J=NIA-1,NIB DUL=UO(L,J)-UO(L,J-1) DUR=UO(L,J+1)-UO(L,J) DUX(L,J)=(SIGN(1.0,DUL)+SIGN(1.0,DUR)) & *ABS(DUL*DUR)/(ABS(DUL)+ABS(DUR)+1.0E-10) ENDDO DUX(L,NIA-2)=0.0 DUX(L,NIB+1)=0.0 ENDDO DO J=NIA-1,NIB WN(J)=1.0/(XO(J+1)-XO(J)) ENDDO DO J=NIA-1,NIB SPEED=0.5*(XO(J)+XO(J+1)-XN(J)-XN(J+1))*WN(J) DO L=1,3 IF(SPEED.GE.0.0) THEN UN(L,J)=UO(L,J)-SPEED*(UO(L,J )+0.5*DUX(L,J ) & -UO(L,J-1)-0.5*DUX(L,J-1)) ELSE UN(L,J)=UO(L,J)-SPEED*(UO(L,J+1)-0.5*DUX(L,J+1) & -UO(L,J )+0.5*DUX(L,J )) ENDIF ENDDO ENDDO C goto 3000 else c c Conservative interpolation c 2000 DO L=1,3 DO J=NIA, NIB DUL=UO(L,J)-UO(L,J-1) DUR=UO(L,J+1)-UO(L,J) DUX(L,J)=0.5*(SIGN(1.D0,DUL)+SIGN(1.D0,DUR)) & *2.0*ABS(DUL*DUR)/(ABS(DUL)+ABS(DUR)+1.0E-10) ENDDO DUX(L,NIA-1)=0.0 DUX(L,NIB+1)=0.0 ENDDO DO J=NIA,NIB XIXO=XO(J+1)-XO(J) XIXN=XN(J+1)-XN(J) SP1=(XO(J+1)-XN(J+1))/XIXN SP2=(XO(J)-XN(J))/XIXN DO L=1,3 UJP1=(UO(L,J+1)-0.5*DUX(L,J+1)) UJP0=UO(L,J )-0.5*DUX(L,J) UJM1=UO(L,J-1)+0.5*DUX(L,J-1) UJM0=UO(L,J )+0.5*DUX(L,J) UN(L,J)=UO(L,J)*XIXO/XIXN-.5*( sp1*(UJP1+UJM0)-SP2*(UJP0+UJM1)) & +.5*(ABS(SP1)*(UJP1-UJM0)-ABS(SP2)*(UJP0-UJM1)) ENDDO ENDDO endif 3000 DO J=NIA,NIB DO L=1,3 UO(L,J)=UN(L,J) ENDDO ENDDO DO J=NTA,NTB XO(J)=XN(J) ENDDO c write(*,*) '===========IT=', IT IF(SUM.GE.1.0E-9.AND.IT.LE.NIT) GOTO 900 END SUBROUTINE PDES(IADAP) PARAMETER(NNX=15000) IMPLICIT REAL*8 (A-H,O-Z) COMMON/COM1/UO(3,NNX),UN(3,NNX),XO(NNX),XN(NNX) COMMON/COM2/NPX,NXI COMMON/EULER/GAMA DIMENSION DUX(3,NNX),UUL(3),UUR(3),FLUX(3,NNX) DIMENSION FXL(3),FXR(3),U01(3,NNX) DIMENSION UP(3,NNX) DATA CFL, TSTOP/0.8,0.15/ NIA=3 NIB=NPX-1 DO J=NIA-2,NIB+2 DO L=1,3 UN(L,J)=UO(L,J) U01(L,J)=UO(L,J) ENDDO ENDDO IT=0 TIME=0.0 100 IF(TIME.GE.TSTOP) GOTO 888 SUM=0.0 SUX=1000.0 DO J=NIA,NIB RHO=UO(1,J) UX=UO(2,J)/RHO PRE=(GAMA-1.)*(UO(3,J)-0.5*RHO*UX*UX) SUM=MAX(SUM,ABS(UX)+SQRT(GAMA*PRE/RHO) ) SUX=MIN(SUX,XO(J+1)-XO(J) ) ENDDO DT=SUX*CFL/SUM IF(TIME+DT.GT.TSTOP) DT=TSTOP-TIME C C First step of third-order RK method C DO J=NIA-2,NIB+2 UP(1,J)=UO(1,J) UP(2,J)=UO(2,J)/UP(1,J) UP(3,J)=(GAMA-1.)*(UO(3,J)-0.5*UP(1,J)*UP(2,J)*UP(2,J)) ENDDO DO L=1, 3 DO J=NIA,NIB DUL=UP(L,J)-UP(L,J-1) DUR=UP(L,J+1)-UP(L,J) DUX(L,J)=Slope_Limiter(DUL,DUR) ENDDO DUX(L,NIA-1)=0.0 DUX(L,NIB+1)=0.0 ENDDO DO J=NIA-1,NIB DO L=1,3 UUL(L)=UP(L,J)+0.5*DUX(L,J) ENDDO c RHOL=UUL(1) c UXL=UUL(2)/UUL(1) c PREL=(GAMA-1.)*(UUL(3)-0.5*RHOL*UXL*UXL) RHOL=UUL(1) UXL=UUL(2) PREL=UUL(3) CL=SQRT(GAMA*PREL/RHOL) EIGVL=abs(UXL)+CL UUL(1)=RHOL UUL(2)=RHOL*UXL UUL(3)=PREL/(GAMA-1.)+0.5*RHOL*UXL*UXL FXL(1)=UUL(1)*UXL FXL(2)=UUL(2)*UXL+PREL FXL(3)=UUL(3)*UXL+UXL*PREL DO L=1,3 UUR(L)=UP(L,J+1)-0.5*DUX(L,J+1) ENDDO c RHOR=UUR(1) c UXR=UUR(2)/UUR(1) c PRER=(GAMA-1.)*(UUR(3)-0.5*RHOR*UXR*UXR) RHOR=UUR(1) UXR=UUR(2) PRER=UUR(3) CR=SQRT(GAMA*PRER/RHOR) EIGVR=abs(UXR)+CR UUR(1)=RHOR UUR(2)=RHOR*UXR UUR(3)=PRER/(GAMA-1.)+0.5*RHOR*UXR*UXR FXR(1)=UUR(1)*UXR FXR(2)=UUR(2)*UXR+PRER FXR(3)=UUR(3)*UXR+UXR*PRER SUM=MAX(EIGVR,EIGVL) do L=1,3 FLUX(L,J)=0.5*(FXR(L)+FXL(L)-SUM*(UUR(L)-UUL(L)) ) enddo ENDDO DO J=NIA,NIB RR=DT/(XO(J+1)-XO(J)) DO L=1,3 U01(L,J)=UO(L,J)-RR*(FLUX(L,J)-FLUX(L,J-1)) ENDDO ENDDO C C Second step of 3rd order RK method C DO J=NIA-2,NIB+2 UP(1,J)=U01(1,J) UP(2,J)=U01(2,J)/UP(1,J) UP(3,J)=(GAMA-1.)*(U01(3,J)-0.5*UP(1,J)*UP(2,J)*UP(2,J)) ENDDO DO L=1, 3 DO J=NIA,NIB DUL=UP(L,J)-UP(L,J-1) DUR=UP(L,J+1)-UP(L,J) DUX(L,J)=Slope_Limiter(DUL,DUR) ENDDO DUX(L,NIA-1)=0.0 DUX(L,NIB+1)=0.0 ENDDO DO J=NIA-1,NIB DO L=1,3 UUL(L)=UP(L,J)+0.5*DUX(L,J) ENDDO c RHOL=UUL(1) c UXL=UUL(2)/UUL(1) c PREL=(GAMA-1.)*(UUL(3)-0.5*RHOL*UXL*UXL) RHOL=UUL(1) UXL=UUL(2) PREL=UUL(3) CL=SQRT(GAMA*PREL/RHOL) EIGVL=abs(UXL)+CL UUL(1)=RHOL UUL(2)=RHOL*UXL UUL(3)=PREL/(GAMA-1.)+0.5*RHOL*UXL*UXL FXL(1)=UUL(1)*UXL FXL(2)=UUL(2)*UXL+PREL FXL(3)=UUL(3)*UXL+UXL*PREL DO L=1,3 UUR(L)=UP(L,J+1)-0.5*DUX(L,J+1) ENDDO c RHOR=UUR(1) c UXR=UUR(2)/UUR(1) c PRER=(GAMA-1.)*(UUR(3)-0.5*RHOR*UXR*UXR) RHOR=UUR(1) UXR=UUR(2) PRER=UUR(3) CR=SQRT(GAMA*PRER/RHOR) EIGVR=abs(UXR)+CR UUR(1)=RHOR UUR(2)=RHOR*UXR UUR(3)=PRER/(GAMA-1.)+0.5*RHOR*UXR*UXR FXR(1)=UUR(1)*UXR FXR(2)=UUR(2)*UXR+PRER FXR(3)=UUR(3)*UXR+UXR*PRER SUM=MAX(EIGVR,EIGVL) do L=1,3 FLUX(L,J)=0.5*(FXR(L)+FXL(L)-SUM*(UUR(L)-UUL(L)) ) enddo ENDDO DO J=NIA,NIB RR=DT/(XO(J+1)-XO(J)) DO L=1,3 UN(L,J)=0.75*UO(L,J)+0.25*(U01(L,J)-RR*(FLUX(L,J)-FLUX(L,J-1))) ENDDO ENDDO C C Third step of 3rd order RK method C DO J=NIA-2,NIB+2 UP(1,J)=UN(1,J) UP(2,J)=UN(2,J)/UP(1,J) UP(3,J)=(GAMA-1.)*(UN(3,J)-0.5*UP(1,J)*UP(2,J)*UP(2,J)) ENDDO DO L=1, 3 DO J=NIA,NIB DUL=UP(L,J)-UP(L,J-1) DUR=UP(L,J+1)-UP(L,J) DUX(L,J)=Slope_Limiter(DUL,DUR) ENDDO DUX(L,NIA-1)=0.0 DUX(L,NIB+1)=0.0 ENDDO DO J=NIA-1,NIB DO L=1,3 UUL(L)=UP(L,J)+0.5*DUX(L,J) ENDDO c RHOL=UUL(1) c UXL=UUL(2)/UUL(1) c PREL=(GAMA-1.)*(UUL(3)-0.5*RHOL*UXL*UXL) RHOL=UUL(1) UXL=UUL(2) PREL=UUL(3) CL=SQRT(GAMA*PREL/RHOL) EIGVL=abs(UXL)+CL UUL(1)=RHOL UUL(2)=RHOL*UXL UUL(3)=PREL/(GAMA-1.)+0.5*RHOL*UXL*UXL FXL(1)=UUL(1)*UXL FXL(2)=UUL(2)*UXL+PREL FXL(3)=UUL(3)*UXL+UXL*PREL DO L=1,3 UUR(L)=UP(L,J+1)-0.5*DUX(L,J+1) ENDDO c RHOR=UUR(1) c UXR=UUR(2)/UUR(1) c PRER=(GAMA-1.)*(UUR(3)-0.5*RHOR*UXR*UXR) RHOR=UUR(1) UXR=UUR(2) PRER=UUR(3) CR=SQRT(GAMA*PRER/RHOR) EIGVR=abs(UXR)+CR UUR(1)=RHOR UUR(2)=RHOR*UXR UUR(3)=PRER/(GAMA-1.)+0.5*RHOR*UXR*UXR FXR(1)=UUR(1)*UXR FXR(2)=UUR(2)*UXR+PRER FXR(3)=UUR(3)*UXR+UXR*PRER SUM=MAX(EIGVR,EIGVL) do L=1,3 FLUX(L,J)=0.5*(FXR(L)+FXL(L)-SUM*(UUR(L)-UUL(L)) ) enddo ENDDO DO J=NIA,NIB RR=DT/(XO(J+1)-XO(J)) DO L=1,3 UO(L,J)=UO(L,J)/3.+2.*(UN(L,J)-RR*(FLUX(L,J)-FLUX(L,J-1)))/3. ENDDO ENDDO DO J=NIA,NIB DO L=1,3 UN(L,J)=UO(L,J) ENDDO ENDDO c IF(IT.LE.100.OR.MOD(IT,50).eq.0) THEN c ITA=ITA+1 c DO J=NIA,NIB c WRITE(13,*) XO(J),time c ENDDO c ENDIF CALL AMMM(time,IADAP) IT=IT+1 TIME=TIME+DT WRITE(*,*) 'TIME,IT=',TIME,IT GOTO 100 888 continue RETURN END FUNCTION Slope_Limiter(DUL,DUR) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT REAL*8 (A-H,O-Z) C van Leer Limiter Slope_Limiter=(SIGN(1.D0,DUL)+SIGN(1.D0,DUR)) & *ABS(DUL*DUR)/(ABS(DUL)+ABS(DUR)+1.0E-10) C Minmod Limiter c Slope_Limiter=0.5*(SIGN(1.D0,DUL)+SIGN(1.D0,DUR)) c & *MIN(ABS(DUL), ABS(DUR) ) END c FUNCTION DABS(X) c implicit real*8 (a-h,o-z) c PI=ASIN(1.0)*2.0 c DABS=X*DERF(X) c DABS=DABS+exp(-x*x)/sqrt(PI) c RETURN c END