solve_pressure_correction Subroutine

   subroutine solve_pressure_correction(mesh, flow, bc, params, rhs, phi, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      real(rk), intent(in) :: rhs(:)
      real(rk), intent(inout) :: phi(:)
      type(solver_stats_t), intent(inout) :: stats

      real(rk) :: rr, rr_new
      real(rk) :: rz, rz_new
      real(rk) :: pap
      real(rk) :: alpha, beta
      integer :: iter
      integer :: c

      call ensure_projection_workspace(mesh)
      call ensure_pressure_operator_cache(mesh, flow, bc)

      associate( &
         r => projection_work%r, &
         z => projection_work%z, &
         pvec => projection_work%pvec, &
         ap => projection_work%ap, &
         local_ap => projection_work%local_ap)

         r = zero
         z = zero
         pvec = zero
         local_ap = zero

         ! The projection driver resets phi to zero before every pressure solve.
         ! Therefore A*phi is exactly zero for the initial PCG residual, so avoid
         ! one unnecessary halo exchange and pressure matvec per timestep.
         local_ap = zero

         do c = flow%first_cell, flow%last_cell
            r(c) = rhs(c)
         end do

         if (.not. pressure_cache%has_dirichlet_pressure) then
            r(1) = zero
            phi(1) = zero
         end if

         do c = flow%first_cell, flow%last_cell
            if (c == 1 .and. .not. pressure_cache%has_dirichlet_pressure) then
               z(c) = zero
            else
               z(c) = r(c) / max(pressure_cache%diag(c), tiny_safe)
            end if
         end do
          
         block
            real(rk) :: res(2)
            call flow_global_two_dots_owned(flow, r, r, r, z, res)
            rr = res(1); rz = res(2)
         end block

         pvec = z
         if (.not. pressure_cache%has_dirichlet_pressure) then
            pvec(1) = zero
         end if

         stats%pressure_iterations = 0
         stats%pressure_residual = sqrt(rr / max(real(mesh%ncells, rk), one))

         do iter = 1, params%pressure_max_iter
            if (stats%pressure_residual <= params%pressure_tol) exit

            ! a. matvec: q = A * p
            call flow_exchange_cell_scalar(flow, pvec)
            call pressure_matvec(mesh, flow, bc, pvec, local_ap)

            ! b. step length: alpha = (r, z) / (p, A*p)
            pap = flow_global_dot_owned(flow, pvec, local_ap)

            if (abs(pap) <= tiny_safe) exit

            alpha = rz / pap

            ! c. update solution and residual: phi = phi + alpha*p, r = r - alpha*q
            do c = flow%first_cell, flow%last_cell
               phi(c) = phi(c) + alpha * pvec(c)
               r(c) = r(c) - alpha * local_ap(c)
            end do

            if (.not. pressure_cache%has_dirichlet_pressure) then
               phi(1) = zero
               r(1) = zero
            end if

            do c = flow%first_cell, flow%last_cell
               if (c == 1 .and. .not. pressure_cache%has_dirichlet_pressure) then
                  z(c) = zero
               else
                  z(c) = r(c) / max(pressure_cache%diag(c), tiny_safe)
               end if
            end do

            block
               real(rk) :: res(2)
               call flow_global_two_dots_owned(flow, r, r, r, z, res)
               rr_new = res(1); rz_new = res(2)
            end block

            stats%pressure_residual = sqrt(rr_new / max(real(mesh%ncells, rk), one))
            stats%pressure_iterations = iter

            if (stats%pressure_residual <= params%pressure_tol) exit

            beta = rz_new / max(rz, tiny_safe)

            do c = flow%first_cell, flow%last_cell
               pvec(c) = z(c) + beta * pvec(c)
            end do
            if (.not. pressure_cache%has_dirichlet_pressure) then
               pvec(1) = zero
            end if

            rr = rr_new
            rz = rz_new
         end do

         if (.not. pressure_cache%has_dirichlet_pressure) then
            phi(1) = zero
         end if

         ! Accumulate pressure-solver iteration diagnostics.
         ! pressure_iterations remains the most recent solve count.
         stats%pressure_solve_count = stats%pressure_solve_count + 1
         stats%pressure_iterations_total = stats%pressure_iterations_total + stats%pressure_iterations
         stats%pressure_iterations_max = max(stats%pressure_iterations_max, stats%pressure_iterations)
         stats%pressure_iterations_avg = real(stats%pressure_iterations_total, rk) / &
                                         max(real(stats%pressure_solve_count, rk), one)

         call flow_exchange_cell_scalar(flow, phi)

      end associate
   end subroutine solve_pressure_correction