I'm currently working with this kind of mesh (AMR) : enter image description here

To compute the gradient at the center of the cell P, I use the following formula : $$\nabla \phi_P = \frac{1}{\Omega} \sum_{faces} \phi_f \vec{S_f}$$ where $\Omega$ is the volume of the cell and $\phi_f$ is the approximation of the flux on the interface between cell P and its neighbor. The simplest choice for computing $\phi_f$ is to take the average of the cell-centred variables : $$\phi_f = \frac{1}{2}(\phi_P + \phi_N) ~~~~(1)$$.

However I am wondering if this approximation of $\phi_f$ can be applied in my example above. I read here that the arithmetic average works only for constant mesh size. For non conform grid like mine, it suggests to take a weighting factor depending on the geometry : $$\phi_f = \alpha \phi_P + (1-\alpha) \phi_N$$ where : $$\alpha = \frac{\vec{r_N}-\vec{r_f}}{\vec{r_N}-\vec{r_P}}$$

That said, it makes sense to me. Then, my question is : do I make a huge mistake by considering an arithmetic average ? If I have correctly understood, the correct approximation would be, in my case ($\alpha = \frac{1}{3}$) : $$ \phi_f = \frac{1}{3} \phi_P + \frac{2}{3} \phi_N~~~~(2)$$ The solver is stable though, but do you think I lose accuracy with (1) instead of (2) ?

  • $\begingroup$ Can you confirm that your grid is orthogonal as pictured and can you say anything about what you're using the gradient information for? There are a lot of methods for generating reconstructions from data and some are more appropriate than others in certain situations. $\endgroup$
    – origimbo
    Sep 1, 2016 at 18:26
  • $\begingroup$ The mesh is orthogonal, the elements are only squares or cubes (no triangle or pyramid) and I apply the 2:1 constraint (no more than 2 small cells for a neighbouring coarse cell, 4 in 3D).The gradient is used to get second order in space in a finite volume scheme (MUSCL reconstruction) $\endgroup$
    – Coriolis
    Sep 1, 2016 at 19:18

2 Answers 2


For an orthogonal grid, using a proper linear interpolation (2) for the two point flux approximation in this case will result in a theoretical second order accuracy in the local equation truncation error in the calculation of the gradient (the first order error error term arriving from face f ends up cancelling with the error term arising from using the same formula with the face above it).

However you also say that you are implementing a MUSCL scheme, implying that you are performing some flux limiting. This implies that your global error will generally be dominated by the necessity for a first order scheme in the vicinity of shocks to damp spurious oscillations. This makes makes a priori analysis difficult. If you're interested enough and have enough time, you could perform a mesh convergence study for the actual problem you're interested in to see how much difference it actually makes.

Editing in some raw mathematics:

Define the point $q$ to be the midpoint of the full righthand face of cell P and label the cell above cell $N$ as cell $M$. Then I claim that using eqn (2) to give the face values for the two sub surfaces in your numerical integration is equivalent to taking

$\phi_q = \frac{1}{3}\phi_P + \frac{2}{3} \frac{\phi_N +\phi_M}{2}$

and integrating over the whole face. Let $\phi$ be a sufficiently smooth exact solution, then applying Taylor series expansions

$\phi(x_P) = \phi(x_q)-\left. h \frac{\partial \phi}{\partial x}\right|_{x_q}+\left.\frac{h^2}{2} \frac{\partial^2 \phi}{\partial x^2}\right|_{x_h}$

$\phi(x_M) = \phi(x_q)+\left.\frac{h}{2}\frac{\partial \phi}{\partial x}\right|_{x_q}+\left.\frac{h}{2}\frac{\partial \phi}{\partial y}\right|_{x_q}+\mathcal{O}(h^2)$

$\phi(x_N) = \phi(x_q)+\left.\frac{h}{2}\frac{\partial \phi}{\partial x}\right|_{x_q}-\left.\frac{h}{2}\frac{\partial \phi}{\partial y}\right|_{x_q}+\mathcal{O}(h^2)$

And simple substitution shows

$\phi_q = \phi (x_q) + \mathcal{O}(h^2)$.

  • $\begingroup$ True. In all cases, I will need a convergence study. What I mean is that from a strictly geometrical point of view, using formulation (1) for the mesh above is the same thing than considering $\phi_N$ located on the right face center and not at the center of cell N, which is obviously false according to my sketch. I'm trying to quantify this error to see if it is absolutely necessary to use (2) instead of (1). $\endgroup$
    – Coriolis
    Sep 2, 2016 at 15:36
  • $\begingroup$ Absolutely necessary in which sense? For point based flux approximations, the default is usually to apply a Taylor series analysis applied to an exact solution. This will show first order terms error terms proportional $h$ times the partial derivatives of $\phi$, the exact solution for both (1) and (2). However for (2) [but not (1)] this term will cancel over the entire right hand face of cell P. giving a second order term term proportional to $h^2$ times the second derivatives of $\phi$. $\endgroup$
    – origimbo
    Sep 2, 2016 at 16:28
  • $\begingroup$ If not cumbersome, may I ask you to write the demonstration (or part of) to get the overall idea? Your comment makes sense to me, I think I can understand with some equations. $\endgroup$
    – Coriolis
    Sep 2, 2016 at 18:38
  • $\begingroup$ This is one of those annoying areas where it's easy to give a flavour in an answer suitable for stack exchange, but to work through all the details from scratch would take a chapter of a textbook (of which I've certain many have been written). I've added some text and equations to the answer to give an idea of where this goes, but this does skip steps like showing the numerical quadrature in the integrals is accurate to second order. $\endgroup$
    – origimbo
    Sep 5, 2016 at 13:43
  • $\begingroup$ Your edit is very helpful. According to that, equation (1) indeed only gives an approximation in $O(h)$. Does this mean that my MUSCL reconstruction is in fact first order instead of second order ? $\endgroup$
    – Coriolis
    Sep 5, 2016 at 14:40

You can use area averaging (where weight is cell area) technique to calculate the variable values at the cell vertices. The variable value at the edge center can be calculated just by taking simple average of the values at the nodes of that edge. We in our group commonly use this technique for gradient calculation using Green-Gauss approach. Refer section 3.4 of this paper for more information. All the best.


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.