maccormack.f90

用maccormack算法求解偏微分方程

program maccormack
implicit none
integer,parameter::Nmax=200
real(kind=8),parameter::r=1.4
real(kind=8),parameter::dt0=0.1
real(kind=8),parameter::time=0.14
real(kind=8),parameter::dx=1.0/Nmax

real(kind=8)::dt,t,tez
integer i,n
character fname*10
real(kind=8)::rou(0:Nmax),u(0:Nmax),p(0:Nmax),e(0:Nmax)
real(kind=8)::u1(0:Nmax,3),u2(0:Nmax,3),u3(0:Nmax,3)
real(kind=8)::F1(0:Nmax),F2(0:Nmax),F3(0:Nmax),c(0:Nmax)
real(kind=8)::x(0:Nmax)

!网格划分
x(0)=0.0
do i=1,Nmax
x(i)=x(i-1)+dx
end do

!初始条件
do i=0,Nmax/2
 rou(i)=0.445
 u(i)=0.698
 p(i)=3.528
end do

do i=Nmax/2+1,Nmax
 rou(i)=0.5
 u(i)=0.0
 p(i)=0.571
end do

tez=0.0
do i=0,Nmax
 c(i)=sqrt(r*p(i)/rou(i))
 tez=max(tez,abs(u(i))+c(i)) 
end do

dt=min(dt0,dx/tez)
t=0.0
t=t+dt

do i=0,Nmax
u1(i,1)=rou(i)
u2(i,1)=rou(i)*u(i)
u3(i,1)=p(i)/(r-1.0)+0.5*rou(i)*(u(i)**2.0)
end do

!maccormack 迭代
do while(t<=time)
!计算中间量迭代步
  print*,n
  n=n+1
  do i=0,Nmax
    F1(i)=rou(i)*u(i)
    F2(i)=p(i)+rou(i)*(u(i)**2.0)
    F3(i)=(r*p(i)/(r-1.0)+0.5*rou(i)*(u(i)**2.0))*u(i)
  end do
 
  do i=0,Nmax-1
    u1(i,2)=u1(i,1)-(dt/dx)*(F1(i+1)-F1(i))
    u2(i,2)=u2(i,1)-(dt/dx)*(F2(i+1)-F2(i))
    u3(i,2)=u3(i,1)-(dt/dx)*(F3(i+1)-F3(i))
  end do
   u1(Nmax,2)=u1(Nmax-1,2)
   u2(Nmax,2)=u2(Nmax-1,2)
   u3(Nmax,2)=u3(Nmax-1,2)
   
  
   do i=0,Nmax
    rou(i)=u1(i,2)
    u(i)=u2(i,2)/u1(i,2)
    p(i)=(r-1.0)*(u3(i,2)-0.5*rou(i)*u(i)**2.0)
   end do
 
  !最后迭代步
  do i=0,Nmax
    F1(i)=rou(i)*u(i)
    F2(i)=p(i)+rou(i)*(u(i)**2.0)
    F3(i)=(r*p(i)/(r-1.0)+0.5*rou(i)*(u(i)**2.0))*u(i)
   end do
   
   do i=1,Nmax
    u1(i,3)=0.5*(u1(i,1)+u1(i,2)-(dt/dx)*(F1(i)-F1(i-1)))
    u2(i,3)=0.5*(u2(i,1)+u2(i,2)-(dt/dx)*(F2(i)-F2(i-1)))
    u3(i,3)=0.5*(u3(i,1)+u3(i,2)-(dt/dx)*(F3(i)-F3(i-1)))
   end do
     
   do i=1,Nmax
    u1(i,1)=u1(i,3)
    u2(i,1)=u2(i,3)
    u3(i,1)=u3(i,3)
   end do
    
    u1(0,1)=u1(1,1)
    u2(0,1)=u2(1,1)
    u3(0,1)=u3(1,1)
   
    tez=0.0
    do i=0,Nmax
      rou(i)=u1(i,1)
      u(i)=u2(i,1)/u1(i,1)
      p(i)=(r-1.0)*(u3(i,1)-0.5*rou(i)*u(i)**2.0)
      c(i)=sqrt(r*p(i)/rou(i))
      tez=max(tez,abs(u(i))+c(i))
    end do
   dt=min(dt0,dx/tez)
   t=t+dt
!   print*,u
 !  pause
end do

 do i=0,Nmax
  rou(i)=u1(i,1)
  u(i)=u2(i,1)/u1(i,1)
  p(i)=(r-1.0)*(u3(i,1)-0.5*rou(i)*u(i)**2.0)
  e(i)=u3(i,1)
 end do
!给定文件名并输出数据
open(2006,file="density1.txt")
  do i=0,Nmax
    write(2006,*)rou(i)
  end do
  
  open(2007,file="u1.txt")
  do i=0,Nmax
    write(2007,*)u(i)
  end do
  open(2008,file="pressure1.txt")
  do i=0,Nmax
    write(2008,*)p(i)
  end do

  open(2009,file="x1.txt")
  do i=0,Nmax
    write(2009,*)x(i)
  end do

end program maccormack