Main Page
src/generic/Generic_CFL.f90


MODULE Generic_CFL

  ! ---- Given the conserved variables, this returns
  ! ---- the correct time step


  INTERFACE Eval_CFL
     MODULE PROCEDURE CFL_line
     MODULE PROCEDURE CFL_area
     MODULE PROCEDURE CFL_volume
  END INTERFACE



CONTAINS


  SUBROUTINE CFL_line(dt,cfl,ux,vx)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms
    USE Interp_coeffs
    USE method_parms

    ! ----  Shared Procedures
    USE Generic_Transport
    USE Generic_EvalGamma

    IMPLICIT NONE

    include 'cles.i'
    
    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
    DOUBLE PRECISION, INTENT(IN) :: dt
    DOUBLE PRECISION, INTENT(OUT) :: cfl

    DOUBLE PRECISION :: csqrd, csnd
    DOUBLE PRECISION :: rho, u,v,w,P
   
    DOUBLE PRECISION :: gamma, g1

    INTEGER :: i
    INTEGER :: bad_pressure
       
    ! ----  compute maximum CFL number

    cfl = 0.d0
    g1 = 0.0d0
    bad_pressure = 0
          
    DO i = 1,nx,1
       ! ---- get the local primative variables      
       
       rho = vx(1,i)
       u = vx(2,i)
       v = vx(3,i)
       w = vx(4,i)
       P = vx(5,i)
       
       CALL GetGamma(ux(:,i), vx(:,i), gamma)
       ! compute average gamma
       g1 = g1 + gamma
       
       ! ---- the square of the speed of sound
       csqrd = gamma*P/rho
       
       IF ( P.lt.0.0000001) THEN
          bad_pressure = 1
          CYCLE
       END IF

       csnd = dsqrt (csqrd)
       cfl = max(cfl,(dabs(u)+csnd)/dx)

    END DO

    IF (useViscous.eq.CLES_TRUE) THEN
       g1 = (g1-nx)/g1
       cfl = max(cfl,2.0d0*maxMu/(dx*dx))
       cfl = max(cfl,2.0d0*g1*maxKappa/(dx*dx))
    END IF
        
    cfl = cfl * dt

    IF (bad_pressure.eq.1) THEN
       cfl = 2.0d0
    END IF
    
    
  END SUBROUTINE CFL_line


  SUBROUTINE CFL_area(dt,cfl,ux,vx)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms
    USE Interp_coeffs
    USE method_parms

    ! ----  Shared Procedures
    USE Generic_EvalGamma
    USE Generic_Transport

    IMPLICIT NONE

    include 'cles.i'

    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
    DOUBLE PRECISION, INTENT(IN) :: dt
    DOUBLE PRECISION, INTENT(OUT) :: cfl

    DOUBLE PRECISION :: csqrd, csnd
    DOUBLE PRECISION :: rho, u,v,w,P,dl
   
    DOUBLE PRECISION :: gamma, g1
     
    INTEGER :: i,j
    INTEGER :: bad_pressure
       
    ! ----  compute maximum CFL number

    cfl = 0.0d0
    g1 = 0.0d0
    bad_pressure = 0

    
    dl = 1.d0/dsqrt(1.d0/(dx*dx) + 1.d0/(dy*dy))
    
    DO j = 1,ny,1
       
       DO i = 1,nx,1

          ! ---- get the local primative variables
          
          rho = vx(1,i,j)
          u = vx(2,i,j)
          v = vx(3,i,j)
          w = vx(4,i,j)
          P = vx(5,i,j)
          
          CALL GetGamma(ux(:,i,j), vx(:,i,j), gamma)
          g1 = g1 + gamma

          ! ---- the square of the speed of sound
          csqrd = gamma*P/rho
          
          IF ( P.lt.0.0000001) THEN
             bad_pressure = 1
             CYCLE
          END IF

          csnd = dsqrt (csqrd)
          
          cfl = max(cfl,(dabs(u)/dx+dabs(v)/dy+csnd/dl))


       ENDDO
       
    ENDDO
    
    

    IF (useViscous.eq.CLES_TRUE) THEN
       g1 = (g1-nx*ny)/g1
       cfl = max(cfl,2.0d0*maxMu*(1.0d0/(dx*dx)+1.0d0/(dy*dy)))
       cfl = max(cfl,2.0d0*g1*maxKappa*(1.0d0/(dx*dx)+1.0d0/(dy*dy)))
    END IF

    cfl = cfl * dt
    
    IF (bad_pressure.eq.1) THEN
       cfl = 2.0d0
    END IF
       
    
  END SUBROUTINE CFL_area


  SUBROUTINE CFL_volume(dt,cfl,ux,vx)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms
    USE Interp_coeffs
    USE method_parms

    ! ----  Shared Procedures
    USE Generic_EvalGamma
    USE Generic_Transport

    IMPLICIT NONE

    include 'cles.i'
    
    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
    DOUBLE PRECISION, INTENT(IN) :: dt
    DOUBLE PRECISION, INTENT(OUT) :: cfl
    
    DOUBLE PRECISION :: csqrd, csnd
    DOUBLE PRECISION :: rho, u,v,w,P,dl

    DOUBLE PRECISION :: gamma, g1

    INTEGER:: i,j,k

    INTEGER :: bad_pressure
       
    call cleslog_log_enter('CFL_volume')

    ! ----  compute maximum CFL number

    cfl = 0.0d0
    g1 = 0.0d0
    bad_pressure = 0

    dl = 1.d0/dsqrt( 1.d0/(dx*dx) + 1.d0/(dy*dy) + 1.d0/(dz*dz))

    DO k = 1,nz,1

       DO j = 1,ny,1

          DO i = 1,nx,1
             ! ---- get the local primative variables
             
             rho = vx(1,i,j,k)
             u = vx(2,i,j,k)
             v = vx(3,i,j,k)
             w = vx(4,i,j,k)
             P = vx(5,i,j,k)

             CALL GetGamma(ux(:,i,j,k), vx(:,i,j,k), gamma)
             g1 = g1 + gamma
             
             ! ---- the square of the speed of sound
             csqrd = gamma*P/rho
             
             IF ( P.lt.0.0000001) THEN
                bad_pressure = 1
                CYCLE
                
             END IF

             csnd = dsqrt (csqrd)
             cfl = max(cfl,(dabs(u)/dx+dabs(v)/dy+dabs(w)/dz+csnd/dl))

             
          ENDDO

       ENDDO

    ENDDO
    
 
    IF (useViscous.eq.CLES_TRUE) THEN
       g1 = (g1-nx*ny*nz)/g1
       cfl = max(cfl,2.0d0*maxMu*(1.0d0/(dx*dx)+1.0d0/(dy*dy)+1.0d0/(dz*dz)))
       cfl = max(cfl,2.0d0*g1*maxKappa*(1.0d0/(dx*dx)+1.0d0/(dy*dy) &
            &    +1.0d0/(dz*dz)))
    END IF

    cfl = cfl * dt
    
    IF (bad_pressure.eq.1) THEN
       cfl = 2.0d0
    END IF
    
    call cleslog_log_exit('CFL_volume')

  END SUBROUTINE CFL_volume


END MODULE Generic_CFL