!==============================================================================
! ThermoFluidCalc — Calculator #32 : FLOP Estimation for CFD
!==============================================================================
! Estimates computational cost of a CFD simulation:
!   Total FLOPs = FLOP_per_cell * N_cells * N_iterations
!   Wall time   = Total_FLOPs / (FLOP_rate * N_procs * efficiency)
!
! Typical FLOP/cell/iter for common solvers:
!   Simple explicit:             50-200
!   Implicit (SIMPLE/PISO):      500-2000
!   Compressible (density-based): 200-1000
!   LES:                         1000-5000
!   DNS:                         5000-20000
!
! Reference : Gupta, §1.5
! Build: gfortran -O2 -o flop_estimation flop_estimation.f90
!==============================================================================
program flop_estimation
  implicit none
  integer, parameter :: dp = selected_real_kind(15,307)
  integer, parameter :: MX = 10000

  integer  :: mode, nprocs, npts, i
  real(dp) :: ncells, niter, flop_cell, flop_rate, eff
  real(dp) :: total_flop, wall_time, speedup
  real(dp) :: nc_min, nc_max, dnc, ncv, wt

  read(*,*) mode

  select case(mode)

  ! MODE 1 : Single estimation
  case(1)
    backspace(5)
    read(*,*) mode, ncells, niter, flop_cell, flop_rate, nprocs, eff
    if(nprocs<1) nprocs=1
    if(eff<=0.0_dp.or.eff>1.0_dp) eff=1.0_dp

    total_flop = flop_cell * ncells * niter
    wall_time  = total_flop / (flop_rate * 1.0e9_dp * real(nprocs,dp) * eff)

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Single Estimation'
    write(*,'(A,ES15.8)') 'NCELLS=', ncells
    write(*,'(A,ES15.8)') 'NITER=', niter
    write(*,'(A,ES15.8)') 'FLOP_CELL=', flop_cell
    write(*,'(A,ES15.8)') 'FLOP_RATE_GFLOPS=', flop_rate
    write(*,'(A,I8)')     'NPROCS=', nprocs
    write(*,'(A,F10.4)')  'EFFICIENCY=', eff
    write(*,'(A,ES15.8)') 'TOTAL_FLOP=', total_flop
    write(*,'(A,ES15.8)') 'WALL_TIME_S=', wall_time
    write(*,'(A,ES15.8)') 'WALL_TIME_H=', wall_time/3600.0_dp
    write(*,'(A,ES15.8)') 'WALL_TIME_D=', wall_time/86400.0_dp

    ! Time breakdown: 1-proc vs multi-proc
    write(*,'(A,ES15.8)') 'SERIAL_TIME_S=', total_flop/(flop_rate*1.0e9_dp)
    write(*,'(A,F12.4)')  'SPEEDUP=', (total_flop/(flop_rate*1.0e9_dp)) / max(wall_time,1.0e-30_dp)

  ! MODE 2 : Mesh size sweep
  case(2)
    backspace(5)
    read(*,*) mode, nc_min, nc_max, niter, flop_cell, flop_rate, nprocs, eff, npts
    if(nprocs<1) nprocs=1; if(npts<2) npts=2; if(npts>MX) npts=MX
    if(eff<=0) eff=1

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Mesh Sweep'
    write(*,'(A,ES15.8)') 'NITER=', niter
    write(*,'(A,ES15.8)') 'FLOP_CELL=', flop_cell
    write(*,'(A,I8)')     'NPROCS=', nprocs

    ! Log-spaced sweep
    write(*,'(A)') 'DATA_START'
    do i = 0, npts-1
      ncv = nc_min * (nc_max/nc_min)**(real(i,dp)/real(npts-1,dp))
      total_flop = flop_cell * ncv * niter
      wt = total_flop / (flop_rate*1.0e9_dp*real(nprocs,dp)*eff)
      write(*,'(ES12.4,A,ES12.4,A,ES12.4)') ncv, ',', total_flop, ',', wt
    end do
    write(*,'(A)') 'DATA_END'

  ! MODE 3 : Processor scaling (for fixed problem)
  case(3)
    backspace(5)
    read(*,*) mode, ncells, niter, flop_cell, flop_rate, npts
    if(npts<2) npts=2; if(npts>MX) npts=MX
    total_flop = flop_cell * ncells * niter

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Processor Scaling'
    write(*,'(A,ES15.8)') 'TOTAL_FLOP=', total_flop

    write(*,'(A)') 'DATA_START'
    do i = 0, npts-1
      nprocs = 2**i
      if(nprocs > 65536) exit
      ! Assume efficiency degrades: E = 1/(1 + 0.05*log2(N))
      eff = 1.0_dp / (1.0_dp + 0.05_dp * log(real(nprocs,dp))/log(2.0_dp))
      wt = total_flop / (flop_rate*1.0e9_dp*real(nprocs,dp)*eff)
      speedup = (total_flop/(flop_rate*1.0e9_dp)) / max(wt,1.0e-30_dp)
      write(*,'(I8,A,F10.4,A,ES12.4,A,F10.4)') nprocs, ',', eff, ',', wt, ',', speedup
    end do
    write(*,'(A)') 'DATA_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (1-3).'; stop
  end select
end program flop_estimation