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