数组越界

IDCY

自己编了一个求激波管问题精确解的代码 但是运行不出来 求各位大佬帮我看看哪里有问题
给下修改意见

IDCY

[Fortran] 纯文本查看 复制代码

查看代码

打印?

001	program shock_tube

002	implicit none

003

004

! 常数定义

005	real, parameter :: gamma = 1.4       ! 气体比热比

006	real, parameter :: L = 10.0           ! 激波管长度

007	real, parameter :: nx = 20         ! 空间网格数

008	real, parameter :: dx = L / nx       ! 空间步长

009	real, parameter :: t_final = 1.6     ! 求解的时间

010	real, parameter :: dt = 0.08       ! 时间步长

011

012

! 变量定义

013	integer :: i, n, tsteps

014	real :: x, t

015	real, dimension(0:nx) :: rho, u, p

016	open(100,file='resulr.dat')

017

! 初始条件

018	do i = 0, nx

019	x = i * dx

020

021	if (x <= 0.5 * L) then

022	rho(i) =0.445d0

023	u(i) = 0.698d0

024	p(i) = 3.528d0

025

else

026	rho(i) = 0.5d0

027	u(i) = 0.0d0

028	p(i) = 0.571d0

029

end if

030

end do

031

032

! 求解时间步进

033

t = 0.0

034	tsteps = t_final / dt

035

036	do n = 1, tsteps

037	! Riemann问题求解

038	do i = 0, nx-1

039	x = (i + 0.5) * dx

040

041	if (x <= 0.5 * L) then

042	rho(i) =0.445d0

043	u(i) = 0.698d0

044	p(i) = 3.528d0

045

else

046	rho(i) = 0.5d0

047	u(i) = 0.0d0

048	p(i) = 0.571d0

049

end if

050

end do

051

052	! 更新各个网格的值

053	do i = 1, nx-2

054	rho(i) = rho(i) - dt / dx * (flux(rho, u, p, i) - flux(rho, u, p, i-1))

055	u(i) = u(i) - dt / dx * (flux(rho*u, u*u + p, u*(u(i)+p(i)/(rho(i)*(gamma-1))), i) - &

056	flux(rho*u, u*u + p, u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))), i-1))

057	p(i) = p(i) - dt / dx * (flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i)+p(i)/(rho(i)*(gamma-1)))), &

058	((rho*u*u+p)*(u(i)+p(i)/(rho(i)*(gamma-1))))/(gamma-1), i) - &

059	flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1)))), &

060	((rho*u*u+p)*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))))/(gamma-1), i-1))

061

end do

062

063

! 更新左边界值

064

i = 0

065	rho(i) = rho(i+1)

066	u(i) = -u(i+1)

067	p(i) = p(i+1)

068

069

! 更新右边界值

070

i = nx-1

071	rho(i) = rho(i-1)

072	u(i) = -u(i-1)

073	p(i) = p(i-1)

074

075

! 更新时间

076	t = t + dt

077

end do

078

079	write(100,*) 'variables=', '"x","rho","u","p"'

080

081

! 输出结果

082	do i = 0, nx

083	x = i * dx

084	write(100,"(4(e15.7))") x, rho(i), u(i), p(i)

085

end do

086	close(100)

087

088

contains

089

090

! 定义通量函数

091	function flux(rho, u, p, i) result(f)

092	real, intent(in) :: rho(:), u(:), p(:)

093	integer, intent(in) :: i

094

real :: f

095

096	f = rho(i) * u(i)

097

098	end function flux

099

100	end program shock_tube

Transpose

[Fortran] 纯文本查看 复制代码

查看代码

打印?

1	function flux(rho, u, p, i) result(f)

2	real, intent(in) :: rho(0:), u(0:), p(0:)

3	integer, intent(in) :: i

4

real :: f

5

6	f = rho(i) * u(i)

7

8	end function flux

这种传递数组的方式如果没有写下界，默认是1，你这里需要改为默认下界是0

这个最好写成integer

[Fortran] 纯文本查看 复制代码

查看代码

打印?

1	integer, parameter :: nx = 20         ! 空间网格数

IDCY

Transpose 发表于 2023-9-8 18:23
[mw_shl_code=fortran,true]  function flux(rho, u, p, i) result(f)
    real, intent(in) :: rho(0:), u ...

十分感谢大佬