solve_pressure_correction Subroutine

   subroutine solve_pressure_correction(mesh, flow, bc, params, transport, 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
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: rhs(:)
      real(rk), intent(inout) :: phi(:)
      type(solver_stats_t), intent(inout) :: stats

      real(rk) :: rr, rr_new, rr_initial
      real(rk) :: rz, rz_new
      real(rk) :: pap
      real(rk) :: alpha, beta
      integer :: iter
      integer :: c
      logical :: use_zero_mean_gauge
      real(rk) :: abs_residual, rel_residual
      logical :: do_periodic_gauge

      call ensure_projection_workspace(mesh)
      call ensure_pressure_operator_cache(mesh, flow, bc)
      if (params%enable_variable_density) then
         call update_variable_pressure_operator_cache(mesh, flow, bc, transport)
      end if
      use_zero_mean_gauge = (.not. pressure_cache%has_dirichlet_pressure) .and. params%enable_variable_density

      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.
         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
            if (params%enable_variable_density) then
               call profiler_start('Pressure_Nullspace_Project')
               call subtract_owned_mean_scalar(flow, r)
               call subtract_owned_mean_scalar(flow, phi)
               call profiler_stop('Pressure_Nullspace_Project')
             else
               r(1) = zero
               phi(1) = zero
             end if
         end if

         call profiler_start('Pressure_Preconditioner')
         do c = flow%first_cell, flow%last_cell
            if (c == 1 .and. .not. pressure_cache%has_dirichlet_pressure .and. &
                .not. params%enable_variable_density) then
               z(c) = zero
            else
               z(c) = r(c) / max(pressure_preconditioner_diag(params, c), tiny_safe)
            end if
         end do
         call profiler_stop('Pressure_Preconditioner')

         if (use_zero_mean_gauge) then
            call profiler_start('Pressure_Nullspace_Project')
            call subtract_owned_mean_scalar(flow, z)
            call profiler_stop('Pressure_Nullspace_Project')
         end if
         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
            if (.not. params%enable_variable_density) then
               pvec(1) = zero
            end if
         end if

         stats%pressure_iterations = 0
         rr_initial = max(rr, tiny_safe)
         abs_residual = sqrt(max(rr, zero) / max(real(mesh%ncells, rk), one))
         rel_residual = sqrt(max(rr, zero) / rr_initial)
         stats%pressure_residual_abs = abs_residual
         stats%pressure_residual = rel_residual

         do iter = 1, params%pressure_max_iter
            if (stats%pressure_residual <= params%pressure_tol .or. &
                stats%pressure_residual_abs <= params%pressure_abs_tol) exit

            call profiler_start('Pressure_Exchange_Pvec')
            call flow_exchange_cell_scalar(flow, pvec)
            call profiler_stop('Pressure_Exchange_Pvec')
            call pressure_matvec(mesh, flow, bc, params, pvec, local_ap)

            call profiler_start('Pressure_Reduce_PAp')
            pap = flow_global_dot_owned(flow, pvec, local_ap)
            call profiler_stop('Pressure_Reduce_PAp')

            if (abs(pap) <= tiny_safe) exit

            alpha = rz / pap

            call profiler_start('Pressure_Vector_Update')
            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
            call profiler_stop('Pressure_Vector_Update')

            if (.not. pressure_cache%has_dirichlet_pressure) then
               if (params%enable_variable_density) then
                  do_periodic_gauge = (mod(iter, pcg_nullspace_projection_interval) == 0)
                  if (do_periodic_gauge) then
                     call profiler_start('Pressure_Nullspace_Project')
                     call subtract_owned_mean_scalar(flow, phi)
                     call subtract_owned_mean_scalar(flow, r)
                     call profiler_stop('Pressure_Nullspace_Project')
                  end if
               else
                  phi(1) = zero
                  r(1) = zero
               end if
            end if

            call profiler_start('Pressure_Preconditioner')
            do c = flow%first_cell, flow%last_cell
               if (c == 1 .and. .not. pressure_cache%has_dirichlet_pressure .and. &
                .not. params%enable_variable_density) then
                  z(c) = zero
               else
                  z(c) = r(c) / max(pressure_preconditioner_diag(params, c), tiny_safe)
               end if
            end do
            call profiler_stop('Pressure_Preconditioner')

            if (use_zero_mean_gauge) then
               call profiler_start('Pressure_Nullspace_Project')
               call subtract_owned_mean_scalar(flow, z)
               call profiler_stop('Pressure_Nullspace_Project')
            end if
            block
               real(rk) :: res(2)
               call profiler_start('Pressure_Reduce_RZ')
               call flow_global_two_dots_owned(flow, r, r, r, z, res)
               call profiler_stop('Pressure_Reduce_RZ')
               rr_new = res(1); rz_new = res(2)
            end block

            abs_residual = sqrt(max(rr_new, zero) / max(real(mesh%ncells, rk), one))
            rel_residual = sqrt(max(rr_new, zero) / rr_initial)
            stats%pressure_residual_abs = abs_residual
            stats%pressure_residual = rel_residual
            stats%pressure_iterations = iter

            if (stats%pressure_residual <= params%pressure_tol .or. &
                stats%pressure_residual_abs <= params%pressure_abs_tol) exit

            beta = rz_new / max(rz, tiny_safe)

            call profiler_start('Pressure_Vector_Update')
            do c = flow%first_cell, flow%last_cell
               pvec(c) = z(c) + beta * pvec(c)
            end do
            call profiler_stop('Pressure_Vector_Update')
            if (.not. pressure_cache%has_dirichlet_pressure) then
               if (.not. params%enable_variable_density) then
                  pvec(1) = zero
               end if
            end if

            rr = rr_new
            rz = rz_new
         end do

         if (.not. pressure_cache%has_dirichlet_pressure) then
            if (params%enable_variable_density) then
               call profiler_start('Pressure_Nullspace_Project')
               call subtract_owned_mean_scalar(flow, phi)
               call profiler_stop('Pressure_Nullspace_Project')
            else
               phi(1) = zero
            end if
         end if

         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 profiler_start('Pressure_Exchange_FinalPhi')
         call flow_exchange_cell_scalar(flow, phi)
         call profiler_stop('Pressure_Exchange_FinalPhi')

      end associate
   end subroutine solve_pressure_correction