数组越界
自己编了一个求激波管问题精确解的代码 但是运行不出来 求各位大佬帮我看看哪里有问题给下修改意见
program shock_tube
implicit none
! 常数定义
real, parameter :: gamma = 1.4 ! 气体比热比
real, parameter :: L = 10.0 ! 激波管长度
real, parameter :: nx = 20 ! 空间网格数
real, parameter :: dx = L / nx ! 空间步长
real, parameter :: t_final = 1.6 ! 求解的时间
real, parameter :: dt = 0.08 ! 时间步长
! 变量定义
integer :: i, n, tsteps
real :: x, t
real, dimension(0:nx) :: rho, u, p
open(100,file='resulr.dat')
! 初始条件
do i = 0, nx
x = i * dx
if (x <= 0.5 * L) then
rho(i) =0.445d0
u(i) = 0.698d0
p(i) = 3.528d0
else
rho(i) = 0.5d0
u(i) = 0.0d0
p(i) = 0.571d0
end if
end do
! 求解时间步进
t = 0.0
tsteps = t_final / dt
do n = 1, tsteps
! Riemann问题求解
do i = 0, nx-1
x = (i + 0.5) * dx
if (x <= 0.5 * L) then
rho(i) =0.445d0
u(i) = 0.698d0
p(i) = 3.528d0
else
rho(i) = 0.5d0
u(i) = 0.0d0
p(i) = 0.571d0
end if
end do
! 更新各个网格的值
do i = 1, nx-2
rho(i) = rho(i) - dt / dx * (flux(rho, u, p, i) - flux(rho, u, p, i-1))
u(i) = u(i) - dt / dx * (flux(rho*u, u*u + p, u*(u(i)+p(i)/(rho(i)*(gamma-1))), i) - &
flux(rho*u, u*u + p, u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))), i-1))
p(i) = p(i) - dt / dx * (flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i)+p(i)/(rho(i)*(gamma-1)))), &
((rho*u*u+p)*(u(i)+p(i)/(rho(i)*(gamma-1))))/(gamma-1), i) - &
flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1)))), &
((rho*u*u+p)*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))))/(gamma-1), i-1))
end do
! 更新左边界值
i = 0
rho(i) = rho(i+1)
u(i) = -u(i+1)
p(i) = p(i+1)
! 更新右边界值
i = nx-1
rho(i) = rho(i-1)
u(i) = -u(i-1)
p(i) = p(i-1)
! 更新时间
t = t + dt
end do
write(100,*) 'variables=', '"x","rho","u","p"'
! 输出结果
do i = 0, nx
x = i * dx
write(100,"(4(e15.7))") x, rho(i), u(i), p(i)
end do
close(100)
contains
! 定义通量函数
function flux(rho, u, p, i) result(f)
real, intent(in) :: rho(:), u(:), p(:)
integer, intent(in) :: i
real :: f
f = rho(i) * u(i)
end function flux
end program shock_tube function flux(rho, u, p, i) result(f)
real, intent(in) :: rho(0:), u(0:), p(0:)
integer, intent(in) :: i
real :: f
f = rho(i) * u(i)
end function flux
这种传递数组的方式如果没有写下界，默认是1，你这里需要改为默认下界是0
这个最好写成integer
integer, parameter :: nx = 20 ! 空间网格数
Transpose 发表于 2023-9-8 18:23
function flux(rho, u, p, i) result(f)
real, intent(in) :: rho(0:), u ...
十分感谢大佬:-lol
页:
[1]