Solving a set of mixed conservative/non-conservative equations with the finite volume method

Question

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?

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.