checkpoint.f90 Source File


Source Code

!> @file checkpoint.f90
!> @brief Checkpoint write and read for the 1D Euler solver (MPI-correct),
!!        built on the dimension-blind `checkpoint_io` shared core (R4).
!!
!! Checkpoints are Fortran unformatted (stream) binary files that capture
!! everything needed to resume a run exactly. The on-disk format is the
!! unified format (currently version 6): the solution is gathered to rank 0
!! and written as a single file over the GLOBAL point count, so a checkpoint
!! is independent of the rank count it was written under. 1D's record list
!! is UNCHANGED since v5 (Q2d bumped only the shared header's version
!! number): a freshly WRITTEN 1D checkpoint differs from its v5 form only
!! in the stamped header version. That does NOT make existing v5 files
!! readable — `read_ckpt_header` enforces strict version equality, so v5
!! files remain rejected. Contents:
!!
!!   - Unified header: magic (42), version (6), dim (1 for this 1D solver).
!!   - Scalars: iter, t, extents (= [n_pt_global]), neq, dt.
!!   - Grid record: grid_type (character(len=32)) and a coordinate
!!     fingerprint — (sum(x), sum(x^2)) over the global node array. On
!!     restart, the fingerprint is compared against the configured grid; a
!!     mismatch (changed grid_type or node positions) is a fatal error.
!!   - The GLOBAL conserved-variable array ub(neq, n_pt_global).
!!   - An optional BDF2 previous-step array bdf2_ub_prev(neq, n_pt_global),
!!     preceded by a presence flag (integer 1 = present, 0 = absent).
!!   - (v5, F2) The Kahan `t_comp` compensation scalar, the halo width `h`
!!     (validated against the current scheme's halo width on read), and the
!!     boundary halo/ghost cells `left(neq,h)`/`right(neq,h)` — see below.
!!
!! Why the halo/ghost cells: for "prescribed" BCs (dirichlet, inflow,
!! supersonic_inlet), `apply_bcs` seeds the boundary ghost cells once at
!! `apply_initial_condition` and never rewrites them afterwards (`no_write`
!! in `boundary_conditions.f90`). But those cells are NOT actually constant
!! across a run: the explicit steppers' whole-array stage arithmetic (see
!! `time_integration.f90`) spans the halo too, so a mathematically-null
!! update still perturbs the ghost value by ~1 ulp on the very first step,
!! and that perturbed value is a stable fixed point for the rest of the run.
!! Pre-v5, checkpoints never persisted this settled ghost, so a restart
!! reloaded the stale IC-time value instead, feeding a 1-ulp-different ghost
!! into the boundary reconstruction stencil on the first post-restart step —
!! a permanent, silent divergence from an uninterrupted run. v5 closes this
!! gap by gathering/scattering the two edge ranks' halo columns alongside
!! the interior array (see `gather_edge_ghosts_to_root`/
!! `scatter_edge_ghosts_from_root` in `solution_gather.F90`), and also
!! persists `t_comp` (previously silently reset to 0 on restart, a much
!! smaller — sub-ulp-in-`t` — but real source of restart non-reproducibility
!! on adaptive-dt runs). See doc/superpowers/2026-07-08-f2-restart-diagnosis.md for the full
!! root-cause analysis.
!!
!! MPI contract:
!!   - `write_checkpoint` gathers every rank's interior slab to rank 0 via
!!     `gather_solution_to_root` and writes the single global file only on
!!     `my_rank() == 0`. Write success is agreed across ranks with
!!     `agree_outcome` (`checkpoint_io`'s collective wrapper over `par_lor`).
!!   - `read_checkpoint` opens the global file READ-ONLY on every rank
!!     (concurrent reads do not race), reads the global arrays with `iostat=`
!!     on every read, validates, then copies out this rank's slab using
!!     `decomp % i_first_global`. Restart success is agreed with
!!     `conclude_read` (one collective `par_lor`) so all ranks proceed or
!!     fail together (no restart deadlock).
!!   - All ranks must have consistent BDF2 allocation state (the same scheme
!!     applied uniformly across ranks); `write_checkpoint` additionally enforces
!!     this defensively via a collective OR (`par_lor`) before the collective
!!     BDF2 gather, so the file flag and the gather stay consistent on all ranks.
!!   - The v5 edge-ghost gather/scatter is a fixed collective position reached
!!     by every rank unconditionally on both write and read (mirroring the
!!     BDF2 gather's positioning), so it cannot desync ranks the way a
!!     rank-divergent early return would.
!!
!! Version compatibility: this is the unified v6 format (magic 42, version 6,
!! dim=1, extents=[n_pt_global]) — 1D's on-disk record list is unchanged from
!! v5 (see the module-level v5 (F2) notes above); Q2d only advanced the
!! shared header's version number, so a freshly WRITTEN 1D checkpoint differs
!! from its v5 form only in that stamped version. This does NOT make
!! existing v5 files readable: `read_ckpt_header`'s version check is strict
!! equality, so v5 — like the legacy v2/v3/v4 formats before it — is
!! rejected cleanly with a clear error (owner-authorized break; no
!! user-generated checkpoints exist). Everything else — the global gather,
!! the BDF2 record, the pointer-file mechanics, and the MPI contract above —
!! is unchanged.
!!
!! File naming: `<base>_N.bin` where N is the (width-agnostic) iteration
!! number.  A companion text file `latest_checkpoint` (written by rank 0)
!! records the path of the most recent checkpoint so that restart drivers
!! need not track the iteration themselves.
!!
!! Typical usage in the driver:
!! @code
!!   use checkpoint, only: write_checkpoint, read_checkpoint
!!   ! --- writing ---
!!   call write_checkpoint(state, cfg%checkpoint_file, t, iter)
!!   ! --- reading ---
!!   call read_checkpoint(state, cfg%restart_file, t, iter)
!! @endcode

module checkpoint
  use precision, only: wp
  use solver_state, only: solver_state_t, neq
  use logger, only: log_info, log_warn
  use mpi_runtime, only: my_rank, parallel_fatal
  use solution_gather, only: gather_solution_to_root, gather_edge_ghosts_to_root, &
                             scatter_edge_ghosts_from_root
  use parallel_reductions, only: par_lor
  use checkpoint_io, only: ckpt_filename, update_latest_pointer, resolve_latest_pointer, &
                           write_ckpt_header, read_ckpt_header, write_grid_record, &
                           read_grid_record_and_check, agree_outcome, conclude_read
  implicit none
  private
  public :: write_checkpoint, read_checkpoint

contains

  ! ---------------------------------------------------------------------------
  !> Coordinate fingerprint of the global node grid: (sum x, sum x^2).
  !! Cheap, sensitive to the node distribution, and bit-reproducible because the
  !! grid is rebuilt deterministically from the same grid_file.
  pure function grid_fingerprint(x) result(fp)
    real(wp), intent(in) :: x(:)
    real(wp) :: fp(2)
    fp(1) = sum(x)
    fp(2) = sum(x**2)
  end function grid_fingerprint

  ! ---------------------------------------------------------------------------
  !> Write a checkpoint file for the current solver state.
  !!
  !! The file is written as an unformatted stream binary.  After a successful
  !! write, the plain text file `latest_checkpoint` is updated to the new
  !! file's path.
  !!
  !! @param[in] state   Current solver state (must have ub allocated).
  !! @param[in] base    Base name (e.g. 'checkpoint'); file = base_NNNNNN.bin.
  !! @param[in] t       Current simulation time [s].
  !! @param[in] iter    Current iteration number.
  !! @param[in] t_comp  Optional Kahan compensation term for `t` (v5, F2);
  !!   absent is written as 0, matching a fresh (non-restarted) run's value.
  subroutine write_checkpoint(state, base, t, iter, is_ok, message, written_file, t_comp)
    type(solver_state_t), intent(in) :: state
    character(len=*), intent(in) :: base
    real(wp), intent(in) :: t
    integer, intent(in) :: iter
    logical, intent(out), optional :: is_ok
    character(len=*), intent(out), optional :: message
    !> Path of the file actually written (e.g. `<base>_<iter>.bin`), set only
    !! on success. Lets callers (e.g. the Cortex WRITE_CHECKPOINT reply) report
    !! the real on-disk path rather than reconstructing it from `base`.
    character(len=*), intent(out), optional :: written_file
    real(wp), intent(in), optional :: t_comp

    character(len=512) :: fname
    integer :: u, info
    integer :: has_bdf2, n_glob_buf
    integer :: h
    logical :: local_ok, lok
    logical :: local_bdf2, any_bdf2, all_bdf2
    character(len=256) :: lmsg
    real(wp) :: t_comp_val
    real(wp), allocatable :: glob(:, :)
    real(wp), allocatable :: edge_left(:, :), edge_right(:, :)

    if (present(is_ok)) is_ok = .true.
    if (present(message)) message = ''
    if (present(written_file)) written_file = ''

    ! Width-agnostic iteration field (I0 never overflows; cf. the old I6.6
    ! which produced 'base_******.bin' for iter >= 1,000,000).
    fname = ckpt_filename(base, iter)

    ! Rank 0 owns the (neq, n_pt_global) gather target; non-root passes a
    ! zero-length buffer (MPI_Gatherv ignores the recv side off-root).
    if (my_rank() == 0) then
      n_glob_buf = state % n_pt_global
    else
      n_glob_buf = 0
    end if
    allocate (glob(neq, n_glob_buf), stat=info)
    ! An allocation failure here (OOM on rank 0's global buffer) is
    ! unrecoverable; abort the whole job collectively rather than return early
    ! and strand the other ranks at the gather below.
    if (info /= 0) call parallel_fatal('checkpoint: cannot allocate global gather buffer')

    ! Gather every rank's interior slab into the global buffer on rank 0.
    call gather_solution_to_root(state % ub(:, 1:state % n_pt), state % n_pt, &
                                 state % decomp, glob)

    ! `gather_solution_to_root` (below) is collective, so all ranks must agree
    ! on whether the BDF2 gather happens. Presence is keyed off the BDF2
    ! integrator having actually been bootstrapped (`bdf2_initialized`), not the
    ! mere allocation of the array. Detect cross-rank inconsistency and fail LOUD
    ! rather than have a rank lacking the array dereference unallocated memory in
    ! the collective gather: `all_bdf2` gates BOTH the gather and the write, so
    ! every rank that enters the BDF2 gather has the array allocated.
    local_bdf2 = state % bdf2_initialized .and. allocated(state % bdf2_ub_prev)
    any_bdf2 = par_lor(local_bdf2)
    all_bdf2 = .not. par_lor(.not. local_bdf2)
    if (any_bdf2 .and. .not. all_bdf2) &
      call parallel_fatal('checkpoint: inconsistent BDF2 state across ranks')
    has_bdf2 = 0
    if (all_bdf2) has_bdf2 = 1

    local_ok = .true.
    if (my_rank() == 0) then
      open (newunit=u, file=trim(fname), status='replace', action='write', &
            form='unformatted', access='stream', iostat=info)
      if (info /= 0) then
        call log_warn('checkpoint: cannot open "'//trim(fname)//'" for writing')
        local_ok = .false.
      else
        ! Header + scalars (global point count, not the per-rank local count)
        call write_ckpt_header(u, 1, iter, t, [state % n_pt_global], neq, state % dt, info)
        if (info /= 0) local_ok = .false.
        if (local_ok) then
          ! Grid record: grid type + coordinate fingerprint, so a restart can
          ! detect a changed mesh instead of silently solving on the wrong grid.
          call write_grid_record(u, state % cfg % grid_type, &
                                 grid_fingerprint(state % mesh % x_global), info)
          ! Primary global solution
          if (info == 0) write (u, iostat=info) glob
          ! Optional BDF2 presence flag
          if (info == 0) write (u, iostat=info) has_bdf2
          if (info /= 0) local_ok = .false.
        end if
      end if
    end if

    ! Optional BDF2 previous-step array: gather it (collective, so every rank
    ! must participate) and write the global array on rank 0.
    if (all_bdf2) then
      call gather_solution_to_root(state % bdf2_ub_prev(:, 1:state % n_pt), &
                                   state % n_pt, state % decomp, glob)
      if (my_rank() == 0 .and. local_ok) then
        write (u, iostat=info) glob
        if (info /= 0) local_ok = .false.
      end if
    end if

    ! v5 (F2): Kahan t_comp + halo width + boundary halo/ghost cells. The
    ! edge-ghost gather is a FIXED collective position reached by every rank
    ! unconditionally (mirroring the BDF2 gather above), not gated on
    ! rank-0-only `local_ok`, so it cannot desync ranks.
    h = state % halo_width
    allocate (edge_left(neq, h), edge_right(neq, h), stat=info)
    if (info /= 0) call parallel_fatal('checkpoint: cannot allocate edge-ghost gather buffer')
    call gather_edge_ghosts_to_root(state, edge_left, edge_right)

    if (my_rank() == 0 .and. local_ok) then
      t_comp_val = 0.0_wp
      if (present(t_comp)) t_comp_val = t_comp
      write (u, iostat=info) t_comp_val
      if (info == 0) write (u, iostat=info) h
      if (info == 0) write (u, iostat=info) edge_left
      if (info == 0) write (u, iostat=info) edge_right
      if (info /= 0) local_ok = .false.
    end if
    deallocate (edge_left, edge_right, stat=info)
    if (info /= 0) local_ok = .false.

    if (my_rank() == 0 .and. local_ok) then
      close (u, iostat=info)
      if (info /= 0) then
        call log_warn('checkpoint: close failed for "'//trim(fname)//'"')
        local_ok = .false.
      end if
    end if

    ! W4: fold the dealloc stat into local_ok, matching the 2D twin
    ! (checkpoint_2d.f90) and the edge-buffer dealloc just above -- the 1D
    ! path previously discarded this one, the lone asymmetry in otherwise
    ! symmetric code.
    deallocate (glob, stat=info)
    if (info /= 0) local_ok = .false.

    ! Agree the outcome across ranks so every rank returns the same status.
    if (agree_outcome(local_ok, 'checkpoint: global write failed on rank 0', is_ok, message)) return

    ! The pointer-file update runs only on rank 0 but its outcome must be agreed
    ! across ranks too: otherwise rank 0 can fail to write `latest_checkpoint`
    ! and return failure while the other ranks return success, reintroducing the
    ! inconsistent-status orphan/deadlock risk. Capture the rank-0 status into
    ! locals and feed it into a second collective reduction every rank reaches.
    local_ok = .true.
    if (my_rank() == 0) then
      call log_info('checkpoint: wrote "'//trim(fname)//'"')
      call update_latest_pointer('latest_checkpoint', fname, 'checkpoint', lok, lmsg)
      if (.not. lok) local_ok = .false.
    end if
    if (agree_outcome(local_ok, 'checkpoint: failed to update latest_checkpoint pointer file', &
                      is_ok, message)) return

    if (present(written_file)) written_file = trim(fname)

  end subroutine write_checkpoint

  ! ---------------------------------------------------------------------------
  !> Read a checkpoint file and restore solver state.
  !!
  !! Validates the magic number, version, grid size, and equation count
  !! against the current state to catch mismatches early.
  !!
  !! @param[in,out] state   Solver state with ub already allocated (from setup_solver).
  !! @param[in]     fname   Path to the checkpoint file (or 'latest_checkpoint').
  !! @param[out]    t       Restored simulation time [s].
  !! @param[out]    iter    Restored iteration count.
  !! @param[out]    t_comp  Optional: restored Kahan compensation term for `t`
  !!   (v5, F2); always read from the file when present, 0 on any failure.
  subroutine read_checkpoint(state, fname, t, iter, is_ok, message, t_comp)
    type(solver_state_t), intent(inout) :: state
    character(len=*), intent(in) :: fname
    real(wp), intent(out) :: t
    integer, intent(out) :: iter
    logical, intent(out), optional :: is_ok
    character(len=*), intent(out), optional :: message
    real(wp), intent(out), optional :: t_comp

    character(len=512) :: actual_file
    integer :: u, info, i0, dstat
    integer :: has_bdf2
    integer :: h_ck
    real(wp) :: dt_ck
    real(wp) :: t_comp_val
    logical :: ok, local_ok, opened
    character(len=256) :: err
    real(wp), allocatable :: glob(:, :)
    real(wp), allocatable :: edge_left(:, :), edge_right(:, :)

    if (present(is_ok)) is_ok = .true.
    if (present(message)) message = ''
    if (present(t_comp)) t_comp = 0.0_wp

    ! Defaults so all ranks return defined values even if a peer rank fails.
    iter = 0
    t = 0.0_wp
    t_comp_val = 0.0_wp
    local_ok = .true.
    opened = .false.
    err = ''

    ! v5 (F2): allocate the edge-ghost buffers UP FRONT, sized by this rank's
    ! own expected halo width (uniform across ranks by construction — the
    ! reconstruction scheme is configured identically everywhere), so that
    ! `scatter_edge_ghosts_from_root` below is always reachable regardless of
    ! whether the sequential `local_ok`-gated reads below succeeded. Calling
    ! it unconditionally (every rank, every time) is what makes it a FIXED
    ! collective position — gating it on a value that could in principle
    ! differ per rank would risk an MPI hang, not just a wrong answer.
    allocate (edge_left(neq, state % halo_width), edge_right(neq, state % halo_width), stat=info)
    if (info /= 0) call parallel_fatal('checkpoint: cannot allocate edge-ghost restart buffer')
    edge_left = 0.0_wp
    edge_right = 0.0_wp

    ! If fname is the pointer file, read the real path from it. Every rank
    ! reads independently (read-only; concurrent reads do not race).
    actual_file = resolve_latest_pointer('latest_checkpoint', fname, 'checkpoint', ok, err)
    if (.not. ok) then
      local_ok = .false.
    else
      ! Every rank opens the single GLOBAL file READ-ONLY.
      open (newunit=u, file=trim(actual_file), status='old', action='read', &
            form='unformatted', access='stream', iostat=info)
      if (info /= 0) then
        err = 'checkpoint: cannot open restart file "'//trim(actual_file)//'"'
        local_ok = .false.
      else
        opened = .true.
      end if
    end if

    ! Header + scalars (the stored extents are the GLOBAL point count).
    if (local_ok) then
      call read_ckpt_header(u, 'checkpoint', 1, [state % n_pt_global], neq, &
                            iter, t, dt_ck, local_ok, err)
      if (local_ok) state % dt = dt_ck
    end if

    ! Grid record: reject a restart whose configured grid differs from the
    ! checkpoint's.
    if (local_ok) then
      call read_grid_record_and_check(u, 'checkpoint', state % cfg % grid_type, &
                                      grid_fingerprint(state % mesh % x_global), local_ok, err)
    end if

    ! Primary global solution: read the full (neq, n_pt_global) array, then
    ! slice out this rank's interior slab via the 1-based global index.
    if (local_ok) then
      allocate (glob(neq, state % n_pt_global), stat=info)
      if (info /= 0) then
        err = 'checkpoint: cannot allocate global restart buffer'
        local_ok = .false.
      end if
    end if
    if (local_ok) then
      read (u, iostat=info) glob
      if (info /= 0) then
        err = 'checkpoint: truncated solution'
        local_ok = .false.
      else
        i0 = state % decomp % i_first_global   ! global index of local cell 1
        state % ub(:, 1:state % n_pt) = glob(:, i0:i0 + state % n_pt - 1)
      end if
    end if

    ! Optional BDF2 previous-step (global), sliced like the primary array.
    if (local_ok) then
      read (u, iostat=info) has_bdf2
      if (info /= 0) then
        err = 'checkpoint: truncated bdf2 flag'
        local_ok = .false.
      else if (has_bdf2 == 1) then
        read (u, iostat=info) glob
        if (info /= 0) then
          err = 'checkpoint: truncated bdf2 array'
          local_ok = .false.
        else
          info = 0
          if (.not. allocated(state % bdf2_ub_prev)) &
            allocate (state % bdf2_ub_prev(neq, state % n_pt), stat=info)
          if (info /= 0) then
            err = 'checkpoint: cannot allocate bdf2 restart array'
            local_ok = .false.
          else
            state % bdf2_ub_prev(:, 1:state % n_pt) = glob(:, i0:i0 + state % n_pt - 1)
            ! BDF2 was bootstrapped when this checkpoint was written; mark it so
            ! the integrator resumes on its normal 2nd-order path rather than the
            ! backward-Euler bootstrap step.
            state % bdf2_initialized = .true.
          end if
        end if
      end if
    end if

    ! v5 (F2): Kahan t_comp, then the halo-width record (validated against
    ! this rank's current reconstruction scheme), then the boundary
    ! halo/ghost cells.
    if (local_ok) then
      read (u, iostat=info) t_comp_val
      if (info /= 0) then
        err = 'checkpoint: truncated t_comp'
        local_ok = .false.
      end if
    end if
    if (local_ok) then
      read (u, iostat=info) h_ck
      if (info /= 0) then
        err = 'checkpoint: truncated halo width'
        local_ok = .false.
      else if (h_ck /= state % halo_width) then
        err = 'checkpoint: halo width mismatch — checkpoint written with a different reconstruction scheme'
        local_ok = .false.
      end if
    end if
    if (local_ok) then
      read (u, iostat=info) edge_left
      if (info == 0) read (u, iostat=info) edge_right
      if (info /= 0) then
        err = 'checkpoint: truncated edge ghost cells'
        local_ok = .false.
      end if
    end if

    ! W4: fold the dealloc stat in (a separate `dstat`, since `info` is reused
    ! by the close immediately below), matching the 2D twin and the edge-buffer
    ! dealloc further down.
    if (allocated(glob)) then
      deallocate (glob, stat=dstat)
      if (dstat /= 0) then
        local_ok = .false.
        err = 'checkpoint: restart buffer deallocation failed'
      end if
    end if
    if (opened) then
      close (u, iostat=info)
      if (info /= 0) &
        call log_warn('checkpoint: close failed for "'//trim(actual_file)//'"')
    end if

    ! Fixed collective position reached by every rank UNCONDITIONALLY (see the
    ! allocation comment above) — delivers the edge ghosts read above (zeros,
    ! harmlessly, if a preceding read/validation failed on this rank; the
    ! collective conclude_read below aborts the whole restart in that case).
    call scatter_edge_ghosts_from_root(state, edge_left, edge_right)
    deallocate (edge_left, edge_right, stat=info)
    if (info /= 0) then
      local_ok = .false.
      err = 'checkpoint: edge-ghost buffer deallocation failed'
    end if

    if (present(t_comp) .and. local_ok) t_comp = t_comp_val

    ! Single collective decision point: every rank proceeds or fails together,
    ! so a failure on one rank cannot leave survivors marching into a collective
    ! the aborted ranks never enter (no restart deadlock).
    if (conclude_read(local_ok, err, 'checkpoint: restart failed on at least one rank', &
                      is_ok, message)) return

    if (my_rank() == 0) call log_info('checkpoint: resumed from "'//trim(actual_file)//'"')

  end subroutine read_checkpoint

end module checkpoint