!> @file mesh_1d.f90
!> @brief 1D grid metric layer: node coordinates and the dx/dξ Jacobian.
!!
!! The 1D solver is node-based finite-difference WENO. On a non-uniform grid the
!! physical coordinate is mapped to a uniform computational grid ξ (Δξ = 1); the
!! residual divides the flux difference by the per-node Jacobian J_i = dx/dξ
!! rather than a scalar dx. This module owns the node coordinates and J.
!!
!! Dependency rule: this module must NOT `use config` (config uses this module).
module mesh_1d
  use precision, only: wp
  use domain_decomposition, only: decomp_t
  implicit none
  private

  public :: mesh_1d_t
  public :: read_node_coords
  public :: build_mesh_global
  public :: build_mesh_local
  public :: build_mesh_cellcentered_1d
  public :: build_mesh_cellcentered_global
  public :: build_mesh_cellcentered_local

  !> Grid metrics for the 1D node-based scheme.
  !!
  !! The nodal fields (`x_global`, `jac_global`, `x_node`, `jac`) are used by
  !! the FDM path.  The cell-centered fields (`x_cell`, `dx_cell`) are used by
  !! the FVM path and are allocated by `build_mesh_cellcentered_1d`; they are
  !! unallocated (and unused) on the FDM path.
  type :: mesh_1d_t
    logical :: uniform = .true.
    integer :: n_global = 0
    real(wp) :: dx_uniform = 0.0_wp   !< constant spacing (uniform path only)
    real(wp) :: h_min = 0.0_wp        !< min node spacing (CFL)
    real(wp) :: h_left = 0.0_wp       !< first interval width (Neumann-gradient BC)
    real(wp) :: h_right = 0.0_wp      !< last  interval width (Neumann-gradient BC)
    real(wp), allocatable :: x_global(:)   !< global node coords (1:n_global)
    real(wp), allocatable :: jac_global(:) !< global dx/dξ (1:n_global)
    real(wp), allocatable :: x_node(:)     !< local node coords (1:n_local)
    real(wp), allocatable :: jac(:)        !< local dx/dξ (1:n_local)
    !> Cell-center coordinates for the FVM path.
    !! Allocated with bounds (1-h:n_cell+h) by build_mesh_cellcentered_1d.
    real(wp), allocatable :: x_cell(:)
    !> Cell widths for the FVM path (uniform: all equal dx).
    !! Allocated with bounds (1-h:n_cell+h) by build_mesh_cellcentered_1d.
    real(wp), allocatable :: dx_cell(:)
    !> Replicated GLOBAL cell-center coordinates for a NON-UNIFORM FVM grid
    !! (1:n_cell).  Populated by build_mesh_cellcentered_global from the face
    !! positions stored in x_global; sliced per rank by
    !! build_mesh_cellcentered_local.  Unallocated on the uniform FVM path,
    !! which recomputes cell centers analytically from x_left/x_right.
    real(wp), allocatable :: x_cell_global(:)
    !> Replicated GLOBAL cell widths for a NON-UNIFORM FVM grid (1:n_cell),
    !! the cell-centered twin of jac_global.  See x_cell_global.
    real(wp), allocatable :: dx_cell_global(:)
  end type mesh_1d_t

contains

  !> Read strictly-increasing node coordinates from a text file.
  !! One value per line; blank lines and lines beginning with '#' are ignored.
  subroutine read_node_coords(path, coords, n, is_ok, message)
    character(len=*), intent(in) :: path
    real(wp), allocatable, intent(out) :: coords(:)
    integer, intent(out) :: n
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: u, info, nline, k, dstat
    real(wp) :: val, prev
    character(len=512) :: line

    is_ok = .false.
    message = ''
    n = 0
    prev = 0.0_wp

    open (newunit=u, file=trim(path), status='old', action='read', iostat=info)
    if (info /= 0) then
      message = 'mesh_1d: cannot open grid_file "'//trim(path)//'"'
      return
    end if

    ! Pass 1: count data lines.
    nline = 0
    do
      read (u, '(A)', iostat=info) line
      if (info /= 0) exit
      line = adjustl(line)
      if (len_trim(line) == 0) cycle
      if (line(1:1) == '#') cycle
      nline = nline + 1
    end do

    if (nline < 2) then
      close (u, iostat=info)
      message = 'mesh_1d: grid_file must contain at least 2 node coordinates'
      return
    end if

    allocate (coords(nline), stat=info)
    if (info /= 0) error stop 'mesh_1d: coords allocation failed'

    ! Pass 2: parse and check strict monotonicity.
    rewind (u)
    k = 0
    do
      read (u, '(A)', iostat=info) line
      if (info /= 0) exit
      line = adjustl(line)
      if (len_trim(line) == 0) cycle
      if (line(1:1) == '#') cycle
      read (line, *, iostat=info) val
      if (info /= 0) then
        close (u, iostat=info)
        deallocate (coords, stat=dstat)
        message = 'mesh_1d: malformed number in grid_file'
        if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)'
        return
      end if
      k = k + 1
      if (k > 1 .and. val <= prev) then
        close (u, iostat=info)
        deallocate (coords, stat=dstat)
        write (message, '(A,I0)') 'mesh_1d: grid_file not strictly increasing at node ', k
        if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)'
        return
      end if
      coords(k) = val
      prev = val
    end do
    close (u, iostat=info)

    if (k /= nline) then
      deallocate (coords, stat=dstat)
      message = 'mesh_1d: grid_file changed during read (line-count mismatch)'
      if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)'
      return
    end if

    n = k
    is_ok = .true.
  end subroutine read_node_coords

  !> Build the replicated global mesh; for file grids, derive n_cell/endpoints.
  subroutine build_mesh_global(mesh, grid_type, grid_file, n_cell, x_left, x_right, is_ok, message)
    type(mesh_1d_t), intent(inout) :: mesh
    character(len=*), intent(in) :: grid_type, grid_file
    integer, intent(inout) :: n_cell
    real(wp), intent(inout) :: x_left, x_right
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: k, n, astat
    real(wp), allocatable :: coords(:)

    is_ok = .true.
    message = ''

    if (trim(grid_type) == 'file') then
      call read_node_coords(grid_file, coords, n, is_ok, message)
      if (.not. is_ok) return
      mesh % uniform = .false.
      mesh % n_global = n
      mesh % dx_uniform = 0.0_wp          ! no single spacing on a non-uniform grid
      call move_alloc(coords, mesh % x_global)
      n_cell = n - 1                      ! file is authoritative
      x_left = mesh % x_global(1)
      x_right = mesh % x_global(n)
    else if (trim(grid_type) == 'uniform') then
      if (n_cell <= 0) then
        is_ok = .false.
        message = 'mesh_1d: uniform grid requires n_cell > 0'
        return
      end if
      mesh % uniform = .true.
      n = n_cell + 1
      mesh % n_global = n
      mesh % dx_uniform = (x_right - x_left) / real(n_cell, wp)
      if (allocated(mesh % x_global)) then
        deallocate (mesh % x_global, stat=astat)
        if (astat /= 0) error stop 'mesh_1d: x_global deallocate failed'
      end if
      allocate (mesh % x_global(n), stat=astat)
      if (astat /= 0) error stop 'mesh_1d: x_global allocation failed'
      do k = 1, n
        mesh % x_global(k) = x_left + mesh % dx_uniform * real(k - 1, wp)
      end do
    else
      ! Reject unknown grid_type rather than silently treating it as uniform.
      ! validate_config also guards this, but build_mesh_global is called
      ! directly by tests that bypass validation.
      is_ok = .false.
      message = "mesh_1d: unknown grid_type '"//trim(grid_type)// &
                "' (expected 'uniform' or 'file')"
      return
    end if

    ! Jacobian J_k = dx/dξ on the uniform ξ grid (Δξ = 1).
    if (allocated(mesh % jac_global)) then
      deallocate (mesh % jac_global, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: jac_global deallocate failed'
    end if
    allocate (mesh % jac_global(n), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: jac_global allocation failed'
    if (mesh % uniform) then
      mesh % jac_global = mesh % dx_uniform     ! exact -> bit-for-bit with /dx
      mesh % h_min = mesh % dx_uniform
      mesh % h_left = mesh % dx_uniform
      mesh % h_right = mesh % dx_uniform
    else
      mesh % jac_global(1) = mesh % x_global(2) - mesh % x_global(1)
      mesh % jac_global(n) = mesh % x_global(n) - mesh % x_global(n - 1)
      do k = 2, n - 1
        mesh % jac_global(k) = 0.5_wp * (mesh % x_global(k + 1) - mesh % x_global(k - 1))
      end do
      mesh % h_left = mesh % x_global(2) - mesh % x_global(1)
      mesh % h_right = mesh % x_global(n) - mesh % x_global(n - 1)
      mesh % h_min = mesh % h_left
      do k = 2, n - 1
        mesh % h_min = min(mesh % h_min, mesh % x_global(k + 1) - mesh % x_global(k))
      end do
    end if
  end subroutine build_mesh_global

  !> Slice the replicated global mesh into this rank's local node arrays.
  subroutine build_mesh_local(mesh, decomp)
    type(mesh_1d_t), intent(inout) :: mesh
    type(decomp_t), intent(in) :: decomp
    integer :: ipt, g, astat

    if (.not. allocated(mesh % x_global) .or. .not. allocated(mesh % jac_global)) &
      error stop 'mesh_1d: build_mesh_local called before build_mesh_global'

    if (allocated(mesh % x_node)) then
      deallocate (mesh % x_node, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: x_node deallocate failed'
    end if
    if (allocated(mesh % jac)) then
      deallocate (mesh % jac, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: jac deallocate failed'
    end if
    allocate (mesh % x_node(decomp % n_local), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: x_node allocation failed'
    allocate (mesh % jac(decomp % n_local), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: jac allocation failed'
    do ipt = 1, decomp % n_local
      g = decomp % i_first_global - 1 + ipt
      mesh % x_node(ipt) = mesh % x_global(g)
      mesh % jac(ipt) = mesh % jac_global(g)
    end do
  end subroutine build_mesh_local

  !> Build cell-centered coordinates and widths for the FVM path.
  !!
  !! For a uniform domain `[x_left, x_right]` divided into `n_cell` cells:
  !! - `dx = (x_right - x_left) / n_cell`
  !! - interior cell centers: `x_cell(i) = x_left + (i - 0.5) * dx`  (i = 1..n_cell)
  !! - ghost centers extrapolated by `dx`: `x_cell(i) = x_left + (i - 0.5) * dx`
  !!   (the same formula applies for halo indices 1-h..0 and n_cell+1..n_cell+h)
  !! - `dx_cell(i) = dx` for all i in `1-halo_width:n_cell+halo_width`
  !!
  !! Both arrays are allocated with lower bound `1 - halo_width`.
  !!
  !! @param mesh        Mesh object to populate (x_cell and dx_cell are set).
  !! @param n_cell      Number of interior cells.
  !! @param x_left      Left domain boundary coordinate.
  !! @param x_right     Right domain boundary coordinate.
  !! @param halo_width  Number of ghost-cell layers on each side (h >= 1).
  subroutine build_mesh_cellcentered_1d(mesh, n_cell, x_left, x_right, halo_width)
    type(mesh_1d_t), intent(inout) :: mesh
    integer, intent(in) :: n_cell
    real(wp), intent(in) :: x_left, x_right
    integer, intent(in) :: halo_width

    integer :: i, lo, hi, astat
    real(wp) :: dx

    lo = 1 - halo_width
    hi = n_cell + halo_width
    dx = (x_right - x_left) / real(n_cell, wp)

    if (allocated(mesh % x_cell)) then
      deallocate (mesh % x_cell, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: x_cell deallocate failed'
    end if
    if (allocated(mesh % dx_cell)) then
      deallocate (mesh % dx_cell, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: dx_cell deallocate failed'
    end if

    allocate (mesh % x_cell(lo:hi), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: x_cell allocation failed'
    allocate (mesh % dx_cell(lo:hi), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: dx_cell allocation failed'

    do i = lo, hi
      mesh % x_cell(i) = x_left + (real(i, wp) - 0.5_wp) * dx
      mesh % dx_cell(i) = dx
    end do
  end subroutine build_mesh_cellcentered_1d

  !> Build the replicated GLOBAL cell-centered mesh for a NON-UNIFORM FVM grid.
  !!
  !! The cell-centered twin of the nodal arrays built by build_mesh_global: where
  !! that routine fills x_global/jac_global, this fills x_cell_global and
  !! dx_cell_global.  For an FVM file grid the grid_file nodes ARE the cell
  !! faces, and build_mesh_global has already loaded those n_cell+1 faces into
  !! `mesh % x_global`.  This routine derives, for each interior cell
  !! i = 1..n_cell:
  !!   - center  x_cell_global(i)  = 0.5 * (face(i) + face(i+1))
  !!   - width   dx_cell_global(i) = face(i+1) - face(i)
  !! and is the authoritative source of non-uniform cell geometry, later sliced
  !! per rank by build_mesh_cellcentered_local.
  !!
  !! It is intentionally NOT called on the uniform FVM path: there the local
  !! builder recomputes cell centers/widths analytically from x_left/x_right so
  !! the uniform result stays bit-for-bit identical to the original formula.
  !!
  !! @param mesh    Mesh object; x_global (faces) must already be populated.
  !! @param n_cell  Number of interior cells (= number of faces - 1).
  subroutine build_mesh_cellcentered_global(mesh, n_cell)
    type(mesh_1d_t), intent(inout) :: mesh
    integer, intent(in) :: n_cell

    integer :: i, astat

    if (.not. allocated(mesh % x_global)) &
      error stop 'mesh_1d: build_mesh_cellcentered_global requires x_global (cell faces)'
    if (size(mesh % x_global) < n_cell + 1) &
      error stop 'mesh_1d: build_mesh_cellcentered_global: x_global has fewer than n_cell+1 faces'

    if (allocated(mesh % x_cell_global)) then
      deallocate (mesh % x_cell_global, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: x_cell_global deallocate failed'
    end if
    if (allocated(mesh % dx_cell_global)) then
      deallocate (mesh % dx_cell_global, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: dx_cell_global deallocate failed'
    end if

    allocate (mesh % x_cell_global(n_cell), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: x_cell_global allocation failed'
    allocate (mesh % dx_cell_global(n_cell), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: dx_cell_global allocation failed'

    do i = 1, n_cell
      mesh % x_cell_global(i) = 0.5_wp * (mesh % x_global(i) + mesh % x_global(i + 1))
      mesh % dx_cell_global(i) = mesh % x_global(i + 1) - mesh % x_global(i)
    end do
  end subroutine build_mesh_cellcentered_global

  !> Build this rank's local cell-centered coordinates and widths (FVM path).
  !!
  !! The MPI twin of `build_mesh_cellcentered_1d`: where that routine builds the
  !! replicated GLOBAL cell-centered mesh, this slices out the per-rank LOCAL
  !! cell-centered mesh from the decomposition, mirroring how `build_mesh_local`
  !! slices the nodal mesh.  For local cell `i` (1..n_local) the 1-based global
  !! cell index is `decomp % i_first_global - 1 + i`; the uniform cell width is
  !! `dx = (x_right - x_left) / decomp % n_global` (the FVM decomposition counts
  !! CELLS, so `n_global` is the global cell count); the cell center is
  !! `x_left + (global_i - 0.5) * dx`.
  !!
  !! UNIFORM grids: the ghost layers `1-h .. 0` and `n_local+1 .. n_local+h` are
  !! filled with the same center formula evaluated at their (possibly
  !! out-of-domain) global cell index, so the global-edge ghosts extrapolate
  !! linearly by `dx`.  At np=1 (`i_first_global == 1`, `n_local == n_global`)
  !! this reduces exactly to
  !! `build_mesh_cellcentered_1d(mesh, n_global, x_left, x_right, halo_width)`.
  !!
  !! NON-UNIFORM grids (`.not. mesh % uniform`, i.e. an FVM file grid): interior
  !! cells are sliced from the replicated global arrays built by
  !! build_mesh_cellcentered_global — `x_cell(i) = x_cell_global(g)`,
  !! `dx_cell(i) = dx_cell_global(g)` with `g = i_first_global - 1 + i`, mirroring
  !! build_mesh_local's nodal slice.  Ghost cells beyond the global domain
  !! (`g < 1` or `g > n_global`) are filled by linear extrapolation of the edge
  !! face positions: the nearest edge cell width is held constant and the center
  !! steps out by that width.  These ghost geometry values are NOT used by the
  !! residual (which only divides interior cells `1..n_local` by their own
  !! dx_cell) nor by the uniform-coefficient reconstruction (which is
  !! geometry-agnostic); they are filled only so the arrays are consistent.
  !!
  !! @param mesh     Mesh object to populate (x_cell and dx_cell are set).
  !! @param decomp   This rank's decomposition (cell-counted for the FVM path).
  !! @param x_left   Left domain boundary coordinate.
  !! @param x_right  Right domain boundary coordinate.
  subroutine build_mesh_cellcentered_local(mesh, decomp, x_left, x_right)
    type(mesh_1d_t), intent(inout) :: mesh
    type(decomp_t), intent(in) :: decomp
    real(wp), intent(in) :: x_left, x_right

    integer :: i, g, lo, hi, h, ng, astat
    real(wp) :: dx

    h = decomp % halo_width
    lo = 1 - h
    hi = decomp % n_local + h
    ng = decomp % n_global

    if (allocated(mesh % x_cell)) then
      deallocate (mesh % x_cell, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: x_cell deallocate failed'
    end if
    if (allocated(mesh % dx_cell)) then
      deallocate (mesh % dx_cell, stat=astat)
      if (astat /= 0) error stop 'mesh_1d: dx_cell deallocate failed'
    end if

    allocate (mesh % x_cell(lo:hi), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: x_cell allocation failed'
    allocate (mesh % dx_cell(lo:hi), stat=astat)
    if (astat /= 0) error stop 'mesh_1d: dx_cell allocation failed'

    if (mesh % uniform) then
      dx = (x_right - x_left) / real(ng, wp)
      do i = lo, hi
        g = decomp % i_first_global - 1 + i
        mesh % x_cell(i) = x_left + (real(g, wp) - 0.5_wp) * dx
        mesh % dx_cell(i) = dx
      end do
    else
      if (.not. allocated(mesh % x_cell_global) .or. .not. allocated(mesh % dx_cell_global)) &
        error stop 'mesh_1d: non-uniform build_mesh_cellcentered_local requires the global cell mesh'
      do i = lo, hi
        g = decomp % i_first_global - 1 + i
        if (g >= 1 .and. g <= ng) then
          mesh % x_cell(i) = mesh % x_cell_global(g)
          mesh % dx_cell(i) = mesh % dx_cell_global(g)
        else if (g < 1) then
          ! Left ghost: extend beyond the first cell by its width.
          mesh % dx_cell(i) = mesh % dx_cell_global(1)
          mesh % x_cell(i) = mesh % x_cell_global(1) - real(1 - g, wp) * mesh % dx_cell_global(1)
        else
          ! Right ghost: extend beyond the last cell by its width.
          mesh % dx_cell(i) = mesh % dx_cell_global(ng)
          mesh % x_cell(i) = mesh % x_cell_global(ng) + real(g - ng, wp) * mesh % dx_cell_global(ng)
        end if
      end do
    end if
  end subroutine build_mesh_cellcentered_local

end module mesh_1d