applications/euler/2d/Cyl/src/srcrx2.f

c
c
c =========================================================
      subroutine src(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
     &     aux,maux,t,dt,ibnd)
c =========================================================
      implicit double precision(a-h,o-z)
c
      dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc, 
     &     1-iby*mbc:maxmy+iby*mbc)
      dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc, 
     &     1-iby*mbc:maxmy+iby*mbc)
c
c     # source terms for cylindrical symmetry in 2d Euler equations
c     # about x-axis (so y=radius)
c
      dimension qstar(4)
      common /param/  gamma,gamma1
c
c     # 2-stage Runge-Kutta method
c
      dt2    = dt/2.d0
      press  = 0.d0
      ndim   = 2

      do 10 i=1-ibx*ibnd*mbc,mx+ibx*ibnd*mbc
         do 10 j=1-iby*ibnd*mbc,my+iby*ibnd*mbc
            y        = aux(1,i,j)
            rho      = q(1,i,j)
            u        = q(2,i,j)/q(1,i,j)
            v        = q(3,i,j)/q(1,i,j)
            press    = gamma1*(q(4,i,j) - 0.5d0*rho*(u**2 + v**2))

            if (y.eq.0.d0) write(6,*) 'y = 0 in src'
            if (rho.eq.0.d0) write(6,*) 'rho = 0 in q'

            qstar(1) = q(1,i,j) - dt2*(ndim-1)/y * q(3,i,j)
            qstar(2) = q(2,i,j) - dt2*(ndim-1)/y * 
     &           (rho*u*v)
            qstar(3) = q(3,i,j) - dt2*(ndim-1)/y * 
     &           (rho*v**2)
            qstar(4) = q(4,i,j) - dt2*(ndim-1)/y * 
     &           v*(q(4,i,j) + press)
c
c        # second stage
c
            rho      = qstar(1)
            u        = qstar(2)/qstar(1)
            v        = qstar(3)/qstar(1)
            press    = gamma1*(qstar(4) - 0.5d0*rho*(u**2 + v**2))
            if (rho.eq.0.d0) write(6,*) 'rho = 0 in qstar'

            q(1,i,j) = q(1,i,j) - dt*(ndim-1)/y * qstar(3)
            q(2,i,j) = q(2,i,j) - dt*(ndim-1)/y * 
     &           (rho*u*v)
            q(3,i,j) = q(3,i,j) - dt*(ndim-1)/y * 
     &           (rho*v**2)
            q(4,i,j) = q(4,i,j) - dt*(ndim-1)/y * 
     &           v*(qstar(4) + press)
 10   continue

      return
      end