pressure_matvec Subroutine

   subroutine pressure_matvec(mesh, flow, bc, params, x, local_ax)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      real(rk), intent(in) :: x(:)
      real(rk), intent(out) :: local_ax(:)

      integer :: c, lf, nb

      call profiler_start('Pressure_Matvec')
      call ensure_pressure_operator_cache(mesh, flow, bc)

      local_ax = zero

      if (params%enable_variable_density) then
         do c = flow%first_cell, flow%last_cell
            local_ax(c) = pressure_cache%var_diag(c) * x(c)

            do lf = 1, mesh%ncell_faces(c)
               nb = pressure_cache%nb(lf, c)
               if (nb <= 0) cycle

               local_ax(c) = local_ax(c) - pressure_cache%var_coeff(lf, c) * x(nb)
            end do
         end do

         call profiler_stop('Pressure_Matvec')
         return
      end if

      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
            local_ax(c) = x(c)
            cycle
         end if

         local_ax(c) = pressure_cache%diag(c) * x(c)

         do lf = 1, mesh%ncell_faces(c)
            nb = pressure_cache%nb(lf, c)
            if (nb <= 0) cycle

            local_ax(c) = local_ax(c) - pressure_cache%coeff(lf, c) * x(nb)
         end do
      end do
      call profiler_stop('Pressure_Matvec')
   end subroutine pressure_matvec