💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
CFL Numerical Stability Criterion
Core Numerical Engine in Fortran 90 • 22 total downloads
cfl_criterion.f90
! =========================================================================
! Source File: 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
Solver Description
Verify Courant-Friedrichs-Lewy (CFL) conditions for explicit hyperbolic solvers in 1D, 2D, and 3D.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 cfl_criterion.f90 -o cfl_criterion
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
cfl_criterion < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! 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
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0
! Parameter 11
0.0
! Parameter 12
e
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0
! Parameter 11
0.0
! Parameter 12
e