dgm(charactvar).f90

间断galerkin方法

!!  discontinuous galerkin method, in one-dimensional case

    module mainvar
    implicit none
    integer,parameter :: m=3, NN=400  !! m for number of basis function, 
    real,parameter :: gama=1.4,small=1.0e-10
    integer :: rk, i_step, step, i,k, n_output  !! i,k should delete for affect
    real :: ro(m,0:NN+1),rou(m,0:NN+1),roe(m,0:NN+1),  &    !! for the reconstrucion coefficient of uh(x,t)
	     u(m,0:NN+1),p(m,0:NN+1),  &
         ro1(m,0:NN),rou1(m,0:NN),roe1(m,0:NN),  &    !! for temporary use in R_K
         roL(0:NN),uL(0:NN),pL(0:NN), roR(0:NN),uR(0:NN),pR(0:NN), &
         rouL(0:NN),roeL(0:NN), rouR(0:NN),roeR(0:NN), &
         ma1(m,NN),ma2(m,NN),ma3(m,NN), mv1(m,NN),mv2(m,NN),mv3(m,NN), &
         F1(0:NN),F2(0:NN),F3(0:NN), FF1(m,NN),FF2(m,NN),FF3(m,NN),FF1p(m,NN),FF2p(m,NN),FF3p(m,NN)
    real :: dt,dx

!    !! basis function, not directly used
!    real :: bs0,bs1,bs2
!    bs0(x)=1
!    bs1(x,xc)=(x-xc)/dx*2.
!    bs2(x,xc)=((x-xc)/dx*2.)**2-0.33333333333333

    end module mainvar

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!  MAIN PROGRAM !!!!!!!!!!!!!!!!!!!!
    program main
    use mainvar
    implicit none

    call initiate
    call integrateOne
    do i_step=1,step
	  write(*,*) "now compute the",i_step,"-th step"
      do rk=1,3
        
        call ghostcell
	    call reconst2
		call boundary
 		call fluxes
!		call flux_LAX
		call integrateThree

        call integrateTwo
               
        call evolution
!		call UpdateNatureVar
      enddo
      if(mod(i_step,n_output)==0)then
        write(*,*) "backup..."
        call output
      endif
    enddo
    end program

!!!!!!!!!! init theory: make the residual least!!!!!!
    subroutine initiate
    use mainvar
    implicit none
	real :: a1,b1,c1,  a2,b2,c2

	open(unit=101,file="init.dat",status="old")
	read(101,*) dt ,step
	read(101,*) a1,b1,c1, a2,b2,c2

    dx=50./NN
!     dt=0.001
!     step=5000
	n_output=100

    do i=1,NN/2
      ro(1,i)=a1
      ro(2,i)=0.
      ro(3,i)=0.

      rou(1,i)=a1*b1
      rou(2,i)=0.
      rou(3,i)=0.

      roe(1,i)=c1/(gama-1.)+0.5*a1*b1*b1
      roe(2,i)=0.
      roe(3,i)=0.
    enddo
    do i=NN/2+1,NN
      ro(1,i)=a2
      ro(2,i)=0.
      ro(3,i)=0.

      rou(1,i)=a2*b2
      rou(2,i)=0.
      rou(3,i)=0.

      roe(1,i)=c2/(gama-1.)+0.5*a2*b2*b2
      roe(2,i)=0.
      roe(3,i)=0.
    enddo
!	call UpdateNatureVar
    end subroutine initiate

!!!!!! ghostcell should aim at the the nature var
    subroutine ghostcell
    use mainvar
    implicit none
    do k=1,m
      ro(k,0)= ro(k,1)
      rou(k,0)= -rou(k,1)
      roe(k,0) = roe(k,1)
	  u(k,0) = -u(k,1)
	  p(k,0) = p(k,1)

      ro(k,NN+1)= ro(k,NN)
      rou(k,NN+1) = -rou(k,NN)
      roe(k,NN+1) = roe(k,NN)
	  u(k,NN+1) = -u(k,NN)
	  p(k,NN+1) = p(k,NN)
    enddo
    end subroutine

    subroutine boundary
    use mainvar
    implicit none
    !! reflect wall condition
    ROL(0) = roR(0)
    UL(0)  = uR(0)
    PL(0)  = pR(0)
	rouL(0)= rouR(0)
	roeL(0)= roeR(0)

    ROR(NN)= roL(NN)
    UR(NN) = uL(NN)
    PR(NN) = pL(NN)
	rouR(NN)= rouL(NN)
	roeR(NN)= roeL(NN)
    end subroutine boundary

    real function MinMod(x,y,z)
    if( x*y>0. .and. y*z>0.) then
	  MinMod=nsgn(x)*min(abs(x),abs(y),abs(z))
	else
	  MinMod=0.0
	endif
	return
	end
	function nsgn(x)
	if(x>0)then
	  nsgn=1;
	else if(x<0)then
	  nsgn=-1
	else
	  nsgn=0;
	endif
	return
	end

	real function dot(x,y)
	real :: x(3),y(3)
	dot=x(1)*y(1)+x(2)*y(2)+x(3)*y(3)
	return
	end

    subroutine reconst
    use mainvar
    implicit none
    real :: ROULt,ROELt,ROURt,ROERt, rop, uup,pp
	real :: UP(3), tro,tu,tp,tH,ta, Lv1(3),Lv2(3),Lv3(3),Rv1(3),Rv2(3),Rv3(3), &
        su1(3),su2(3),su3(3),suL(3),suR(3), v1(3),v2(3),v3(3),vL(3),vR(3),vp(3), &
        tuL(3),tUR(3),   testv(3)

	real,external :: minmod, dot,  tvb,fi1

    do i=1,NN
       !! TVD limiter
       ro(2,i)= tvb(dx*dx, ro(2,i), fi1(ro(2,i), 1./2*fi1( ro(1,i+1)-ro(1,i)+(ro(3,i)-ro(3,i+1))/3., &
                                                           ro(1,i)-ro(1,i-1)+(ro(3,i-1)-ro(3,i))/3.) ))
       ro(3,i)= tvb(dx*dx*dx,ro(3,i),fi1(ro(3,i),1./4*fi1( ro(2,i+1)-ro(2,i), ro(2,i)-ro(2,i-1)) ))

       rou(2,i)= tvb(dx*dx, rou(2,i), fi1(rou(2,i), 1./2*fi1( rou(1,i+1)-rou(1,i)+(rou(3,i)-rou(3,i+1))/3., &
                                                              rou(1,i)-rou(1,i-1)+(rou(3,i-1)-rou(3,i))/3.) ))
       rou(3,i)= tvb(dx*dx*dx,rou(3,i),fi1(rou(3,i),1./4*fi1( rou(2,i+1)-rou(2,i), rou(2,i)-rou(2,i-1)) ))

       roe(2,i)= tvb(dx*dx, roe(2,i), fi1(roe(2,i), 1./2*fi1( roe(1,i+1)-roe(1,i)+(roe(3,i)-roe(3,i+1))/3., &
                                                              roe(1,i)-roe(1,i-1)+(roe(3,i-1)-roe(3,i))/3.) ))
       roe(3,i)= tvb(dx*dx*dx,roe(3,i),fi1(roe(3,i),1./4*fi1( roe(2,i+1)-roe(2,i), roe(2,i)-roe(2,i-1)) ))

        roL(i)=ro(1,i)  +ro(2,i) +2./3*ro(3,i)
        rouLt =rou(1,i) +rou(2,i)+2./3*rou(3,i)
        rouL(i)=rouLt
        roeLt =roe(1,i) +roe(2,i)+2./3*roe(3,i)
        roeL(i)=roeLt
        uL(i) =rouLt/roL(i)
        pL(i) =(roeLt-0.5*roL(i)*uL(i)*uL(i))*(gama-1.)
  
        roR(i-1)=ro(1,i)  -ro(2,i) +2./3*ro(3,i)
        rouRt   =rou(1,i) -rou(2,i)+2./3*rou(3,i)
        rouR(i-1)=rouRt
        roeRt   =roe(1,i) -roe(2,i)+2./3*roe(3,i)
        roeR(i-1)=roeRt
        uR(i-1) = rouRt/roR(i-1)
        pR(i-1) =(roeRt-0.5*roR(i-1)*uR(i-1)*uR(i-1))*(gama-1.)
      enddo
      end subroutine reconst

	  real function tvb(hl,a,b)
	  use mainvar
      real :: ml
	  ml=0.
	  if(abs(a)<= ml*hl)then
	    tvb=a
	  else
	    tvb=b
	  endif
	  return
	  end
	  real function fi1(a,b)
      integer,external :: nsgn
	  if(a*b<=0)then
	    fi1=0.
	  else
	    fi1=nsgn(a)*min(abs(a),abs(b))
	  endif
	  return
	  end

!!!!!  flux: lax-friedrichs 
    subroutine flux_LAX
    use mainvar
    implicit none
    real :: ro_L,u_L,p_L, ro_R,u_R,p_R, al,ar,DL,DR,F1L,F2L,f3l,f1r,f2r,f3r, rom(3), &
            rou_L,roe_L, rou_R,roe_R

    do i=0,NN
      ro_L= max(roL(i),small)
      u_L = uL(i)
      p_L = max(pL(i),small)
      rou_L=rouL(i)
      roe_L=roeL(i)
      ro_R= max(roR(i),small)
      u_R = uR(i)
      p_R = max(pR(i),small)
      rou_R=rouR(i)
      roe_R=roeR(i)

      al=sqrt(gama*p_L/ro_L)
      ar=sqrt(gama*p_R/ro_R)
      rom(1)=max(abs(u_L-al), abs(u_R-ar))
      rom(2)=max(abs(u_L), abs(u_R))
      rom(3)=max(abs(u_L+al), abs(u_R+ar))

      F1L= ro_L*u_L
      F2L= ro_L*u_L*u_L + p_L
      F3L= u_L*p_L*gama/(gama-1.) + 0.5*ro_L*u_L*u_L*u_L

      F1R= ro_R*u_R
      F2R= ro_R*u_R*u_R + p_R
      F3R= u_R*p_R*gama/(gama-1.) + 0.5*ro_R*u_R*u_R*u_R

      F1(i)=0.5*(F1L+F1R - rom(1)*(ro_R-ro_L))
      F2(i)=0.5*(F2L+F2R - rom(2)*(rou_R-rou_L))
      F3(i)=0.5*(F3L+F3R - rom(3)*(roe_R-roe_L))
!	  if(max(abs(f1(i)),abs(f2(i)),abs(f3(i))) > 100. )then
!	    write(*,*) "error"
!		stop
!	  endif
    enddo
    end subroutine flux_LAX

    subroutine integrateOne
    use mainvar
    implicit none
    do i=1,NN
      !! 再求逆
      Ma1(1:m,i)=(/ 1./dx,3./dx,45./4./dx /)
      Ma2(1:m,i)=(/ 1./dx,3./dx,45./4./dx /)
      Ma3(1:m,i)=(/ 1./dx,3./dx,45./4./dx /)
    enddo
    end subroutine

    subroutine integrateTwo
    use mainvar
    implicit none
    real :: a0,a1,a2, b0,b1,b2, c0,c1,c2
    real,external :: fun1,fun2,fun3,fun4, GaussInteg
    common /f123/ a0,a1,a2, b0,b1,b2, c0,c1,c2
    do i=1,NN
      a0=ro(1,i)
      a1=ro(2,i)
      a2=ro(3,i)
      b0=rou(1,i)
      b1=rou(2,i)
      b2=rou(3,i)
      c0=roe(1,i)
      c1=roe(2,i)
      c2=roe(3,i)