一维黎曼问题求助
求助:我的计算结果,压力算的不对,我检查代码好多次,都不知道哪错了。这是我根据网上的代码略微修改,基本没改,求大神帮助。
! MAIN (Explicit method)(我的代码)! MacCormack_1D_Rimenn problem
!===============================================================
PROGRAM MacCormack_1D_Rimenn
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL(KIND=8), DIMENSION(NMAX) :: X
REAL(KIND=8), DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB, EV
REAL(KIND=8) :: rou, uu, p, Te, a, M, h
REAL(KIND=8) :: DT, t, TT=0.4 !声明时间离散变量
INTEGER :: n
REAL(KIND=8) :: rou1, u1, p1, rou2, u2, p2
REAL(KIND=8) :: DX
INTEGER :: NX, IX !声明空间离散变量
REAL(KIND=8) :: sigma, PI = 4.0d0*DATAN(1.0d0), GAMA = 1.4, R = 278.0
REAL(KIND=8) :: maxvel = 1e-10, vel, eta=0.25, theta
INTEGER :: I, K
WRITE (*, *) "Enter NX"
!READ(*, *) NX
! DT = adjustable
NX = 1000 ! 网格数
IX = NX+1
DX = 2.0d0/NX ! 计算域(0,2)
DO I = 1, IX
x(i) = dble((I-1)*DX)
end do ! 网格
!-------------------------------------------------------
! Shock tube initialtion condition
!-------------------------------------------------------
rou1=1.0
u1=0.0
p1=1.0
rou2=0.125
u2=0.0
p2=0.1
DO I = 1, (IX-1)/2
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
END DO
DO I = (IX-1)/2+1, IX
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2
END DO
t = 0
! March
! time step(capture signal---stability)
10 DO I = 2, IX-1
uu = U(i,1)/U(i,0)
p = (GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
vel = dsqrt(gama*p/U(i,0))+dabs(uu)
if (vel .gt. maxvel) then
maxvel = vel
end if
END DO
DT = 0.8*DX/maxvel
sigma = DT/DX
!---------------------------------------------------------
! Computation
t = t + DT
write(*,*)'T=',T,'dt=',dt
DO I = 2,IX-1
DO K = 0,2
!开关函数
theta=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0)))/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)
!人工粘性修正
EV(i,k) = U(i,k)+0.5*eta*theta*(U(i+1,k)-2*U(i,k)+U(i-1,k))
END DO
END DO
DO K = 0,2
DO I = 2,IX-1
U(i,k) = EV(i,k)
END DO
END DO
! Predictor
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i)
do i=1,IX
uu=U(i,1)/U(i,0)
p=(GAMA-1)*U(i,2)-0.5*U(i,1)*U(i,1)/U(i,0)
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
DO I = 1, IX-1
DO K = 0, 2
UB(I,K) = U(I,K) - sigma*( E(I+1,K) - E(I,K) )
END DO
END DO
! Corrector
do i=1,IX-1
uu=UB(i,1)/UB(i,0)
p=(GAMA-1)*UB(i,2)-0.5*UB(i,1)*UB(i,1)/UB(i,0)
EB(i,0)=UB(i,1)
EB(i,1)=UB(i,0)*uu*uu+p
EB(i,2)=(UB(i,2)+p)*uu
end do
DO I = 2, IX-1
DO k = 0,2
U(i,k)=0.5*(U(i,k)+UB(i,k))-0.5*sigma*(EB(i,k)-EB(i-1,k))
END DO
END DO
! Boundary condition
DO K = 0, 2
U(1,k) = U(2,k)
U(IX,k) = U(IX-1,k)
END DO
!-------------------------------------------------------
if (t .lt. 0.4) goto 10
!-------------------------------------------------------
! Plot
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)
write(1,81)i*dx,rou,uu,p,U(i,2)
END DO
close(1)
81 format(D20.10,D20.10,D20.10,D20.10,D20.10)
END PROGRAM MacCormack_1D_Rimenn
这是网上的代码
! MacCormack1D.for
!------------------------------------------------------------------------------------------------------------
!利用差分格式求解一维激波管问题(语言版本)(网上的代码)
!--------------------------------------------------------------------------------------------------------------*/
program MacCormack1D
implicit double precision (a-h,o-z)
parameter (M=1000)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:M+1,0:2),Uf(0:M+1,0:2)
dimension Ef(0:M+1,0:2)
!气体常数
GAMA=1.4
PI=3.1415926
!网格数
J=M
!计算区域
dL=2.0
!总时间
TT=0.45
!时间步长因子
Sf=0.8
call Init(U,dx)
T=0
1 dt=CFL(U,dx)
T=T+dt
write(*,*)'T=',T,'dt=',dt
call MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
call bound(U,dx)
if(T.lt.TT)goto 1
call Output(U,dx)
end
!-------------------------------------------------------
!计算时间步长
!入口: U, 当前物理量,dx, 网格宽度;
!返回: 时间步长。
!-------------------------------------------------------
double precision function CFL(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
dmaxvel=1e-10
do 10 i=1,J
uu=U(i,1)/U(i,0)
p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
vel=dsqrt(GAMA*p/U(i,0))+dabs(uu)
if(vel.gt.dmaxvel)dmaxvel=vel
10 continue
CFL=Sf*dx/dmaxvel
end
!-------------------------------------------------------
!初始化
!入口: 无;
!出口: U, 已经给定的初始值,dx,网格宽度。
!-------------------------------------------------------
subroutine Init(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
!初始条件
rou1=1.0
u1=0
v1=0
p1=1.0
rou2=0.125
u2=0
v2=0
p2=0.1
dx=dL/J
do 20 i=0,J/2
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
20 continue
do 21 i=J/2+1,J+1
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2
21 continue
end
!-------------------------------------------------------
!边界条件
!入口: dx,网格宽度;
!出口: U, 已经给定边界。
!-------------------------------------------------------
subroutine bound(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
!左边界
do 30 k=0,2
U(0,k)=U(1,k)
30 continue
!右边界
do 31 k=0,2
U(J+1,k)=U(J,k)
31 continue
end
!-------------------------------------------------------
!根据U计算E
!入口: U,当前U矢量;
!出口: E,计算得到的E矢量,
! U、E定义见Euler方程组。
!-------------------------------------------------------
subroutine U2E(U,E,is,in)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2),E(0:J+1,0:2)
do 40 i=is,in
uu=U(i,1)/U(i,0)
p=(GAMA-1)*(U(i,2)-0.5*U(i,1)*U(i,1)/U(i,0))
E(i,0)=U(i,1)
E(i,1)=U(i,0)*uu*uu+p
E(i,2)=(U(i,2)+p)*uu
40 continue
end
!-------------------------------------------------------
!一维差分格式求解器
!入口: U, 上一时刻U矢量,
! Uf、Ef,临时变量,
! dx,网格宽度,dt,,时间步长;
!出口: U, 计算得到得当前时刻U矢量。
!-------------------------------------------------------
subroutine MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2),Uf(0:J+1,0:2)
dimension Ef(0:J+1,0:2)
r=dt/dx
dnu=0.25
do 60 i=1,J
do 60 k=0,2
!开关函数
q=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0)))/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)
!人工黏性项
Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k))
60 continue
do 61 k=0,2
do 61 i=1,J
U(i,k)=Ef(i,k)
61 continue
call U2E(U,Ef,0,J+1)
do 63 i=0,J
do 63 k=0,2
!U(n+1/2)(i+1/2)
Uf(i,k)=U(i,k)-r*(Ef(i+1,k)-Ef(i,k))
63 continue
!E(n+1/2)(i+1/2)
call U2E(Uf,Ef,0,J)
do 64 i=1,J
do 64 k=0,2
!U(n+1)(i)
U(i,k)=0.5*(U(i,k)+Uf(i,k))-0.5*r*(Ef(i,k)-Ef(i-1,k))
64 continue
end
!-------------------------------------------------------
!输出结果, 用数据格式画图
!入口: U, 当前时刻U矢量,
! dx,网格宽度;
!出口: 无。
!-------------------------------------------------------
subroutine Output(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
open(1,file='result.txt',status='unknown')
do 80 i=0,J+1
rou=U(i,0)
uu=U(i,1)/rou
p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
write(1,81)i*dx,rou,uu,p,U(i,2)
80 continue
close(1)
81 format(D20.10,D20.10,D20.10,D20.10,D20.10)
end 你需要写清楚,改动了哪里? 改动了所有的变量声明,并且取消了所有的子程序 li913 发表于 2020-6-13 10:41
你需要写清楚,改动了哪里?
您好,我改动了所有变量声明,并取消了子程序,可以帮我看看吗。很的感谢 我的建议是先从原始程序开始,确认原始程序没问题,把你改的地方一点一点加进去,比如先改初始化部分,编译执行看有没有问题,然后改其它的。
原始程序可能有些不严谨的地方,你可以修改得更严谨一些,但是不建议改动下标范围、变量名称,使用别人的程序最好适应别人的思路,小程序好改,大程序就会很难办。
最好先改不容易错的地方,多做备份。 亲,我测试结果,两个代码的结果是一样的,都有台阶。 带子程序的我的代码,和我的一整个主程序都算的不对
! MAIN (Explicit method)
! MacCormack_1D_Rimenn problem
!===============================================================
PROGRAM MacCormack_1D_Rimenn
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL(KIND=8), DIMENSION(NMAX) :: X
REAL(KIND=8), DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB
REAL(KIND=8) :: uu, p
REAL(KIND=8) :: DT, t, TT=0.4 !声明时间离散变量
INTEGER :: n
REAL(KIND=8) :: DX
INTEGER :: NX, IX !声明空间离散变量
REAL(KIND=8) :: sigma, PI = 4.0d0*DATAN(1.0d0), GAMA = 1.4, R = 278.0
REAL(KIND=8) :: maxvel = 1e-10, vel
INTEGER :: I, K
WRITE (*, *) "Enter NX"
READ(*, *) NX
! DT = adjustable
! NX = 1000
IX = NX+1
DX = 2.0d0/NX ! 计算域(0,2)
DO I = 1, IX
x(i) = dble((I-1)*DX)
end do ! 网格
! Initial Condition
call Init(IX, U)
t = 0
! March
! time step(capture signal---stability)
10 DO I = 2, IX-1
uu = U(i,1)/U(i,0)
p = (GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
vel = dsqrt(gama*p/U(i,0))+dabs(uu)
if (vel .gt. maxvel) then
maxvel = vel
end if
END DO
DT = 0.8*DX/maxvel
sigma = DT/DX
! Computation
t = t + DT
call MacCormack_1D_Solver(U,sigma,IX)
if (t .lt. 0.4) goto 10
! Plot
call Results(IX, U, DX)
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(IIX, U)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
INTEGER, INTENT(IN) :: IIX
REAL(KIND=8), INTENT(OUT), DIMENSION(NMAX,0:2) :: U
REAL(KIND=8), DIMENSION(NMAX) :: X
REAL(KIND=8) :: GAMA=1.4, rou1, u1, p1, rou2, u2, p2
REAL(KIND=8) :: dx
INTEGER :: i
rou1=1.0
u1=0.0
p1=1.0
rou2=0.125
u2=0.0
p2=0.1
DO I = 1, IIX
if ( i .le. (IIX-1)/2 ) 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(KIND=8), DIMENSION(NMAX) :: X
INTEGER, INTENT(IN) :: ist,ie
REAL(KIND=8), INTENT(IN), DIMENSION(NMAX,0:2) :: U
REAL(KIND=8), INTENT(OUT), DIMENSION(NMAX,0:2) :: E
REAL(KIND=8) :: uu, p
REAL(KIND=8) :: gama=1.4
INTEGER :: i
do i=ist,ie
uu=U(i,1)/U(i,0)
p=(GAMA-1)*U(i,2)-0.5*U(i,1)*U(i,1)/U(i,0)
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,sigma,IX)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL(KIND=8), INTENT(IN) :: sigma
INTEGER, INTENT(IN) :: IX
REAL(KIND=8), INTENT(INOUT), DIMENSION(NMAX,0:2) :: U
REAL(KIND=8), DIMENSION(NMAX) :: X
REAL(KIND=8), DIMENSION(NMAX,0:2) :: UB, UBB, E, EB,EV
INTEGER :: I, K
REAL(KIND=8) :: eta=0.25, theta
DO I = 2,IX-1
DO K = 0,2
!开关函数
theta=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0)))/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)
!人工粘性修正
EV(i,k) = U(i,k)+0.5*eta*theta*(U(i+1,k)-2*U(i,k)+U(i-1,k))
END DO
END DO
DO K = 0,2
DO I = 2,IX-1
U(i,k) = EV(i,k)
END DO
END DO
! Predictor
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i)
call U2E(U,E,1,IX)
DO I = 1, IX-1
DO K = 0, 2
UB(I,K) = U(I,K) - sigma*( E(I+1,K) - E(I,K) )
END DO
END DO
! Corrector
call U2E(UB,EB,1,IX-1)
DO I = 2, IX-1
DO k = 0,2
!UBB(I,K) = UB(I,K) - sigma * (EB(I,K) - EB(I-1,K))
U(i,k)=0.5*(U(i,k)+UB(i,k))-0.5*sigma*(EB(i,k)-EB(i-1,k))
END DO
END DO
! Updating
!DO k = 0,2
!DO I = 2, IX-1
! U(I,k) =0.5d0*( U(I,k) + UBB(I,k) )
!END DO
!END DO
! Boundary condition
DO K = 0, 2
U(1,k) = U(2,k)
U(IX,k) = U(IX-1,k)
END DO
END Subroutine MacCormack_1D_Solver
!-------------------------------------------------------
! Plot Data
! Input: U, 上一时刻U矢量
! Export: , 计算得到得当前时刻原始变量( conserved variable )
!-------------------------------------------------------
SUBROUTINE Results(IX, U, DX)
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
INTEGER, INTENT(IN) :: IX
REAL(KIND=8), INTENT(IN) :: DX
REAL(KIND=8), INTENT(IN), DIMENSION(NMAX,0:2) :: U
REAL(KIND=8), DIMENSION(NMAX) :: X
REAL(KIND=8) :: rou, uu, p, Te,a,M,h
INTEGER :: I
REAL(KIND=8) :: GAMA = 1.4, R = 278.0
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) (I-1)*DX,rou,uu,p,Te,a,M,h
81 FORMAT(8(F10.5, 1x))
End DO
close(1)
END SUBROUTINE Results
necrohan 发表于 2020-6-15 12:38
我的建议是先从原始程序开始,确认原始程序没问题,把你改的地方一点一点加进去,比如先改初始化部分,编译 ...
谢谢,我就先按您的方法找找问题在哪,加油 li913 发表于 2020-6-15 20:37
亲,我测试结果,两个代码的结果是一样的,都有台阶。
不一样的,仔细看第三个变量的值,我的代码多了一个台阶
页:
[1]
2