Do we really have a fourth order scheme?
 

Hi,
I'd like to find out an answer to some doubts regarding the cubic scheme.
The aim is to evaluate in 1D fvc::div(phi, T) with phi=constant flux .
This should give the gradient in 1D. For T I'm using a Gaussian centered
at 0.5 and 0<=x<=1. The Gaussian is narrow enough to fit to the domain.
OpenFOAM manual claims that the cubic scheme is fourth order.
I'm not getting this decay rate but more likely second order decay with
slightly lower error (also dispersion in advection problem is a bit reduced in my experience).
In 1D Cartesian grid (uniform) a finite volume scheme can be viewed
as a finite difference scheme and thereby the fourth order decay could be
obtained here. I have some doubts that the cubic is not really fourth order or
then I'm doing something wrong.

My settings in fvSchemes are:

divSchemes
{
default Gauss cubic;
}

gradSchemes
{
default Gauss cubic;
}

interpolationSchemes
{
default cubic;
}

snGradSchemes
{
default cubic corrected;
}

fluxRequired
{
default no;
p ;
}


Any comments? I know that it is possible to construct a "good" fourth order
scheme "easily" on unstructured grids by using face neighbor cell gradients
i.e. essentially linear interpolation + correction based on the gradients. Is
this what the cubic is aimed at?
I'd thereby expect that constructing a fourth order scheme in OF is very possible
and I'm wondering why we are not seeing the expected behavior.

Best,
Ville

Hi,
your question has a trick ... FD and FV are not the same, also on uniform grids, if you go to accuracy greater than 2.
The difference is in the resolved variable. For example:

FD -> d /dt (phi) + Div (f) = 0

FV -> d /dt Int [V] (phi) dx + Int [SV] n.(f) dS = 0

You can see that in FV you solve for:

phi_bar = (1/|V|) Int [V] (phi) dx

which is analytically a second order approximation of phi (for a point centered around V).

In conclusion, while in FD if you discretize space derivative to fourth order accuracy then is correct to see a fourth order slope in the error decay, in FV when you use a cubic reconstruction for the flux (fourth order accurate FV) is correct to see a second order slope in the error because you are comparing phi_bar to phi. If you want to see in FV a real fourth order slope you have to average the exact solution.
This is because you compute the discretization error, not directly the local truncation error of the scheme

Note that some FV-based ENO scheme are built in such a way to be equivalent to a fourth order FD derivative and you see an effective fourth order convergence.

This issue is quite misleaded in the literature, but for example the book of Ferziger and Peric said clearly the difference between FD and FV, as well as some infos are in the book of LeVeque. I have also more references on this issue.

I hope to be of some help ....

The OPs question is correctly posed. He is putting a smooth function and computing its derivative. He is not solving a PDE at this stage.

He is doing this in 1d on uniform grids.

He should see the theoretical convergence rates and there is no issue of FD versus FV in this case.

Actually it is not easy to construct stable and efficient higher than second order schemes on unstructured grids with FVM.

Hi,
that is interesting. So what you say is that a "more right" way of formulating
the problem would be something like this:

FV-method does not give us an approximation for dg/dx (g=Gaussian here)
but more likely to dg_ave/dx . And to compare correctly means having to smooth
the Gaussian a posteriori - assuming that the order is higher than 2.

I'm not completely eating this yet but let me try it out;)

-Ville

Yes, you got it ;) if you have some like this

[f(i+1/2) - f(i-1/2)]/dx (1)

in a FV formulation this means that you reconstruct the fluxes in i+1/2 and i-1/2 also with high order reconstruction (e.g. a cubic interpolation) but (1) is a discrete approximation of the surface integration that balances d/dt(phi_bar). If want, you can check the convergence error on the exact flux and you will recover a fourth order accuracy.

Contact me if you need more referrences.



Conversely

Another try
 

It is unbelievable but it indeed seems that it might be important to really
explicitly consider the filtering. I repeated the test and stated that
we are not seeking for approximation for dg/dx but instead for volume
averaged quantity dg_ave/dx = (1/dx) * \int dg/dx which has the
explicit "solution" for cell i i.e. analytically dg_ave/dx = 1/dx * (g(i+1/2) - g(i-1/2) )
g=Gaussian function.
What we note is that by comparing the finite volume approximation to this
analytical expression
we get quite a different outcome and we note the slope 4 for the cubic scheme.
Although at the finest resolution the decay rate settles down a bit deviating from
the 4th order decay.
Quite a miracle and I'm still not eating this completely ;)
Btw, I'd really appreciate if you could send me some more references on the
topic. Here or to my email ville.vuorinen at aalto.fi

This really makes me think a bit further on the FV-method. The volume averaging
as such does not filter enough in high-Reynolds number flow when central differences
are applied because the diagonal is missing from the stencil (well, on uniform grids).

Best,
Ville

I stand corrected. The averaged derivative can indeed be shown to be fourth order accurate by doing a truncation error analysis, at least in the interior. I dont know what openfoam would do at boundaries. It may not be fourth order accurate at the boundaries.

I am happy you found the trick :) FV and FD are actually different things!

If you are interested:

- G. De Stefano, F.M. Denaro, G. Riccardi High Order Filtering Formulation for Control Volumes Simulation of Turbulent Flows, Int. J. Num. Methods in Fluids, 37, 7, 2001.
- P.Iannelli, F.M. Denaro, G. De Stefano,A Deconvolution-based Fourth Order Finite Volume Method for Incompressible Flows On Non-Uniform Grids Int. J. Num. Methods in Fluids, 43, 4, 2003, 431-462.
- F.M. Denaro, What does Finite Volume-based implicit filtering really resolve in large-eddy simulations? J. Comp.Phys., 230, 10, May 2011, 3849-3883, DOI 10.1016/j.jcp.2011.02.011

I absolutely agree that very few authors stress the differences between FV and FD methods, and in particular almost no one underlines that FV methods introduce two sources of approximation - integration and interpolation. I had read this only in the book by Ferziger & Peric. Which Leveque's book are you making reference to?

• Finite Volume Methods for Hyperbolic Problems, Cambridge University Press (2002). ISBN 0-521-00924-3[4]
• Finite Difference Methods for Ordinary and Partial Differential Equations, Steady State and Time Dependent Problems, SIAM (2007). ISBN 978-0-89871-629-0

Note that the accuracy order of a FV scheme is defined on the basis of the local truncation error of the flux reconstruction, which differs much from the accuracy order of a FD scheme, based on the local truncation error of the derivatives

Lots of reading for tomorrow.

Does anyone have further info on what the cubic scheme in OF is actually
is based on?

I assume that what OF is doing for cubic in the background is to use the explicitly evaluated (second order) gradients and the flux would become something like

0.5*( T_i + T_i+1) + 0.5*(grad(T_i) dot (Rface - Ri) + grad(T_{i+1}) dot (Rface -R_{i+1} )

i.e. uses the gradient/Jacobian to correct the flux on the cell face center.
In 1D central (finite) difference framework (2nd order) this strategy would lead
eventually to fourth order approximate for the derivative which is seen easily.
The source for cubic.H could be something like that (since gaussGrad is used
in a mystical way) but the correction is encrypted into the syntax.

-V

Thanks. I know pretty well the first one (excellent book :)), but didn't remember where such topics were discussed in the text. For those interested, section 8.1 should be a good start.

One useful reference could be

http://www.tandfonline.com/doi/abs/1...07790500314947

in which the vairous schemes are catalogued with the same OpenFOAM's logic. I'm investigating OF numerical schemes too, so let's keep each other updated in case we find out something :)

Francesco

Just an update on convergence rate
 

To avoid any confusion I attach a final picture of L2-error. The error cleanly
goes down with slope=4 (earlier it didn't, purely a matter of boundary effects
which are avoided by making the Gaussian more narrow so that it really
fits into the domain).
-V

good to see that OF works!
As Francesco said, the matter of distinguish FV from FD is quite unexplored, I tried in some papers to give inlights since also in the book of Leveque is necessary to "extrapolate" the concept.
This matter become relevant both for peaople working with Euler flows and thos working on LES (FV are a way to filter)

Regards
Filippo