Gradient computation in AMR framework (Green-Gauss theorem)

Question

I'm currently working with this kind of mesh (AMR) :

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) ?

Answer 1

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)$.

Answer 2

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.