vtf-logo

src/equations/ip1eustd.f

c -----------------------------------------------------
c Predefined internal physical boundary conditions
c for Euler equations in WENO solver
c -----------------------------------------------------

c Transformation of vector of conserved quantities
c into primitives (rho,u,0,0,p,s1,s2,dc)

c =====================================================
      SUBROUTINE it1eu(mx,meqn,q,qt)
c =====================================================
      
      IMPLICIT NONE
      
      INTEGER mx, meqn
      DOUBLE PRECISION  q(meqn,mx)
      DOUBLE PRECISION  qt(meqn,mx)
      
c     ---- Local variables
      INTEGER i, m, nvars, ierr

      call cles_getiparam('nvars', nvars, ierr)

      DO i = 1, mx 
         ! rho
         qt(1,i) = q(1,i)
         ! u, v, w
         do m=2, nvars
            qt(m,i) = q(m,i)/q(1,i)
         enddo
         ! p
         call cles_eqstate(q(1,i),meqn,qt(1,i),nvars,1,0)
         ! temperature
         qt(nvars+1,i) = q(nvars+1,i)
         ! dcflag
         qt(nvars+2,i) = 0.0
         ! all others
         DO m=nvars+3, meqn
            qt(m,i) = q(m,i)
         END Do
      END DO
      
      RETURN
      END

c -----------------------------------------------------
c Construction of reflective boundary conditions from
c mirrored primitive values and application in
c conservative form in local patch in 1D
c -----------------------------------------------------

c =====================================================
      SUBROUTINE ip1eurfl(q,mx,lb,ub,meqn,nc,idx, 
     $     qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================

      IMPLICIT NONE
     
      INTEGER mx, meqn, maux, nc, idx(1,nc), lb(1), ub(1)
      DOUBLE PRECISION xc(1,nc), 
     $     phi(nc), vn(1,nc), auex(maux,nc), dx(1), time
      DOUBLE PRECISION q(meqn, mx)
      DOUBLE PRECISION qex(meqn,nc)
      
c     ---- Local variables
      INTEGER i, n, m, stride, getindx, nvars, useViscous, ierr
      DOUBLE PRECISION u, ul
      
      call cles_getiparam('nvars', nvars, ierr)
      call cles_getiparam('useviscous', useViscous, ierr)

      stride = (ub(1) - lb(1))/(mx-1)
      
      DO n = 1, nc
         i = getindx(idx(1,n), lb(1), stride)
         
         u = -qex(2,n)
c     ---- Add boundary velocities if available
         if (maux.ge.1) u = u + auex(1,n)
         u = 2.d0*u
         
c      ---- Invert entire velocity vector for Navier-Stokes
         IF (useViscous.eq.1) THEN
            qex(2,n) = qex(2,n) + u
c      ---- Invert only normal velocity vector for Euler
         ELSE
            ul = u*vn(1,n)
            qex(2,n) = qex(2,n) + ul*vn(1,n) 
         ENDIF

         q(1,i) = qex(1,n)
         do m=2, nvars
            q(m,i) = qex(m,n)*qex(1,n)
         enddo
         call cles_inveqst(q(1,i),meqn,qex(1,n),nvars,1,0)
         ! temperature
         q(nvars+1,i) = qex(nvars+1,n) 
         do m=nvars+3, meqn  ! skip dcflag
            q(m,i) = qex(m,n)
         enddo
         
      END DO
      
      RETURN
      END
      
c -----------------------------------------------------
c Injection of conservative extrapolated values in local patch
c -----------------------------------------------------

c =====================================================
      SUBROUTINE ip1euex(q,mx,lb,ub,meqn,nc,idx, 
     $     qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
      
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) mx, meqn, maux, nc, idx(1,nc), lb(1), ub(1)
      DOUBLE PRECISION xc(1,nc), 
     $     phi(nc), vn(1,nc), auex(maux,nc), dx(1), time
      DOUBLE PRECISION q(meqn, mx)
      DOUBLE PRECISION qex(meqn,nc)
      
c     ---- Local variables
      INTEGER i, n, m, stride, getindx, nvars, ierr
      DOUBLE PRECISION rho, u, p

      call cles_getiparam('nvars', nvars, ierr)

      stride = (ub(1) - lb(1))/(mx-1)

      DO n = 1, nc
         i = getindx(idx(1,n), lb(1), stride)

         ! rho
         q(1,i) = qex(1,n)
         ! rho (u,v,w)
         do m=2, nvars
            q(m,i) = qex(m,n)*qex(1,n)
         enddo
         ! E 
         call cles_inveqst(q(1,i),meqn,qex(1,n),nvars,1,0)
         ! temperature
         q(nvars+1,i) = qex(nvars+1,n) 
         do m=nvars+3, meqn  ! skip dcflag
            q(m,i) = qex(m,n)
         enddo
      END DO
      
      RETURN
      END

<