[Onim]

Thanks for visiting this thread.

I am working on coding numerical algorithm for 3D Euler equation with general coordinate. Since the flow is compressible, I need to utilize upwinding.

After I implement AUSM+ and MUSCL for upwinding/flux limiting, I tried to solve Sod shock tube problem, and got an errorneous result. Rarefaction wave and shock wave are well refined, but contact discontinuity is oscillating.

If I use 1st order upwind difference for numerical flux, I can get a proper solution. I tried to figure out any bug in MUSCL code, but it was unsuccessful.

In this case, which part of the code should be reviewed?

[duri]

You said it is 1st order scheme but using MUSCL. This is quite confusing. 1st order scheme is monotone preserving no need to use techniques like MUSCL to preserve monotonicity.
Basic first order AUSM+ is sufficient enough to produce 1st order solution for sod problem.
Could you provide the details of face flux estimation from cell averaged data that would help in fixing this issue.

[Onim]

I attach some MUSCL code for your information.

#define indcn i,j,k
#define indup i-ish,j-jsh,k-ksh
#define inddn i+ish,j+jsh,k+ksh

C     For left flux
      DO k=0,nz; DO j=0,ny; DO i=0,nx

C     Obtain gradient ratio r
c
c     Definition of r:
c           r_f= (phi_dnstream - phi_center)
c               /(phi_center - phi_upstream)
      rrho = slopeFrac(rho(indup),rho(indcn),rho(inddn))
      rpx  = slopeFrac(cpx(indup),cpx(indcn),cpx(inddn))

C     Obtain flux limiter function
      psi_rho =fluxlimit(rrho)
      psi_px  =fluxlimit(rpx)
C     Apply flux limiter and obtain final vector
      phi_left(indcn,1)   =rho(indcn)
     +                 + 0.5d0*psi_rho*(rho(indcn)-rho(indup))
      p_left(indcn,2)     =cpx(indcn)
     +                 + 0.5d0*psi_px*(cpx(indcn)-cpx(indup))
c     index 2 means that pressure in xi direction affects
c        quantity rho*u
      ENDDO;ENDDO;ENDDO

#undef indup
#undef indcn
#undef inddn


#define indup i+ish,j+jsh,k+ksh
#define indcn i,j,k
#define inddn i-ish,j-jsh,k-ksh

C     For right flux
      DO k=1,nz+1; DO j=1,ny+1; DO i=1,nx+1

      rrho = slopeFrac(rho(indup),rho(indcn),rho(inddn))


C     Obtain flux limiter function
      psirho =fluxlimit(rrho)

C     Apply flux limiter and obtain final vector
      phi_right(inddn,1)   =rho(indcn)
     +                 + 0.5d0*psirho*(rho(indcn)-rho(indup))

      ENDDO;ENDDO;ENDDO

      FUNCTION slopeFrac(upwind,center,downwind) RESULT(frac)

      nom   = downwind-center
      denom = center-upwind
c     make sure division by 0 does not happen
      if(abs(nom).lt.1d-14)then
        nom = 0.0d0
        denom = 1.0d0
      elseif(nom.gt.1d-14.and.abs(denom).lt.1d-14)then
        nom = 1d14
        denom = 1d0
      elseif(nom.lt.-1d-14.and.abs(denom).lt.1d-14)then
        nom = -1d14
        denom = 1d0
      endif
      frac = nom/denom

      FUNCTION fluxlimit(var) RESULT(fl)

      fl=MAX(0.0d0,MIN(2.0d0*var,0.5d0*(var+1.0d0),2.0d0))

[Onim]

Thank your for your reply, duri.

Well, I want to make things clear, I used both 1st order and 2nd order differencing(simple upwind scheme and Fromm's scheme) and 1st order case was successfule while 2nd order case was unsuccessful.

[duri]

It is not quite clear how you are applying limiter. It seems you are limiting the slope but comments says flux limiter. I hope you are using slope limiter. Stencil size for slope limiter will be 5, but your code says it is only 3. on left face slope need to be limited from i-2,i-1 and i-1,i on the right face it should be from i,i+1 and i+1,i+2.

[Onim]

You`re right. Because I omitted definition of indices, it is still confusing. I`ll add their definition.
By the way, I don`t know what is slope limiter. I thought flux limiter gets slope rate as its argument and the result limits numerical flux. Are slope limiter and flux limiter different?

[Onim]

I added some definitions to the code.

[Onim]

I have made my program work, and I would like to post the reason for oscillation.

At the first time, when I "extrapolate" physical quantities from the cell center(?) to cell interface, I just calculated for each of them independently. That made these physical quantities be inconsistent to each other. Especially for the case of contact discontinuity, where some quantities are continuous while others are not, so on the interface some quantities are extrapolated at high order accuracy but others are at first order . That was resource of unphysical phenomenom. After I reconstruct all physical quantities on the interface from selected primitive variables which are extrapolated with flux limiter, the problem has gone away.

Additionally, I understood the difference between slope limiter and flux limiter, but the difference does not matter that much. We may find some discussion here. http://www.cfd-online.com/Forums/mai...x-limiter.html