🔬
Solver Purpose & Physical Scope
Verify Courant-Friedrichs-Lewy (CFL) conditions for explicit hyperbolic solvers in 1D, 2D, and 3D.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 1,072 times
📄 Source File:
cfl_criterion.f90
📁 calcul/CFD /
cfl_criterion.f90
program cfl_criterion
implicit none
double precision :: u,v,w,c_s,dx,dy,dz,dt_in,alph,nu_v,gam
integer :: scheme, ios, ndim, i, npts
double precision :: Vmag, Mach
double precision :: cfl_conv, cfl_acou, cfl_tot
double precision :: d_th, d_vi
double precision :: inv_dx2, cfl_max_scheme
double precision :: dt_cfl, dt_diff, dt_visc, dt_comb, dt_use
double precision :: Re_cell, Pe_cell, combined
character(len=30) :: regime, sch_name, stab_type
double precision :: tc, dc, vc, cc, st_val
read(*,*,iostat=ios) u; if(ios/=0)then;write(*,*)'ERROR: Invalid u.';stop;end if
read(*,*,iostat=ios) v; if(ios/=0)then;write(*,*)'ERROR: Invalid v.';stop;end if
read(*,*,iostat=ios) w; if(ios/=0)then;write(*,*)'ERROR: Invalid w.';stop;end if
read(*,*,iostat=ios) c_s; if(ios/=0)then;write(*,*)'ERROR: Invalid c_sound.';stop;end if
read(*,*,iostat=ios) dx; if(ios/=0)then;write(*,*)'ERROR: Invalid dx.';stop;end if
read(*,*,iostat=ios) dy; if(ios/=0)then;write(*,*)'ERROR: Invalid dy.';stop;end if
read(*,*,iostat=ios) dz; if(ios/=0)then;write(*,*)'ERROR: Invalid dz.';stop;end if
read(*,*,iostat=ios) dt_in; if(ios/=0)then;write(*,*)'ERROR: Invalid dt.';stop;end if
read(*,*,iostat=ios) alph; if(ios/=0)then;write(*,*)'ERROR: Invalid alpha.';stop;end if
read(*,*,iostat=ios) nu_v; if(ios/=0)then;write(*,*)'ERROR: Invalid nu.';stop;end if
read(*,*,iostat=ios) gam; if(ios/=0)then;write(*,*)'ERROR: Invalid gamma.';stop;end if
read(*,*,iostat=ios) scheme; if(ios/=0)then;write(*,*)'ERROR: Invalid scheme.';stop;end if
if(dx<=0d0)then;write(*,*)'ERROR: dx must be > 0.';stop;end if
! Determine dimensionality
ndim = 1
if (dy > 0d0) ndim = 2
if (dy > 0d0 .and. dz > 0d0) ndim = 3
if (dy <= 0d0) dy = 1d30
if (dz <= 0d0) dz = 1d30
Vmag = sqrt(u**2 + v**2 + w**2)
if (c_s > 0d0) then
Mach = Vmag / c_s
else
Mach = 0d0
c_s = 0d0
end if
! Scheme info
select case(scheme)
case(1)
sch_name = 'Euler Explicit'
cfl_max_scheme = 1.0d0
stab_type = 'Conditionally Stable'
case(2)
sch_name = 'Leapfrog'
cfl_max_scheme = 1.0d0
stab_type = 'Conditionally Stable'
case(3)
sch_name = 'Runge-Kutta 4'
cfl_max_scheme = 2.83d0
stab_type = 'Conditionally Stable'
case(4)
sch_name = 'Implicit (1st order)'
cfl_max_scheme = 1d10
stab_type = 'Unconditionally Stable'
case default
sch_name = 'Euler Explicit'
cfl_max_scheme = 1.0d0
stab_type = 'Conditionally Stable'
end select
! Inverse grid spacing squared sum
inv_dx2 = 1d0/dx**2
if (ndim >= 2) inv_dx2 = inv_dx2 + 1d0/dy**2
if (ndim >= 3) inv_dx2 = inv_dx2 + 1d0/dz**2
! Max stable time steps
! CFL-based
if (Vmag + c_s > 0d0) then
dt_cfl = cfl_max_scheme * dx / (Vmag + c_s)
else
dt_cfl = 1d30
end if
! Diffusion-based
if (alph > 0d0) then
dt_diff = 0.5d0 / (alph * inv_dx2)
else
dt_diff = 1d30
end if
! Viscous-based
if (nu_v > 0d0) then
dt_visc = 0.5d0 / (nu_v * inv_dx2)
else
dt_visc = 1d30
end if
dt_comb = min(dt_cfl, dt_diff, dt_visc)
! Determine dt to use
if (dt_in > 0d0) then
dt_use = dt_in
else
dt_use = dt_comb * 0.9d0
end if
! CFL numbers
cfl_conv = abs(u)*dt_use/dx
if (ndim >= 2) cfl_conv = cfl_conv + abs(v)*dt_use/dy
if (ndim >= 3) cfl_conv = cfl_conv + abs(w)*dt_use/dz
if (c_s > 0d0) then
cfl_acou = (abs(u)+c_s)*dt_use/dx
if (ndim >= 2) cfl_acou = cfl_acou + (abs(v)+c_s)*dt_use/dy
if (ndim >= 3) cfl_acou = cfl_acou + (abs(w)+c_s)*dt_use/dz
else
cfl_acou = cfl_conv
end if
cfl_tot = cfl_acou
! Diffusion numbers
d_th = alph * dt_use * inv_dx2
d_vi = nu_v * dt_use * inv_dx2
! Combined
combined = cfl_conv + 2d0*max(d_th, d_vi)
! Cell numbers
if (nu_v > 0d0) then
Re_cell = abs(u)*dx/nu_v
else
Re_cell = 0d0
end if
if (alph > 0d0) then
Pe_cell = abs(u)*dx/alph
else
Pe_cell = 0d0
end if
! Stability status
if (scheme == 4) then
regime = 'STABLE (Implicit)'
else if (cfl_tot <= cfl_max_scheme .and. d_th <= 0.5d0 .and. d_vi <= 0.5d0) then
regime = 'STABLE'
else
regime = 'UNSTABLE'
end if
! ---- Output ------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' CFL (COURANT-FRIEDRICHS-LEWY) CRITERION CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
write(*,'(A,F14.6,A)') ' u (x-velocity) = ', u, ' m/s'
write(*,'(A,F14.6,A)') ' v (y-velocity) = ', v, ' m/s'
write(*,'(A,F14.6,A)') ' w (z-velocity) = ', w, ' m/s'
write(*,'(A,F14.4,A)') ' c (speed of sound) = ', c_s, ' m/s'
write(*,'(A,F14.6,A)') ' |V| (magnitude) = ', Vmag, ' m/s'
write(*,'(A,F14.6)') ' Mach number = ', Mach
write(*,*)
write(*,'(A)') '--- GRID INFORMATION ----------------------------------------'
write(*,'(A,ES14.6,A)') ' dx = ', dx, ' m'
if (ndim >= 2) write(*,'(A,ES14.6,A)') ' dy = ', dy, ' m'
if (ndim >= 3) write(*,'(A,ES14.6,A)') ' dz = ', dz, ' m'
write(*,'(A,I2,A)') ' Dimensions = ', ndim, 'D'
if (ndim >= 2) write(*,'(A,F14.4)') ' Aspect ratio dy/dx = ', dy/dx
write(*,*)
write(*,'(A)') '--- SCHEME INFORMATION --------------------------------------'
write(*,'(A,A)') ' Scheme = ', trim(sch_name)
write(*,'(A,F14.4)') ' CFL_max (scheme) = ', min(cfl_max_scheme, 999d0)
write(*,'(A,A)') ' Stability type = ', trim(stab_type)
write(*,*)
write(*,'(A)') '--- TIME STEP -----------------------------------------------'
write(*,'(A,ES14.6,A)') ' dt (used) = ', dt_use, ' s'
write(*,'(A,ES14.6,A)') ' dt_max (CFL limit) = ', dt_cfl, ' s'
write(*,'(A,ES14.6,A)') ' dt_max (diffusion) = ', dt_diff, ' s'
write(*,'(A,ES14.6,A)') ' dt_max (viscous) = ', dt_visc, ' s'
write(*,'(A,ES14.6,A)') ' dt_max (combined) = ', dt_comb, ' s'
write(*,*)
write(*,'(A)') '--- CFL NUMBERS ---------------------------------------------'
write(*,'(A,F14.6)') ' CFL (convective) = ', cfl_conv
write(*,'(A,F14.6)') ' CFL (acoustic) = ', cfl_acou
write(*,'(A,F14.6)') ' CFL_max (scheme limit) = ', min(cfl_max_scheme, 999d0)
write(*,'(A,F14.2,A)') ' CFL margin = ', (cfl_max_scheme - cfl_tot)/cfl_max_scheme*100d0, ' %'
write(*,*)
write(*,'(A)') '--- DIFFUSION NUMBERS (Fourier) -----------------------------'
write(*,'(A,F14.6)') ' d_thermal = alpha*dt/dx2 = ', d_th
write(*,'(A,F14.6)') ' d_viscous = nu*dt/dx2 = ', d_vi
write(*,'(A,F14.6)') ' Limit (explicit) = 0.500000'
write(*,*)
write(*,'(A)') '--- COMBINED STABILITY --------------------------------------'
write(*,'(A,F14.6)') ' CFL + 2*d_max = ', combined
write(*,'(A,F14.6)') ' Limit = 1.000000'
write(*,*)
write(*,'(A)') '--- CELL NUMBERS --------------------------------------------'
write(*,'(A,F14.4)') ' Re_cell = |u|*dx/nu = ', Re_cell
write(*,'(A,F14.4)') ' Pe_cell = |u|*dx/alpha = ', Pe_cell
if (Re_cell > 2d0) then
write(*,'(A)') ' WARNING: Re_cell > 2, may need upwinding'
end if
write(*,*)
write(*,'(A)') '--- STABILITY STATUS ----------------------------------------'
write(*,'(A,A)') ' STATUS = ', trim(regime)
write(*,*)
! ---- Profile sweep -----------------------------------------------
write(*,'(A)') '--- PROFILE vs dt ------------------------------------------'
write(*,'(A)') ' dt CFL_conv CFL_acou d_thermal d_viscous CFL+2d Status'
write(*,'(A)') ' --------------------------------------------------------------------------------------'
npts = 40
do i = 1, npts
tc = dt_comb/10d0 + (dt_comb*3d0 - dt_comb/10d0)*dble(i-1)/dble(npts-1)
if (tc <= 0d0) tc = 1d-15
cc = abs(u)*tc/dx
if (ndim >= 2) cc = cc + abs(v)*tc/dy
if (ndim >= 3) cc = cc + abs(w)*tc/dz
if (c_s > 0d0) then
st_val = (abs(u)+c_s)*tc/dx
if (ndim >= 2) st_val = st_val + (abs(v)+c_s)*tc/dy
if (ndim >= 3) st_val = st_val + (abs(w)+c_s)*tc/dz
else
st_val = cc
end if
dc = alph * tc * inv_dx2
vc = nu_v * tc * inv_dx2
write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6,2X,A)') &
tc, cc, st_val, dc, vc, cc+2d0*max(dc,vc), &
merge('OK ','UNSTAB', (st_val<=cfl_max_scheme .and. dc<=0.5d0 .and. vc<=0.5d0))
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' CFL = |u|*dt/dx + |v|*dt/dy + |w|*dt/dz (Eq. 2.86)'
write(*,'(A)') ' CFL_acoustic = (|u|+c)*dt/dx'
write(*,'(A)') ' d = alpha*dt/dx^2 (diffusion/Fourier number)'
write(*,'(A)') ' Combined: CFL + 2d <= 1'
write(*,'(A)') ' dt_cfl = CFL_max * dx / (|u|+c)'
write(*,'(A)') ' dt_diff = 0.5 / (alpha * sum(1/dxi^2))'
write(*,'(A)') ' Re_cell = |u|*dx/nu, Pe_cell = |u|*dx/alpha'
write(*,'(A)') '============================================================'
end program cfl_criterion
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 cfl_criterion.f90 -o cfl_criterion
2. Execution with input.txt redirection:
cfl_criterion < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 e
Parameter Description:
u [m/s]\nv [m/s]\nw [m/s]\nc_i\nx [m]\ny [m]\nz [m]\nt [s]\nal_i\nu [m/s]\nga_i\nsc_i