# How to approximate flux (with gradient) when using finite volumes?

In finite volume method one is using cell averages. In nonlinear conservation laws discontinuities can be created in the solution process. How to compute the flux when the flux contains a gradient and the Taylor approximation is not valid due to discontinuities?

-

Non-degenerate diffusion will prevent discontinuities in the true solution, but those may be poorly resolved. I will assume that your question is about how to evaluate diffusive fluxes of the form $f(u,\nabla u,p) = \kappa(u,\nabla u,p)\nabla u$ where $\kappa$ is a $d\times d$ SPD tensor and $p$ represents some spatially distributed parameters that are cell-centered and discontinuous. There are a few approaches.

### Two-point fluxes

Take the centroids $c_0$ and $c_1$ of the two cells and compute $\nabla u\cdot(c_1 - c_0) = u(c_1) - u(c_0)$. If the mesh is orthogonal (with face quadrature point on line segment connecting cell centroids and normal with aligned normal) and the coefficients are constant and isotropic, this is a consistent second order accurate flux (i.e., superconvergent at centroids and face centers). Unfortunately, it loses accuracy and has stability problems on general meshes and for variable coefficients. Two-point fluxes are often used because they are simple and compact, but it is a fragile methodology since the grids regularity and constant coefficient assumptions are so frequently violated.

The underlying principle is to reconstruct a gradient on each cell, typically with a weighted least squares procedure, and use that gradient evaluated at quadrature points (face centers) to define a flux. This procedure works simply for nonlinear fluxes. There are many variants that you can find in the literature, but you may want to look at these notes for background and the most commonly used variants. The downside of this approach is that the stencil is larger than strictly necessary and the resulting linear systems are generally not symmetric or M-matrices (thus allowing oscillations).

### Reinterpretation as mixed finite element methods

Mixed finite element methods explicitly discretize the flux functions and work with inherently discontinuous (perhaps tensor-valued) coefficients and solution variables. The downside of mixed FEM is that there are more degrees of freedom and the resulting system of equations is a saddle point problem (which makes it more difficult to solve). After a special choice of quadrature, however, the "flux" degrees of freedom can be eliminated resulting in a (usually symmetric) positive definite system for the cell centered variables. An example of this is detailed in Wheeler and Yotov (2006) which uses an underlying BDM-1 space (linear discontinuous flux functions on edges, constant pressure functions on cells) and a cell quadrature at cell corners which makes the mass matrix block diagonal (with one block per corner). After eliminating the fluxes, the resulting system is SPD and has a minimal stencil (only connecting cells that share a corner). It is second order superconvergent for tensor coefficients and general grids, though the fluxes are only first order accurate for irregular grids.

Note that in addition to the relationship between finite volume and mixed finite element, there is a duality between mixed finite element methods and mimetic finite difference methods. In all cases I know of, it is possible to choose basis functions and quadratures to make the methods coincide, though one or the other may be more natural for certain types of analysis or for implementation in a given software package.

-
+1 For the notes, thanks! :) –  tmaric Jan 5 '13 at 21:11
My flux contains also a convective part, but I was interested in how to treat the diffusive part. Thanks for your answer so far. Do you think it's possible to use slope limiters as gradient reconstruction? The idea was arising as I was reading about high-resolution schemes today which instead of using a constant cell value use a linear interpolation which is slope limited to limit oscillations. –  delta5 Jan 5 '13 at 21:47
@delta5 Inviscid fluxes are computed using a Riemann solver of some sort. In general, you should not try to use the "same procedure" for the diffusive fluxes. (That is, keep the Riemann solver and the diffusive flux routines separate rather than piling them into one big jumble.) Depending on the form of $\kappa(u, \nabla u, p)$, you may need to use a limited reconstruction when computing $\kappa$ (note that typically there is only one value of $\kappa$ per cell). –  Jed Brown Jan 6 '13 at 1:15
So I propose to use a Riemann solver for the convective part and compute the diffusive flux with a two point flux and add those fluxes together? My equations are like that: $\frac{\partial d}{\partial t} = - \nabla \cdot ( \vec e(d) d - f \nabla d)$ –  delta5 Jan 6 '13 at 10:42
@delta5 Good, that is the standard approach. If $f$ is variable or tensor-valued, or if your mesh is irregular, you may want to consider replacing the two-point flux with one of the more sophisticated methods I mentioned. –  Jed Brown Jan 6 '13 at 16:07