pres_psi.f.html
Source file:   pres_psi.f
Directory:   /var/www/html/clawpack/links/shockvacuum10/shockvacuum
Converted:   Mon Aug 16 2010 at 14:13:31   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 

      subroutine pres_psi(q,pres,cs,psidx,mx,maxmx,mbc,meqn)
c
c     On input: q = cell averages 
c     Return:
c       pres = pressure in each cell
c       cs = sound speed in each cell
c       psidx = dx*averaged source term across interface at left of cell
c     The psi*dx source term will be subtracted from the pressure difference
c     in the Riemann solver to incorporate source terms into the f-wave
c     splitting.  
c
c     The formulas used here are designed to make the method well-balanced
c     for zero-velocity equilibria. 

c     For polytropic and Euler the formulas are from [RJL, 2009, to appear].

      implicit double precision (a-h,o-z)
      dimension   q(1-mbc:maxmx+mbc, meqn)
      dimension   pres(1-mbc:maxmx+mbc)
      dimension   cs(1-mbc:maxmx+mbc)
      dimension   psidx(1-mbc:maxmx+mbc)
      double precision npoly, nstar
      common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos
      common /comxt/ dtcom,dxcom,tcom

c     # parameters set in setrun.py

c     # NOTE: ieos=1 is now Euler.  For shallow water use ieos=2 with npoly=1.

      if (ieos.eq.1) then
c        # ideal gas:  
         psifactor = gamma1*grav*dxcom / (gamma*sK)
         gratio = gamma/gamma1
         do i=1-mbc,mx+mbc
            if (q(i,1).gt.0.d0) then
               pres(i) = gamma1 * (q(i,3) - 0.5d0*q(i,2)**2 / q(i,1))
               cs(i) = dsqrt(gamma*pres(i)/q(i,1))
              else
               pres(i) = 0.d0
               cs(i) = 0.d0
              endif
              if (cs(i).ne.cs(i)) then
                 write(6,*) '*** NaN in computing cs at i = ',i
                 stop 
                 endif
c           # Source term assuming equilibrium has constant entropy,
c           # so p=sK*rho**gamma as in polytropic case 

c           # Try choosing sK for each i as pres / rho**gamma:
            sK = pres(i) / q(i,1)**gamma
            psifactor = gamma1*grav*dxcom / (gamma*sK)

            if (i.gt.1 .and. i.le.mx) then
                psidx(i) = -0.5d0*sK*(-q(i,1)**gamma  
     &                       - (q(i-1,1)**gamma1 - psifactor)**gratio
     &                       + (q(i,1)**gamma1 + psifactor)**gratio
     &                       + q(i-1,1)**gamma)
              else
                psidx(i) = 0.d0
              if (psidx(i).ne.psidx(i)) then
                 write(6,*) '*** NaN in computing psidx at i = ',i
                 stop 
                 endif
              endif
            enddo

      else if (ieos.eq.2) then
c        # polytropic:
         psifactor = gamma1*grav*dxcom / (gamma*sK)
         gratio = gamma/gamma1
         do i=1-mbc,mx+mbc
            pres(i) = sK*q(i,1)**gamma
            cs(i) = dsqrt(sK*gamma*q(i,1)**gamma1)
            if (i.gt.1 .and. i.le.mx) then
                psidx(i) = -0.5d0*sK*(-q(i,1)**gamma  
     &                       - (q(i-1,1)**gamma1 - psifactor)**gratio
     &                       + (q(i,1)**gamma1 + psifactor)**gratio
     &                       + q(i-1,1)**gamma)
              else
                psidx(i) = 0.d0
              endif
            enddo
c
        endif

      return
      end