bw.f

四十三种差分格式源代码,这只是计算用的

           endif
           a = sqrt(gamma*p/u1(i))
           rho(i) = max(abs(u+a),abs(u-a))
           rhomax = max(rhomax,rho(i))
 70     continue

        tt = tt + delta_t
        if(tstep.eq.1.and.lambda*rhomax.gt.1.) then
          write(*,*) 'WARNING: CFL condition violated'
          write(*,*) '# time steps = ', j
          write(*,*) 'Maximum CFL number = ', lambda*rhomax
        elseif(tstep.eq.2.and.tt.lt.t) then
          delta_t = cfl*delta_x/rhomax
          if(tt+delta_t.ge.t) then
            delta_t = t - tt
          endif
          lambda = delta_t/delta_x
        elseif(tstep.eq.2.and.tt.ge.t) then
          goto 200
        endif

 100  continue
      
 200  write(*,*) 'Calculation complete.'
      write(*,*) 'Final time = ', t 
      if(tstep.eq.2) write(*,*) '# iterations =  ', j
      if(tstep.eq.1) write(*,*) 'Final CFL number = ' , lambda*rhomax

      open (unit = 13, file = 'pressure.out')
      open (unit = 14, file = 'velocity.out')
      open (unit = 15, file = 'sound.out')
      open (unit = 16, file = 'density.out')
      open (unit = 17, file = 'entropy.out')
      open (unit = 18, file = 'mach.out')
      open (unit = 19, file = 'massflux.out')
      open (unit = 20, file = 'spectral.out')

c...Report results.
      do 110, i=1,N+1
        x = aa + (bb-aa)*real(i-1)/real(N)
        p = (gamma-1.)*(u3(i) - .5*u2(i)*u2(i)/u1(i))
        a = sqrt(gamma*p/u1(i))
        u = u2(i)/u1(i)
        write(13,*) x, p
        write(14,*) x, u
        write(15,*) x, a
        write(16,*) x, u1(i)
        write(17,*) x, cu*log(p)-cp*log(u1(i))
        write(18,*) x, u/a
        write(19,*) x, u2(i)
        write(20,*) x, lambda*rho(i)
 110  continue
        
      close(unit=13)
      close(unit=14)
      close(unit=15)
      close(unit=16)
      close(unit=17)
      close(unit=18)
      close(unit=19)
      close(unit=20)
       
      stop
      end

      subroutine sw(rl,rul,retl,rr,rur,retr,fp1,fp2,fp3,fm1,fm2,fm3)

c...Steger-Warming flux vector splitting.

      real*8 gamma, rl, rul, retl, rr, rur, retr
      real*8 rhol, rhor, ul, ur, pl, pr, hl, hr, al, ar
      real*8 lambda1p, lambda2p, lambda3p, lambda1m, lambda2m, lambda3m
      real*8 fp1, fp2, fp3, fm1, fm2, fm3

      parameter(gamma=1.4)

c...Convert conservative variables to primitive variables.
      rhol = rl
      rhor = rr
      ul = rul/rhol
      ur = rur/rhor
      pl = (gamma-1.)*(retl - .5*rul*rul/rhol)
      pr = (gamma-1.)*(retr - .5*rur*rur/rhor)
      hl = (retl+pl)/rhol
      hr = (retr+pr)/rhor
      al = sqrt(gamma*pl/rhol)
      ar = sqrt(gamma*pr/rhor)

c...Compute wave speed splitting.

      lambda1p = max(0.,ul)
      lambda2p = max(0.,ul+al)
      lambda3p = max(0.,ul-al)

      lambda1m = min(0.,ur)
      lambda2m = min(0.,ur+ar)
      lambda3m = min(0.,ur-ar)

c...Compute flux splitting.
      fp1 = .5*rhol*(2.*(gamma-1.)*lambda1p+lambda2p+lambda3p)/gamma
      fp2 = .5*rhol*(2.*(gamma-1.)*ul*lambda1p+(ul+al)*lambda2p
     *                + (ul-al)*lambda3p)/gamma
      fp3 = .5*(gamma-1.)*rhol*lambda1p*ul*ul/gamma
     *       + .5*rhol*((hl+al*ul)*lambda2p+(hl-al*ul)*lambda3p)/gamma

      fm1 = .5*rhor*(2.*(gamma-1.)*lambda1m+lambda2m+lambda3m)/gamma
      fm2 = .5*rhor*(2.*(gamma-1.)*ur*lambda1m+(ur+ar)*lambda2m
     *                + (ur-ar)*lambda3m)/gamma
      fm3 = .5*(gamma-1.)*rhor*lambda1m*ur*ur/gamma
     *       + .5*rhor*((hr+ar*ur)*lambda2m+(hr-ar*ur)*lambda3m)/gamma

      return
      end

      subroutine vl(rl,rul,retl,rr,rur,retr,fp1,fp2,fp3,fm1,fm2,fm3)

c...Van Leer flux vector splitting.

      real*8 gamma, rl, rul, retl, rr, rur, retr
      real*8 rhol, rhor, ul, ur, pl, pr, hl, hr, al, ar, Ml, Mr
      real*8 Mp, Mm, tp, tm, fp1, fp2, fp3, fm1, fm2, fm3

      parameter(gamma=1.4)

c...Convert conservative variables to primitive variables.
      rhol = rl
      rhor = rr
      ul = rul/rhol
      ur = rur/rhor
      pl = (gamma-1.)*(retl - .5*rul*rul/rhol)
      pr = (gamma-1.)*(retr - .5*rur*rur/rhor)
      hl = (retl+pl)/rhol
      hr = (retr+pr)/rhor
      al = sqrt(gamma*pl/rhol)
      ar = sqrt(gamma*pr/rhor)
      Ml = ul/al
      Mr = ur/ar

c...Compute positive flux splitting.

      if(Ml.le.-1.) then
        fp1 = 0.
        fp2 = 0.
        fp3 = 0.
      elseif(Ml.lt.1.) then
        Mp = .25*(Ml+1.)*(Ml+1.)
        tp = (gamma-1.)*ul+2.*al
        fp1 = rhol*al*Mp
        fp2 = fp1*tp/gamma
        fp3 = .5*fp1*tp*tp/(gamma*gamma-1.)
      else
        fp1 = rul
        fp2 = rul*ul + pl
        fp3 = rhol*hl*ul
      endif

c...Compute negative flux splitting.

      if(Mr.le.-1.) then
        fm1 = rur
        fm2 = rur*ur + pr
        fm3 = rhor*hr*ur
      elseif(Mr.lt.1.) then
        Mm = -.25*(Mr-1.)*(Mr-1.)
        tm = (gamma-1.)*ur-2.*ar
        fm1 = rhor*ar*Mm
        fm2 = fm1*tm/gamma
        fm3 = .5*fm1*tm*tm/(gamma*gamma-1.)
      else
        fm1 = 0.
        fm2 = 0.
        fm3 = 0.
      endif

      return
      end

      subroutine zb(rl,rul,retl,rr,rur,retr,fp1,fp2,fp3,fm1,fm2,fm3)

c...Zha-Bilgen flux vector splitting.

      real*8 gamma, rl, rul, retl, rr, rur, retr, fp1, fp2, fp3
      real*8 rhol, rhor, ul, ur, pl, pr, hl, hr, al, ar, Ml, Mr
      real*8 pp, pm, pap, pam, fm1, fm2, fm3

      parameter(gamma=1.4)

c...Convert conservative variables to primitive variables.
      rhol = rl
      rhor = rr
      ul = rul/rhol
      ur = rur/rhor
      pl = (gamma-1.)*(retl - .5*rul*rul/rhol)
      pr = (gamma-1.)*(retr - .5*rur*rur/rhor)
      hl = (retl+pl)/rhol
      hr = (retr+pr)/rhor
      al = sqrt(gamma*pl/rhol)
      ar = sqrt(gamma*pr/rhor)
      Ml = ul/al
      Mr = ur/ar

c...Compute positive splitting of p.
      if(Ml.le.-1.) then
        pp = 0.
        pap = 0.
      elseif(Ml.lt.1.) then
        pp = .5*(1.+Ml)*pl
        pap = .5*(ul+al)*pl
      else
        pp = pl
        pap = pl*ul
      endif

c...Compute negative splitting of M and p.
      if(Mr.le.-1.) then
        pm = pr
        pam = pr*ur
      elseif(Mr.lt.1.) then
        pm = .5*(1.-Mr)*pr
        pam =.5*(ur-ar)*pr
      else
        pm = 0.
        pam = 0.
      endif

c...Compute conserative numerical fluxes.

      fp1 = max(0.,ul)*rhol
      fm1 = min(0.,ur)*rhor
      fp2 = max(0.,ul)*rul + pp
      fm2 = min(0.,ur)*rur + pm
      fp3 = max(0.,ul)*retl + pap
      fm3 = min(0.,ur)*retr + pam

      return
      end