|
子程序是个数据提取功能,希望x一列,然后不同时间点下u的值一列(筛选时间点)
问题:1.这个问题中,当我选用不同的DT、DX时,比如DT=0.1,DX=40,会运行出错, 子程序里read(9, *) it, t
read(9, *) (tempRes(I), I = 1, IX)提示,end of file
2.取dt=0.01,dx=0.04运行成功时,当我的筛选时间为t=1.0、2.0、3.0都可以正确的将txt中的行数据转化为列数据,但是当我取t=0.1、0.2等时就会出现不匹配的数据,而且在t=4.0时就全部变成了0。
我附上需要转化的数据文件,希望大神帮我看看,真的是困扰好多天了,感谢感谢!超感谢!
! MAIN (Explicit method)
! Two-step predictor-corrector scheme of MacCormack for NON-linear Euler Equationn
!===============================================================
PROGRAM PC
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
REAL*8, DIMENSION(NMAX) :: U, UB, UBB, E, EB
REAL*8 :: DT, t !声明时间离散变量
INTEGER :: n, NLIMIT = 400
REAL*8 :: DX, x
INTEGER :: NX, IX !声明空间离散变量
REAL*8 :: sigma, PI = 4.0d0*DATAN(1.0d0)
INTEGER :: I
REAL*8 :: fun
WRITE (*, *) "Enter DT and NX"
READ (*, *) DT,NX
! DT = 0.01
! NX = 50
open(unit=1, file="MAC-BOut.txt")
IX = NX+1
DX = 4.0d0/NX
sigma = DT/DX
write(1, 310) DX, DT, sigma
310 FORMAT("DX =",F6.3," DT =",F7.5," sigma = 1.0*dt/dx =",F6.3)
! Initial COndition
DO I = 1, IX
x = dble((I-1)*DX)
if ( x .le. 2.0d0 ) then
U(I) = 1
else
U(I) = 0.0d0
end if
END DO
n = 0
t = dble(n*DT)
open(unit=2, file="MAC-BSum.txt")
WRITE(1,*) "Time = ", t
WRITE (1,300)(U(I), I = 1, IX)
write(2, 300) ((I-1)*DX, I = 1, IX )
write(2, 330) n, t
330 FORMAT(I6, f17.7)
write(2, 300) (U(I), I = 1, IX)
! March
10 n = n + 1
t = dble(n*DT)
! Boundary condition
U(1) = 0.0d0
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i) (5.3.3a)
! Predictor (5.3.3a)
DO I = 1, NX
E(I) = fun(U(I))
UB(I) = U(I) - sigma * ( E(I+1) - E(I) )
END DO
! Numerical Boundary Condition
UB(IX) = 2.0d0*UB(IX-1) - UB(IX-2)
! Corrector (5.3.3b)
DO I = 2, IX
EB(I) = fun(UB(I))
UBB(I) = U(I) - sigma * (EB(I) - EB(I-1))
END DO
! Boundary condition
UBB(1) = 0.0d0
! Updating (5.3.3c)
DO I = 1, IX
U(I) = 0.5d0*( UB(I) + UBB(I) )
END DO
! Numerical boundary condition (extrapolation)
! U(IX) = 2.0d0*U(IX-1) - U(IX-2)
WRITE(1,370) n, t
WRITE (1,300)(U(I), I = 1, IX)
write(2, 330) n, t
write(2,300) (U(I), I = 1, IX)
370 FORMAT("n = ",I6," Time = ", f10.4)
if ( n .lt. NLIMIT ) goto 10
close(1)
close(2)
call Results(IX, NLIMIT, DX, DT)
WRITE(*, *) "Numerical Solution is in MAC-BOut.txt"
WRITE(*, *) "Summary file is in MAC-BSum.txt"
write(*, *)
WRITE(*, *) "Calculations are successfully completed. "
WRITE(*, *) "Hit any key to close DOS window!"
300 FORMAT(10(F10.5, 1x))
STOP
END PROGRAM PC
SUBROUTINE Results(IX, n, dx, dt)
! to write results in column for plot purpose
IMPLICIT NONE
INTEGER, PARAMETER :: NMAX = 1001
INTEGER :: IX, n, ic, k, i ,j, it
REAL*8 :: t, dx, dt
REAL*8 :: Res(6,NMAX), tempRes(NMAX)
open(unit=9, file="MAC-BSum.txt")
rewind 9
ic = 1
read(9,*) (Res(ic, K), k=1, IX )
do j = 0, n
read(9, *) it, t
read(9, *) (tempRes(I), I = 1, IX)
if (t == 0.0 .or. &
t .eq. 0.2 .or. &
t .eq. 1.5 .or. &
t .eq. 2.0 .or. &
t .eq. 3.0) then
ic = ic + 1
do k = 1, IX
Res(ic,K) = tempRes(K)
end do
end if
end do
close(9)
open(unit=10, file="MAC-BSumecplot-1.txt")
rewind 10
write(10,31) DX, DT, DT/DX
31 FORMAT(21x,"DX =",F6.3," DT =",F7.5," tau = dt/dx =",F6.3)
write(10, 32)
32 FORMAT(8x,"x ",13x,"t=0 ",8x,"t=1.0 ",7x,"t=2.0 ",8x, &
"t=3.0 ",7x,"t=4.0")
DO k = 1, IX
write(10, 33) (Res(ic,k), ic =1, 6)
END DO
33 FORMAT(6(f14.7, 1x))
close(10)
RETURN
END
REAL*8 FUNCTION fun(X)
IMPLICIT NONE
REAL*8 :: X
! Burger's Equation
fun = 0.5*X*X
END FUNCTION
|
|