the program is about solving continuity [w(1)], momentum [w(2)] and energy [w(3)] equations using FVM method and Lax Wendroff scheme for shock tube. please provide suggestions in equations also. and give suggestion how to eliminate the warnings. the line truncated is also one more error i aced in long equations.
enter code here
program oneD_heatEqn
implicit none ! Variables must be declared
this program aout shock tube solving the continuity, momentum and energy equations using FVM and Lax Wendroff Scheme.
integer,parameter:: n = 100 ! number of points in the grid
real:: tmax = 0.1! time of the simulation
real:: dt !time step
real:: c !
real::h
real:: time = 0
real(kind=4),dimension(n):: w1
real(kind=4),dimension(n):: w2
real(kind=4),dimension(n):: w3
real(kind=4),dimension(n):: w4
real(kind=4),dimension(n):: w5
real(kind=4),dimension(n):: u
real(kind=4),dimension(n):: w1new
real(kind=4),dimension(n):: w2new
real(kind=4),dimension(n):: w3new
real(kind=4),dimension(n):: w4new
real(kind=4),dimension(n):: w5new
real(kind=4),dimension(n):: A
real(kind=4),dimension(n):: A1
real(kind=4),dimension(n):: A2
real(kind=4),dimension(n):: E
real(kind=4),dimension(n):: Enew
real(kind=4),dimension(n):: p
real(kind=4),dimension(n):: pnew
real(kind=4),dimension(n):: T
real(kind=4),dimension(n):: unew
real(kind=4),dimension(n):: Tnew
integer::i
h = 1/real(n) ! Be careful to operate with variables of the same type.
dt=h/(340)
c=(340*dt)/h
!Set Boundary Conditions
w1 (1) =1.346
w2 (1)=0
w3 (1)=387849.9
e (1)=1000*288.15
p (1)=1113250
w1 (n)=1.225
w2 (n)=0
w3 (n)=(101325*1000*288.15)/(287*288.15)
e (n) =1000*288.15
p (n)=101325
!Set Initial Conditions
w1(2:50) = 1.346
w2 (2:50)=0
w3 (2:50)=387849.9
e(2:50)=1000*288.15
p(2:50)=111325
w1(51:N-1)=1.225
w2 (51:N-1)=0
w3 (51:n-1)=(101325*1000*288.15)/(287*288.15)
e (51:n-1) =1000*288.15
p (51:n-1)=101325
T(1:n)=288.15
u(1:n)=0
do i=1,n
w4(i)= w2(i)*((w3(i)/w1(i))+(p(i)/w1(i)))
w5(i)=(((w2(i))**2/w1(i))+p(i))
end do
write(*,*) "c=",c, 'dt=',dt
write(*,*) "w4=",w4
write(*,*) "w5=",w5
OPEN(unit=1,file="datau=0final .dat")
c1: do while ( time < tmax) !Note that you can numerate/name the loop (easier to know where you are)
w1new(1)= w1(1)
w4new(1)=w4(1)
w3new(1)=w3(1)
w5new(1)=w5(1)
w2new(1)=w2(1)
unew(1)=u(1)
enew(1)=e(1)
pnew(1)=p(1)
Tnew(1)=T(1)
L2: do i = 2, n-1
!continuity equation
A(i-(1/2))=(w2(i-1)+w2(i))/2
A(i+(1/2))=(w2(i)+w2(i+1))/2
w1new(i)=w1(i)-((dt/2*h)*(w2(i+1)-w2(i-1)))+(dt**3)/(2*h**2)*((A(i+(1/2))*(w2(i+1)-w2(i)))-(A(i-(1/2))*(w2(i)-w2(i-1))))
!energy equation
w4new(i)= w2(i)*((w3(i)/w1(i))+(p(i)/w1(i)))
A1(i-(1/2))=(w4(i-1)+w4(i))/2
A1(i+(1/2))=(w4(i)+w4(i+1))/2
w3new(i)=w3(i)-((dt/2*h)*(w4(i+1)-w4(i-1)))+(dt**3)/(2*h**2)*((A1(i+(1/2))*(w4(i+1)-w4(i)))-(A1(i-(1/2))*(w4(i)-w4(i-1))))
!momentum equation
w5new(i)=(((w2(i))**2/w1(i))+p(i))
A2(i-(1/2))=(w5(i-1)+w5(i))/2
A2(i+(1/2))=(w5(i)+w5(i+1))/2
w2new(i)=w2(i)-((dt/2*h)*(w5(i+1)-w5(i-1)))+(dt**3)/(2*h**2)*((A2(i+(1/2))*(w5(i+1)-w5(i)))-(A2(i-(1/2))*(w5(i)-w5(i-1))))
unew(i)=w2(i)/w1(i)
!other equations
Enew(i)=(w3(i)/w1(i))-((1/2)*(w2(i)/w1(i))**2)
pnew(i)=0.4*w1(i)*e(i)
Tnew(i)=p(i)/(w1(i)*287)
end do L2
w1new(n)= w1(n)
w4new(n)=w4(n)
w3new(n)=w3(n)
w5new(n)=w5(n)
w2new(n)=w2(n)
Unew(n)=u(n)
Enew(n)=E(n)
pnew(n)=p(n)
Tnew(n)=T(n)
w1(1:N)=w1new(1:N) ! Overwrite old value by new values
w4(1:N)=W4new(1:N)
w3(1:N)=w3new(1:N)
w5(1:N)=W5new(1:N)
w2(1:N)=w2new(1:N)
u(1:N)=unew(1:N)
e(1:N)=enew(1:N)
p(1:N)=pnew(1:N)
T(1:N)=Tnew(1:N)
time = time + dt
if (mod(int(time/Dt),int(tmax/Dt)/100)==0) then
write(*,fmt="(F10.3,1X,I5,'% completed',A)",advance='yes') time, int(100*time/tmax), char(13) ! Overwrite the line
end if
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', w1(1:N) ! Do output but keep only 100 profils
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', w2(1:N)
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', w3(1:N)
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', e(1:N)
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', p(1:N)
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', u(1:N)
if (mod(int(time/dt),int(tmax/dt)/100)==0) write(1,*) 'node(i)=1,N', T(1:N)
end do c1
write(*,fmt="(1X,I5,'% completed',A)") int(100*time/tmax), char(13)
call system('gnuplot -p output.plt')
close(1)
write(*,*) "Program finished successfully"
end program oneD_heatEqn
