I want to solve this set of 2D advection-diffusion equations of this form in spherical coordinates:

$$ \frac{\partial f}{\partial t}=-\mathbf{u}\cdot\nabla f+\eta\nabla^2f+\eta_1(\mathbf{e}_1\cdot\nabla)(\mathbf{e}_1\cdot\nabla)f+\eta_2(\mathbf{e}_2\cdot\nabla)(\mathbf{e}_2\cdot\nabla)f+s_f(g) $$

$$ \frac{\partial g}{\partial t}=-\nabla\cdot(\mathbf{u}g)+\nabla\cdot(\eta\nabla g)+\nabla\cdot(\eta_1\mathbf{e}_1\mathbf{e}_1\cdot\nabla g)+\nabla\cdot(\eta_2\mathbf{e}_2\mathbf{e}_2\cdot\nabla g)+\nabla\cdot h(f)+s_g(g) $$ where the $\mathbf{e}$ are unit vectors not parallel to the unit vectors in spherical coordinates. Most of the terms depend in a nonlinear manner on $f$ and or $g$, so an explicit algorithm is necessary.

I was thinking of finite volumes because most of the terms in the equation for $g$ are in conservative form, but for the equation for $f$, none of them are. I mean, I guess I can do it anyway, so my question is this: is there a way to take maximal advantage of the finite volume method here? If both equations are not in conservative form, I suspect I will loose a lot of the advantages of having the other equation in such a form? I thought that it could be advantageous to get rid of one derivative order for f by writing (which also works for the directional gradients because $\nabla\cdot\mathbf{e}=0$): $$ \eta\nabla^2f=\nabla\cdot(\eta\nabla f)-\nabla\eta\cdot\nabla f. $$ Or is the finite volume even the preferred method here?

  • $\begingroup$ What do you know about the advective field? Is there something we should know about the behaviour of $\mathbf{u}(r \to 0)$ ? $\endgroup$
    – MPIchael
    Commented Feb 10, 2022 at 8:19
  • $\begingroup$ The domain is an annulus, and the $\mathbf{u}\rightarrow 0$ there. It does also depend on $f$ and $g$ $\endgroup$ Commented Feb 10, 2022 at 13:42

1 Answer 1


Whether the finite-volume method is preferred for your problem depends on what requirements you have for the solution. For example if geometric flexibility and local mesh refinement are important then it is hard to beat finite elements. Or, if you need very high accuracy order, in a simple domain, then the choice would probably be spectral or spectral elements. But if enforcing the conservation law for the quantity $g$ is important for you, then finite volumes would be the choice. For finite volume implementation, what you'd need to do is to write your equation for $g$ as

$ \partial_t g = -\nabla \cdot \vec{G} + s_g, $

where $\vec{G}$ is the sum of those terms under the divergence operator in the right-hand side of the PDE for $g$.

Then, using consistent calculation of those fluxes $\vec{G}$ though the cell faces, you'd end up with a scheme conserving $g$ to the machine accuracy, e.g., in the steady state the flux $\vec{G}$ integrated over the domain (or any subdomain) boundary will be equal to the source term $s_g$ integrated over the domain volume.

For the other PDE, you are not going to have a similar conservation law there, because the equation is not in the conservative form. All those terms on the right-hand side of that equation can be viewed just as a source term for $f$, which can be perfectly well implemented in a finite volume method.

  • $\begingroup$ I don’t think I’m going to get into local mesh refinements. What I want, is a good balance between accuracy, speed and stability. The problem is highly non-linear and time-depend, so I guess I need a rather refined time resolution. I already intend, if possible, to have a grid equally spaced in $\cos{\theta}$ in the $\theta$-direction and Chebyshev polynomials in the $r$-direction. I think finite volumes could help at least with the $g$-equation with stability by being conservative. It could also be good from a physical perspective, since I am interested in the ‘recycling’ of $g$. $\endgroup$ Commented Feb 10, 2022 at 13:51
  • $\begingroup$ Or I assume I can use different schemes for each equation? $\endgroup$ Commented Feb 10, 2022 at 13:58
  • $\begingroup$ Conservative formulation is needed for enforcing physical conservation laws. If there is a conservation law in underlying PDEs then any consistent spatial discretization method will reproduce those laws in the limit of fine grid resolution. But a conservative finite volume method would enforce those laws as an algebraic identity for $any$ grid resolution. Numerical stability is something independent of conservation properties, there is no direct connection. But sometimes enforcing conservation laws helps keep the numerical solution "physical". $\endgroup$ Commented Feb 10, 2022 at 16:35
  • 2
    $\begingroup$ If the solution blows up developing large positive and negative spikes, that may still be consistent with the conservation law in the system because those spikes would average out. So, no - conservation properties of the solution is not something that would prevent it from blowing up. What you need is stability analysis to prevent blowing up, that's separate from conservation properties. $\endgroup$ Commented Feb 10, 2022 at 18:51
  • 1
    $\begingroup$ @BitterDecoction You might be thinking of the shock-capturing schemes that are sometimes used with finite volume methods (e.g., WENO, MUSCL with flux limiters). In these cases, you can prevent non-physical oscillatory behaviors that otherwise still satisfy the conservation law. Still, if that's your goal then there are finite-difference WENO varients. In the end, the principal use of finite volume schemes is to enforce conservation laws, as MaximUmansky says. $\endgroup$
    – user20857
    Commented Mar 12, 2022 at 11:12

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.