!==============================================================================
! ThermoFluidCalc — Calculator #22 : Numerical Péclet Number
!==============================================================================
! Physics : The cell Péclet number governs the ratio of convective to
!           diffusive transport within a single FVM cell:
!
!             Pe = ρ·u·L / Γ   =   u·L / α
!
!           where
!             u = characteristic velocity              (m/s)
!             L = cell length / face-to-face distance  (m)
!             α = thermal (or mass) diffusivity        (m²/s)
!             Γ = diffusion coefficient                (kg/(m·s) or W/(m·K))
!             ρ = density                              (kg/m³)
!
!           When |Pe| > 2 the Central Differencing Scheme (CDS) becomes
!           unbounded and alternative schemes (UDS, Hybrid, Power-Law,
!           QUICK) must be used.
!
! Reference : Patankar, "Numerical Heat Transfer and Fluid Flow", Ch. 5
!             Versteeg & Malalasekera, "An Introduction to CFD", Ch. 5
!
! Build:
!   gfortran -O2 -o peclet_number peclet_number.f90
!
! Modes:
!   1 = Single cell evaluation
!   2 = 1-D mesh sweep  (Pe vs cell size)
!   3 = Scheme comparison profiles
!
! Input (stdin):
!   Mode 1:  1  u  L  alpha  rho  gamma
!            (if alpha > 0 use it, else compute alpha = gamma/rho)
!   Mode 2:  2  u  alpha  rho  gamma  L_start  L_end  ncells
!   Mode 3:  3  Pe_value  npts
!==============================================================================
program peclet_number
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: MAX_CELLS = 5000

  integer  :: mode, ncells, npts, i
  real(dp) :: u, L, alpha, rho, gamma_coeff
  real(dp) :: Pe, Pe_val
  real(dp) :: L_start, L_end, dL, Lc
  real(dp) :: x, dx, phi_exact, phi_cds, phi_uds
  real(dp) :: phi_hybrid, phi_power, phi_quick
  real(dp) :: Pe_abs, w_power
  character(len=60) :: diagnosis, recommendation

  ! ── Read mode ─────────────────────────────────────────────────────────────
  read(*,*) mode

  select case (mode)

  ! ========================================================================
  ! MODE 1 : Single cell
  ! ========================================================================
  case (1)
    backspace(5)
    read(*,*) mode, u, L, alpha, rho, gamma_coeff

    ! Compute alpha if not provided
    if (alpha <= 0.0_dp) then
      if (rho > 0.0_dp .and. gamma_coeff > 0.0_dp) then
        alpha = gamma_coeff / rho
      else
        write(*,'(A)') 'ERROR=Must provide either alpha>0 or both rho>0 and gamma>0.'
        stop
      end if
    end if

    if (L <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Cell size L must be positive.'
      stop
    end if

    Pe = u * L / alpha

    ! Diagnosis
    Pe_abs = abs(Pe)
    if (Pe_abs < 0.1_dp) then
      diagnosis      = 'Pure diffusion'
      recommendation = 'CDS is ideal'
    else if (Pe_abs < 2.0_dp) then
      diagnosis      = 'Diffusion-dominated'
      recommendation = 'CDS is safe and 2nd-order accurate'
    else if (Pe_abs < 2.5_dp) then
      diagnosis      = 'Transition zone'
      recommendation = 'Consider Hybrid or QUICK scheme'
    else if (Pe_abs < 10.0_dp) then
      diagnosis      = 'Convection-dominated'
      recommendation = 'Use UDS, Hybrid, Power-Law, or QUICK'
    else
      diagnosis      = 'Strongly convective'
      recommendation = 'UDS or Power-Law; refine mesh to lower Pe'
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Single Cell'
    write(*,'(A,ES15.8)')   'VELOCITY=', u
    write(*,'(A,ES15.8)')   'CELL_SIZE=', L
    write(*,'(A,ES15.8)')   'ALPHA=', alpha
    write(*,'(A,ES15.8)')   'RHO=', rho
    write(*,'(A,ES15.8)')   'GAMMA=', gamma_coeff
    write(*,'(A,F15.6)')    'PE=', Pe
    write(*,'(A,F15.6)')    'PE_ABS=', Pe_abs
    write(*,'(A,A)')        'DIAGNOSIS=', trim(diagnosis)
    write(*,'(A,A)')        'RECOMMENDATION=', trim(recommendation)

  ! ========================================================================
  ! MODE 2 : 1-D mesh sweep
  ! ========================================================================
  case (2)
    backspace(5)
    read(*,*) mode, u, alpha, rho, gamma_coeff, L_start, L_end, ncells

    if (alpha <= 0.0_dp) then
      if (rho > 0.0_dp .and. gamma_coeff > 0.0_dp) then
        alpha = gamma_coeff / rho
      else
        write(*,'(A)') 'ERROR=Must provide either alpha>0 or both rho>0 and gamma>0.'
        stop
      end if
    end if

    if (ncells < 2) ncells = 2
    if (ncells > MAX_CELLS) ncells = MAX_CELLS
    if (L_start <= 0.0_dp .or. L_end <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Cell sizes must be positive.'
      stop
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=1-D Mesh Sweep'
    write(*,'(A,ES15.8)')   'VELOCITY=', u
    write(*,'(A,ES15.8)')   'ALPHA=', alpha
    write(*,'(A,ES15.8)')   'L_START=', L_start
    write(*,'(A,ES15.8)')   'L_END=', L_end
    write(*,'(A,I5)')       'NCELLS=', ncells

    dL = (L_end - L_start) / real(ncells - 1, dp)

    write(*,'(A)') 'DATA_START'
    do i = 0, ncells - 1
      Lc = L_start + real(i, dp) * dL
      Pe = u * Lc / alpha
      write(*,'(I5,A,ES15.8,A,F15.6)') i+1, ',', Lc, ',', Pe
    end do
    write(*,'(A)') 'DATA_END'

    ! Count flagged cells
    write(*,'(A,I5)') 'FLAGGED=', count_flagged(u, alpha, L_start, dL, ncells)

  ! ========================================================================
  ! MODE 3 : Scheme comparison profiles
  ! ========================================================================
  case (3)
    backspace(5)
    read(*,*) mode, Pe_val, npts

    if (npts < 10) npts = 10
    if (npts > MAX_CELLS) npts = MAX_CELLS

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Scheme Comparison'
    write(*,'(A,F15.6)')   'PE=', Pe_val
    write(*,'(A,I5)')      'NPTS=', npts

    Pe_abs = abs(Pe_val)
    dx = 1.0_dp / real(npts, dp)

    write(*,'(A)') 'DATA_START'
    do i = 0, npts
      x = real(i, dp) * dx

      ! Exact solution: phi = (exp(Pe*x) - 1) / (exp(Pe) - 1)
      if (abs(Pe_val) < 1.0e-12_dp) then
        phi_exact = x    ! linear when Pe→0
      else
        phi_exact = (exp(Pe_val * x) - 1.0_dp) / (exp(Pe_val) - 1.0_dp)
      end if

      ! CDS approximation (2nd order, central):
      ! For a uniform 1D grid the CDS discretisation gives the exact
      ! solution of the *discretised* equations.  The well-known
      ! wiggles appear when those discrete values overshoot/undershoot.
      ! We model the CDS discrete profile on npts cells:
      !   a_P phi_P = a_W phi_W + a_E phi_E
      !   with F = rho*u*A, D = Gamma*A/dx, Pe_cell = F/D
      ! For illustration we solve the tri-diagonal system directly.
      ! (computed below after the loop in a dedicated pass)

      ! For the per-point output we store exact only; CDS etc. come from
      ! the tri-diagonal solves.  Placeholder -999 will be replaced.
      write(*,'(F12.8,A,ES15.8,A,A,A,A,A,A,A,A)') x, ',', phi_exact, &
        ',', '0', ',', '0', ',', '0', ',', '0'
    end do
    write(*,'(A)') 'DATA_END'

    ! Now output the full discrete profiles via TDMA
    call output_discrete_profiles(Pe_val, npts)

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

contains

  !------------------------------------------------------------------------
  integer function count_flagged(u, alpha, L0, dL, n)
    real(dp), intent(in) :: u, alpha, L0, dL
    integer, intent(in)  :: n
    integer :: j
    real(dp) :: Lj, Pej
    count_flagged = 0
    do j = 0, n-1
      Lj  = L0 + real(j, dp) * dL
      Pej = abs(u * Lj / alpha)
      if (Pej > 2.0_dp) count_flagged = count_flagged + 1
    end do
  end function

  !------------------------------------------------------------------------
  ! Solve 1-D convection-diffusion with different schemes via TDMA
  ! Domain [0,1], phi(0)=0, phi(1)=1, uniform grid, npts cells
  !------------------------------------------------------------------------
  subroutine output_discrete_profiles(Pe_global, n)
    real(dp), intent(in) :: Pe_global
    integer, intent(in)  :: n
    real(dp) :: dx_l, F, D, Pe_cell
    real(dp) :: aW, aE, aP, Su, Sp
    real(dp), allocatable :: phi_cds_a(:), phi_uds_a(:)
    real(dp), allocatable :: phi_hyb_a(:), phi_pow_a(:), phi_quick_a(:)
    real(dp), allocatable :: a(:), b(:), c(:), d_rhs(:), sol(:)
    integer :: j, scheme
    real(dp) :: aWq, aEq, aWW

    allocate(phi_cds_a(0:n), phi_uds_a(0:n))
    allocate(phi_hyb_a(0:n), phi_pow_a(0:n), phi_quick_a(0:n))
    allocate(a(n), b(n), c(n), d_rhs(n), sol(n))

    dx_l = 1.0_dp / real(n, dp)
    ! F = rho*u (per unit area), D = Gamma/dx
    ! Pe_cell = F*dx/Gamma = Pe_global * dx_l  (since Pe_global = u*L_total/alpha)
    ! But here L_total = 1, so Pe_cell = Pe_global / n * n ... = Pe_global * dx_l
    ! Actually Pe_global = u * 1 / alpha, Pe_cell = u * dx_l / alpha = Pe_global * dx_l
    F = Pe_global   ! normalised: F/D_ref
    D = 1.0_dp / dx_l
    Pe_cell = F / D  ! = Pe_global * dx_l

    ! ── Loop over 5 schemes ──────────────────────────────────────────────
    do scheme = 1, 5

      ! Build TDMA coefficients for interior cells j=1..n
      do j = 1, n
        select case (scheme)
        case (1) ! CDS
          aW = D + F / 2.0_dp
          aE = D - F / 2.0_dp
        case (2) ! UDS
          aW = D + max(F, 0.0_dp)
          aE = D + max(-F, 0.0_dp)
        case (3) ! Hybrid
          aW = max(F, D + F/2.0_dp, 0.0_dp)
          aE = max(-F, D - F/2.0_dp, 0.0_dp)
        case (4) ! Power-Law
          w_power = max(0.0_dp, (1.0_dp - 0.1_dp*abs(Pe_cell))**5)
          aW = D * w_power + max(F, 0.0_dp)
          aE = D * w_power + max(-F, 0.0_dp)
        case (5) ! QUICK (deferred correction simplified as UDS base)
          aW = D + max(F, 0.0_dp)
          aE = D + max(-F, 0.0_dp)
        end select

        Su = 0.0_dp
        Sp = 0.0_dp

        ! Boundary conditions
        if (j == 1) then
          ! West boundary: phi = 0
          aP = aW + aE + (aW) - Sp  ! extra aW from boundary
          Su = Su + (2.0_dp * D + F) * 0.0_dp  ! phi_boundary = 0
          aW = 0.0_dp
          aP = aE + (2.0_dp * D + F)
        else if (j == n) then
          ! East boundary: phi = 1
          aP = aW + aE + aE - Sp
          Su = Su + (2.0_dp * D - F) * 1.0_dp
          aE = 0.0_dp
          aP = aW + (2.0_dp * D - F)
        else
          aP = aW + aE - Sp
        end if

        ! TDMA arrays (a=sub, b=diag, c=super, d=rhs)
        a(j) = -aW
        b(j) = aP
        c(j) = -aE
        d_rhs(j) = Su
      end do

      ! ── TDMA solve ────────────────────────────────────────────────────
      ! Forward sweep
      do j = 2, n
        if (abs(b(j-1)) < 1.0e-30_dp) then
          b(j-1) = 1.0e-30_dp
        end if
        c(j-1) = c(j-1) / b(j-1)
        d_rhs(j-1) = d_rhs(j-1) / b(j-1)
        b(j) = b(j) - a(j) * c(j-1)
        d_rhs(j) = d_rhs(j) - a(j) * d_rhs(j-1)
      end do
      if (abs(b(n)) < 1.0e-30_dp) b(n) = 1.0e-30_dp
      sol(n) = d_rhs(n) / b(n)
      do j = n-1, 1, -1
        sol(j) = d_rhs(j) / b(j) - c(j) * sol(j+1)
      end do

      ! Store — cell centers at (j-0.5)*dx, plus boundaries
      select case (scheme)
      case (1)
        phi_cds_a(0) = 0.0_dp; phi_cds_a(n) = 1.0_dp
        do j = 1, n; phi_cds_a(j) = sol(j); end do  ! approximate: store at j index
      case (2)
        phi_uds_a(0) = 0.0_dp; phi_uds_a(n) = 1.0_dp
        do j = 1, n; phi_uds_a(j) = sol(j); end do
      case (3)
        phi_hyb_a(0) = 0.0_dp; phi_hyb_a(n) = 1.0_dp
        do j = 1, n; phi_hyb_a(j) = sol(j); end do
      case (4)
        phi_pow_a(0) = 0.0_dp; phi_pow_a(n) = 1.0_dp
        do j = 1, n; phi_pow_a(j) = sol(j); end do
      case (5)
        phi_quick_a(0) = 0.0_dp; phi_quick_a(n) = 1.0_dp
        do j = 1, n; phi_quick_a(j) = sol(j); end do
      end select

    end do  ! scheme loop

    ! ── Output discrete profiles ────────────────────────────────────────
    write(*,'(A)') 'PROFILES_START'
    do j = 0, n
      x = real(j, dp) * dx_l
      ! Exact
      if (abs(Pe_global) < 1.0e-12_dp) then
        phi_exact = x
      else
        phi_exact = (exp(Pe_global * x) - 1.0_dp) / (exp(Pe_global) - 1.0_dp)
      end if
      write(*,'(F12.8,5(A,ES15.8))') x, &
        ',', phi_exact, &
        ',', phi_cds_a(j), &
        ',', phi_uds_a(j), &
        ',', phi_hyb_a(j), &
        ',', phi_pow_a(j)
    end do
    write(*,'(A)') 'PROFILES_END'

    deallocate(phi_cds_a, phi_uds_a, phi_hyb_a, phi_pow_a, phi_quick_a)
    deallocate(a, b, c, d_rhs, sol)

  end subroutine output_discrete_profiles

end program peclet_number