How do we correctly define the flux in a finite volume method applied to Poisson's equation where we have a piecewise constant material? Specifically, say we have the equation \begin{align*} -\nabla\cdot(k\nabla u) =& f\textrm{ on }\Omega\\ u =& 0\textrm{ on }\partial\Omega \end{align*} where $\Omega$ is some kind of easy, square domain in $\mathbb{R}^2$. Then, let's mesh $\Omega$ using square elements where one of the cells is $K$. Then, on the interior, we can integrate over $K$, which gives $$ \int_{K}-\nabla\cdot(k\nabla u) = \int_{K}f $$ After using the divergence theorem, we have $$ \int_{\partial K}-k(\nabla u\cdot n) = \int_{K}f $$ Now, we have to discretize the quantity on the left, which can be divided into four separate integrals. Assume that center of cell $K$ is $(x_i,y_j)$ and let $u_{ij}=u(x_i,y_j)$ and $k_{ij}=k(x_i,y_j)$. Then, we can discretize $\nabla u \cdot n$ into four different pieces, depending on what side of the square we're on, as \begin{align*} &\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\\ &\frac{u_{ij}-u_{{i+1},j}}{\Delta x}\\ &\frac{u_{ij}-u_{i,{j-1}}}{\Delta y}\\ &\frac{u_{ij}-u_{i,{j+1}}}{\Delta y} \end{align*} What I'm not clear about is what to do about $k$ if we assume that $k$ is piecewise constant and defined on the same mesh as $u$. We could just multiply the fluxes by $k_{ij}$. For example, discretize $-k(\nabla u\cdot n)$ along the left face as $$ k_{ij}\left(\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\right) $$ However, in this case, the flux on the other side would be $$ k_{{i-1},j}\left(\frac{u_{{i-1},j}-u_{i,j}}{\Delta x}\right) $$ so they wouldn't match and be conservative because they're using two different material values. The best I can think of is to average the material and use $$ \left(\frac{k_{ij}+k_{{i-1},j}}{2}\right)\left(\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\right) $$ but I'm not sure if this is correct or the implications of this choice.
As such, how do we discretize the flux when we have a piecewise constant material, $k$, defined on the same mesh as our material $u$.
