๐ป 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.
FVM Solver Time Step Estimator
Core Numerical Engine in Fortran 90 โข 31 total downloads
fvm_timestep.f90
! =========================================================================
! Source File: fvm_timestep.f90
! =========================================================================
!==============================================================================
! ThermoFluidCalc โ Calculator #20 : FVM Time-Step Estimator
!==============================================================================
! Physics : In the Finite-Volume Method the time step must satisfy
!
! Convective limit : ฮt_conv โค V / ฮฃ_f |v_{n,f} ยท S_f|
! Diffusive limit : ฮt_diff โค V / (2 ยท ฮฑ ยท ฮฃ_f (S_f / d_f))
! Combined : ฮt = min(ฮt_conv, ฮt_diff) ร safety
!
! where
! V = cell volume (mยณ or mยฒ in 2-D)
! S_f = face area (mยฒ or m in 2-D)
! v_{n,f} = velocity component normal to face f (m/s)
! d_f = distance from cell centroid to neighbour centroid through f (m)
! ฮฑ = diffusivity (mยฒ/s)
! safety = user-specified safety factor (0 < safety โค 1)
!
! Modes:
! 1 = Quick single-face : ฮt = V / (|v|ยทS) ร safety
! 2 = Convective (multi) : ฮt = V / ฮฃ|v_nยทS_f| ร safety
! 3 = Diffusive (multi) : ฮt = V / (2ยทฮฑยทฮฃ(S_f/d_f)) ร safety
! 4 = Combined (multi) : min of modes 2 and 3
!
! Reference : Versteeg & Malalasekera, "An Introduction to CFD โ The FVM"
! Moukalled, Mangani & Darwiche, "The FVM in CFD", ยง7
!
! Build:
! gfortran -O2 -o fvm_timestep fvm_timestep.f90
!
! Input (stdin):
! Line 1 : mode safety
! Mode 1 line 2 : V v S
! Mode 2-4 line 2 : V nfaces alpha
! then nfaces lines : S_f v_nf d_f
!
! Output (stdout): structured key=value + per-face data
!==============================================================================
program fvm_timestep
implicit none
! --- Double precision kind -------------------------------------------------
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: MAX_FACES = 20
! --- Variables -------------------------------------------------------------
integer :: mode, nfaces, i
real(dp) :: safety, V, v_single, S_single
real(dp) :: alpha
real(dp) :: Sf(MAX_FACES), vnf(MAX_FACES), df(MAX_FACES)
real(dp) :: flux_f, sum_flux, sum_diff_ratio
real(dp) :: dt_conv, dt_diff, dt_combined
character(len=40) :: mode_name, limiting
! --- Read mode and safety --------------------------------------------------
read(*,*) mode, safety
! Clamp safety
if (safety <= 0.0_dp) safety = 0.1_dp
if (safety > 1.0_dp) safety = 1.0_dp
! โโ Mode 1 : Quick single-face โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
if (mode == 1) then
read(*,*) V, v_single, S_single
if (V <= 0.0_dp .or. S_single <= 0.0_dp) then
write(*,'(A)') 'ERROR=Volume and face area must be positive.'
stop
end if
mode_name = 'Quick Single-Face'
sum_flux = abs(v_single) * S_single
dt_conv = 0.0_dp
dt_diff = 0.0_dp
dt_combined = 0.0_dp
sum_diff_ratio = 0.0_dp
if (sum_flux > 0.0_dp) then
dt_conv = V / sum_flux * safety
end if
write(*,'(A,I2)') 'MODE=', mode
write(*,'(A,A)') 'MODE_NAME=', trim(mode_name)
write(*,'(A,ES15.8)') 'CELL_VOLUME=', V
write(*,'(A,ES15.8)') 'SAFETY=', safety
write(*,'(A,I2)') 'NFACES=', 1
write(*,'(A,ES15.8,A,ES15.8,A,ES15.8)') &
'FACE_1: S_f=', S_single, ', v_n=', v_single, ', flux=', sum_flux
write(*,'(A,ES15.8)') 'SUM_FLUX=', sum_flux
write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv
write(*,'(A)') 'DT_DIFF=N/A'
write(*,'(A,ES15.8)') 'DT_RESULT=', dt_conv
write(*,'(A)') 'LIMITING=convective'
! โโ Modes 2-4 : Multi-face โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
else if (mode >= 2 .and. mode <= 4) then
read(*,*) V, nfaces, alpha
if (V <= 0.0_dp) then
write(*,'(A)') 'ERROR=Cell volume must be positive.'
stop
end if
if (nfaces < 1 .or. nfaces > MAX_FACES) then
write(*,'(A,I3,A)') 'ERROR=Number of faces must be 1-', MAX_FACES, '.'
stop
end if
! Read face data
do i = 1, nfaces
read(*,*) Sf(i), vnf(i), df(i)
end do
! Compute sums
sum_flux = 0.0_dp
sum_diff_ratio = 0.0_dp
! Print header
select case (mode)
case (2); mode_name = 'Convective (Multi-Face)'
case (3); mode_name = 'Diffusive (Multi-Face)'
case (4); mode_name = 'Combined (Multi-Face)'
end select
write(*,'(A,I2)') 'MODE=', mode
write(*,'(A,A)') 'MODE_NAME=', trim(mode_name)
write(*,'(A,ES15.8)') 'CELL_VOLUME=', V
write(*,'(A,ES15.8)') 'ALPHA=', alpha
write(*,'(A,ES15.8)') 'SAFETY=', safety
write(*,'(A,I3)') 'NFACES=', nfaces
write(*,'(A)') 'FACES_START'
do i = 1, nfaces
flux_f = abs(vnf(i)) * Sf(i)
sum_flux = sum_flux + flux_f
if (df(i) > 0.0_dp) then
sum_diff_ratio = sum_diff_ratio + Sf(i) / df(i)
end if
write(*,'(I3,A,ES15.8,A,ES15.8,A,ES15.8,A,ES15.8)') &
i, ',', Sf(i), ',', vnf(i), ',', df(i), ',', flux_f
end do
write(*,'(A)') 'FACES_END'
write(*,'(A,ES15.8)') 'SUM_FLUX=', sum_flux
write(*,'(A,ES15.8)') 'SUM_DIFF_RATIO=', sum_diff_ratio
! Compute ฮt values
dt_conv = 0.0_dp
dt_diff = 0.0_dp
if (sum_flux > 0.0_dp) then
dt_conv = V / sum_flux * safety
end if
if (alpha > 0.0_dp .and. sum_diff_ratio > 0.0_dp) then
dt_diff = V / (2.0_dp * alpha * sum_diff_ratio) * safety
end if
select case (mode)
case (2)
write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv
write(*,'(A)') 'DT_DIFF=N/A'
dt_combined = dt_conv
limiting = 'convective'
case (3)
write(*,'(A)') 'DT_CONV=N/A'
write(*,'(A,ES15.8)') 'DT_DIFF=', dt_diff
dt_combined = dt_diff
limiting = 'diffusive'
case (4)
write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv
write(*,'(A,ES15.8)') 'DT_DIFF=', dt_diff
if (dt_conv > 0.0_dp .and. dt_diff > 0.0_dp) then
dt_combined = min(dt_conv, dt_diff)
if (abs(dt_conv - dt_diff) < 1.0e-15_dp) then
limiting = 'equal'
else if (dt_conv < dt_diff) then
limiting = 'convective'
else
limiting = 'diffusive'
end if
else if (dt_conv > 0.0_dp) then
dt_combined = dt_conv
limiting = 'convective'
else if (dt_diff > 0.0_dp) then
dt_combined = dt_diff
limiting = 'diffusive'
else
dt_combined = 0.0_dp
limiting = 'none'
end if
end select
write(*,'(A,ES15.8)') 'DT_RESULT=', dt_combined
write(*,'(A,A)') 'LIMITING=', trim(limiting)
else
write(*,'(A)') 'ERROR=Invalid mode (must be 1-4).'
stop
end if
end program fvm_timestep
Solver Description
Determine maximum allowable time step based on local convection and diffusion limits in Finite Volume grids.
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 fvm_timestep.f90 -o fvm_timestep
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
fvm_timestep < input.txt
๐ฅ Downloads & Local Files
Preview of the required input file (input.txt):
! Calculation Mode (1=Quick Single-Face, 2=Multi-Face Convective, 3=Multi-Face Diffusive, 4=Multi-Face Combined)\nSafety factor (CFL)\nCell volume [m3]\nConvective velocity [m/s]\nFace area [m2]
1
! Parameter 2
0.9
! Parameter 3
1.0
! Parameter 4
2.0
! Parameter 5
0.5
1
! Parameter 2
0.9
! Parameter 3
1.0
! Parameter 4
2.0
! Parameter 5
0.5