I am working on a problem in which I need to numerically solve the continuity equation in an axisymmetric coordinate system (i.e. cylindrical with no $\phi$ dependence). For concreteness, I will use my problem, which deals with an electron number density function $n_e$ and electron velocity function $\vec{u_e}$, but of course this relates to any other continuous medium as well. I came up with two possible schemes: one by directly finite differencing the continuity equation in axisymmetric coordinates, and another by deriving the continuity equation directly. The two schemes end up being subtly different, with the second, as far as I can tell, being preferable. In axisymmetric coordinates, the continuity equation is:
$$
\frac{\partial{n_e}}{\partial{t}}=-\nabla\cdot \left[ {n_e \vec{u_e}} \right]=-\frac{\partial}{\partial{z}}\left[ n_e u_{ez} \right] – \frac{1}{r} \frac{\partial}{\partial{r}} \left[r n_e u_{er} \right]
$$
Where $n_e$ is the electron density and $\vec{u_e}$ is the electron velocity function. Using a space centered finite difference method, we can approximate the r-derivative term as:
$$
\frac{\Delta_r n_{j,l}}{\Delta t}=- \frac{1}{r} \frac{\partial}{\partial{r}} \left[r n_e u_{er} \right] \approx -\frac{1}{r_{j,l}}\frac{1}{2\Delta r} \left[r_{j,l+1} n_{j,l+1}^n u_{rj,l+1}-r_{j,l-1} n_{j,l-1}^n u_{rj,l-1}\right]
$$
While this is a valid scheme, it turns out the r-derivative term does not quite represent the continuity equation as accurately as possible. Lets now consider deriving the continuity equation directly in axisymmetric coordinates (see e.g. Schey 50-52 for a derivation in Cartesian coordinates). In axisymmetric coordinates, our volume elements are essentially concentric hollow cylinders. If we use a five point stencil, then we know the electron density and velocity at each face, as well as relevant information about the density at the point we are updating. Then, the width of the concentric cylinder is $r_{j,l+1}-r_{j,l-1}$ and height is $2\Delta z$.
In $z$, the flux across the bottom and top boundaries can be treated in the normal way because the areas at each face are the same, and this term comes out identical to the $z$ term when we directly finite differenced the continuity equation. However, this is not the case in the $r$ direction. In this case, the flux across the boundary at $r_{j,l-1}$ is:
$$\Phi_{1j,l}=A_{j,l-1}n_{j,l-1}u_{rj,l-1}$$
While that at $r_{j,l+1}$ is:
$$\Phi_{2j,l}=-A_{j,l+1}n_{j,l+1}u_{rj,l+1}$$
The change in the electron density (due to the r-direction) is then the sum of these two fluxes divided by our volume element. From geometry, we have that $A_{j,l-1}=4 \Delta z \pi r_{j,l-1}$ and $A_{j,l+1}=4 \Delta z \pi r_{j,l+1}$ while our volume is $V=2 \Delta z \pi (r_{j,l+1}^2-r_{j,l-1}^2)$. Putting this all together, we find that the r-component of the change in $n_e$ is:
$$\frac{\Delta_r n_{j,l}}{\Delta t} = -\frac{2}{r_{j,l+1}^2-r_{j,l-1}^2} (r_{j,l+1}n_{j,l+1}u_{rj,l+1}-r_{j,l-1}n_{j,l-1}u_{rj,l-1})$$
Which is different from what we found by directly finite differencing the continuity equation. It is clear where this difference comes from: if we take this change in the limit as $\Delta r \to 0$, we can take our volume to be $2 \pi r dr$. Then, when we finite difference this derivative, we are essentially still using this (incorrect) volume, so we get our first expression. When we use the correct volume expression for our finite concentric cylinder, we get the second expression.
It seems to me this adjustment is strictly more accurate than the first expression with no real drawback. Is there a reason to use the first expression? Also, this adjustment naturally arises for the continuity equation because it is obviously about tracking fluxes across surfaces, but I assume it would be a useful adjustment when finite differencing other differential equations in non-Cartesian coordinate systems. What is this adjustment called in computational physics and how is it used in other cases?
Best Answer
Note that, $$\frac{2}{r_{i+1}^2-r_{i-1}^2}=\frac{2}{(r_{i+1}-r_{i-1})\cdot(r_{i+1}+r_{i-1})}=\frac{2}{\Delta r\cdot(r_{i+1}+r_{i-1})}=\frac{1}{\Delta r\cdot\frac{1}{2}(r_{i+1}+r_{i-1})}$$ Then also note that $\frac{1}{2}(r_{i+1}+r_{i-1})=r_i$. Thus, $$\frac{2}{r_{i+1}^2-r_{i-1}^2}\equiv\frac{1}{\Delta r\cdot r_i}.$$ So your usage of the volume and area has only proven the original case is correct.
That said, some hydrodynamics codes I have worked with have opted to express the divergence term using the area and volume approach, $$\nabla\cdot\mathbf{F}\approx\frac{A_{i+1/2}R_{i+1/2}F_{i+1/2}-A_{i-1/2}R_{i-1/2}F_{i-1/2}}{R_i\Delta\mathcal{V}_{ijk}}+\cdots$$ where, in polar coordinates, $\Delta\mathcal{V}=R_i\Delta R_i\Delta\phi_j\Delta z_k$ is the volume. This was done, mind you, to allow the code to be written for arbitrary coordinate systems (namely Cartesian, polar & spherical) without adding/removing functions and/or co-factors.