FVM Solver Time Step Estimator

!==============================================================================
! 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