I'm currently trying to implement my own FVM code in cpp, but when I try to calculate the laplacian of a test function, given by \begin{align}\phi_0=\sin(2\pi x)\sin(2\pi y),\end{align} I get unphysical boundary values as seen in the figure. I used a Dirichlet boundary condition with 0 value on every boundary.
My code follows the methods given in "The Finite Volume Method in Computational Fluid Dynamics" by Moukalled, F. et al. according to which the laplacian (diffusion in this case) of a field is discretized for a cell using \begin{align}\int_V\nabla\cdot(\Gamma\nabla\phi)dV\approx\sum_{f_i=1}^{n_{f,i}}\Gamma_{f_i}\frac{\mid\textbf{S}_f\mid}{\mid\textbf{d}_n\mid}(\phi_N-\phi_C)+\sum_{f_b=1}^{n_{f,b}}\Gamma_{f_b}\frac{\mid\textbf{S}_f\mid}{\mid\textbf{d}_n\mid}(\phi_b-\phi_C),\end{align} where $\phi_C$ is the field value at the cell center, $\phi_N$ is the field value at the neighbour face(s), $\phi_b$, is the field value at the boundary, $\mid\textbf{S}_f\mid$ is the magnitude of the surface normal of the face, $\mid\textbf{d}_n\mid $ is the distance between the cell centers of the owner cell and the neighbour cells.
The above grid represents a test case with equal distances between cell centers ($\delta x$ and $\delta y$) and equal grid spacing ($\Delta x$ and $\Delta y$). So, as an example, the discretization of element 0 in the above grid is given by the algebraic relation \begin{align}a_C\phi_C+\sum a_N\phi_N=b_C , \end{align} where \begin{align}\sum a_N\phi_N=a_{Nn}\phi_n+a_{Ne}\phi_e+a_{Ns}\phi_s+a_{Nw}\phi_w=\Gamma\Bigg(\frac{\Delta x}{\delta y}\phi_n+\frac{\Delta y}{\delta x}\phi_e+0\phi_s+0\phi_w\Bigg)=\Gamma\phi_3+\Gamma\phi_1 ,\end{align} and
\begin{align}a_C = a_{Cn}+a_{Ce}+a_{Cs}+a_{Cw}=\Gamma \Bigg[-\frac{\Delta x}{\delta y}+ \Bigg(-\frac{\Delta y}{\delta x}\Bigg)+\Bigg(-\frac{\Delta x}{\frac{1}{2}\delta y}\Bigg)+\Bigg(-\frac{\Delta y}{\frac{1}{2}\delta x}\Bigg)\Bigg]=-6\Gamma\end{align}
$b_c$ is a zero vector due to the boundary conditions. Repeating this for all cells we get the following matrix of coefficients
\begin{align}\textbf{A}=\Gamma\begin{bmatrix}-6&1&0&1&0&0&0&0&0\\1&-5&1&0&1&0&0&0\\0&1&-6&0&0&1&0&0&0\\1&0&0&-5&1&0&1&0&0\\0&1&0&1&-4&1&0&1&0\\0&0&1&0&1&-5&0&0&1\\0&0&0&1&0&0&-6&1&0\\0&0&0&0&1&0&1&-5&1\\0&0&0&0&0&1&0&1&-6\end{bmatrix}\end{align}
The coefficients at the boundaries are larger since the cells centers are closer to the faces than the interior point(s). This is the exact same matrix I get from my code but the resulting field is incorrect as seen in the figure below
I don't think there is anything wrong with my code because the resulting coefficient matrix for the test case with grid size 3x3 agrees with my hand calculations. I don't want to use ghost cells either, since I want to have the possibility of generalizing the grid. Maybe some of you can shed some light on my problem.
I guess my main question is that is this a valid way of implementing the FVM or have I missed something? I wasn't sure what to include in this post so it might be lacking, but please do ask for more information in that case.
Also, I know that I haven't been thorough in explaining all the nomenclature but as most of it is pretty standard notation I think you'll manage.
