Let $\rho$ and $p$ be 2D scalar fields. How do I construct a symmetric Matrix $A$ so that is satisfies
\begin{equation} Ap = \nabla \cdot (\rho\nabla p) \end{equation}
in finite differences? I am using a staggered grid, where $p$ is stored in the center of the cells, and $\rho$ is stored in both the center and the edges.
EDIT1:
So I tried the method suggested by User7530, and what I get are these coefficients:
$p_{i,j}$: $-\frac{4\rho_{i,j}}{h^2}$
$p_{i+1,j}$: $\frac{\rho_{i,j}}{h^2} + \frac{\rho_{i+1/2,j} - \rho_{i-1/2,j}}{2h^2}$
$p_{i-1,j}$: $\frac{\rho_{i,j}}{h^2} - \frac{\rho_{i+1/2,j} - \rho_{i-1/2,j}}{2h^2}$
$p_{i,j+1}$: $\frac{\rho_{i,j}}{h^2} + \frac{\rho_{i,j+1/2} - \rho_{i,j-1/2}}{2h^2}$
$p_{i,j-1}$: $\frac{\rho_{i,j}}{h^2} - \frac{\rho_{i,j+1/2} - \rho_{i,j-1/2}}{2h^2}$
The problem now is that, since I want the matrix to be symmetric, shouldn't the coefficient for $p_{i+1,j}$ be the same as for $p_{i-1,j}$?
There is not one single way of discretizing this equation with finite differences; most approaches however will begin by applying the product rule to get
$$\nabla \cdot (\rho \nabla p) = \nabla \rho \cdot \nabla p + \rho \Delta p.$$
On a regular grid, now you can discretize the first term using e.g. central differences and the second term using e.g. the five-point stencil.
EDIT: If you have $\rho$ on the edges then you can get a symmetric expression by directly discretizing the original expression, as suggested by Winther:
$$\frac{\rho_{i+1/2,j}(p_{i+1,j} - p_{i,j}) + \rho_{i,j+1/2}(p_{i,j+1}-p_{i,j}) + \rho_{i-1/2,j}(p_{i-1,j}-p_{i,j})+\rho_{i,j-1/2}(p_{i,j-1/2}-p_{i,j})}{h^2}$$
Note that this expression agrees with the one you would get by finite volumes (draw a "dual" cell centered at $p_{i,j}$ and use Stokes's Theorem).