I have part of a problem that is described by the momentum conservation equation:

$\frac{\partial \rho}{\partial t} + \frac{1}{\sin\theta} \frac{\partial}{\partial \theta}(\rho u \sin \theta) =0$

Where $u=f(\theta)$ and $ \rho = f(\theta,t) $ (constant velocity).

Naively one could apply one of the solutions listed here. The problem at hand is best described in spherical polar coordinates (thin spherical shell) and those solutions in Cartesian. must I make some kind of coordinate transformation before discretizing this equation or can I discretize it directly?

Secondly is there any reason one should first expand the derivative in $\theta$ and then attempt to discretize?

As a note - I have done several of the above and have obtained solutions that do not seem to be consistent (physically a couple seem to make sense). I am interested in if there is a proper coordinate transformation that should be made or if any of the previously mentioned methods will suffice.


I define flux as:

$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} \sin{\theta_{i}} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1} \sin{\theta_{i+1}} $

$\Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} \sin{\theta_{i-1}} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i} \sin{\theta_{i}}$

I'm guessing the 'proper' way to define the flux would be to evaluate $\sin{\theta}$ at the $\pm\frac{1}{2}$ 'cell boundaries' rather than at the cell center. This would be more in line with the definition of flux.

One final question - what should I do at the boundaries ($\theta = 0 $ is particularly a problem and I just avoid that point altogether).

  • 2
    $\begingroup$ I don't see $h$ in your equation. Please correct sin\theta ($sin\theta$) in \sin\theta ($\sin\theta$): proper math formatting matters! $\endgroup$
    – Stefano M
    Jul 19, 2012 at 7:14

1 Answer 1


There is no reason you have to transform the coordinates before discretizing the equation. However, since you are using spherical coordinates, you will end up with a non-linear system rather than the linear system that would appear in cartesian coordinates. If you want to work with a linear system, you will have to do variable transformations (much like the creation of $\Phi$ that appears in the first answer of your linked question).

As far as expanding the derivative in $\theta$, I would recommend not to, unless you must for the above transformations. The equation in its current state is in the conservative form, with all $u$'s being coupled with $\rho$. When you expand the derivative you will lose this coupling. The two solutions will converge to each other with finer and finer meshes, however there is a non-zero difference between them. Also, someone else with more experience may be able to elaborate on the advantage/disadvantages of the two forms, since even though it is called the conservative form of the equation, finite difference methods are, by their nature, non-conservative on some level.

  • $\begingroup$ Thanks for the response - I have investigated possible solutions to this problem but have yet to find a solution that seems reasonable. I have the BC's down - but no benchmarks to test my scheme against. $\endgroup$
    – Marm0t
    Jul 23, 2012 at 13:11
  • $\begingroup$ what was your thinking behind the formulation of $\Phi$? It seems that you are choosing your index on $\rho\sin(\theta)$ based on the sign of $u$. Is this for stability? $\endgroup$ Jul 23, 2012 at 17:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.