output_2d.f90 Source File

!> @file output_2d.f90
!> @brief Structured 2D field writer: columns x, y, rho, u, v, p.
module output_2d
  use precision, only: wp, safe_vel
  use solver_state_2d, only: solver_state_2d_t
  use option_registry, only: method_fdm
  implicit none
  private
  public :: write_result_2d

contains

  !> Write this rank's interior cells as "x y rho u v p" (6ES24.12E3), i
  !! fastest, with a blank line after each constant-j row (gnuplot/ParaView
  !! structured grid). This is the per-rank/test-only writer; the production
  !! gathered-field writer is `solution_gather_2d.write_global_field`. The two
  !! share the same ASCII layout and column format (6ES24.12E3) — this routine
  !! writes per-rank local cells using each cell's global-index offset, while
  !! `write_global_field` writes the assembled global field on rank 0.
  subroutine write_result_2d(state, filename, is_ok, message)
    type(solver_state_2d_t), intent(in) :: state
    character(len=*), intent(in) :: filename
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message
    integer :: u, ios, i, j
    real(wp) :: x, y, rho, vu, vv, p, gam

    is_ok = .false.
    message = ''
    gam = state % cfg % gam
    open (newunit=u, file=filename, status='replace', action='write', iostat=ios)
    if (ios /= 0) then
      message = 'output_2d: cannot open "'//trim(filename)//'"'
      return
    end if

    do j = 1, state % ny_local
      do i = 1, state % nx_local
        if (trim(state % blocks(1) % method) == method_fdm) then
          ! FDM: x_node/y_node are GLOBAL arrays (build_mesh_2d_local is a no-op for
          ! uniform meshes), so index with the GLOBAL node index, mirroring the FDM
          ! path in initial_conditions_2d.f90.  Node g = ix_first_global + i - 1 lies
          ! at x_left + (g-1)*dx.  At np=1 (ix_first_global=1) this is bit-identical
          ! to x_node(i,j).
          x = state % cfg % x_left + real(state % decomp_2d % ix_first_global + i - 2, wp) * state % dx
          y = state % cfg % y_left + real(state % decomp_2d % iy_first_global + j - 2, wp) * state % dy
        else if (state % mesh % uniform) then
          x = state % cfg % x_left + (real(state % decomp_2d % ix_first_global + i - 1, wp) - 0.5_wp) * state % dx
          y = state % cfg % y_left + (real(state % decomp_2d % iy_first_global + j - 1, wp) - 0.5_wp) * state % dy
        else
          x = state % mesh % xc(i, j)
          y = state % mesh % yc(i, j)
        end if
        rho = state % ub(1, i, j)
        ! safe_vel emits 0 velocity (no division) for a vacuum/non-physical cell
        ! (rho <= 0) so output cannot raise a SIGFPE; a no-op for valid rho > 0.
        vu = safe_vel(state % ub(2, i, j), rho)
        vv = safe_vel(state % ub(3, i, j), rho)
        p = (gam - 1.0_wp) * (state % ub(4, i, j) - 0.5_wp * rho * (vu * vu + vv * vv))
        ! ES24.12E3: explicit 3-digit exponent so gfortran always emits 'E±NNN'
        ! (even for |exp| >= 100), matching the production write_global_field
        ! format so np.loadtxt parses both writers' output identically.
        write (u, '(6ES24.12E3)', iostat=ios) x, y, rho, vu, vv, p
        if (ios /= 0) then
          message = 'output_2d: write error'
          close (u, iostat=ios)
          return
        end if
      end do
      write (u, '(a)', iostat=ios) ''   ! blank line between rows
      if (ios /= 0) then
        message = 'output_2d: write error'
        close (u, iostat=ios)
        return
      end if
    end do

    close (u, iostat=ios)
    if (ios /= 0) then
      message = 'output_2d: close failed for "'//trim(filename)//'"'
      return
    end if
    is_ok = .true.
  end subroutine write_result_2d

end module output_2d