# A good 2D finite difference for the continuity equation

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?

• What you have written is an upwind scheme which is very robust. But are you solving this coupled to some other equations ? What is the time scheme, time step ? What does the velocity field look like ? – cpraveen Jun 20 '18 at 12:24
• Hey. So I'm using a Runge-Kutta 4th order time integrator, with a time step that's $\Delta t = \Delta x/4$. The density is indeed coupled to another set of equations (Maxwell's equations) that I'm solving using a 4th order in space finite difference scheme (I'm comparing with an steady-state analytical solution, and if I plug in the correct steady-state density, the Maxwell solver works well). The velocity is a function of the electric and magnetic field, $\propto (EB-E-B)/(E^2)$. – grizzlyjoker Jun 20 '18 at 15:03
• where $E > 0$ always – grizzlyjoker Jun 20 '18 at 19:06
• Does your time step satisfy CFL condition ? Make sure that $\Delta t \le \frac{1}{\max(|u_x|/\Delta x + |u_y|/\Delta y)}$ where the max is taken over the whole grid. Also your spatial scheme is first order accurate, so you dont need RK4 in time. Moreover, the first order upwind scheme is too dissipative. – cpraveen Jun 21 '18 at 4:11
• My time step does indeed satisfy the CFL condition. Ok, I'll try getting a higher order upwind scheme to work. – grizzlyjoker Jun 21 '18 at 7:07