dgm(charactvar).f90

间断galerkin方法

      mv1(1:m,i)=(/0., 2*rou(1,i), 4./3*rou(2,i)/)

      mv2(1,i)=0.
      mv2(2,i)=GaussInteg(-1.,1.,fun1, 1.0 )
      mv2(3,i)=GaussInteg(-1.,1.,fun2, 1.0 )

      mv3(1,i)=0.
      mv3(2,i)=GaussInteg(-1.,1.,fun3, 1.0 )
      mv3(3,i)=GaussInteg(-1.,1.,fun4, 1.0 )
    enddo
    end subroutine

    real function fun1(x)
    use mainvar
    implicit none
    real :: a0,a1,a2, b0,b1,b2, c0,c1,c2,x
    common /f123/ a0,a1,a2, b0,b1,b2, c0,c1,c2
    fun1=((3 - gama)*(b0 + b1*x + b2*(-0.3333333333333333 + x ** 2)) **  &
          2)/(2.*(a0 + a1*x + a2*(-0.3333333333333333 + x ** 2))) + (-1 +  &
        gama)*(c0 + c1*x + c2*(-0.3333333333333333 + x ** 2))
    return
    end
    real function fun2(x)
    use mainvar
    implicit none
    real :: a0,a1,a2, b0,b1,b2, c0,c1,c2,x
    common /f123/ a0,a1,a2, b0,b1,b2, c0,c1,c2
    fun2=2*x*(((3 - gama)*(b0 + b1*x + b2*(-0.3333333333333333 + x ** 2)) **  &
              2)/(2.*(a0 + a1*x + a2*(-0.3333333333333333 + x ** 2))) + (-1 +  &
            gama)*(c0 + c1*x + c2*(-0.3333333333333333 + x ** 2)))
    return
    end
    real function fun3(x)
    use mainvar
    implicit none
    real :: a0,a1,a2, b0,b1,b2, c0,c1,c2,x
    common /f123/ a0,a1,a2, b0,b1,b2, c0,c1,c2
    fun3=((b0 + b1*x + b2*(-0.3333333333333333 + x ** 2))*((1 - gama)*(b0 + b1*x +  &
	     b2*(-0.3333333333333333 + x ** 2)) ** 2 + 2*gama*(a0 + a1*x + a2*(-0.3333333333333333 + x ** 2))*(c0 + c1*x + &
        c2*(-0.3333333333333333 + x ** 2))))/(2.*(a0 + a1*x + a2*(-0.3333333333333333 + x ** 2)) ** 2)
    return
    end
    real function fun4(x)
    use mainvar
    implicit none
    real :: a0,a1,a2, b0,b1,b2, c0,c1,c2,x
    common /f123/ a0,a1,a2, b0,b1,b2, c0,c1,c2
    fun4=(x*(b0 + b1*x + b2*(-0.3333333333333333 + x ** 2))*((1 - &
         gama)*(b0 + b1*x + b2*(-0.3333333333333333 + x ** 2)) ** 2 +  &
         2*gama*(a0 + a1*x + a2*(-0.3333333333333333 + x ** 2))*(c0 + c1*x +  &
         c2*(-0.3333333333333333 + x ** 2))))/(a0 + a1*x + &
         a2*(-0.3333333333333333 + x ** 2)) ** 2
    return
    end

    subroutine IntegrateThree  !! namely FF1=f1*M^(-1)*(0,1/2,....),  FF2, FF3
    use mainvar
    implicit none
    real :: v1(3),v2(3)
    v1=(/1., 1.,0.66666666667/)
    v2=(/1.,-1.,0.66666666667/)
    do i=1,NN
      !! in fact, (/ Ma1(1,1)*v1(1), Ma1(2,2)*v1(2), .. /) is OK enough
      FF1(1:m,i)= f1(i)  *(/v1(1),v1(2),v1(3)/)
      FF1p(1:m,i)=f1(i-1)*(/v2(1),v2(2),v2(3)/)

      FF2(1:m,i) =f2(i)  *(/v1(1),v1(2),v1(3)/)
      FF2p(1:m,i)=f2(i-1)*(/v2(1),v2(2),v2(3)/)

      FF3(1:m,i) =f3(i)  *(/v1(1),v1(2),v1(3)/)
      FF3p(1:m,i)=f3(i-1)*(/v2(1),v2(2),v2(3)/)
    enddo
    end subroutine IntegrateThree

!!!!!!!!!! evolution !!!!!!!!!!!!!!
    subroutine evolution
    use mainvar
    implicit none
	real :: temp

    if(rk==1)then
     do i=1,NN
     do k=1,m
      ro1(k,i) = ro(k,i)  !! for temporary storage
      rou1(k,i)= rou(k,i)
      roe1(k,i)= roe(k,i)
!! ??? mv1=M^(-1)*mv1, 
      ro(k,i) = ro(k,i) + ma1(k,i)*dt*( mv1(k,i)-FF1(k,i)+FF1p(k,i) )
      rou(k,i)= rou(k,i)+ ma2(k,i)*dt*( mv2(k,i)-FF2(k,i)+FF2p(k,i) )
      roe(k,i)= roe(k,i)+ ma3(k,i)*dt*( mv3(k,i)-FF3(k,i)+FF3p(k,i) )
	  temp=max(abs(ro(k,i)),abs(rou(k,i)),abs(roe(k,i)))
!	  if(temp>1000. )then
!	    write(*,*) "error here in evo!"
!		stop
!	  endif
     enddo
     enddo
    else if(rk==2)then
     do i=1,NN
     do k=1,m
      ro(k,i)  = 0.75*ro1(k,i) +0.25*ro(k,i) +0.25*ma1(k,i)*dt*(mv1(k,i)-FF1(k,i)+FF1p(k,i))
      rou(k,i) = 0.75*rou1(k,i)+0.25*rou(k,i)+0.25*ma2(k,i)*dt*(mv2(k,i)-FF2(k,i)+FF2p(k,i))
      roe(k,i) = 0.75*roe1(k,i)+0.25*roe(k,i)+0.25*ma3(k,i)*dt*(mv3(k,i)-FF3(k,i)+FF3p(k,i))
     enddo
     enddo
    else ! rk==3
     do i=1,NN
     do k=1,m
      ro(k,i)  = 0.3333333333*ro1(k,i) +0.6666666667*ro(k,i) +0.66666666667*ma1(k,i)*dt*(mv1(k,i)-FF1(k,i)+FF1p(k,i))
      rou(k,i) = 0.3333333333*rou1(k,i)+0.6666666667*rou(k,i)+0.66666666667*ma2(k,i)*dt*(mv2(k,i)-FF2(k,i)+FF2p(k,i))
      roe(k,i) = 0.3333333333*roe1(k,i)+0.6666666667*roe(k,i)+0.66666666667*ma3(k,i)*dt*(mv3(k,i)-FF3(k,i)+FF3p(k,i))
     enddo
     enddo
    endif

    end subroutine evolution

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine output
    use mainvar
    character*20 outfile
	k = i_step/n_output
	if(k<10) then
      write(outfile,'(a,i1,a)') '00',k,'.dat'
    else if(k<100) then
      write(outfile,'(a,i2,a)') '0',k,'.dat'
    else
      write(outfile,'(i3,a)') k,'.dat'
    end if
    open(unit=103,file=outfile)
    write(103,*) "title= shocktube"
    write(103,*) 'variables="x","ro","u","p"'
    write(103,*) 'zone i=',NN," f= POINT"
    do i=1,NN
      xx=dx*(i-0.5)
      tro=ro(1,i)-ro(3,i)/3.
      tu =(rou(1,i)-rou(3,i)/3.)/tro
      troe=roe(1,i)-roe(3,i)/3.
      tp =(troe-0.5*tro*tu*tu)*(gama-1.)
      write(103,*) xx, tro, tu, tp
    enddo
    close(103)
    end subroutine output


!!!!!!!!!!!!!!  Assist function !!!!!!!
    REAL FUNCTION GaussInteg(A,B,F,EPS)
	REAL :: T(5),C(5)
	REAL :: A,B,F,G,S,P,H,AA,BB,W,X,Q
    T=(/-0.9061798459,-0.5384693101,0.0, &
               0.5384693101,0.9061798459/)
	C=(/0.2369268851,0.4786286705,0.5688888889, &
               0.4786286705,0.2369268851/)
	M=5
	S=(B-A)*0.001
	P=0.0
10	H=(B-A)/M
	G=0.0
	DO 30 I=1,M
	  AA=A+(I-1)*H
	  BB=A+I*H
	  W=0.0
	  DO 20 J=1,5
	    X=((BB-AA)*T(J)+(BB+AA))/2.0
	    W=W+F(X)*C(J)
20	  CONTINUE
	  G=G+W
30	CONTINUE
	G=G*H/2.0
	Q=ABS(G-P)/(1.0+ABS(G))
	IF ((Q.GE.EPS).AND.(ABS(H).GT.ABS(S))) THEN
	  P=G
	  M=M+1
	  GOTO 10
	END IF
	GaussInteg=G
	RETURN
	END


!!!!!!!!!! FLUEX:  HLL !!!!!!!!!!!!!!!!
    subroutine fluxes
    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)

    do i=0,NN
      ro_L= max(roL(i),small)
      u_L = uL(i)
      p_L = max(pL(i),small)
      ro_R= max(roR(i),small)
      u_R = uR(i)
      p_R = max(pR(i),small)

      !! use HLL
      al=sqrt(gama*p_L/ro_L)
      ar=sqrt(gama*p_R/ro_R)
      DL=min(u_L-al, u_R-ar)
      DR=max(u_L+al, 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

      !! compute the flux
	  if(DL>0.) then
	    F1(i)=F1L
		F2(i)=F2L
		F3(i)=F3L
      else if(DR<0) then
	    F1(i)=F1R
		F2(i)=F2R
		F3(i)=F3R
	  else
        F1(i) = ( DR*F1L-DL*F1R+DR*DL*(ro_R-ro_L) )/(DR-DL)
        F2(i) = ( DR*F2L-DL*F2R+DR*DL*(u_R - u_L) )/(DR-DL)
        F3(i) = ( DR*F3L-DL*F3R+DR*DL*(p_R - p_L) )/(DR-DL)
	  endif
    enddo
    end subroutine fluxes


     subroutine reconst2
     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), v1(3),v2(3),v3(3), &
         tuL(3),tUR(3)
 	 real,external :: minmod, dot
 
     do i=1,NN
        UP=(/0.5*(ro(1,i)+ro(1,i+1)),0.5*(rou(1,i)+rou(1,i+1)),0.5*(roe(1,i)+roe(1,i+1))/)
        tro=Up(1)
        tu =Up(2)/Up(1) 
        tp =(Up(3)-0.5*tro*tu*tu)*(gama-1.)
        tH=gama/(gama-1.)*tp/tro+ 0.5*tu*tu
        ta = sqrt(gama*tp/tro)
        !! notice: Lv 横向量，Rv列向量
        Lv1=(/tu*tu/2.+tu*ta/(gama-1.), -tu-ta/(gama-1.), 1./) *(gama-1.)/2./ta/ta
        Lv2=(/2*ta*ta/(gama-1.)-tu*tu, 2*tu, -2. /)*(gama-1.)/2./ta/ta
        Lv3=(/tu*tu/2-tu*ta/(gama-1.), -tu+ta/(gama-1.), 1./) *(gama-1.)/2./ta/ta
        Rv1=(/1.,1.,1./)
        Rv2=(/tu-ta,tu,tu+ta/)
        Rv3=(/tH-tu*ta,0.5*tu*tu,tH+tu*ta/)

        su1=(/ro(3,i), rou(3,i),roe(3,i)/)
        su2=(/ro(2,i+1)-ro(2,i), rou(2,i+1)-rou(2,i), roe(2,i+1)-roe(2,i)/) 
        su3=(/ro(2,i)-ro(2,i-1), rou(2,i)-rou(2,i-1), roe(2,i)-roe(2,i-1)/) 
        v1 =(/dot(Lv1,su1), dot(Lv2,su1), dot(Lv3,su1)/)  
        v2 =(/dot(Lv1,su2), dot(Lv2,su2), dot(Lv3,su2)/)  
        v3 =(/dot(Lv1,su3), dot(Lv2,su3), dot(Lv3,su3)/)    
        do k=1,3
          v1(k)=minmod(v1(k), 0.25*v2(k), 0.25*v3(k))
        enddo
        tUR=(/ dot(Rv1,v1), dot(Rv2,v1), dot(Rv3,v1) /)  !! ??????
        ro(3,i) = tUR(1)
        rou(3,i)= tUR(2)
        roe(3,i)= tUR(3)
  
        su1=(/ro(2,i), rou(2,i),roe(2,i)/)
        su2=(/ro(1,i+1)-ro(1,i), rou(1,i+1)-rou(1,i), roe(1,i+1)-roe(1,i)/) 
        su3=(/ro(1,i)-ro(1,i-1), rou(1,i)-rou(1,i-1), roe(1,i)-roe(1,i-1)/) 
        v1=(/dot(Lv1,su1), dot(Lv2,su1), dot(Lv3,su1)/)  
        v2=(/dot(Lv1,su2), dot(Lv2,su2), dot(Lv3,su2)/)  
        v3=(/dot(Lv1,su3), dot(Lv2,su3), dot(Lv3,su3)/) 
        do k=1,3
          v1(k)=minmod(v1(k), 0.5*v2(k), 0.5*v3(k))
        enddo
        tUL=(/ dot(Rv1,v1), dot(Rv2,v1), dot(Rv3,v1) /)  !! ??????
        ro(2,i) = tUL(1)
        rou(2,i)= tUL(2)
        roe(2,i)= tUL(3)
        
        roL(i) = ro(1,i) +ro(2,i) +2./3*ro(3,i)
        rouL(i)= rou(1,i)+rou(2,i)+2./3*rou(3,i)
        roeL(i)= roe(1,i)+roe(2,i)+2./3*roe(3,i)
        uL(i)  = rouL(i)/roL(i)
        pL(i)  = (roeL(i)-0.5*roL(i)*uL(i)*uL(i))*(gama-1.)
 
        roR(i-1) = ro(1,i) -ro(2,i) +2./3*ro(3,i) 
        rouR(i-1)= rou(1,i)-rou(2,i)+2./3*rou(3,i)
        roeR(i-1)= roe(1,i)-roe(2,i)+2./3*roe(3,i)
        uR(i-1)  = rouR(i-1)/roR(i-1)
        pR(i-1)  = (roeR(i-1)-0.5*roR(i-1)*uR(i-1)*uR(i-1))*(gama-1.)
       enddo
       end subroutine reconst2