How could I go about solving the continuity equation below in 2D?
$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$$
I saw that a similar question was posted here: A good finite difference for the continuity equation, but that only discussed the 1-D case.
In it, one of the answers defines incoming, and outgoing fluxes, as $$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}\\ \Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$$
With these, one could use the finite difference scheme
$$\dfrac{\partial(\rho)}{\partial t}=-\dfrac{\partial(\rho u)}{\partial x} = -\dfrac{\partial \Phi}{\partial x} \approx -\dfrac{\Delta \Phi}{\Delta x} \approx -\dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$$
The problem is that when I use a naive generalization of this, namely defining
$$\Phi_{i+1/2,j;x} = \dfrac{u_{i+1/2,j;x}+|u_{i+1/2,j;x}|}{2}\rho_{i,j} + \dfrac{u_{i+1/2,j;x}-|u_{i+1/2,j;x}|}{2}\rho_{i+1,j}\\ \Phi_{i,j+1/2;y} = \dfrac{u_{i,j+1/2;y}+|u_{i,j+1/2;y}|}{2}\rho_{i,j} + \dfrac{u_{i,j+1/2;y}-|u_{i,j+1/2;y}|}{2}\rho_{i,j+1}$$
and taking
$$\dfrac{\partial (\rho u_x)}{\partial x} \approx \frac{\Phi_{i+1/2,j;x}-\Phi_{i-1/2,j;x}}{\Delta x}$$ and $$\dfrac{\partial (\rho u_y)}{\partial y} \approx \frac{\Phi_{i,j+1/2;y}-\Phi_{i,j-1/2;y}}{\Delta y}$$ my density diverges. What scheme could I use to solve this equation properly in 2D?
