Projection method FVM poisson part, adding source term

Question

The idea of the method is to decompose the Navier-Stokes equation into the solenoidal and irrotational parts.

$$\frac{\partial u}{\partial t}+u(\nabla \cdot u)=-\frac{1}{\rho}\nabla p+\nabla ^2 u$$ Discretizing the transient term and initializing the intermediate velocity $u^*$.
$$\frac{(u^{n+1}-u^*) + (u^*-u^n)}{\Delta t}+u(\nabla \cdot u)=-\frac{1}{\rho}\nabla p+\nabla ^2 u$$ Decomposing the equation

1. $$\frac{(u^*-u^n)}{\Delta t}+u(\nabla \cdot u)=+\nabla ^2 u$$
2. $$\frac{(u^{n+1}-u^*)}{\Delta t}=-\frac{1}{\rho}\nabla p$$ expressing the $u^{n+1}$ from (2) and plugging into continuity equation leads to
3. $$\nabla^2p=\frac{\rho}{\Delta t}(\nabla \cdot u^{n+1})$$

Algorithm:

1. Solve equation (1) to find $u^*$
2. Solve equation (3) to find the $p$
3. Solve equation (2) to find $u^{n+1}$

I was struggling with the curved no-slip boundaries in step 1 of the algorithm described here. Keep the approach of nullifying the boundary cells for now.

The new issue arose when solving the Poisson equation from step 2 of the algorithm. The discretization of Poisson equation: $$\nabla\cdot(\nabla p)=\frac{\rho}{\Delta t}(\nabla \cdot u^{n+1})$$

Integrating over the finite volume $$\int\nabla\cdot(\nabla p)dV=\frac{\rho}{\Delta t}\int(\nabla \cdot u^{n+1})dV$$

Applying Gauss divergence theorem $$\int(\nabla p)\cdot \hat ndS=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{n+1})dS$$

Decomposing the left part into normal and tangential flux:

$$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{n+1})dS$$

$\delta_f$ is the normal distance between the cell and its neighbour through the face $f$
$p_{a(f)},{p_{b(f)}}$ are the nodal value of $p$ at the face $f$
$A_f$ is the length of face $f$
$l_f$ is the vector pointing from centre of the cell into its neighbours centre
$\hat{t_f}$ is the tangential unit vector

RHS discretization of the source term:

$$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{*})dS$$ $$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\sum_{f=1}^{N_{f,O}} (u^*\cdot\hat n)A_f$$

Rearranging (considering triangular mesh):

$$\bigg( \frac{A_1}{\delta_1}+\frac{A_2}{\delta_2}+\frac{A_3}{\delta_3}\bigg)p_O - \frac{A_1}{\delta_1}p_1 - \frac{A_2}{\delta_2}p_2 - \frac{A_3}{\delta_3}p_3=$$ $$=-\frac{\rho}{\Delta t}\sum_{f=1}^{N_{f,O}} (u^*\cdot\hat n)A_f + \bigg( \frac{p_a-p_b}{\delta_1|t_1|}(\hat t_1 \cdot l_1)\bigg)A_1+\bigg( \frac{p_b-p_c}{\delta_2|t_2|}(\hat t_2 \cdot l_2)\bigg)A_2 + \bigg( \frac{p_c-p_a}{\delta_3|t_3|}(\hat t_3 \cdot l_3)\bigg)A_3$$

The RHS of the equation is named SCSKEW

Used Gauss Seidel algorithm

Running the code without the source term gives perfect results (left boundary 1, right Neumann condition = 0): Adding the source term (set $u=(1,0)$ at the left faces):

What could be wrong with the discretization or the implementation of it into the code?
As it is expected to slightly stretch the pressure field in the direction that it is pointing to, not move the pressure and increase it by 16 times from what is currently being obtained.
(the array is named source and is added to the scskew as they both are on the RHS of the discretized equation)? (running the code on refined mesh, playing with the values of $\Delta t$ did not affect the result)

Can paste or email the full code and input file if needed.

Update:

After application of Pressure boundary conditions below:
Left:$p=1$
Right:$p=0$
TopBot:$\frac{\partial p}{\partial n}=0$
The results seem to look nicer, however, diverge, issue might be in the vector field applied improperly

Answer 1

The smooth solution turned out to have BC's applied in the following way:

Walls and inlet: $\frac{\partial p}{\partial n}=0$
Outlet: $p=0$

Actually thought that we need only one value of P to pin, not really outlet P = 0; can take any point, the reason behind is that infinite number of answers are imposed by all Neuman conditions, need one point to fix the solution. But for some reason this is not the case.