ensure_pressure_operator_cache Subroutine

   subroutine ensure_pressure_operator_cache(mesh, flow, bc)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      integer :: ierr, c, lf, f, nb
      real(rk) :: dist, coeff
      integer :: local_dirichlet_flag, global_dirichlet_flag

      if (pressure_cache%initialized) then
         if (pressure_cache%ncells == mesh%ncells) return
         call finalize_flow_projection_workspace()
      end if

      pressure_cache%ncells = mesh%ncells
      pressure_cache%max_faces = size(mesh%cell_faces, 1)

      allocate(pressure_cache%nb(pressure_cache%max_faces, mesh%ncells))
      allocate(pressure_cache%coeff(pressure_cache%max_faces, mesh%ncells))
      allocate(pressure_cache%diag(mesh%ncells))
      allocate(pressure_cache%var_coeff(pressure_cache%max_faces, mesh%ncells))
      allocate(pressure_cache%var_diag(mesh%ncells))

      pressure_cache%nb = 0
      pressure_cache%coeff = zero
      pressure_cache%diag = zero
      pressure_cache%var_coeff = zero
      pressure_cache%var_diag = zero
      pressure_cache%has_neumann_outlet = .false.

      local_dirichlet_flag = 0

      do c = 1, mesh%ncells
         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)
            nb = face_effective_neighbor(mesh, bc, f, c)

            if (nb <= 0) then
               if (patch_type_for_face(mesh, bc, f) == bc_neumann) then
                  pressure_cache%has_neumann_outlet = .true.
               end if

               if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
                  dist = face_normal_distance(mesh, bc, f, c, 0)
                  coeff = mesh%faces(f)%area / dist
                  pressure_cache%diag(c) = pressure_cache%diag(c) + coeff
                  local_dirichlet_flag = 1
               end if
               cycle
            end if

            dist = face_normal_distance(mesh, bc, f, c, nb)
            coeff = mesh%faces(f)%area / dist

            pressure_cache%nb(lf, c) = nb
            pressure_cache%coeff(lf, c) = coeff
            pressure_cache%diag(c) = pressure_cache%diag(c) + coeff
         end do

         if (pressure_cache%diag(c) <= tiny_safe) then
            call fatal_error('numerics', 'non-positive cached pressure diagonal')
         end if
      end do

      call MPI_Allreduce(local_dirichlet_flag, global_dirichlet_flag, 1, MPI_INTEGER, MPI_MAX, flow%comm, ierr)
      pressure_cache%has_dirichlet_pressure = (global_dirichlet_flag > 0)

      pressure_cache%initialized = .true.
   end subroutine ensure_pressure_operator_cache