|
一维黎曼问题的欧拉求解代码,利用MacCormack格式计算,最后输出原始变量密度、压力、速度。我自己编写后,利用gfortran编译(code-block集成环境),弹出一个串口:
出现了一个问题,导致MacCormack_1D_Rimenn.exe无法运行,请关闭程序
请问这是咋回事?是写的代码有问题?还是软件存在问题啊,我截了图,希望大神帮我解答!感谢感谢感感谢!
! MAIN (Explicit method)
! MacCormack_1D_Rimenn problem
!===============================================================
PROGRAM MacCormack_1D_Rimenn
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: X
REAL*8, DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB
REAL*8 :: DT, t, TT=0.01 !声明时间离散变量
INTEGER :: n
REAL*8 :: DX
INTEGER :: NX, IX !声明空间离散变量
REAL*8 :: sigma, PI = 4.0d0*DATAN(1.0d0), GAMA = 1.4, R = 278
INTEGER :: I, K
WRITE (*, *) "Enter DT and NX"
READ (*, *) DT,NX
! DT = 0.1
! NX = 50
IX = NX+1
DX = 20d0/NX
sigma = DT/DX
DO I = 1, IX
x(i) = dble((I-1)*DX)
end do
! Initial COndition
call Init(U)
n = 0
t = dble(n*DT)
! March
10 n = n + 1
t = dble(n*DT)
call MacCormack_1D_Solver(U,UB,UBB,E,EB,dx,dt)
if ( t .lt. 0.01 ) goto 10
! Plot
call Results(IX, U)
WRITE(*, *) "Numerical Solution is in MACCORMACK-1D-Rimenn.txt"
WRITE(*, *) "Calculations are successfully completed. "
WRITE(*, *) "Hit any key to close DOS window!"
STOP
END PROGRAM MacCormack_1D_Rimenn
!-------------------------------------------------------
! Shock tube initialtion condition
!-------------------------------------------------------
subroutine Init(U)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: X
REAL*8, DIMENSION(NMAX,0:2) :: U
REAL*8 :: GAMA, rou1, u1, p1, rou2, u2, p2
REAL*8 :: dx
INTEGER :: i, IX
rou1=1.0
u1=0
p1=100
rou2=0.125
u2=0
p2=10
DO I = 1, IX
x(i) = dble((I-1)*DX)
if ( x(i) .le. 0 ) then
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
else
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=0.5*rou2*u2*u2 + p2/(GAMA-1)
end if
END DO
END subroutine Init
!-------------------------------------------------------
! 根据U计算E
! Input:U, 当前U矢量
! Export: E,计算得到的E矢量( flux term )
!-------------------------------------------------------
Subroutine U2E(U,E,ist,ie)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: X
REAL*8, DIMENSION(NMAX,0:2) :: U, E
REAL*8 :: uu, p
REAL*8 :: gama
INTEGER :: i, ist, ie
do i=ist,ie
uu=U(i,1)/U(i,0)
p=(GAMA-1)*U(i,2)
E(i,0)=U(i,1)
E(i,1)=U(i,0)*uu*uu+p
E(i,2)=(U(i,2)+p)*uu
end do
end Subroutine U2E
!-------------------------------------------------------
! 1D_MacCormack差分格式求解器
! Input: U, 上一时刻U矢量
! Export: U, 计算得到得当前时刻U矢量( conserved variable )
!-------------------------------------------------------
Subroutine MacCormack_1D_Solver(U,UB,UBB,E,EB)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: X
REAL*8, DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB
INTEGER :: I, K, IX, NX
REAL*8 :: sigma
call U2E(U,E,1,IX)
! Predictor
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i)
DO K = 0, 2
DO I = 1, NX
UB(I,K) = U(I,K) - sigma * ( E(I+1,K) - E(I,K) )
END DO
END DO
! Numerical Boundary Condition
DO K = 0, 2
UB(IX,K) = UB(IX-1,K)
END DO
! Corrector (5.3.3b)
call U2E(UB,EB,1,IX)
DO k = 0,2
DO I = 2, IX
UBB(I,K) = UB(I,K) - sigma * (EB(I,K) - EB(I-1,K))
END DO
END DO
! Boundary condition
DO K = 0, 2
UBB(1,k) = UBB(2,k)
END DO
! Updating
DO k = 0,2
DO I = 1, IX
U(I,k) = 0.5d0*( U(I,k) + UBB(I,k) )
END DO
END DO
END Subroutine MacCormack_1D_Solver
!-------------------------------------------------------
! Plot Data
! Input: U, 上一时刻U矢量
! Export: , 计算得到得当前时刻原始变量( conserved variable )
!-------------------------------------------------------
SUBROUTINE Results(IX, U)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: X
REAL*8, DIMENSION(NMAX,0:2) :: U
REAL*8 :: rou, uu, p, Te,a,M,h
INTEGER :: I, IX
REAL*8 :: gama, R
open(1,file='MacCormack_1D_Rimenn.txt',status='unknown')
DO i=1,IX
rou=U(i,0)
uu=U(i,1)/rou
p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
Te=p/(rou*R)
a=sqrt( gama*R*Te )
M=uu/a
h=gama*R*Te/(gama-1)
write(1,81) x(i),rou,uu,p,Te,a,M,h
81 FORMAT(8(F10.5, 1x))
End DO
close(1)
END SUBROUTINE Results
|
|